Actual source code: plexrefregular.c

  1: #include <petsc/private/dmplextransformimpl.h>

  3: /* Regular Refinement of Hybrid Meshes

  5:    We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
  6:    to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
  7:    transformations, such as changing from one type of cell to another, as simplex to hex.

  9:    To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
 10:    and the types of the new cells.

 12:    We need the group multiplication table for group actions from the dihedral group for each cell type.

 14:    We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
 15:    we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
 16:    (face, orient) pairs for each subcell.
 17: */

 19: /*@
 20:   DMPlexRefineRegularGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell

 22:   Input Parameters:
 23: + cr - The DMPlexCellRefiner object
 24: - ct - The cell type

 26:   Output Parameters:
 27: + Nf   - The number of faces for this cell type
 28: . v0   - The translation of the first vertex for each face
 29: . J    - The Jacobian for each face (map from original cell to subcell)
 30: . invJ - The inverse Jacobian for each face
 31: - detJ - The determinant of the Jacobian for each face

 33:   Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
 34: $         v0 + j x_face = x_cell
 35: $    invj (x_cell - v0) = x_face

 37:   Level: developer

 39: .seealso: `DMPlexCellRefinerGetAffineTransforms()`, `Create()`
 40: @*/
 41: PetscErrorCode DMPlexRefineRegularGetAffineFaceTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
 42: {
 43:   /*
 44:    2
 45:    |\
 46:    | \
 47:    |  \
 48:    |   \
 49:    |    \
 50:    |     \
 51:    |      \
 52:    2       1
 53:    |        \
 54:    |         \
 55:    |          \
 56:    0---0-------1
 57:    v0[Nf][dc]:       3 x 2
 58:    J[Nf][df][dc]:    3 x 1 x 2
 59:    invJ[Nf][dc][df]: 3 x 2 x 1
 60:    detJ[Nf]:         3
 61:    */
 62:   static PetscReal tri_v0[]   = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
 63:   static PetscReal tri_J[]    = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
 64:   static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0};
 65:   static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0};
 66:   /*
 67:    3---------2---------2
 68:    |                   |
 69:    |                   |
 70:    |                   |
 71:    3                   1
 72:    |                   |
 73:    |                   |
 74:    |                   |
 75:    0---------0---------1

 77:    v0[Nf][dc]:       4 x 2
 78:    J[Nf][df][dc]:    4 x 1 x 2
 79:    invJ[Nf][dc][df]: 4 x 2 x 1
 80:    detJ[Nf]:         4
 81:    */
 82:   static PetscReal quad_v0[]   = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 - 1.0, 0.0};
 83:   static PetscReal quad_J[]    = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
 84:   static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
 85:   static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0};

 87:   switch (ct) {
 88:   case DM_POLYTOPE_TRIANGLE:
 89:     if (Nf) *Nf = 3;
 90:     if (v0) *v0 = tri_v0;
 91:     if (J) *J = tri_J;
 92:     if (invJ) *invJ = tri_invJ;
 93:     if (detJ) *detJ = tri_detJ;
 94:     break;
 95:   case DM_POLYTOPE_QUADRILATERAL:
 96:     if (Nf) *Nf = 4;
 97:     if (v0) *v0 = quad_v0;
 98:     if (J) *J = quad_J;
 99:     if (invJ) *invJ = quad_invJ;
100:     if (detJ) *detJ = quad_detJ;
101:     break;
102:   default:
103:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
104:   }
105:   return 0;
106: }

108: /*@
109:   DMPlexRefineRegularGetAffineTransforms - Gets the affine map from the reference cell to each subcell

111:   Input Parameters:
112: + cr - The DMPlexCellRefiner object
113: - ct - The cell type

115:   Output Parameters:
116: + Nc   - The number of subcells produced from this cell type
117: . v0   - The translation of the first vertex for each subcell
118: . J    - The Jacobian for each subcell (map from reference cell to subcell)
119: - invJ - The inverse Jacobian for each subcell

121:   Level: developer

123: .seealso: `DMPlexRefineRegularGetAffineFaceTransforms()`, `DMPLEXREFINEREGULAR`
124: @*/
125: PetscErrorCode DMPlexRefineRegularGetAffineTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
126: {
127:   /*
128:    2
129:    |\
130:    | \
131:    |  \
132:    |   \
133:    | C  \
134:    |     \
135:    |      \
136:    2---1---1
137:    |\  D  / \
138:    | 2   0   \
139:    |A \ /  B  \
140:    0---0-------1
141:    */
142:   static PetscReal tri_v0[]   = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
143:   static PetscReal tri_J[]    = {0.5, 0.0,  0.0, 0.5,

145:                                  0.5, 0.0,  0.0, 0.5,

147:                                  0.5, 0.0,  0.0, 0.5,

149:                                  0.0, -0.5, 0.5, 0.5};
150:   static PetscReal tri_invJ[] = {2.0, 0.0, 0.0,  2.0,

152:                                  2.0, 0.0, 0.0,  2.0,

154:                                  2.0, 0.0, 0.0,  2.0,

156:                                  2.0, 2.0, -2.0, 0.0};
157:   /*
158:      3---------2---------2
159:      |         |         |
160:      |    D    2    C    |
161:      |         |         |
162:      3----3----0----1----1
163:      |         |         |
164:      |    A    0    B    |
165:      |         |         |
166:      0---------0---------1
167:      */
168:   static PetscReal quad_v0[]   = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
169:   static PetscReal quad_J[]    = {0.5, 0.0, 0.0, 0.5,

171:                                   0.5, 0.0, 0.0, 0.5,

173:                                   0.5, 0.0, 0.0, 0.5,

175:                                   0.5, 0.0, 0.0, 0.5};
176:   static PetscReal quad_invJ[] = {2.0, 0.0, 0.0, 2.0,

178:                                   2.0, 0.0, 0.0, 2.0,

180:                                   2.0, 0.0, 0.0, 2.0,

182:                                   2.0, 0.0, 0.0, 2.0};
183:   /*
184:      Bottom (viewed from top)    Top
185:      1---------2---------2       7---------2---------6
186:      |         |         |       |         |         |
187:      |    B    2    C    |       |    H    2    G    |
188:      |         |         |       |         |         |
189:      3----3----0----1----1       3----3----0----1----1
190:      |         |         |       |         |         |
191:      |    A    0    D    |       |    E    0    F    |
192:      |         |         |       |         |         |
193:      0---------0---------3       4---------0---------5
194:      */
195:   static PetscReal hex_v0[]   = {-1.0, -1.0, -1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
196:   static PetscReal hex_J[]    = {0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

198:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

200:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

202:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

204:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

206:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

208:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

210:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
211:   static PetscReal hex_invJ[] = {2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

213:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

215:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

217:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

219:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

221:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

223:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

225:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0};

227:   switch (ct) {
228:   case DM_POLYTOPE_TRIANGLE:
229:     if (Nc) *Nc = 4;
230:     if (v0) *v0 = tri_v0;
231:     if (J) *J = tri_J;
232:     if (invJ) *invJ = tri_invJ;
233:     break;
234:   case DM_POLYTOPE_QUADRILATERAL:
235:     if (Nc) *Nc = 4;
236:     if (v0) *v0 = quad_v0;
237:     if (J) *J = quad_J;
238:     if (invJ) *invJ = quad_invJ;
239:     break;
240:   case DM_POLYTOPE_HEXAHEDRON:
241:     if (Nc) *Nc = 8;
242:     if (v0) *v0 = hex_v0;
243:     if (J) *J = hex_J;
244:     if (invJ) *invJ = hex_invJ;
245:     break;
246:   default:
247:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
248:   }
249:   return 0;
250: }

252: #if 0
253: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
254: {
255:   static PetscReal seg_v[]  = {-1.0,  0.0,  1.0};
256:   static PetscReal tri_v[]  = {-1.0, -1.0,  1.0, -1.0,  -1.0, 1.0,  0.0, -1.0,  0.0, 0.0,  -1.0, 0.0};
257:   static PetscReal quad_v[] = {-1.0, -1.0,  1.0, -1.0,   1.0, 1.0,  -1.0, 1.0,  0.0, -1.0,  1.0, 0.0,   0.0, 1.0,  -1.0, 0.0,  0.0, 0.0};
258:   static PetscReal tet_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
259:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,  -1.0,  1.0, -1.0,
260:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,  -1.0,  0.0,  0.0,  -1.0, -1.0,  1.0};
261:   static PetscReal hex_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
262:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,   1.0,  0.0, -1.0,
263:                                -1.0,  1.0, -1.0,   0.0,  1.0, -1.0,   1.0,  1.0, -1.0,
264:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,   1.0, -1.0,  0.0,
265:                                -1.0,  0.0,  0.0,   0.0,  0.0,  0.0,   1.0,  0.0,  0.0,
266:                                -1.0,  1.0,  0.0,   0.0,  1.0,  0.0,   1.0,  1.0,  0.0,
267:                                -1.0, -1.0,  1.0,   0.0, -1.0,  1.0,   1.0, -1.0,  1.0,
268:                                -1.0,  0.0,  1.0,   0.0,  0.0,  1.0,   1.0,  0.0,  1.0,
269:                                -1.0,  1.0,  1.0,   0.0,  1.0,  1.0,   1.0,  1.0,  1.0};

271:   switch (ct) {
272:     case DM_POLYTOPE_SEGMENT:       *Nv =  3; *subcellV = seg_v;  break;
273:     case DM_POLYTOPE_TRIANGLE:      *Nv =  6; *subcellV = tri_v;  break;
274:     case DM_POLYTOPE_QUADRILATERAL: *Nv =  9; *subcellV = quad_v; break;
275:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 10; *subcellV = tet_v;  break;
276:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 27; *subcellV = hex_v;  break;
277:     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
278:   }
279:   return 0;
280: }

