Commits

Matt Knepley  committed 99dec3a

DMPlex: All quadrilateral geometry tests working
- Had to pass in the coordinate size to the projection routine to make it work

  • Participants
  • Parent commits c68fa55

Comments (0)

Files changed (3)

File src/dm/impls/plex/examples/tests/ex8.c

   /* Create reference triangle */
   dim  = 2;
   ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) dm, "triangle");CHKERRQ(ierr);
   ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
   ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
   {
       DM idm;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) idm, "triangle");CHKERRQ(ierr);
       ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
       ierr = DMDestroy(&dm);CHKERRQ(ierr);
       dm   = idm;
     ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
     ierr = PetscRandomDestroy(&ang);CHKERRQ(ierr);
   }
-
   /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */
   dim = 3;
   {
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "TestQuadrilateral"
+PetscErrorCode TestQuadrilateral(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
+{
+  DM             dm;
+  PetscRandom    r, ang, ang2;
+  PetscInt       dim, t;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  /* Create reference quadrilateral */
+  dim  = 2;
+  ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);
+  ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
+  {
+    PetscInt    numPoints[2]        = {4, 1};
+    PetscInt    coneSize[5]         = {4, 0, 0, 0, 0};
+    PetscInt    cones[4]            = {1, 2, 3, 4};
+    PetscInt    coneOrientations[4] = {0, 0, 0, 0};
+    PetscScalar vertexCoords[8]     = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
+
+    ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+    if (interpolate) {
+      DM idm;
+
+      ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) idm, "quadrilateral");CHKERRQ(ierr);
+      ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
+      ierr = DMDestroy(&dm);CHKERRQ(ierr);
+      dm   = idm;
+    }
+    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+  }
+  /* Check reference geometry: determinant is scaled by reference volume (2.0) */
+  {
+    PetscReal v0Ex[2]       = {-1.0, -1.0};
+    PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
+    PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
+    PetscReal detJEx        = 1.0;
+    PetscReal centroidEx[2] = {0.0, 0.0};
+    PetscReal normalEx[2]   = {0.0, 0.0};
+    PetscReal volEx         = 4.0;
+
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
+  }
+  /* Check random quadrilaterals: rotate, scale, then translate */
+  if (transform) {
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(r, 0.0, 10.0);CHKERRQ(ierr);
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(ang, 0.0, 2*PETSC_PI);CHKERRQ(ierr);
+    for (t = 0; t < 100; ++t) {
+      PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}, trans[2];
+      PetscReal   v0Ex[2]         = {-1.0, -1.0};
+      PetscReal   JEx[4]          = {1.0, 0.0, 0.0, 1.0}, R[4], rot[2], rotM[4];
+      PetscReal   invJEx[4]       = {1.0, 0.0, 0.0, 1.0};
+      PetscReal   detJEx          = 1.0, scale, phi;
+      PetscReal   centroidEx[2]   = {0.0, 0.0};
+      PetscReal   normalEx[2]     = {0.0, 0.0};
+      PetscReal   volEx           = 4.0;
+      PetscInt    d, e, f, p;
+
+      ierr = PetscRandomGetValueReal(r, &scale);CHKERRQ(ierr);
+      ierr = PetscRandomGetValueReal(ang, &phi);CHKERRQ(ierr);
+      R[0] = cos(phi); R[1] = -sin(phi);
+      R[2] = sin(phi); R[3] =  cos(phi);
+      for (p = 0; p < 4; ++p) {
+        for (d = 0; d < dim; ++d) {
+          for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+            rot[d] += R[d*dim+e] * vertexCoords[p*dim+e];
+          }
+        }
+        for (d = 0; d < dim; ++d) vertexCoords[p*dim+d] = rot[d];
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+          rot[d] += R[d*dim+e] * centroidEx[e];
+        }
+      }
+      for (d = 0; d < dim; ++d) centroidEx[d] = rot[d];
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          for (f = 0, rotM[d*dim+e] = 0.0; f < dim; ++f) {
+            rotM[d*dim+e] += R[d*dim+f] * JEx[f*dim+e];
+          }
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          JEx[d*dim+e] = rotM[d*dim+e];
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          for (f = 0, rotM[d*dim+e] = 0.0; f < dim; ++f) {
+            rotM[d*dim+e] += invJEx[d*dim+f] * R[e*dim+f];
+          }
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          invJEx[d*dim+e] = rotM[d*dim+e];
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        ierr = PetscRandomGetValueReal(r, &trans[d]);CHKERRQ(ierr);
+        for (p = 0; p < 4; ++p) {
+          vertexCoords[p*dim+d] *= scale;
+          vertexCoords[p*dim+d] += trans[d];
+        }
+        v0Ex[d] = vertexCoords[d];
+        for (e = 0; e < dim; ++e) {
+          JEx[d*dim+e]    *= scale;
+          invJEx[d*dim+e] /= scale;
+        }
+        detJEx *= scale;
+        centroidEx[d] *= scale;
+        centroidEx[d] += trans[d];
+        volEx *= scale;
+      }
+      ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
+      ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+      if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
+    }
+    ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
+    ierr = PetscRandomDestroy(&ang);CHKERRQ(ierr);
+  }
+  /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (4.0) */
+  dim = 3;
+  {
+    PetscScalar vertexCoords[12] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0};
+    PetscReal v0Ex[3]            = {-1.0, -1.0, 0.0};
+    PetscReal JEx[9]             = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+    PetscReal invJEx[9]          = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+    PetscReal detJEx             = 1.0;
+    PetscReal centroidEx[3]      = {0.0, 0.0, 0.0};
+    PetscReal normalEx[3]        = {0.0, 0.0, 1.0};
+    PetscReal volEx              = 4.0;
+
+    ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
+  }
+  /* Check random quadrilaterals: scale, translate, then rotate */
+  if (transform) {
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(r, 0.0, 10.0);CHKERRQ(ierr);
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(ang, 0.0, 2*PETSC_PI);CHKERRQ(ierr);
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang2);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(ang2);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(ang2, 0.0, PETSC_PI);CHKERRQ(ierr);
+    for (t = 0; t < 100; ++t) {
+      PetscScalar vertexCoords[12] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0}, trans[3];
+      PetscReal   v0Ex[3]          = {-1.0, -1.0, 0.0};
+      PetscReal   JEx[9]           = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}, R[9], rot[3], rotM[9];
+      PetscReal   invJEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+      PetscReal   detJEx           = 1.0, scale, phi, theta, psi = 0.0;
+      PetscReal   centroidEx[3]    = {0.0, 0.0, 0.0};
+      PetscReal   normalEx[3]      = {0.0, 0.0, 1.0};
+      PetscReal   volEx            = 4.0;
+      PetscInt    d, e, f, p;
+
+      ierr = PetscRandomGetValueReal(r, &scale);CHKERRQ(ierr);
+      ierr = PetscRandomGetValueReal(ang, &phi);CHKERRQ(ierr);
+      ierr = PetscRandomGetValueReal(ang2, &theta);CHKERRQ(ierr);
+      for (d = 0; d < dim; ++d) {
+        ierr = PetscRandomGetValueReal(r, &trans[d]);CHKERRQ(ierr);
+        for (p = 0; p < 4; ++p) {
+          vertexCoords[p*dim+d] *= scale;
+          vertexCoords[p*dim+d] += trans[d];
+        }
+        centroidEx[d] *= scale;
+        centroidEx[d] += trans[d];
+        for (e = 0; e < dim-1; ++e) {
+          JEx[d*dim+e]    *= scale;
+          invJEx[d*dim+e] /= scale;
+        }
+        if (d < dim-1) {
+          detJEx *= scale;
+          volEx  *= scale;
+        }
+      }
+      R[0] = cos(theta)*cos(psi); R[1] = sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi); R[2] = sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi);
+      R[3] = cos(theta)*sin(psi); R[4] = cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi); R[5] = cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi);
+      R[6] = -sin(theta);         R[7] = sin(phi)*cos(theta);                              R[8] = cos(phi)*cos(theta);
+      for (p = 0; p < 4; ++p) {
+        for (d = 0; d < dim; ++d) {
+          for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+            rot[d] += R[d*dim+e] * vertexCoords[p*dim+e];
+          }
+        }
+        for (d = 0; d < dim; ++d) vertexCoords[p*dim+d] = rot[d];
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+          rot[d] += R[d*dim+e] * centroidEx[e];
+        }
+      }
+      for (d = 0; d < dim; ++d) centroidEx[d] = rot[d];
+      for (d = 0; d < dim; ++d) {
+        for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+          rot[d] += R[d*dim+e] * normalEx[e];
+        }
+      }
+      for (d = 0; d < dim; ++d) normalEx[d] = rot[d];
+      for (d = 0; d < dim; ++d) {
+        v0Ex[d] = vertexCoords[d];
+        for (e = 0; e < dim; ++e) {
+          for (f = 0, rotM[d*dim+e] = 0.0; f < dim; ++f) {
+            rotM[d*dim+e] += R[d*dim+f] * JEx[f*dim+e];
+          }
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          JEx[d*dim+e] = rotM[d*dim+e];
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          for (f = 0, rotM[d*dim+e] = 0.0; f < dim; ++f) {
+            rotM[d*dim+e] += invJEx[d*dim+f] * R[e*dim+f];
+          }
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          invJEx[d*dim+e] = rotM[d*dim+e];
+        }
+      }
+      ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
+      ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+      if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
+    }
+    ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
+    ierr = PetscRandomDestroy(&ang);CHKERRQ(ierr);
+    ierr = PetscRandomDestroy(&ang2);CHKERRQ(ierr);
+  }
+  /* Cleanup */
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "TestTetrahedron"
 PetscErrorCode TestTetrahedron(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
 {
   /* Create reference tetrahedron */
   dim  = 3;
   ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) dm, "tetrahedron");CHKERRQ(ierr);
   ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
   ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
   {
       DM idm;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) idm, "tetrahedron");CHKERRQ(ierr);
       ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
       ierr = DMDestroy(&dm);CHKERRQ(ierr);
       dm   = idm;
   /* Create reference hexahedron */
   dim  = 3;
   ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) dm, "hexahedron");CHKERRQ(ierr);
   ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
   ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
   {
       DM idm;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) idm, "hexahedron");CHKERRQ(ierr);
       ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
       ierr = DMDestroy(&dm);CHKERRQ(ierr);
       dm   = idm;
   if (transform) {ierr = PetscPrintf(PETSC_COMM_WORLD, "Using random transforms");CHKERRQ(ierr);}
   ierr = TestTriangle(PETSC_COMM_SELF, PETSC_FALSE, transform);CHKERRQ(ierr);
   ierr = TestTriangle(PETSC_COMM_SELF, PETSC_TRUE,  transform);CHKERRQ(ierr);