282: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
283: {
284:   static PetscInt seg_v[]  = {0, 1, 1, 2};
285:   static PetscInt tri_v[]  = {0, 3, 5,  3, 1, 4,  5, 4, 2,  3, 4, 5};
286:   static PetscInt quad_v[] = {0, 4, 8, 7,  4, 1, 5, 8,  8, 5, 2, 6,  7, 8, 6, 3};
287:   static PetscInt tet_v[]  = {0, 3, 1, 6,  3, 2, 4, 8,  1, 4, 5, 7,  6, 8, 7, 9,
288:                               1, 6, 3, 7,  8, 4, 3, 7,  7, 3, 1, 4,  7, 3, 8, 6};
289:   static PetscInt hex_v[]  = {0,  3,  4,  1,  9, 10, 13, 12,   3,  6,  7,  4, 12, 13, 16, 15,   4,  7,  8,  5, 13, 14, 17, 16,   1,  4 , 5 , 2, 10, 11, 14, 13,
290:                               9, 12, 13, 10, 18, 19, 22, 21,  10, 13, 14, 11, 19, 20, 23, 22,  13, 16, 17, 14, 22, 23, 26, 25,  12, 15, 16, 13, 21, 22, 25, 24};

293:   switch (ct) {
294:     case DM_POLYTOPE_SEGMENT:       *Nv = 2; *subcellV = &seg_v[r*(*Nv)];  break;
295:     case DM_POLYTOPE_TRIANGLE:      *Nv = 3; *subcellV = &tri_v[r*(*Nv)];  break;
296:     case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
297:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 4; *subcellV = &tet_v[r*(*Nv)];  break;
298:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 8; *subcellV = &hex_v[r*(*Nv)];  break;
299:     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
300:   }
301:   return 0;
302: }
303: #endif

305: static PetscErrorCode DMPlexTransformView_Regular(DMPlexTransform tr, PetscViewer viewer)
306: {
307:   PetscBool isascii;

311:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
312:   if (isascii) {
313:     const char *name;

315:     PetscObjectGetName((PetscObject)tr, &name);
316:     PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : "");
317:   } else {
318:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
319:   }
320:   return 0;
321: }

323: static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr)
324: {
325:   return 0;
326: }

328: static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr)
329: {
330:   DMPlexRefine_Regular *f = (DMPlexRefine_Regular *)tr->data;

332:   PetscFree(f);
333:   return 0;
334: }

336: PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
337: {
338:   static PetscInt seg_seg[]   = {1, -1, 0, -1, 0, 0, 1, 0};
339:   static PetscInt tri_seg[]   = {2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
340:   static PetscInt tri_tri[]   = {1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2};
341:   static PetscInt quad_seg[]  = {1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 1, 0, 2, 0};
342:   static PetscInt quad_quad[] = {2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0, -3, 3, -3, 2, -3, 0, -2, 3, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 3, 1, 0, 1, 2, 2, 3, 2, 0, 2, 1, 2, 3, 3, 0, 3, 1, 3, 2, 3};
343:   static PetscInt tseg_seg[]  = {0, -1, 0, 0, 0, 0, 0, -1};
344:   static PetscInt tseg_tseg[] = {1, -2, 0, -2, 1, -1, 0, -1, 0, 0, 1, 0, 0, 1, 1, 1};
345:   static PetscInt tet_seg[]   = {0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0};
346:   static PetscInt tet_tri[] = {2, -1, 3, -1, 1, -3, 0, -2, 6, 1, 7, -3, 5, 2, 4, -3, 3, -1, 1, -1, 2, -3, 0, -1, 4, 0,  5, 0,  6, 0,  7, 0,  1, -1, 2, -1, 3, -3, 0, -3, 4, 0,  5, 0,  6, 0,  7, 0, 3, -2, 2, -3, 0, -1, 1, -1, 7, -3, 6, 1, 4, 2,  5, -3,
347:                                2, -3, 0, -2, 3, -1, 1, -3, 4, 0, 5, 0,  6, 0, 7, 0,  0, -2, 3, -2, 2, -2, 1, -2, 4, 0,  5, 0,  6, 0,  7, 0,  0, -1, 1, -2, 3, -2, 2, -2, 7, 1,  6, -3, 5, -3, 4, 2, 1, -2, 3, -3, 0, -3, 2, -1, 4, 0,  5, 0, 6, 0,  7, 0,
348:                                3, -3, 0, -1, 1, -1, 2, -3, 4, 0, 5, 0,  6, 0, 7, 0,  1, -3, 0, -3, 2, -1, 3, -3, 6, -3, 7, 1,  4, -3, 5, 2,  0, -3, 2, -2, 1, -2, 3, -2, 4, 0,  5, 0,  6, 0,  7, 0, 2, -2, 1, -3, 0, -2, 3, -1, 4, 0,  5, 0, 6, 0,  7, 0,
349:                                0, 0,  1, 0,  2, 0,  3, 0,  4, 0, 5, 0,  6, 0, 7, 0,  1, 0,  2, 2,  0, 1,  3, 1,  4, 0,  5, 0,  6, 0,  7, 0,  2, 2,  0, 0,  1, 1,  3, 2,  4, 0,  5, 0,  6, 0,  7, 0, 1, 2,  0, 1,  3, 1,  2, 2,  5, 0,  4, 0, 7, -1, 6, -1,
350:                                0, 1,  3, 0,  1, 0,  2, 0,  4, 0, 5, 0,  6, 0, 7, 0,  3, 0,  1, 2,  0, 2,  2, 1,  4, 0,  5, 0,  6, 0,  7, 0,  2, 0,  3, 2,  0, 0,  1, 1,  4, -2, 5, -2, 7, 0,  6, 0, 3, 2,  0, 2,  2, 1,  1, 2,  4, 0,  5, 0, 6, 0,  7, 0,
351:                                0, 2,  2, 0,  3, 0,  1, 0,  4, 0, 5, 0,  6, 0, 7, 0,  3, 1,  2, 1,  1, 2,  0, 2,  5, -2, 4, -2, 6, -1, 7, -1, 2, 1,  1, 1,  3, 2,  0, 0,  4, 0,  5, 0,  6, 0,  7, 0, 1, 1,  3, 1,  2, 2,  0, 1,  4, 0,  5, 0, 6, 0,  7, 0};
352:   static PetscInt tet_tet[] = {2, -12, 3, -12, 1, -12, 0, -12, 6, -9, 7, -9, 5, -12, 4, -12, 3, -11, 1, -11, 2, -11, 0, -11, 4, 0,  5, 0,  6, 0,  7, 0,  1, -10, 2, -10, 3, -10, 0, -10, 4, 0,  5, 0,  6, 0,  7, 0,  3, -9, 2, -9, 0, -9, 1, -9,
353:                                7, -9,  6, -9,  4, -12, 5, -12, 2, -8, 0, -8, 3, -8,  1, -8,  4, 0,   5, 0,   6, 0,   7, 0,   0, -7, 3, -7, 2, -7, 1, -7, 4, 0,   5, 0,   6, 0,   7, 0,   0, -6, 1, -6, 3, -6, 2, -6, 4, -3, 5, -3, 7, -6, 6, -6,
354:                                1, -5,  3, -5,  0, -5,  2, -5,  4, 0,  5, 0,  6, 0,   7, 0,   3, -4,  0, -4,  1, -4,  2, -4,  4, 0,  5, 0,  6, 0,  7, 0,  1, -3,  0, -3,  2, -3,  3, -3,  5, -3, 4, -3, 6, -6, 7, -6, 0, -2, 2, -2, 1, -2, 3, -2,
355:                                4, 0,   5, 0,   6, 0,   7, 0,   2, -1, 1, -1, 0, -1,  3, -1,  4, 0,   5, 0,   6, 0,   7, 0,   0, 0,  1, 0,  2, 0,  3, 0,  4, 0,   5, 0,   6, 0,   7, 0,   1, 1,  2, 1,  0, 1,  3, 1,  4, 0,  5, 0,  6, 0,  7, 0,
356:                                2, 2,   0, 2,   1, 2,   3, 2,   4, 0,  5, 0,  6, 0,   7, 0,   1, 3,   0, 3,   3, 3,   2, 3,   5, 0,  4, 0,  7, 0,  6, 0,  0, 4,   3, 4,   1, 4,   2, 4,   4, 0,  5, 0,  6, 0,  7, 0,  3, 5,  1, 5,  0, 5,  2, 5,
357:                                4, 0,   5, 0,   6, 0,   7, 0,   2, 6,  3, 6,  0, 6,   1, 6,   6, 6,   7, 6,   4, 6,   5, 6,   3, 7,  0, 7,  2, 7,  1, 7,  4, 0,   5, 0,   6, 0,   7, 0,   0, 8,  2, 8,  3, 8,  1, 8,  4, 0,  5, 0,  6, 0,  7, 0,
358:                                3, 9,   2, 9,   1, 9,   0, 9,   7, 6,  6, 6,  5, 6,   4, 6,   2, 10,  1, 10,  3, 10,  0, 10,  4, 0,  5, 0,  6, 0,  7, 0,  1, 11,  3, 11,  2, 11,  0, 11,  4, 0,  5, 0,  6, 0,  7, 0};
359:   static PetscInt hex_seg[] = {2, 0, 3, 0, 4, 0, 5, 0, 1, 0, 0, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 2, 0, 3, 0,
360:                                2, 0, 3, 0, 5, 0, 4, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 4, 0, 5, 0,
361:                                1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 2, 0, 3, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0,
362:                                1, 0, 0, 0, 3, 0, 2, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 5, 0, 4, 0,
363:                                0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 3, 0, 2, 0, 2, 0, 3, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 4, 0, 5, 0,
364:                                3, 0, 2, 0, 0, 0, 1, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 5, 0, 4, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0,
365:                                2, 0, 3, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 2, 0, 3, 0, 5, 0, 4, 0, 1, 0, 0, 0,
366:                                4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 1, 0, 0, 0, 5, 0, 4, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 2, 0, 3, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0};
367:   static PetscInt hex_quad[]    = {7,  -2, 4,  -2, 5,  -2, 6,  -2, 8,  -3, 11, -3, 10, -3, 9,  -3, 3,  1,  2,  1,  1,  1,  0,  1,  8,  -2, 9,  -2, 10, -2, 11, -2, 3,  -4, 0,  -4, 1,  -4, 2,  -4, 7,  0,  4,  0,  5,  0,  6,  0,  9,  1,  8,  1,  11,
368:                                    1,  10, 1,  0,  3,  3,  3,  2,  3,  1,  3,  5,  2,  6,  2,  7,  2,  4,  2,  6,  3,  5,  3,  4,  3,  7,  3,  10, -1, 9,  -1, 8,  -1, 11, -1, 2,  -4, 3,  -4, 0,  -4, 1,  -4, 4,  1,  7,  1,  6,  1,  5,  1,  11, 2,
369:                                    8,  2,  9,  2,  10, 2,  1,  3,  0,  3,  3,  3,  2,  3,  10, -4, 11, -4, 8,  -4, 9,  -4, 2,  1,  1,  1,  0,  1,  3,  1,  6,  -1, 5,  -1, 4,  -1, 7,  -1, 5,  -4, 6,  -4, 7,  -4, 4,  -4, 9,  0,  10, 0,  11, 0,  8,
370:                                    0,  0,  -2, 1,  -2, 2,  -2, 3,  -2, 11, 3,  10, 3,  9,  3,  8,  3,  1,  -2, 2,  -2, 3,  -2, 0,  -2, 4,  -3, 7,  -3, 6,  -3, 5,  -3, 11, -1, 8,  -1, 9,  -1, 10, -1, 7,  3,  4,  3,  5,  3,  6,  3,  2,  2,  1,  2,
371:                                    0,  2,  3,  2,  10, 2,  9,  2,  8,  2,  11, 2,  5,  1,  6,  1,  7,  1,  4,  1,  1,  -1, 2,  -1, 3,  -1, 0,  -1, 5,  2,  4,  2,  7,  2,  6,  2,  1,  2,  0,  2,  3,  2,  2,  2,  10, -4, 9,  -4, 8,  -4, 11, -4, 4,
372:                                    -3, 5,  -3, 6,  -3, 7,  -3, 0,  -3, 1,  -3, 2,  -3, 3,  -3, 8,  -2, 11, -2, 10, -2, 9,  -2, 3,  1,  0,  1,  1,  1,  2,  1,  9,  -4, 8,  -4, 11, -4, 10, -4, 6,  3,  7,  3,  4,  3,  5,  3,  1,  3,  2,  3,  3,  3,
373:                                    0,  3,  10, 1,  11, 1,  8,  1,  9,  1,  5,  -4, 4,  -4, 7,  -4, 6,  -4, 8,  0,  11, 0,  10, 0,  9,  0,  4,  -4, 7,  -4, 6,  -4, 5,  -4, 0,  0,  3,  0,  2,  0,  1,  0,  0,  0,  1,  0,  2,  0,  3,  0,  5,  -1, 4,
374:                                    -1, 7,  -1, 6,  -1, 9,  -3, 8,  -3, 11, -3, 10, -3, 9,  -3, 10, -3, 11, -3, 8,  -3, 6,  -2, 5,  -2, 4,  -2, 7,  -2, 3,  -3, 0,  -3, 1,  -3, 2,  -3, 7,  0,  6,  0,  5,  0,  4,  0,  2,  -1, 3,  -1, 0,  -1, 1,  -1,
375:                                    11, 3,  8,  3,  9,  3,  10, 3,  2,  2,  3,  2,  0,  2,  1,  2,  6,  2,  7,  2,  4,  2,  5,  2,  10, 2,  11, 2,  8,  2,  9,  2,  6,  -1, 7,  -1, 4,  -1, 5,  -1, 3,  0,  2,  0,  1,  0,  0,  0,  9,  1,  10, 1,  11,
376:                                    1,  8,  1,  2,  -4, 1,  -4, 0,  -4, 3,  -4, 11, -2, 10, -2, 9,  -2, 8,  -2, 7,  -2, 6,  -2, 5,  -2, 4,  -2, 1,  -1, 0,  -1, 3,  -1, 2,  -1, 4,  0,  5,  0,  6,  0,  7,  0,  11, -1, 10, -1, 9,  -1, 8,  -1, 0,  -2,
377:                                    3,  -2, 2,  -2, 1,  -2, 8,  3,  9,  3,  10, 3,  11, 3,  4,  1,  5,  1,  6,  1,  7,  1,  3,  -3, 2,  -3, 1,  -3, 0,  -3, 7,  -3, 6,  -3, 5,  -3, 4,  -3, 8,  0,  9,  0,  10, 0,  11, 0,  0,  0,  1,  0,  2,  0,  3,
378:                                    0,  4,  0,  5,  0,  6,  0,  7,  0,  8,  0,  9,  0,  10, 0,  11, 0,  1,  3,  2,  3,  3,  3,  0,  3,  11, -2, 10, -2, 9,  -2, 8,  -2, 4,  1,  5,  1,  6,  1,  7,  1,  2,  2,  3,  2,  0,  2,  1,  2,  7,  -3, 6,  -3,
379:                                    5,  -3, 4,  -3, 11, -1, 10, -1, 9,  -1, 8,  -1, 3,  1,  0,  1,  1,  1,  2,  1,  8,  3,  9,  3,  10, 3,  11, 3,  7,  -2, 6,  -2, 5,  -2, 4,  -2, 5,  2,  4,  2,  7,  2,  6,  2,  0,  -3, 1,  -3, 2,  -3, 3,  -3, 9,
380:                                    1,  10, 1,  11, 1,  8,  1,  1,  -1, 0,  -1, 3,  -1, 2,  -1, 5,  -1, 4,  -1, 7,  -1, 6,  -1, 10, 2,  11, 2,  8,  2,  9,  2,  4,  -3, 5,  -3, 6,  -3, 7,  -3, 1,  2,  0,  2,  3,  2,  2,  2,  11, 3,  8,  3,  9,  3,
381:                                    10, 3,  8,  0,  11, 0,  10, 0,  9,  0,  7,  3,  4,  3,  5,  3,  6,  3,  3,  -3, 0,  -3, 1,  -3, 2,  -3, 3,  -3, 2,  -3, 1,  -3, 0,  -3, 6,  2,  7,  2,  4,  2,  5,  2,  9,  -3, 8,  -3, 11, -3, 10, -3, 9,  -3, 10,
382:                                    -3, 11, -3, 8,  -3, 5,  1,  6,  1,  7,  1,  4,  1,  0,  0,  3,  0,  2,  0,  1,  0,  0,  -2, 3,  -2, 2,  -2, 1,  -2, 9,  -4, 8,  -4, 11, -4, 10, -4, 5,  -4, 4,  -4, 7,  -4, 6,  -4, 2,  -4, 1,  -4, 0,  -4, 3,  -4,
383:                                    10, 1,  11, 1,  8,  1,  9,  1,  6,  3,  7,  3,  4,  3,  5,  3,  7,  0,  6,  0,  5,  0,  4,  0,  3,  0,  2,  0,  1,  0,  0,  0,  8,  -2, 11, -2, 10, -2, 9,  -2, 6,  -1, 7,  -1, 4,  -1, 5,  -1, 2,  -1, 3,  -1, 0,
384:                                    -1, 1,  -1, 10, -4, 9,  -4, 8,  -4, 11, -4, 11, -1, 8,  -1, 9,  -1, 10, -1, 4,  -4, 7,  -4, 6,  -4, 5,  -4, 1,  -1, 2,  -1, 3,  -1, 0,  -1, 10, 2,  9,  2,  8,  2,  11, 2,  6,  -2, 5,  -2, 4,  -2, 7,  -2, 2,  2,
385:                                    1,  2,  0,  2,  3,  2,  8,  -2, 9,  -2, 10, -2, 11, -2, 0,  3,  3,  3,  2,  3,  1,  3,  4,  -3, 7,  -3, 6,  -3, 5,  -3, 4,  1,  7,  1,  6,  1,  5,  1,  8,  -3, 11, -3, 10, -3, 9,  -3, 0,  -2, 1,  -2, 2,  -2, 3,
386:                                    -2, 9,  1,  8,  1,  11, 1,  10, 1,  3,  -4, 0,  -4, 1,  -4, 2,  -4, 6,  -1, 5,  -1, 4,  -1, 7,  -1, 5,  -4, 6,  -4, 7,  -4, 4,  -4, 10, -1, 9,  -1, 8,  -1, 11, -1, 1,  3,  0,  3,  3,  3,  2,  3,  7,  -2, 4,  -2,
387:                                    5,  -2, 6,  -2, 11, 2,  8,  2,  9,  2,  10, 2,  2,  -4, 3,  -4, 0,  -4, 1,  -4, 10, -4, 11, -4, 8,  -4, 9,  -4, 1,  -2, 2,  -2, 3,  -2, 0,  -2, 5,  2,  6,  2,  7,  2,  4,  2,  11, 3,  10, 3,  9,  3,  8,  3,  2,
388:                                    1,  1,  1,  0,  1,  3,  1,  7,  0,  4,  0,  5,  0,  6,  0,  6,  3,  5,  3,  4,  3,  7,  3,  9,  0,  10, 0,  11, 0,  8,  0,  3,  1,  2,  1,  1,  1,  0,  1};
389:   static PetscInt hex_hex[]     = {3, -24, 0, -24, 4, -24, 5, -24, 2, -24, 6, -24, 7, -24, 1, -24, 3, -23, 5, -23, 6, -23, 2, -23, 0, -23, 1, -23, 7, -23, 4, -23, 4, -22, 0, -22, 1, -22, 7, -22, 5, -22, 6, -22, 2, -22, 3, -22, 6, -21, 7, -21,
390:                                    1, -21, 2, -21, 5, -21, 3, -21, 0, -21, 4, -21, 1, -20, 2, -20, 6, -20, 7, -20, 0, -20, 4, -20, 5, -20, 3, -20, 6, -19, 2, -19, 3, -19, 5, -19, 7, -19, 4, -19, 0, -19, 1, -19, 4, -18, 5, -18, 3, -18, 0, -18,
391:                                    7, -18, 1, -18, 2, -18, 6, -18, 1, -17, 7, -17, 4, -17, 0, -17, 2, -17, 3, -17, 5, -17, 6, -17, 2, -16, 3, -16, 5, -16, 6, -16, 1, -16, 7, -16, 4, -16, 0, -16, 7, -15, 4, -15, 0, -15, 1, -15, 6, -15, 2, -15,
392:                                    3, -15, 5, -15, 7, -14, 1, -14, 2, -14, 6, -14, 4, -14, 5, -14, 3, -14, 0, -14, 0, -13, 4, -13, 5, -13, 3, -13, 1, -13, 2, -13, 6, -13, 7, -13, 5, -12, 4, -12, 7, -12, 6, -12, 3, -12, 2, -12, 1, -12, 0, -12,
393:                                    7, -11, 6, -11, 5, -11, 4, -11, 1, -11, 0, -11, 3, -11, 2, -11, 0, -10, 1, -10, 7, -10, 4, -10, 3, -10, 5, -10, 6, -10, 2, -10, 4, -9,  7, -9,  6, -9,  5, -9,  0, -9,  3, -9,  2, -9,  1, -9,  5, -8,  6, -8,
394:                                    2, -8,  3, -8,  4, -8,  0, -8,  1, -8,  7, -8,  2, -7,  6, -7,  7, -7,  1, -7,  3, -7,  0, -7,  4, -7,  5, -7,  6, -6,  5, -6,  4, -6,  7, -6,  2, -6,  1, -6,  0, -6,  3, -6,  5, -5,  3, -5,  0, -5,  4, -5,
395:                                    6, -5,  7, -5,  1, -5,  2, -5,  2, -4,  1, -4,  0, -4,  3, -4,  6, -4,  5, -4,  4, -4,  7, -4,  1, -3,  0, -3,  3, -3,  2, -3,  7, -3,  6, -3,  5, -3,  4, -3,  0, -2,  3, -2,  2, -2,  1, -2,  4, -2,  7, -2,
396:                                    6, -2,  5, -2,  3, -1,  2, -1,  1, -1,  0, -1,  5, -1,  4, -1,  7, -1,  6, -1,  0, 0,   1, 0,   2, 0,   3, 0,   4, 0,   5, 0,   6, 0,   7, 0,   1, 1,   2, 1,   3, 1,   0, 1,   7, 1,   4, 1,   5, 1,   6, 1,
397:                                    2, 2,   3, 2,   0, 2,   1, 2,   6, 2,   7, 2,   4, 2,   5, 2,   3, 3,   0, 3,   1, 3,   2, 3,   5, 3,   6, 3,   7, 3,   4, 3,   4, 4,   0, 4,   3, 4,   5, 4,   7, 4,   6, 4,   2, 4,   1, 4,   7, 5,   4, 5,
398:                                    5, 5,   6, 5,   1, 5,   2, 5,   3, 5,   0, 5,   1, 6,   7, 6,   6, 6,   2, 6,   0, 6,   3, 6,   5, 6,   4, 6,   3, 7,   2, 7,   6, 7,   5, 7,   0, 7,   4, 7,   7, 7,   1, 7,   5, 8,   6, 8,   7, 8,   4, 8,
399:                                    3, 8,   0, 8,   1, 8,   2, 8,   4, 9,   7, 9,   1, 9,   0, 9,   5, 9,   3, 9,   2, 9,   6, 9,   4, 10,  5, 10,  6, 10,  7, 10,  0, 10,  1, 10,  2, 10,  3, 10,  6, 11,  7, 11,  4, 11,  5, 11,  2, 11,  3, 11,
400:                                    0, 11,  1, 11,  3, 12,  5, 12,  4, 12,  0, 12,  2, 12,  1, 12,  7, 12,  6, 12,  6, 13,  2, 13,  1, 13,  7, 13,  5, 13,  4, 13,  0, 13,  3, 13,  1, 14,  0, 14,  4, 14,  7, 14,  2, 14,  6, 14,  5, 14,  3, 14,
401:                                    6, 15,  5, 15,  3, 15,  2, 15,  7, 15,  1, 15,  0, 15,  4, 15,  0, 16,  4, 16,  7, 16,  1, 16,  3, 16,  2, 16,  6, 16,  5, 16,  0, 17,  3, 17,  5, 17,  4, 17,  1, 17,  7, 17,  6, 17,  2, 17,  5, 18,  3, 18,
402:                                    2, 18,  6, 18,  4, 18,  7, 18,  1, 18,  0, 18,  7, 19,  6, 19,  2, 19,  1, 19,  4, 19,  0, 19,  3, 19,  5, 19,  2, 20,  1, 20,  7, 20,  6, 20,  3, 20,  5, 20,  4, 20,  0, 20,  7, 21,  1, 21,  0, 21,  4, 21,
403:                                    6, 21,  5, 21,  3, 21,  2, 21,  2, 22,  6, 22,  5, 22,  3, 22,  1, 22,  0, 22,  4, 22,  7, 22,  5, 23,  4, 23,  0, 23,  3, 23,  6, 23,  2, 23,  1, 23,  7, 23};
404:   static PetscInt trip_seg[]    = {1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 0, -1, 2, -1, 1, -1, 1, -1, 0, -1, 2, -1, 2, -1, 1, -1, 0, -1,
405:                                    0, 0, 1, 0, 2, 0, 2, 0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1};
406:   static PetscInt trip_tri[]    = {1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2, 0, 0, 1, 0, 2, 0, 3, 0, 2, -1, 1, -1, 0, -1, 3, -3, 0, -2, 2, -2, 1, -2, 3, -1, 1, -3, 0, -3, 2, -3, 3, -2,
407:                                    0, 0, 1, 0, 2, 0, 3, 0, 2, 2, 0, 2, 1, 2, 3, 2, 1, 1, 2, 1, 0, 1, 3, 1, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3};
408:   static PetscInt trip_quad[]   = {4, -1, 5, -1, 3, -1, 1, -1, 2, -1, 0, -1, 5, -1, 3, -1, 4, -1, 2, -1, 0, -1, 1, -1, 3, -1, 4, -1, 5, -1, 0, -1, 1, -1, 2, -1, 0, -3, 2, -3, 1, -3, 3, -3, 5, -3, 4, -3,
409:                                    1, -3, 0, -3, 2, -3, 4, -3, 3, -3, 5, -3, 2, -3, 1, -3, 0, -3, 5, -3, 4, -3, 3, -3, 0, 0,  1, 0,  2, 0,  3, 0,  4, 0,  5, 0,  2, 0,  0, 0,  1, 0,  5, 0,  3, 0,  4, 0,
410:                                    1, 0,  2, 0,  0, 0,  4, 0,  5, 0,  3, 0,  5, 2,  4, 2,  3, 2,  2, 2,  1, 2,  0, 2,  4, 2,  3, 2,  5, 2,  1, 2,  0, 2,  2, 2,  3, 2,  5, 2,  4, 2,  0, 2,  2, 2,  1, 2};
411:   static PetscInt trip_trip[]   = {5, -6, 6, -6, 4, -6, 7, -6, 1, -6, 2, -6, 0, -6, 3, -6, 6, -5, 4, -5, 5, -5, 7, -5, 2, -5, 0, -5, 1, -5, 3, -5, 4, -4, 5, -4, 6, -4, 7, -4, 0, -4, 1, -4, 2, -4, 3, -4,
412:                                    2, -3, 1, -3, 0, -3, 3, -1, 6, -3, 5, -3, 4, -3, 7, -1, 0, -2, 2, -2, 1, -2, 3, -3, 4, -2, 6, -2, 5, -2, 7, -3, 1, -1, 0, -1, 2, -1, 3, -2, 5, -1, 4, -1, 6, -1, 7, -2,
413:                                    0, 0,  1, 0,  2, 0,  3, 0,  4, 0,  5, 0,  6, 0,  7, 0,  2, 1,  0, 1,  1, 1,  3, 1,  6, 1,  4, 1,  5, 1,  7, 1,  1, 2,  2, 2,  0, 2,  3, 2,  5, 2,  6, 2,  4, 2,  7, 2,
414:                                    5, 3,  4, 3,  6, 3,  7, 4,  1, 3,  0, 3,  2, 3,  3, 4,  4, 4,  6, 4,  5, 4,  7, 5,  0, 4,  2, 4,  1, 4,  3, 5,  6, 5,  5, 5,  4, 5,  7, 3,  2, 5,  1, 5,  0, 5,  3, 3};
415:   static PetscInt ttri_tseg[]   = {2, -2, 1, -2, 0, -2, 1, -2, 0, -2, 2, -2, 0, -2, 2, -2, 1, -2, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1,
416:                                    0, 0,  1, 0,  2, 0,  1, 0,  2, 0,  0, 0,  2, 0,  0, 0,  1, 0,  0, 1,  1, 1,  2, 1,  1, 1,  2, 1,  0, 1,  2, 1,  0, 1,  1, 1};
417:   static PetscInt ttri_ttri[]   = {1, -6, 0, -6, 2, -6, 3, -5, 0, -5, 2, -5, 1, -5, 3, -4, 2, -4, 1, -4, 0, -4, 3, -6, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3,
418:                                    0, 0,  1, 0,  2, 0,  3, 0,  1, 1,  2, 1,  0, 1,  3, 1,  2, 2,  0, 2,  1, 2,  3, 2,  0, 3,  1, 3,  2, 3,  3, 3,  1, 4,  2, 4,  0, 4,  3, 4,  2, 5,  0, 5,  1, 5,  3, 5};
419:   static PetscInt tquad_tvert[] = {0, -1, 0, -1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, -1, 0, -1};
420:   static PetscInt tquad_tseg[]  = {1, 1, 0, 1, 3, 1, 2, 1, 0, 1, 3, 1, 2, 1, 1, 1, 3, 1, 2, 1, 1, 1, 0, 1, 2, 1, 1, 1, 0, 1, 3, 1, 1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0,
421:                                    0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 1, 0, 2, 0, 0, 1, 1, 1, 2, 1, 3, 1, 1, 1, 2, 1, 3, 1, 0, 1, 2, 1, 3, 1, 0, 1, 1, 1, 3, 1, 0, 1, 1, 1, 2, 1};
422:   static PetscInt tquad_tquad[] = {2,  -8, 1,  -8, 0,  -8, 3,  -8, 1,  -7, 0,  -7, 3,  -7, 2,  -7, 0,  -6, 3,  -6, 2,  -6, 1, -6, 3, -5, 2, -5, 1, -5, 0, -5, 2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0,
423:                                    -3, 3,  -3, 2,  -3, 0,  -2, 3,  -2, 2,  -2, 1,  -2, 3,  -1, 2,  -1, 1,  -1, 0,  -1, 0,  0, 1,  0, 2,  0, 3,  0, 1,  1, 2,  1, 3,  1, 0,  1, 2,  2, 3,  2, 0,  2,
424:                                    1,  2,  3,  3,  0,  3,  1,  3,  2,  3,  0,  4,  1,  4,  2,  4,  3,  4,  1,  5,  2,  5,  3, 5,  0, 5,  2, 6,  3, 6,  0, 6,  1, 6,  3, 7,  0, 7,  1, 7,  2, 7};

427:   *rnew = r;
428:   *onew = o;
429:   if (!so) return 0;
430:   switch (sct) {
431:   case DM_POLYTOPE_POINT:
432:     break;
433:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
434:     *onew = so < 0 ? -(o + 1) : o;
435:     break;
436:   case DM_POLYTOPE_SEGMENT:
437:     switch (tct) {
438:     case DM_POLYTOPE_POINT:
439:       break;
440:     case DM_POLYTOPE_SEGMENT:
441:       *rnew = seg_seg[(so + 1) * 4 + r * 2];
442:       *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so + 1) * 4 + r * 2 + 1]);
443:       break;
444:     default:
445:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
446:     }
447:     break;
448:   case DM_POLYTOPE_TRIANGLE:
449:     switch (tct) {
450:     case DM_POLYTOPE_SEGMENT:
451:       *rnew = tri_seg[(so + 3) * 6 + r * 2];
452:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
453:       break;
454:     case DM_POLYTOPE_TRIANGLE:
455:       *rnew = tri_tri[(so + 3) * 8 + r * 2];
456:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 8 + r * 2 + 1]);
457:       break;
458:     default:
459:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
460:     }
461:     break;
462:   case DM_POLYTOPE_QUADRILATERAL:
463:     switch (tct) {
464:     case DM_POLYTOPE_POINT:
465:       break;
466:     case DM_POLYTOPE_SEGMENT:
467:       *rnew = quad_seg[(so + 4) * 8 + r * 2];
468:       *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]);
469:       break;
470:     case DM_POLYTOPE_QUADRILATERAL:
471:       *rnew = quad_quad[(so + 4) * 8 + r * 2];
472:       *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so + 4) * 8 + r * 2 + 1]);
473:       break;
474:     default:
475:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
476:     }
477:     break;
478:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
479:     switch (tct) {
480:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
481:       *rnew = tseg_seg[(so + 2) * 2 + r * 2];
482:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so + 2) * 2 + r * 2 + 1]);
483:       break;
484:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
485:       *rnew = tseg_tseg[(so + 2) * 4 + r * 2];
486:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so + 2) * 4 + r * 2 + 1]);
487:       break;
488:     default:
489:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
490:     }
491:     break;
492:   case DM_POLYTOPE_TETRAHEDRON:
493:     switch (tct) {
494:     case DM_POLYTOPE_SEGMENT:
495:       *rnew = tet_seg[(so + 12) * 2 + r * 2];
496:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 2 + r * 2 + 1]);
497:       break;
498:     case DM_POLYTOPE_TRIANGLE:
499:       *rnew = tet_tri[(so + 12) * 16 + r * 2];
500:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 16 + r * 2 + 1]);
501:       break;
502:     case DM_POLYTOPE_TETRAHEDRON:
503:       *rnew = tet_tet[(so + 12) * 16 + r * 2];
504:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 16 + r * 2 + 1]);
505:       break;
506:     default:
507:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
508:     }
509:     break;
510:   case DM_POLYTOPE_HEXAHEDRON:
511:     switch (tct) {
512:     case DM_POLYTOPE_POINT:
513:       break;
514:     case DM_POLYTOPE_SEGMENT:
515:       *rnew = hex_seg[(so + 24) * 12 + r * 2];
516:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so + 24) * 12 + r * 2 + 1]);
517:       break;
518:     case DM_POLYTOPE_QUADRILATERAL:
519:       *rnew = hex_quad[(so + 24) * 24 + r * 2];
520:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so + 24) * 24 + r * 2 + 1]);
521:       break;
522:     case DM_POLYTOPE_HEXAHEDRON:
523:       *rnew = hex_hex[(so + 24) * 16 + r * 2];
524:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so + 24) * 16 + r * 2 + 1]);
525:       break;
526:     default:
527:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
528:     }
529:     break;
530:   case DM_POLYTOPE_TRI_PRISM:
531:     switch (tct) {
532:     case DM_POLYTOPE_SEGMENT:
533:       *rnew = trip_seg[(so + 6) * 6 + r * 2];
534:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so + 6) * 6 + r * 2 + 1]);
535:       break;
536:     case DM_POLYTOPE_TRIANGLE:
537:       *rnew = trip_tri[(so + 6) * 8 + r * 2];
538:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so + 6) * 8 + r * 2 + 1]);
539:       break;
540:     case DM_POLYTOPE_QUADRILATERAL:
541:       *rnew = trip_quad[(so + 6) * 12 + r * 2];
542:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so + 6) * 12 + r * 2 + 1]);
543:       break;
544:     case DM_POLYTOPE_TRI_PRISM:
545:       *rnew = trip_trip[(so + 6) * 16 + r * 2];
546:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so + 6) * 16 + r * 2 + 1]);
547:       break;
548:     default:
549:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
550:     }
551:     break;
552:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
553:     switch (tct) {
554:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
555:       *rnew = ttri_tseg[(so + 6) * 6 + r * 2];
556:       *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so + 6) * 6 + r * 2 + 1]);
557:       break;
558:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
559:       *rnew = ttri_ttri[(so + 6) * 8 + r * 2];
560:       *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so + 6) * 8 + r * 2 + 1]);
561:       break;
562:     default:
563:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
564:     }
565:     break;
566:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
567:     switch (tct) {
568:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
569:       *rnew = tquad_tvert[(so + 8) * 2 + r * 2];
570:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so + 8) * 2 + r * 2 + 1]);
571:       break;
572:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
573:       *rnew = tquad_tseg[(so + 8) * 8 + r * 2];
574:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so + 8) * 8 + r * 2 + 1]);
575:       break;
576:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
577:       *rnew = tquad_tquad[(so + 8) * 8 + r * 2];
578:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so + 8) * 8 + r * 2 + 1]);
579:       break;
580:     default:
581:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
582:     }
583:     break;
584:   default:
585:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
586:   }
587:   return 0;
588: }