+  ierr = TestQuadrilateral(PETSC_COMM_SELF, PETSC_FALSE, transform);CHKERRQ(ierr);
+  ierr = TestQuadrilateral(PETSC_COMM_SELF, PETSC_TRUE,  transform);CHKERRQ(ierr);
   ierr = TestTetrahedron(PETSC_COMM_SELF, PETSC_FALSE, transform);CHKERRQ(ierr);
   ierr = TestTetrahedron(PETSC_COMM_SELF, PETSC_TRUE,  transform);CHKERRQ(ierr);
   ierr = TestHexahedron(PETSC_COMM_SELF, PETSC_FALSE, transform);CHKERRQ(ierr);

File src/dm/impls/plex/examples/tests/output/ex8_0.out

-Mesh 'DM_0x84000000_0':
+Mesh 'triangle':
 Max sizes cone: 3 support: 1
 orientation is missing
 cap --> base:
   (   1) dim  2 offset   0 -1 -1
   (   2) dim  2 offset   2 1 -1
   (   3) dim  2 offset   4 -1 1
-Mesh 'DM_0x84000000_0':
+Mesh 'triangle':
 Max sizes cone: 3 support: 1
 orientation is missing
 cap --> base:
   (   1) dim  3 offset   0 -1 -1 0
   (   2) dim  3 offset   3 1 -1 0
   (   3) dim  3 offset   6 -1 1 0
-Mesh 'DM_0x84000000_0':
+Mesh 'triangle':
 Max sizes cone: 3 support: 1
 orientation is missing
 cap --> base:
   (   1) dim  3 offset   0 0 -1 -1
   (   2) dim  3 offset   3 0 1 -1
   (   3) dim  3 offset   6 0 -1 1
-Mesh 'DM_0x84000000_1':
+Mesh 'triangle':
 Max sizes cone: 3 support: 2
 orientation is missing
 cap --> base:
   (   1) dim  2 offset   0 -1 -1
   (   2) dim  2 offset   2 1 -1
   (   3) dim  2 offset   4 -1 1
-Mesh 'DM_0x84000000_1':
+Mesh 'triangle':
 Max sizes cone: 3 support: 2
 orientation is missing
 cap --> base:
   (   1) dim  3 offset   0 -1 -1 0
   (   2) dim  3 offset   3 1 -1 0
   (   3) dim  3 offset   6 -1 1 0
-Mesh 'DM_0x84000000_1':
+Mesh 'triangle':
 Max sizes cone: 3 support: 2
 orientation is missing
 cap --> base:
   (   1) dim  3 offset   0 0 -1 -1
   (   2) dim  3 offset   3 0 1 -1
   (   3) dim  3 offset   6 0 -1 1