590: PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
591: {
592:   /* All vertices remain in the refined mesh */
593:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
594:   static PetscInt       vertexS[] = {1};
595:   static PetscInt       vertexC[] = {0};
596:   static PetscInt       vertexO[] = {0};
597:   /* Split all edges with a new vertex, making two new 2 edges
598:      0--0--0--1--1
599:   */
600:   static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
601:   static PetscInt       segS[] = {1, 2};
602:   static PetscInt       segC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
603:   static PetscInt       segO[] = {0, 0, 0, 0};
604:   /* Do not split tensor edges */
605:   static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
606:   static PetscInt       tvertS[] = {1};
607:   static PetscInt       tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
608:   static PetscInt       tvertO[] = {0, 0};
609:   /* Add 3 edges inside every triangle, making 4 new triangles.
610:    2
611:    |\
612:    | \
613:    |  \
614:    0   1
615:    | C  \
616:    |     \
617:    |      \
618:    2---1---1
619:    |\  D  / \
620:    1 2   0   0
621:    |A \ /  B  \
622:    0-0-0---1---1
623:   */
624:   static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
625:   static PetscInt       triS[] = {3, 4};
626:   static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
627:   static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0};
628:   /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
629:      3----1----2----0----2
630:      |         |         |
631:      0    D    2    C    1
632:      |         |         |
633:      3----3----0----1----1
634:      |         |         |
635:      1    A    0    B    0
636:      |         |         |
637:      0----0----0----1----1
638:   */
639:   static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
640:   static PetscInt       quadS[] = {1, 4, 4};
641:   static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
642:   static PetscInt quadO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, -1, 0, 0};
643:   /* Add 1 edge inside every tensor quad, making 2 new tensor quads
644:      2----2----1----3----3
645:      |         |         |
646:      |         |         |
647:      |         |         |
648:      4    A    6    B    5
649:      |         |         |
650:      |         |         |
651:      |         |         |
652:      0----0----0----1----1
653:   */
654:   static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
655:   static PetscInt       tsegS[] = {1, 2};
656:   static PetscInt tsegC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
657:   static PetscInt tsegO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
658:   /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
659:      The vertices of our reference tet are [(-1, -1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, 1)], which we call [v0, v1, v2, v3]. The first
660:      three edges are [v0, v1], [v1, v2], [v2, v0] called e0, e1, and e2, and then three edges to the top point [v0, v3], [v1, v3], [v2, v3]
661:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
662:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
663:      The first four tets just cut off the corners, using the replica notation for new vertices,
664:        [v0,      (e0, 0), (e2, 0), (e3, 0)]
665:        [(e0, 0), v1,      (e1, 0), (e4, 0)]
666:        [(e2, 0), (e1, 0), v2,      (e5, 0)]
667:        [(e3, 0), (e4, 0), (e5, 0), v3     ]
668:      The next four tets match a vertex to the newly created faces from cutting off those first tets.
669:        [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
670:        [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
671:        [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
672:        [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
673:      We can see that a new edge is introduced in the cell [(e0, 0), (e5, 0)] which we call (-1, 0). The first four faces created are
674:        [(e2, 0), (e0, 0), (e3, 0)]
675:        [(e0, 0), (e1, 0), (e4, 0)]
676:        [(e2, 0), (e5, 0), (e1, 0)]
677:        [(e3, 0), (e4, 0), (e5, 0)]
678:      The next four, from the second group of tets, are
679:        [(e2, 0), (e0, 0), (e5, 0)]
680:        [(e4, 0), (e0, 0), (e5, 0)]
681:        [(e0, 0), (e1, 0), (e5, 0)]
682:        [(e5, 0), (e3, 0), (e0, 0)]
683:      I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
684:    */
685:   static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
686:   static PetscInt       tetS[] = {1, 8, 8};
687:   static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
688:   static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -1, -1, 1, 0, 0, -1, -1, -3, 2, -1, 0, -1, 1};
689:   /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
690:      The vertices of our reference hex are (-1, -1, -1), (-1, 1, -1), (1, 1, -1), (1, -1, -1), (-1, -1, 1), (1, -1, 1), (1, 1, 1), (-1, 1, 1) which we call [v0, v1, v2, v3, v4, v5, v6, v7]. The fours edges around the bottom [v0, v1], [v1, v2], [v2, v3], [v3, v0] are [e0, e1, e2, e3], and likewise around the top [v4, v5], [v5, v6], [v6, v7], [v7, v4] are [e4, e5, e6, e7]. Finally [v0, v4], [v1, v7], [v2, v6], [v3, v5] are [e9, e10, e11, e8]. The faces of a hex, given in DMPlexGetRawFaces_Internal(), oriented with outward normals, are
691:        [v0, v1, v2, v3] f0 bottom
692:        [v4, v5, v6, v7] f1 top
693:        [v0, v3, v5, v4] f2 front
694:        [v2, v1, v7, v6] f3 back
695:        [v3, v2, v6, v5] f4 right
696:        [v0, v4, v7, v1] f5 left
697:      The eight hexes are divided into four on the bottom, and four on the top,
698:        [v0,      (e0, 0),  (f0, 0),  (e3, 0),  (e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
699:        [(e0, 0), v1,       (e1, 0),  (f0, 0),  (f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
700:        [(f0, 0), (e1, 0),  v2,       (e2, 0),  (c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
701:        [(e3, 0), (f0, 0),  (e2, 0),  v3,       (f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
702:        [(e9, 0), (f5, 0),  (c0, 0),  (f2, 0),  v4,      (e4, 0),  (f1, 0),  (e7, 0)]
703:        [(f2, 0), (c0, 0),  (f4, 0),  (e8, 0),  (e4, 0), v5,       (e5, 0),  (f1, 0)]
704:        [(c0, 0), (f3, 0),  (e11, 0), (f4, 0),  (f1, 0), (e5, 0),  v6,       (e6, 0)]
705:        [(f5, 0), (e10, 0), (f3, 0),  (c0, 0),  (e7, 0), (f1, 0),  (e6, 0),  v7]
706:      The 6 internal edges will go from the faces to the central vertex. The 12 internal faces can be divided into groups of 4 by the plane on which they sit. First the faces on the x-y plane are,
707:        [(e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
708:        [(f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
709:        [(c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
710:        [(f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
711:      and on the x-z plane,
712:        [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
713:        [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
714:        [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
715:        [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
716:      and on the y-z plane,
717:        [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
718:        [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
719:        [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
720:        [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
721:   */
722:   static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
723:   static PetscInt       hexS[] = {1, 6, 12, 8};
724:   static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 5, 3, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 1, 5, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 1, 5, 2};
725:   static PetscInt hexO[] = {0, 0, 0, 0,  0,  0, 0, 0, 0, 0, 0,  0, 0, 0, -1, -1, 0,  -1, -1, 0, -1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0,  -1, -1, -1, 0,  0, 0,  -1, -1, 0, 0, 0, -1, -1, 0,  0, -1, -1, -1, 0, 0,  -1, -1, -1,
726:                             0, 0, 0, -1, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, -3, 0,  -2, 0,  0,  0, -3, 0,  0, 0, 0,  0, 0, 0,  0,  0, -2, 0,  0,  0,  -2, 0, -2, 0,  0,  0, 0, 0, -2, 0,  -3, 0, 0,  0,  -2, 0, -3, 0,  -2, 0};
727:   /* Add 3 quads inside every triangular prism, making 4 new prisms. */
728:   static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
729:   static PetscInt       tripS[] = {3, 4, 6, 8};
730:   static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5};
731:   static PetscInt tripO[] = {0, 0, 0, 0, 0,  0, 0, -1, -1, -1, 0,  -1, -1, -1, 0, 0, 0, 0, -1, 0, -1, -1, -1, 0, -1, -1, -1, 0, -1, -1, 0,  -1, -1, 0,  0, -1, -1, 0, 0, -1, -1,
732:                              0, 0, 0, 0, -3, 0, 0, 0,  0,  0,  -3, 0,  0,  -3, 0, 0, 2, 0, 0,  0, 0,  -2, 0,  0, -3, 0,  -2, 0, 0,  0,  -3, -2, 0,  -3, 0, 0,  -2, 0, 0, 0,  0};
733:   /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
734:       2
735:       |\
736:       | \
737:       |  \
738:       0---1

740:       2

742:       0   1

744:       2
745:       |\
746:       | \
747:       |  \
748:       0---1
749:   */
750:   static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
751:   static PetscInt       ttriS[] = {3, 4};
752:   static PetscInt ttriC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2};
753:   static PetscInt ttriO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0};
754:   /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
755:   static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
756:   static PetscInt       tquadS[] = {1, 4, 4};
757:   static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 5, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
758:   static PetscInt tquadO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0};
759:   /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */
760:   static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID};
761:   static PetscInt       tpyrS[] = {4, 12, 1, 4, 6};
762:   static PetscInt       tpyrC[] = {
763:     DM_POLYTOPE_POINT,
764:     2,
765:     1,
766:     1,
767:     0,
768:     DM_POLYTOPE_POINT,
769:     1,
770:     0,
771:     0,
772:     DM_POLYTOPE_POINT,
773:     2,
774:     2,
775:     1,
776:     0,
777:     DM_POLYTOPE_POINT,
778:     1,
779:     0,
780:     0,
781:     DM_POLYTOPE_POINT,
782:     2,
783:     3,
784:     1,
785:     0,
786:     DM_POLYTOPE_POINT,
787:     1,
788:     0,
789:     0,
790:     DM_POLYTOPE_POINT,
791:     2,
792:     4,
793:     1,
794:     0,
795:     DM_POLYTOPE_POINT,
796:     1,
797:     0,
798:     0,
799:     /* These four triangle face out of the bottom pyramid */
800:     DM_POLYTOPE_SEGMENT,
801:     1,
802:     1,
803:     1,
804:     DM_POLYTOPE_SEGMENT,
805:     0,
806:     3,
807:     DM_POLYTOPE_SEGMENT,
808:     0,
809:     0,
810:     DM_POLYTOPE_SEGMENT,
811:     1,
812:     2,
813:     1,
814:     DM_POLYTOPE_SEGMENT,
815:     0,
816:     0,
817:     DM_POLYTOPE_SEGMENT,
818:     0,
819:     1,
820:     DM_POLYTOPE_SEGMENT,
821:     1,
822:     3,
823:     1,
824:     DM_POLYTOPE_SEGMENT,
825:     0,
826:     1,
827:     DM_POLYTOPE_SEGMENT,
828:     0,
829:     2,
830:     DM_POLYTOPE_SEGMENT,
831:     1,
832:     4,
833:     1,
834:     DM_POLYTOPE_SEGMENT,
835:     0,
836:     2,
837:     DM_POLYTOPE_SEGMENT,
838:     0,
839:     3,
840:     /* These eight triangles face out of the corner pyramids */
841:     DM_POLYTOPE_SEGMENT,
842:     1,
843:     0,
844:     3,
845:     DM_POLYTOPE_SEGMENT,
846:     0,
847:     3,
848:     DM_POLYTOPE_SEGMENT,
849:     1,
850:     1,
851:     2,
852:     DM_POLYTOPE_SEGMENT,
853:     1,
854:     0,
855:     2,
856:     DM_POLYTOPE_SEGMENT,
857:     0,
858:     0,
859:     DM_POLYTOPE_SEGMENT,
860:     1,
861:     2,
862:     2,
863:     DM_POLYTOPE_SEGMENT,
864:     1,
865:     0,
866:     1,
867:     DM_POLYTOPE_SEGMENT,
868:     0,
869:     1,
870:     DM_POLYTOPE_SEGMENT,
871:     1,
872:     3,
873:     2,
874:     DM_POLYTOPE_SEGMENT,
875:     1,
876:     0,
877:     0,
878:     DM_POLYTOPE_SEGMENT,
879:     0,
880:     2,
881:     DM_POLYTOPE_SEGMENT,
882:     1,
883:     4,
884:     2,
885:     DM_POLYTOPE_SEGMENT,
886:     1,
887:     0,
888:     3,
889:     DM_POLYTOPE_SEGMENT,
890:     1,
891:     1,
892:     0,
893:     DM_POLYTOPE_SEGMENT,
894:     0,
895:     0,
896:     DM_POLYTOPE_SEGMENT,
897:     1,
898:     0,
899:     2,
900:     DM_POLYTOPE_SEGMENT,
901:     1,
902:     2,
903:     0,
904:     DM_POLYTOPE_SEGMENT,
905:     0,
906:     1,
907:     DM_POLYTOPE_SEGMENT,
908:     1,
909:     0,
910:     1,
911:     DM_POLYTOPE_SEGMENT,
912:     1,
913:     3,
914:     0,
915:     DM_POLYTOPE_SEGMENT,
916:     0,
917:     2,
918:     DM_POLYTOPE_SEGMENT,
919:     1,
920:     0,
921:     0,
922:     DM_POLYTOPE_SEGMENT,
923:     1,
924:     4,
925:     0,
926:     DM_POLYTOPE_SEGMENT,
927:     0,
928:     3,
929:     /* This quad faces out of the bottom pyramid */
930:     DM_POLYTOPE_SEGMENT,
931:     1,
932:     1,
933:     1,
934:     DM_POLYTOPE_SEGMENT,
935:     1,
936:     2,
937:     1,
938:     DM_POLYTOPE_SEGMENT,
939:     1,
940:     3,
941:     1,
942:     DM_POLYTOPE_SEGMENT,
943:     1,
944:     4,
945:     1,
946:     /* The bottom face of each tet is on the triangular face */
947:     DM_POLYTOPE_TRIANGLE,
948:     1,
949:     1,
950:     3,
951:     DM_POLYTOPE_TRIANGLE,
952:     0,
953:     8,
954:     DM_POLYTOPE_TRIANGLE,
955:     0,
956:     4,
957:     DM_POLYTOPE_TRIANGLE,
958:     0,
959:     0,
960:     DM_POLYTOPE_TRIANGLE,
961:     1,
962:     2,
963:     3,
964:     DM_POLYTOPE_TRIANGLE,
965:     0,
966:     9,
967:     DM_POLYTOPE_TRIANGLE,
968:     0,
969:     5,
970:     DM_POLYTOPE_TRIANGLE,
971:     0,
972:     1,
973:     DM_POLYTOPE_TRIANGLE,
974:     1,
975:     3,
976:     3,
977:     DM_POLYTOPE_TRIANGLE,
978:     0,
979:     10,
980:     DM_POLYTOPE_TRIANGLE,
981:     0,
982:     6,
983:     DM_POLYTOPE_TRIANGLE,
984:     0,
985:     2,
986:     DM_POLYTOPE_TRIANGLE,
987:     1,
988:     4,
989:     3,
990:     DM_POLYTOPE_TRIANGLE,
991:     0,
992:     11,
993:     DM_POLYTOPE_TRIANGLE,
994:     0,
995:     7,
996:     DM_POLYTOPE_TRIANGLE,
997:     0,
998:     3,
999:     /* The front face of all pyramids is toward the front */
1000:     DM_POLYTOPE_QUADRILATERAL,
1001:     1,
1002:     0,
1003:     0,
1004:     DM_POLYTOPE_TRIANGLE,
1005:     1,
1006:     1,
1007:     0,
1008:     DM_POLYTOPE_TRIANGLE,
1009:     0,
1010:     4,
1011:     DM_POLYTOPE_TRIANGLE,
1012:     0,
1013:     11,
1014:     DM_POLYTOPE_TRIANGLE,
1015:     1,
1016:     4,
1017:     1,
1018:     DM_POLYTOPE_QUADRILATERAL,
1019:     1,
1020:     0,
1021:     3,
1022:     DM_POLYTOPE_TRIANGLE,
1023:     1,
1024:     1,
1025:     1,
1026:     DM_POLYTOPE_TRIANGLE,
1027:     1,
1028:     2,
1029:     0,
1030:     DM_POLYTOPE_TRIANGLE,
1031:     0,
1032:     5,
1033:     DM_POLYTOPE_TRIANGLE,
1034:     0,
1035:     8,
1036:     DM_POLYTOPE_QUADRILATERAL,
1037:     1,
1038:     0,
1039:     2,
1040:     DM_POLYTOPE_TRIANGLE,
1041:     0,
1042:     9,
1043:     DM_POLYTOPE_TRIANGLE,
1044:     1,
1045:     2,
1046:     1,
1047:     DM_POLYTOPE_TRIANGLE,
1048:     1,
1049:     3,
1050:     0,
1051:     DM_POLYTOPE_TRIANGLE,
1052:     0,
1053:     6,
1054:     DM_POLYTOPE_QUADRILATERAL,
1055:     1,
1056:     0,
1057:     1,
1058:     DM_POLYTOPE_TRIANGLE,
1059:     0,
1060:     7,
1061:     DM_POLYTOPE_TRIANGLE,
1062:     0,
1063:     10,
1064:     DM_POLYTOPE_TRIANGLE,
1065:     1,
1066:     3,
1067:     1,
1068:     DM_POLYTOPE_TRIANGLE,
1069:     1,
1070:     4,
1071:     0,
1072:     DM_POLYTOPE_QUADRILATERAL,
1073:     0,
1074:     0,
1075:     DM_POLYTOPE_TRIANGLE,
1076:     1,
1077:     1,
1078:     2,
1079:     DM_POLYTOPE_TRIANGLE,
1080:     1,
1081:     2,
1082:     2,
1083:     DM_POLYTOPE_TRIANGLE,
1084:     1,
1085:     3,
1086:     2,
1087:     DM_POLYTOPE_TRIANGLE,
1088:     1,
1089:     4,
1090:     2,
1091:     DM_POLYTOPE_QUADRILATERAL,
1092:     0,
1093:     0,
1094:     DM_POLYTOPE_TRIANGLE,
1095:     0,
1096:     0,
1097:     DM_POLYTOPE_TRIANGLE,
1098:     0,
1099:     3,
1100:     DM_POLYTOPE_TRIANGLE,
1101:     0,
1102:     2,
1103:     DM_POLYTOPE_TRIANGLE,
1104:     0,
1105:     1,
1106:   };
1107:   static PetscInt tpyrO[] = {0,  0, 0,  0,  0,  0, 0,  0,  0,  0, -1, 0,  0,  -1, 0,  0,  -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0,  -1, 0, 0, -1, 0, 0, -1, -1, -1,
1108:                              -1, 0, -3, -2, -3, 0, -3, -2, -3, 0, -3, -2, -3, 0,  -3, -2, -3, 0, 0, 0,  0, 0,  0, 0, 0,  0, 0, 0,  0, 0, 0,  0, 0,  0, 0, 0,  0, -2, 0,  0, 0, 0,  1, 0, 0,  0,  0};

1110:   if (rt) *rt = 0;
1111:   switch (source) {
1112:   case DM_POLYTOPE_POINT:
1113:     *Nt     = 1;
1114:     *target = vertexT;
1115:     *size   = vertexS;
1116:     *cone   = vertexC;
1117:     *ornt   = vertexO;
1118:     break;
1119:   case DM_POLYTOPE_SEGMENT:
1120:     *Nt     = 2;
1121:     *target = segT;
1122:     *size   = segS;
1123:     *cone   = segC;
1124:     *ornt   = segO;
1125:     break;
1126:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1127:     *Nt     = 1;
1128:     *target = tvertT;
1129:     *size   = tvertS;
1130:     *cone   = tvertC;
1131:     *ornt   = tvertO;
1132:     break;
1133:   case DM_POLYTOPE_TRIANGLE:
1134:     *Nt     = 2;
1135:     *target = triT;
1136:     *size   = triS;
1137:     *cone   = triC;
1138:     *ornt   = triO;
1139:     break;
1140:   case DM_POLYTOPE_QUADRILATERAL:
1141:     *Nt     = 3;
1142:     *target = quadT;
1143:     *size   = quadS;
1144:     *cone   = quadC;
1145:     *ornt   = quadO;
1146:     break;
1147:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1148:     *Nt     = 2;
1149:     *target = tsegT;
1150:     *size   = tsegS;
1151:     *cone   = tsegC;
1152:     *ornt   = tsegO;
1153:     break;
1154:   case DM_POLYTOPE_TETRAHEDRON:
1155:     *Nt     = 3;
1156:     *target = tetT;
1157:     *size   = tetS;
1158:     *cone   = tetC;
1159:     *ornt   = tetO;
1160:     break;
1161:   case DM_POLYTOPE_HEXAHEDRON:
1162:     *Nt     = 4;
1163:     *target = hexT;
1164:     *size   = hexS;
1165:     *cone   = hexC;
1166:     *ornt   = hexO;
1167:     break;
1168:   case DM_POLYTOPE_TRI_PRISM:
1169:     *Nt     = 4;
1170:     *target = tripT;
1171:     *size   = tripS;
1172:     *cone   = tripC;
1173:     *ornt   = tripO;
1174:     break;
1175:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1176:     *Nt     = 2;
1177:     *target = ttriT;
1178:     *size   = ttriS;
1179:     *cone   = ttriC;
1180:     *ornt   = ttriO;
1181:     break;
1182:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1183:     *Nt     = 3;
1184:     *target = tquadT;
1185:     *size   = tquadS;
1186:     *cone   = tquadC;
1187:     *ornt   = tquadO;
1188:     break;
1189:   case DM_POLYTOPE_PYRAMID:
1190:     *Nt     = 5;
1191:     *target = tpyrT;
1192:     *size   = tpyrS;
1193:     *cone   = tpyrC;
1194:     *ornt   = tpyrO;
1195:     break;
1196:   default:
1197:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1198:   }
1199:   return 0;
1200: }

1202: static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr)
1203: {
1204:   tr->ops->view                  = DMPlexTransformView_Regular;
1205:   tr->ops->setup                 = DMPlexTransformSetUp_Regular;
1206:   tr->ops->destroy               = DMPlexTransformDestroy_Regular;
1207:   tr->ops->celltransform         = DMPlexTransformCellRefine_Regular;
1208:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular;
1209:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
1210:   return 0;
1211: }

1213: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr)
1214: {
1215:   DMPlexRefine_Regular *f;

1218:   PetscNew(&f);
1219:   tr->data = f;

1221:   DMPlexTransformInitialize_Regular(tr);
1222:   return 0;
1223: }