-Mesh 'DM_0x84000000_2':
+Mesh 'quadrilateral':
+Max sizes cone: 4 support: 1
+orientation is missing
+cap --> base:
+[0]: 1 ----> 0
+[0]: 2 ----> 0
+[0]: 3 ----> 0
+[0]: 4 ----> 0
+base <-- cap:
+[0]: 0 <---- 1 (0)
+[0]: 0 <---- 2 (0)
+[0]: 0 <---- 3 (0)
+[0]: 0 <---- 4 (0)
+coordinates with 1 fields
+  field 0 with 2 components
+Process 0:
+  (   1) dim  2 offset   0 -1 -1
+  (   2) dim  2 offset   2 1 -1
+  (   3) dim  2 offset   4 1 1
+  (   4) dim  2 offset   6 -1 1
+Mesh 'quadrilateral':
+Max sizes cone: 4 support: 1
+orientation is missing
+cap --> base:
+[0]: 1 ----> 0
+[0]: 2 ----> 0
+[0]: 3 ----> 0
+[0]: 4 ----> 0
+base <-- cap:
+[0]: 0 <---- 1 (0)
+[0]: 0 <---- 2 (0)
+[0]: 0 <---- 3 (0)
+[0]: 0 <---- 4 (0)
+coordinates with 1 fields
+  field 0 with 3 components
+Process 0:
+  (   1) dim  3 offset   0 -1 -1 0
+  (   2) dim  3 offset   3 1 -1 0
+  (   3) dim  3 offset   6 1 1 0
+  (   4) dim  3 offset   9 -1 1 0
+Mesh 'quadrilateral':
+Max sizes cone: 4 support: 2
+orientation is missing
+cap --> base:
+[0]: 1 ----> 5
+[0]: 1 ----> 8
+[0]: 2 ----> 5
+[0]: 2 ----> 6
+[0]: 3 ----> 6
+[0]: 3 ----> 7
+[0]: 4 ----> 7
+[0]: 4 ----> 8
+[0]: 5 ----> 0
+[0]: 6 ----> 0
+[0]: 7 ----> 0
+[0]: 8 ----> 0
+base <-- cap:
+[0]: 0 <---- 5 (0)
+[0]: 0 <---- 6 (0)
+[0]: 0 <---- 7 (0)
+[0]: 0 <---- 8 (0)
+[0]: 5 <---- 1 (0)
+[0]: 5 <---- 2 (0)
+[0]: 6 <---- 2 (0)
+[0]: 6 <---- 3 (0)
+[0]: 7 <---- 3 (0)
+[0]: 7 <---- 4 (0)
+[0]: 8 <---- 4 (0)
+[0]: 8 <---- 1 (0)
+coordinates with 1 fields
+  field 0 with 2 components
+Process 0:
+  (   1) dim  2 offset   0 -1 -1
+  (   2) dim  2 offset   2 1 -1
+  (   3) dim  2 offset   4 1 1
+  (   4) dim  2 offset   6 -1 1
+Mesh 'quadrilateral':
+Max sizes cone: 4 support: 2
+orientation is missing
+cap --> base:
+[0]: 1 ----> 5
+[0]: 1 ----> 8
+[0]: 2 ----> 5
+[0]: 2 ----> 6
+[0]: 3 ----> 6
+[0]: 3 ----> 7
+[0]: 4 ----> 7
+[0]: 4 ----> 8
+[0]: 5 ----> 0
+[0]: 6 ----> 0
+[0]: 7 ----> 0
+[0]: 8 ----> 0
+base <-- cap:
+[0]: 0 <---- 5 (0)
+[0]: 0 <---- 6 (0)
+[0]: 0 <---- 7 (0)
+[0]: 0 <---- 8 (0)
+[0]: 5 <---- 1 (0)
+[0]: 5 <---- 2 (0)
+[0]: 6 <---- 2 (0)
+[0]: 6 <---- 3 (0)
+[0]: 7 <---- 3 (0)
+[0]: 7 <---- 4 (0)
+[0]: 8 <---- 4 (0)
+[0]: 8 <---- 1 (0)
+coordinates with 1 fields
+  field 0 with 3 components
+Process 0:
+  (   1) dim  3 offset   0 -1 -1 0
+  (   2) dim  3 offset   3 1 -1 0
+  (   3) dim  3 offset   6 1 1 0
+  (   4) dim  3 offset   9 -1 1 0
+Mesh 'tetrahedron':
 Max sizes cone: 4 support: 1
 orientation is missing
 cap --> base:
   (   2) dim  3 offset   3 -1 1 -1
   (   3) dim  3 offset   6 1 -1 -1
   (   4) dim  3 offset   9 -1 -1 1
-Mesh 'DM_0x84000000_3':
+Mesh 'tetrahedron':
 Max sizes cone: 4 support: 3
 orientation is missing
 cap --> base:
   (   2) dim  3 offset   3 -1 1 -1
   (   3) dim  3 offset   6 1 -1 -1
   (   4) dim  3 offset   9 -1 -1 1
-Mesh 'DM_0x84000000_4':
+Mesh 'hexahedron':
 Max sizes cone: 8 support: 1
 orientation is missing
 cap --> base:
   (   6) dim  3 offset  15 1 -1 1
   (   7) dim  3 offset  18 1 1 1
   (   8) dim  3 offset  21 -1 1 1
-Mesh 'DM_0x84000000_5':
+Mesh 'hexahedron':
 Max sizes cone: 6 support: 3
 orientation is missing
 cap --> base:

File src/dm/impls/plex/plexgeometry.c

 /*
   DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
 */
-static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscScalar coords[], PetscReal R[])
+static PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
 {
   PetscReal      x1[3],  x2[3], n[3], norm;
-  PetscReal      x1p[3], x2p[3];
+  PetscReal      x1p[3], x2p[3], xnp[3];
   PetscReal      sqrtz, alpha;
   const PetscInt dim = 3;
-  PetscInt       d, e;
+  PetscInt       d, e, p;
 
   PetscFunctionBegin;
   /* 0) Calculate normal vector */
   sqrtz = sqrt(1.0 - n[2]*n[2]);
   /* Check for n = z */
   if (sqrtz < 1.0e-10) {
-    coords[0] = 0.0;
-    coords[1] = 0.0;
     if (n[2] < 0.0) {
-      coords[2] = x2[0];
-      coords[3] = x2[1];
-      coords[4] = x1[0];
-      coords[5] = x1[1];
+      if (coordSize > 9) {
+        coords[2] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
+        coords[3] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
+        coords[4] = x2[0];
+        coords[5] = x2[1];
+        coords[6] = x1[0];
+        coords[7] = x1[1];
+      } else {
+        coords[2] = x2[0];
+        coords[3] = x2[1];
+        coords[4] = x1[0];
+        coords[5] = x1[1];
+      }
       R[0] = 1.0; R[1] = 0.0; R[2] = 0.0;
       R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
       R[6] = 0.0; R[7] = 0.0; R[8] = -1.0;
     } else {
+      for (p = 3; p < coordSize/3; ++p) {
+        coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
+        coords[p*2+1] = PetscRealPart(coords[p*dim+1] - coords[0*dim+1]);
+      }
       coords[2] = x1[0];
       coords[3] = x1[1];
       coords[4] = x2[0];
       R[3] = 0.0; R[4] = 1.0; R[5] = 0.0;
       R[6] = 0.0; R[7] = 0.0; R[8] = 1.0;
     }
+    coords[0] = 0.0;
+    coords[1] = 0.0;
     PetscFunctionReturn(0);
   }
   alpha = 1.0/sqrtz;
   if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
   if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
   /* 2) Project to (x, y) */
+  for (p = 3; p < coordSize/3; ++p) {
+    for (d = 0; d < dim; ++d) {
+      xnp[d] = 0.0;
+      for (e = 0; e < dim; ++e) {
+        xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
+      }
+      if (d < dim-1) coords[p*2+d] = xnp[d];
+    }
+  }
   coords[0] = 0.0;
   coords[1] = 0.0;
   coords[2] = x1p[0];
     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
 
     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
-    ierr = DMPlexComputeProjection3Dto2D_Internal(coords, R);CHKERRQ(ierr);
+    ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr);
     if (J)    {
       const PetscInt pdim = 2;
 
   PetscSection   coordSection;
   Vec            coordinates;
   PetscScalar   *coords;
-  const PetscInt dim = 2;
-  PetscInt       d, f;
+  PetscInt       numCoords, d, f, g;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-  ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
+  ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr);
   *detJ = 0.0;
-  if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
-  if (J)    {
-    for (d = 0; d < dim; d++) {
-      for (f = 0; f < dim; f++) {
-        J[d*dim+f] = 0.5*(PetscRealPart(coords[(f*2+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
+  if (numCoords == 12) {
+    const PetscInt dim = 3;
+    PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
+
+    if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
+    ierr = DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);CHKERRQ(ierr);
+    if (J)    {
+      const PetscInt pdim = 2;
+
+      for (d = 0; d < pdim; d++) {
+        J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
+        J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
       }
+      PetscLogFlops(8.0);
+      Det3D_Internal(detJ, J0);
+      for (d = 0; d < dim; d++) {
+        for (f = 0; f < dim; f++) {
+          J[d*dim+f] = 0.0;
+          for (g = 0; g < dim; g++) {
+            J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
+          }
+        }
+      }
+      PetscLogFlops(18.0);
     }
-    PetscLogFlops(8.0);
-    Det2D_Internal(detJ, J);
-  }
-  if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
-  ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);CHKERRQ(ierr);
+    if (invJ) {Invert3D_Internal(invJ, J, *detJ);}
+  } else if (numCoords == 8) {
+    const PetscInt dim = 2;
+
+    if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
+    if (J)    {
+      for (d = 0; d < dim; d++) {
+        J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
+        J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
+      }
+      PetscLogFlops(8.0);
+      Det2D_Internal(detJ, J);
+    }
+    if (invJ) {Invert2D_Internal(invJ, J, *detJ);}
+  } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %d != 6", numCoords);
+  ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   Vec            coordinates;
   PetscScalar   *coords;
   const PetscInt dim = 3;
-  PetscInt       d, f;
+  PetscInt       d;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
       for (d = 0; d < dim; ++d) normal[d] = 0.0;
     }
   }
-  if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D_Internal(coords, R);CHKERRQ(ierr);}
+  if (dim == 3) {ierr = DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);CHKERRQ(ierr);}
   for (p = 0; p < numCorners; ++p) {
     /* Need to do this copy to get types right */
     for (d = 0; d < tdim; ++d) {