Actual source code: plexrefregular.c

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

3: /*
4:    Regular Refinement of Hybrid Meshes

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

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

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

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

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

23:   Input Parameters:
24: + tr - The `DMPlexTransform` object
25: - ct - The cell type

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

34:   Level: developer

36:   Note:
37:   The Jacobian and inverse Jacobian will be rectangular, and the inverse is really a generalized inverse.
38: .vb
39:     v0 + j x_face = x_cell
40:     invj (x_cell - v0) = x_face
41: .ve

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

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

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

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

116:   Input Parameters:
117: + tr - The `DMPlexTransform` object
118: - ct - The cell type

120:   Output Parameters:
121: + Nc   - The number of subcells produced from this cell type
122: . v0   - The translation of the first vertex for each subcell
123: . J    - The Jacobian for each subcell (map from reference cell to subcell)
124: - invJ - The inverse Jacobian for each subcell

126:   Level: developer

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

150:                                  0.5, 0.0,  0.0, 0.5,

152:                                  0.5, 0.0,  0.0, 0.5,

154:                                  0.0, -0.5, 0.5, 0.5};
155:   static PetscReal tri_invJ[] = {2.0, 0.0, 0.0,  2.0,

157:                                  2.0, 0.0, 0.0,  2.0,

159:                                  2.0, 0.0, 0.0,  2.0,

161:                                  2.0, 2.0, -2.0, 0.0};
162:   static PetscReal tet_v0[]   = {-1.0, -1.0, -1.0, -1.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
163:   static PetscReal tet_J[]    = {0.5,  0.0,  0.0,  0.0, 0.5, 0.0,  0.0,  0.0,  0.5,

165:                                  0.5,  0.0,  0.0,  0.0, 0.5, 0.0,  0.0,  0.0,  0.5,

167:                                  0.5,  0.0,  0.0,  0.0, 0.5, 0.0,  0.0,  0.0,  0.5,

169:                                  0.5,  0.0,  0.0,  0.0, 0.5, 0.0,  0.0,  0.0,  0.5,

171:                                  -0.5, -0.5, 0.0,  0.5, 0.0, 0.0,  0.0,  0.5,  0.5,

173:                                  0.0,  0.5,  0.5,  0.0, 0.0, -0.5, -0.5, -0.5, 0.0,

175:                                  -0.5, 0.0,  0.0,  0.5, 0.0, 0.5,  -0.5, -0.5, -0.5,

177:                                  -0.5, -0.5, -0.5, 0.5, 0.5, 0.0,  0.0,  -0.5, 0.0};
178:   static PetscReal tet_invJ[] = {2.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  2.0,

180:                                  2.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  2.0,

182:                                  2.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  2.0,

184:                                  2.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  2.0,

186:                                  0.0,  2.0,  0.0,  -2.0, -2.0, 0.0,  2.0,  2.0,  2.0,

188:                                  -2.0, -2.0, -2.0, 2.0,  2.0,  0.0,  0.0,  -2.0, 0.0,

190:                                  -2.0, 0.0,  0.0,  0.0,  -2.0, -2.0, 2.0,  2.0,  0.0,

192:                                  0.0,  2.0,  2.0,  0.0,  0.0,  -2.0, -2.0, -2.0, 0.0};
193:   /*
194:      3---------2---------2
195:      |         |         |
196:      |    D    2    C    |
197:      |         |         |
198:      3----3----0----1----1
199:      |         |         |
200:      |    A    0    B    |
201:      |         |         |
202:      0---------0---------1
203:      */
204:   static PetscReal quad_v0[]   = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
205:   static PetscReal quad_J[]    = {0.5, 0.0, 0.0, 0.5,

207:                                   0.5, 0.0, 0.0, 0.5,

209:                                   0.5, 0.0, 0.0, 0.5,

211:                                   0.5, 0.0, 0.0, 0.5};
212:   static PetscReal quad_invJ[] = {2.0, 0.0, 0.0, 2.0,

214:                                   2.0, 0.0, 0.0, 2.0,

216:                                   2.0, 0.0, 0.0, 2.0,

218:                                   2.0, 0.0, 0.0, 2.0};
219:   /*
220:      Bottom (viewed from top)    Top
221:      1---------2---------2       7---------2---------6
222:      |         |         |       |         |         |
223:      |    B    2    C    |       |    H    2    G    |
224:      |         |         |       |         |         |
225:      3----3----0----1----1       3----3----0----1----1
226:      |         |         |       |         |         |
227:      |    A    0    D    |       |    E    0    F    |
228:      |         |         |       |         |         |
229:      0---------0---------3       4---------0---------5
230:      */
231:   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};
232:   static PetscReal hex_J[]    = {0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

234:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

236:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

238:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

240:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

242:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

244:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,

246:                                  0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
247:   static PetscReal hex_invJ[] = {2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

249:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

251:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

253:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

255:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

257:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

259:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,

261:                                  2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0};

263:   PetscFunctionBegin;
264:   switch (ct) {
265:   case DM_POLYTOPE_TRIANGLE:
266:     if (Nc) *Nc = 4;
267:     if (v0) *v0 = tri_v0;
268:     if (J) *J = tri_J;
269:     if (invJ) *invJ = tri_invJ;
270:     break;
271:   case DM_POLYTOPE_QUADRILATERAL:
272:     if (Nc) *Nc = 4;
273:     if (v0) *v0 = quad_v0;
274:     if (J) *J = quad_J;
275:     if (invJ) *invJ = quad_invJ;
276:     break;
277:   case DM_POLYTOPE_TETRAHEDRON:
278:     if (Nc) *Nc = 8;
279:     if (v0) *v0 = tet_v0;
280:     if (J) *J = tet_J;
281:     if (invJ) *invJ = tet_invJ;
282:     break;
283:   case DM_POLYTOPE_HEXAHEDRON:
284:     if (Nc) *Nc = 8;
285:     if (v0) *v0 = hex_v0;
286:     if (J) *J = hex_J;
287:     if (invJ) *invJ = hex_invJ;
288:     break;
289:   default:
290:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
291:   }
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: #if 0
296: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
297: {
298:   static PetscReal seg_v[]  = {-1.0,  0.0,  1.0};
299:   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};
300:   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};
301:   static PetscReal tet_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
302:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,  -1.0,  1.0, -1.0,
303:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,  -1.0,  0.0,  0.0,  -1.0, -1.0,  1.0};
304:   static PetscReal hex_v[]  = {-1.0, -1.0, -1.0,   0.0, -1.0, -1.0,   1.0, -1.0, -1.0,
305:                                -1.0,  0.0, -1.0,   0.0,  0.0, -1.0,   1.0,  0.0, -1.0,
306:                                -1.0,  1.0, -1.0,   0.0,  1.0, -1.0,   1.0,  1.0, -1.0,
307:                                -1.0, -1.0,  0.0,   0.0, -1.0,  0.0,   1.0, -1.0,  0.0,
308:                                -1.0,  0.0,  0.0,   0.0,  0.0,  0.0,   1.0,  0.0,  0.0,
309:                                -1.0,  1.0,  0.0,   0.0,  1.0,  0.0,   1.0,  1.0,  0.0,
310:                                -1.0, -1.0,  1.0,   0.0, -1.0,  1.0,   1.0, -1.0,  1.0,
311:                                -1.0,  0.0,  1.0,   0.0,  0.0,  1.0,   1.0,  0.0,  1.0,
312:                                -1.0,  1.0,  1.0,   0.0,  1.0,  1.0,   1.0,  1.0,  1.0};

314:   PetscFunctionBegin;
315:   switch (ct) {
316:     case DM_POLYTOPE_SEGMENT:       *Nv =  3; *subcellV = seg_v;  break;
317:     case DM_POLYTOPE_TRIANGLE:      *Nv =  6; *subcellV = tri_v;  break;
318:     case DM_POLYTOPE_QUADRILATERAL: *Nv =  9; *subcellV = quad_v; break;
319:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 10; *subcellV = tet_v;  break;
320:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 27; *subcellV = hex_v;  break;
321:     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
322:   }
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
327: {
328:   static PetscInt seg_v[]  = {0, 1, 1, 2};
329:   static PetscInt tri_v[]  = {0, 3, 5,  3, 1, 4,  5, 4, 2,  3, 4, 5};
330:   static PetscInt quad_v[] = {0, 4, 8, 7,  4, 1, 5, 8,  8, 5, 2, 6,  7, 8, 6, 3};
331:   static PetscInt tet_v[]  = {0, 3, 1, 6,  3, 2, 4, 8,  1, 4, 5, 7,  6, 8, 7, 9,
332:                               1, 6, 3, 7,  8, 4, 3, 7,  7, 3, 1, 4,  7, 3, 8, 6};
333:   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,
334:                               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};

336:   PetscFunctionBegin;
337:   PetscCheck(ct == rct,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
338:   switch (ct) {
339:     case DM_POLYTOPE_SEGMENT:       *Nv = 2; *subcellV = &seg_v[r*(*Nv)];  break;
340:     case DM_POLYTOPE_TRIANGLE:      *Nv = 3; *subcellV = &tri_v[r*(*Nv)];  break;
341:     case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
342:     case DM_POLYTOPE_TETRAHEDRON:   *Nv = 4; *subcellV = &tet_v[r*(*Nv)];  break;
343:     case DM_POLYTOPE_HEXAHEDRON:    *Nv = 8; *subcellV = &hex_v[r*(*Nv)];  break;
344:     default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
345:   }
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }
348: #endif

350: static PetscErrorCode DMPlexTransformView_Regular(DMPlexTransform tr, PetscViewer viewer)
351: {
352:   PetscBool isascii;

354:   PetscFunctionBegin;
357:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
358:   if (isascii) {
359:     const char *name;

361:     PetscCall(PetscObjectGetName((PetscObject)tr, &name));
362:     PetscCall(PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : ""));
363:   } else {
364:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
365:   }
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr)
370: {
371:   PetscFunctionBegin;
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr)
376: {
377:   DMPlexRefine_Regular *f = (DMPlexRefine_Regular *)tr->data;

379:   PetscFunctionBegin;
380:   PetscCall(PetscFree(f));
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
385: {
386:   static PetscInt seg_seg[]   = {1, -1, 0, -1, 0, 0, 1, 0};
387:   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};
388:   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};
389:   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};
390:   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};
391:   static PetscInt tseg_seg[]  = {0, -1, 0, 0, 0, 0, 0, -1};
392:   static PetscInt tseg_tseg[] = {1, -2, 0, -2, 1, -1, 0, -1, 0, 0, 1, 0, 0, 1, 1, 1};
393:   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};
394:   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,
395:                                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,
396:                                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,
397:                                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,
398:                                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,
399:                                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};
400:   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,
401:                                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,
402:                                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,
403:                                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,
404:                                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,
405:                                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,
406:                                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};
407:   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,
408:                                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,
409:                                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,
410:                                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,
411:                                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,
412:                                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,
413:                                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,
414:                                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};
415:   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,
416:                                    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,
417:                                    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,
418:                                    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,
419:                                    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,
420:                                    -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,
421:                                    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,
422:                                    -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,
423:                                    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,
424:                                    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,
425:                                    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,
426:                                    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,
427:                                    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,
428:                                    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,
429:                                    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,
430:                                    -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,
431:                                    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,
432:                                    -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,
433:                                    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,
434:                                    -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,
435:                                    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,
436:                                    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};
437:   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,
438:                                    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,
439:                                    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,
440:                                    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,
441:                                    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,
442:                                    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,
443:                                    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,
444:                                    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,
445:                                    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,
446:                                    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,
447:                                    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,
448:                                    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,
449:                                    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,
450:                                    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,
451:                                    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};
452:   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,
453:                                    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};
454:   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,
455:                                    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};
456:   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,
457:                                    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,
458:                                    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};
459:   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,
460:                                    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,
461:                                    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,
462:                                    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};
463:   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,
464:                                    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};
465:   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,
466:                                    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};
467:   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};
468:   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,
469:                                    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};
470:   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,
471:                                    -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,
472:                                    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};

474:   PetscFunctionBeginHot;
475:   *rnew = r;
476:   *onew = o;
477:   if (!so) PetscFunctionReturn(PETSC_SUCCESS);
478:   switch (sct) {
479:   case DM_POLYTOPE_POINT:
480:     break;
481:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
482:     *onew = so < 0 ? -(o + 1) : o;
483:     break;
484:   case DM_POLYTOPE_SEGMENT:
485:     switch (tct) {
486:     case DM_POLYTOPE_POINT:
487:       break;
488:     case DM_POLYTOPE_SEGMENT:
489:       *rnew = seg_seg[(so + 1) * 4 + r * 2];
490:       *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so + 1) * 4 + r * 2 + 1]);
491:       break;
492:     default:
493:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
494:     }
495:     break;
496:   case DM_POLYTOPE_TRIANGLE:
497:     switch (tct) {
498:     case DM_POLYTOPE_SEGMENT:
499:       *rnew = tri_seg[(so + 3) * 6 + r * 2];
500:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
501:       break;
502:     case DM_POLYTOPE_TRIANGLE:
503:       *rnew = tri_tri[(so + 3) * 8 + r * 2];
504:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 8 + 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_QUADRILATERAL:
511:     switch (tct) {
512:     case DM_POLYTOPE_POINT:
513:       break;
514:     case DM_POLYTOPE_SEGMENT:
515:       *rnew = quad_seg[(so + 4) * 8 + r * 2];
516:       *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]);
517:       break;
518:     case DM_POLYTOPE_QUADRILATERAL:
519:       *rnew = quad_quad[(so + 4) * 8 + r * 2];
520:       *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so + 4) * 8 + r * 2 + 1]);
521:       break;
522:     default:
523:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
524:     }
525:     break;
526:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
527:     switch (tct) {
528:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
529:       *rnew = tseg_seg[(so + 2) * 2 + r * 2];
530:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so + 2) * 2 + r * 2 + 1]);
531:       break;
532:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
533:       *rnew = tseg_tseg[(so + 2) * 4 + r * 2];
534:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so + 2) * 4 + r * 2 + 1]);
535:       break;
536:     default:
537:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
538:     }
539:     break;
540:   case DM_POLYTOPE_TETRAHEDRON:
541:     switch (tct) {
542:     case DM_POLYTOPE_SEGMENT:
543:       *rnew = tet_seg[(so + 12) * 2 + r * 2];
544:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 2 + r * 2 + 1]);
545:       break;
546:     case DM_POLYTOPE_TRIANGLE:
547:       *rnew = tet_tri[(so + 12) * 16 + r * 2];
548:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 16 + r * 2 + 1]);
549:       break;
550:     case DM_POLYTOPE_TETRAHEDRON:
551:       *rnew = tet_tet[(so + 12) * 16 + r * 2];
552:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 16 + r * 2 + 1]);
553:       break;
554:     default:
555:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
556:     }
557:     break;
558:   case DM_POLYTOPE_HEXAHEDRON:
559:     switch (tct) {
560:     case DM_POLYTOPE_POINT:
561:       break;
562:     case DM_POLYTOPE_SEGMENT:
563:       *rnew = hex_seg[(so + 24) * 12 + r * 2];
564:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so + 24) * 12 + r * 2 + 1]);
565:       break;
566:     case DM_POLYTOPE_QUADRILATERAL:
567:       *rnew = hex_quad[(so + 24) * 24 + r * 2];
568:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so + 24) * 24 + r * 2 + 1]);
569:       break;
570:     case DM_POLYTOPE_HEXAHEDRON:
571:       *rnew = hex_hex[(so + 24) * 16 + r * 2];
572:       *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so + 24) * 16 + r * 2 + 1]);
573:       break;
574:     default:
575:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
576:     }
577:     break;
578:   case DM_POLYTOPE_TRI_PRISM:
579:     switch (tct) {
580:     case DM_POLYTOPE_SEGMENT:
581:       *rnew = trip_seg[(so + 6) * 6 + r * 2];
582:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so + 6) * 6 + r * 2 + 1]);
583:       break;
584:     case DM_POLYTOPE_TRIANGLE:
585:       *rnew = trip_tri[(so + 6) * 8 + r * 2];
586:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so + 6) * 8 + r * 2 + 1]);
587:       break;
588:     case DM_POLYTOPE_QUADRILATERAL:
589:       *rnew = trip_quad[(so + 6) * 12 + r * 2];
590:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so + 6) * 12 + r * 2 + 1]);
591:       break;
592:     case DM_POLYTOPE_TRI_PRISM:
593:       *rnew = trip_trip[(so + 6) * 16 + r * 2];
594:       *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so + 6) * 16 + r * 2 + 1]);
595:       break;
596:     default:
597:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
598:     }
599:     break;
600:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
601:     switch (tct) {
602:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
603:       *rnew = ttri_tseg[(so + 6) * 6 + r * 2];
604:       *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so + 6) * 6 + r * 2 + 1]);
605:       break;
606:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
607:       *rnew = ttri_ttri[(so + 6) * 8 + r * 2];
608:       *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so + 6) * 8 + r * 2 + 1]);
609:       break;
610:     default:
611:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
612:     }
613:     break;
614:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
615:     switch (tct) {
616:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
617:       *rnew = tquad_tvert[(so + 8) * 2 + r * 2];
618:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so + 8) * 2 + r * 2 + 1]);
619:       break;
620:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
621:       *rnew = tquad_tseg[(so + 8) * 8 + r * 2];
622:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so + 8) * 8 + r * 2 + 1]);
623:       break;
624:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
625:       *rnew = tquad_tquad[(so + 8) * 8 + r * 2];
626:       *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so + 8) * 8 + r * 2 + 1]);
627:       break;
628:     default:
629:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
630:     }
631:     break;
632:   default:
633:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
634:   }
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }

638: PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
639: {
640:   /* All vertices remain in the refined mesh */
641:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
642:   static PetscInt       vertexS[] = {1};
643:   static PetscInt       vertexC[] = {0};
644:   static PetscInt       vertexO[] = {0};
645:   /* Split all edges with a new vertex, making two new 2 edges
646:      0--0--0--1--1
647:   */
648:   static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
649:   static PetscInt       segS[] = {1, 2};
650:   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};
651:   static PetscInt       segO[] = {0, 0, 0, 0};
652:   /* Do not split tensor edges */
653:   static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
654:   static PetscInt       tvertS[] = {1};
655:   static PetscInt       tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
656:   static PetscInt       tvertO[] = {0, 0};
657:   /* Add 3 edges inside every triangle, making 4 new triangles.
658:    2
659:    |\
660:    | \
661:    |  \
662:    0   1
663:    | C  \
664:    |     \
665:    |      \
666:    2---1---1
667:    |\  D  / \
668:    1 2   0   0
669:    |A \ /  B  \
670:    0-0-0---1---1
671:   */
672:   static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
673:   static PetscInt       triS[] = {3, 4};
674:   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};
675:   static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0};
676:   /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
677:      3----1----2----0----2
678:      |         |         |
679:      0    D    2    C    1
680:      |         |         |
681:      3----3----0----1----1
682:      |         |         |
683:      1    A    0    B    0
684:      |         |         |
685:      0----0----0----1----1
686:   */
687:   static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
688:   static PetscInt       quadS[] = {1, 4, 4};
689:   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};
690:   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};
691:   /* Add 1 edge inside every tensor quad, making 2 new tensor quads
692:      2----2----1----3----3
693:      |         |         |
694:      |         |         |
695:      |         |         |
696:      4    A    6    B    5
697:      |         |         |
698:      |         |         |
699:      |         |         |
700:      0----0----0----1----1
701:   */
702:   static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
703:   static PetscInt       tsegS[] = {1, 2};
704:   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};
705:   static PetscInt tsegO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
706:   /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
707:      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
708:      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]
709:      called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
710:        [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
711:      The first four tets just cut off the corners, using the replica notation for new vertices,
712:        [v0,      (e0, 0), (e2, 0), (e3, 0)]
713:        [(e0, 0), v1,      (e1, 0), (e4, 0)]
714:        [(e2, 0), (e1, 0), v2,      (e5, 0)]
715:        [(e3, 0), (e4, 0), (e5, 0), v3     ]
716:      The next four tets match a vertex to the newly created faces from cutting off those first tets.
717:        [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
718:        [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
719:        [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
720:        [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
721:      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
722:        [(e2, 0), (e0, 0), (e3, 0)]
723:        [(e0, 0), (e1, 0), (e4, 0)]
724:        [(e2, 0), (e5, 0), (e1, 0)]
725:        [(e3, 0), (e4, 0), (e5, 0)]
726:      The next four, from the second group of tets, are
727:        [(e2, 0), (e0, 0), (e5, 0)]
728:        [(e4, 0), (e0, 0), (e5, 0)]
729:        [(e0, 0), (e1, 0), (e5, 0)]
730:        [(e5, 0), (e3, 0), (e0, 0)]
731:      I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
732:    */
733:   static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
734:   static PetscInt       tetS[] = {1, 8, 8};
735:   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};
736:   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};
737:   /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
738:      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
739:        [v0, v1, v2, v3] f0 bottom
740:        [v4, v5, v6, v7] f1 top
741:        [v0, v3, v5, v4] f2 front
742:        [v2, v1, v7, v6] f3 back
743:        [v3, v2, v6, v5] f4 right
744:        [v0, v4, v7, v1] f5 left
745:      The eight hexes are divided into four on the bottom, and four on the top,
746:        [v0,      (e0, 0),  (f0, 0),  (e3, 0),  (e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
747:        [(e0, 0), v1,       (e1, 0),  (f0, 0),  (f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
748:        [(f0, 0), (e1, 0),  v2,       (e2, 0),  (c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
749:        [(e3, 0), (f0, 0),  (e2, 0),  v3,       (f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
750:        [(e9, 0), (f5, 0),  (c0, 0),  (f2, 0),  v4,      (e4, 0),  (f1, 0),  (e7, 0)]
751:        [(f2, 0), (c0, 0),  (f4, 0),  (e8, 0),  (e4, 0), v5,       (e5, 0),  (f1, 0)]
752:        [(c0, 0), (f3, 0),  (e11, 0), (f4, 0),  (f1, 0), (e5, 0),  v6,       (e6, 0)]
753:        [(f5, 0), (e10, 0), (f3, 0),  (c0, 0),  (e7, 0), (f1, 0),  (e6, 0),  v7]
754:      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,
755:        [(e9, 0), (f2, 0),  (c0, 0),  (f5, 0)]
756:        [(f5, 0), (c0, 0),  (f3, 0),  (e10, 0)]
757:        [(c0, 0), (f4, 0),  (e11, 0), (f3, 0)]
758:        [(f2, 0), (e8, 0),  (f4, 0),  (c0, 0)]
759:      and on the x-z plane,
760:        [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
761:        [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
762:        [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
763:        [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
764:      and on the y-z plane,
765:        [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
766:        [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
767:        [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
768:        [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
769:   */
770:   static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
771:   static PetscInt       hexS[] = {1, 6, 12, 8};
772:   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};
773:   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,
774:                             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};
775:   /* Add 3 quads inside every triangular prism, making 4 new prisms. */
776:   static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
777:   static PetscInt       tripS[] = {3, 4, 6, 8};
778:   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};
779:   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,
780:                              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};
781:   /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
782:       2
783:       |\
784:       | \
785:       |  \
786:       0---1

788:       2

790:       0   1

792:       2
793:       |\
794:       | \
795:       |  \
796:       0---1
797:   */
798:   static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
799:   static PetscInt       ttriS[] = {3, 4};
800:   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};
801:   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};
802:   /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
803:   static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
804:   static PetscInt       tquadS[] = {1, 4, 4};
805:   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};
806:   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};
807:   /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */
808:   static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID};
809:   static PetscInt       tpyrS[] = {4, 12, 1, 4, 6};
810:   static PetscInt       tpyrC[] = {
811:     DM_POLYTOPE_POINT,
812:     2,
813:     1,
814:     1,
815:     0,
816:     DM_POLYTOPE_POINT,
817:     1,
818:     0,
819:     0,
820:     DM_POLYTOPE_POINT,
821:     2,
822:     2,
823:     1,
824:     0,
825:     DM_POLYTOPE_POINT,
826:     1,
827:     0,
828:     0,
829:     DM_POLYTOPE_POINT,
830:     2,
831:     3,
832:     1,
833:     0,
834:     DM_POLYTOPE_POINT,
835:     1,
836:     0,
837:     0,
838:     DM_POLYTOPE_POINT,
839:     2,
840:     4,
841:     1,
842:     0,
843:     DM_POLYTOPE_POINT,
844:     1,
845:     0,
846:     0,
847:     /* These four triangle face out of the bottom pyramid */
848:     DM_POLYTOPE_SEGMENT,
849:     1,
850:     1,
851:     1,
852:     DM_POLYTOPE_SEGMENT,
853:     0,
854:     3,
855:     DM_POLYTOPE_SEGMENT,
856:     0,
857:     0,
858:     DM_POLYTOPE_SEGMENT,
859:     1,
860:     2,
861:     1,
862:     DM_POLYTOPE_SEGMENT,
863:     0,
864:     0,
865:     DM_POLYTOPE_SEGMENT,
866:     0,
867:     1,
868:     DM_POLYTOPE_SEGMENT,
869:     1,
870:     3,
871:     1,
872:     DM_POLYTOPE_SEGMENT,
873:     0,
874:     1,
875:     DM_POLYTOPE_SEGMENT,
876:     0,
877:     2,
878:     DM_POLYTOPE_SEGMENT,
879:     1,
880:     4,
881:     1,
882:     DM_POLYTOPE_SEGMENT,
883:     0,
884:     2,
885:     DM_POLYTOPE_SEGMENT,
886:     0,
887:     3,
888:     /* These eight triangles face out of the corner pyramids */
889:     DM_POLYTOPE_SEGMENT,
890:     1,
891:     0,
892:     3,
893:     DM_POLYTOPE_SEGMENT,
894:     0,
895:     3,
896:     DM_POLYTOPE_SEGMENT,
897:     1,
898:     1,
899:     2,
900:     DM_POLYTOPE_SEGMENT,
901:     1,
902:     0,
903:     2,
904:     DM_POLYTOPE_SEGMENT,
905:     0,
906:     0,
907:     DM_POLYTOPE_SEGMENT,
908:     1,
909:     2,
910:     2,
911:     DM_POLYTOPE_SEGMENT,
912:     1,
913:     0,
914:     1,
915:     DM_POLYTOPE_SEGMENT,
916:     0,
917:     1,
918:     DM_POLYTOPE_SEGMENT,
919:     1,
920:     3,
921:     2,
922:     DM_POLYTOPE_SEGMENT,
923:     1,
924:     0,
925:     0,
926:     DM_POLYTOPE_SEGMENT,
927:     0,
928:     2,
929:     DM_POLYTOPE_SEGMENT,
930:     1,
931:     4,
932:     2,
933:     DM_POLYTOPE_SEGMENT,
934:     1,
935:     0,
936:     3,
937:     DM_POLYTOPE_SEGMENT,
938:     1,
939:     1,
940:     0,
941:     DM_POLYTOPE_SEGMENT,
942:     0,
943:     0,
944:     DM_POLYTOPE_SEGMENT,
945:     1,
946:     0,
947:     2,
948:     DM_POLYTOPE_SEGMENT,
949:     1,
950:     2,
951:     0,
952:     DM_POLYTOPE_SEGMENT,
953:     0,
954:     1,
955:     DM_POLYTOPE_SEGMENT,
956:     1,
957:     0,
958:     1,
959:     DM_POLYTOPE_SEGMENT,
960:     1,
961:     3,
962:     0,
963:     DM_POLYTOPE_SEGMENT,
964:     0,
965:     2,
966:     DM_POLYTOPE_SEGMENT,
967:     1,
968:     0,
969:     0,
970:     DM_POLYTOPE_SEGMENT,
971:     1,
972:     4,
973:     0,
974:     DM_POLYTOPE_SEGMENT,
975:     0,
976:     3,
977:     /* This quad faces out of the bottom pyramid */
978:     DM_POLYTOPE_SEGMENT,
979:     1,
980:     1,
981:     1,
982:     DM_POLYTOPE_SEGMENT,
983:     1,
984:     2,
985:     1,
986:     DM_POLYTOPE_SEGMENT,
987:     1,
988:     3,
989:     1,
990:     DM_POLYTOPE_SEGMENT,
991:     1,
992:     4,
993:     1,
994:     /* The bottom face of each tet is on the triangular face */
995:     DM_POLYTOPE_TRIANGLE,
996:     1,
997:     1,
998:     3,
999:     DM_POLYTOPE_TRIANGLE,
1000:     0,
1001:     8,
1002:     DM_POLYTOPE_TRIANGLE,
1003:     0,
1004:     4,
1005:     DM_POLYTOPE_TRIANGLE,
1006:     0,
1007:     0,
1008:     DM_POLYTOPE_TRIANGLE,
1009:     1,
1010:     2,
1011:     3,
1012:     DM_POLYTOPE_TRIANGLE,
1013:     0,
1014:     9,
1015:     DM_POLYTOPE_TRIANGLE,
1016:     0,
1017:     5,
1018:     DM_POLYTOPE_TRIANGLE,
1019:     0,
1020:     1,
1021:     DM_POLYTOPE_TRIANGLE,
1022:     1,
1023:     3,
1024:     3,
1025:     DM_POLYTOPE_TRIANGLE,
1026:     0,
1027:     10,
1028:     DM_POLYTOPE_TRIANGLE,
1029:     0,
1030:     6,
1031:     DM_POLYTOPE_TRIANGLE,
1032:     0,
1033:     2,
1034:     DM_POLYTOPE_TRIANGLE,
1035:     1,
1036:     4,
1037:     3,
1038:     DM_POLYTOPE_TRIANGLE,
1039:     0,
1040:     11,
1041:     DM_POLYTOPE_TRIANGLE,
1042:     0,
1043:     7,
1044:     DM_POLYTOPE_TRIANGLE,
1045:     0,
1046:     3,
1047:     /* The front face of all pyramids is toward the front */
1048:     DM_POLYTOPE_QUADRILATERAL,
1049:     1,
1050:     0,
1051:     0,
1052:     DM_POLYTOPE_TRIANGLE,
1053:     1,
1054:     1,
1055:     0,
1056:     DM_POLYTOPE_TRIANGLE,
1057:     0,
1058:     4,
1059:     DM_POLYTOPE_TRIANGLE,
1060:     0,
1061:     11,
1062:     DM_POLYTOPE_TRIANGLE,
1063:     1,
1064:     4,
1065:     1,
1066:     DM_POLYTOPE_QUADRILATERAL,
1067:     1,
1068:     0,
1069:     3,
1070:     DM_POLYTOPE_TRIANGLE,
1071:     1,
1072:     1,
1073:     1,
1074:     DM_POLYTOPE_TRIANGLE,
1075:     1,
1076:     2,
1077:     0,
1078:     DM_POLYTOPE_TRIANGLE,
1079:     0,
1080:     5,
1081:     DM_POLYTOPE_TRIANGLE,
1082:     0,
1083:     8,
1084:     DM_POLYTOPE_QUADRILATERAL,
1085:     1,
1086:     0,
1087:     2,
1088:     DM_POLYTOPE_TRIANGLE,
1089:     0,
1090:     9,
1091:     DM_POLYTOPE_TRIANGLE,
1092:     1,
1093:     2,
1094:     1,
1095:     DM_POLYTOPE_TRIANGLE,
1096:     1,
1097:     3,
1098:     0,
1099:     DM_POLYTOPE_TRIANGLE,
1100:     0,
1101:     6,
1102:     DM_POLYTOPE_QUADRILATERAL,
1103:     1,
1104:     0,
1105:     1,
1106:     DM_POLYTOPE_TRIANGLE,
1107:     0,
1108:     7,
1109:     DM_POLYTOPE_TRIANGLE,
1110:     0,
1111:     10,
1112:     DM_POLYTOPE_TRIANGLE,
1113:     1,
1114:     3,
1115:     1,
1116:     DM_POLYTOPE_TRIANGLE,
1117:     1,
1118:     4,
1119:     0,
1120:     DM_POLYTOPE_QUADRILATERAL,
1121:     0,
1122:     0,
1123:     DM_POLYTOPE_TRIANGLE,
1124:     1,
1125:     1,
1126:     2,
1127:     DM_POLYTOPE_TRIANGLE,
1128:     1,
1129:     2,
1130:     2,
1131:     DM_POLYTOPE_TRIANGLE,
1132:     1,
1133:     3,
1134:     2,
1135:     DM_POLYTOPE_TRIANGLE,
1136:     1,
1137:     4,
1138:     2,
1139:     DM_POLYTOPE_QUADRILATERAL,
1140:     0,
1141:     0,
1142:     DM_POLYTOPE_TRIANGLE,
1143:     0,
1144:     0,
1145:     DM_POLYTOPE_TRIANGLE,
1146:     0,
1147:     3,
1148:     DM_POLYTOPE_TRIANGLE,
1149:     0,
1150:     2,
1151:     DM_POLYTOPE_TRIANGLE,
1152:     0,
1153:     1,
1154:   };
1155:   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,
1156:                              -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};

1158:   PetscFunctionBegin;
1159:   if (rt) *rt = 0;
1160:   switch (source) {
1161:   case DM_POLYTOPE_POINT:
1162:     *Nt     = 1;
1163:     *target = vertexT;
1164:     *size   = vertexS;
1165:     *cone   = vertexC;
1166:     *ornt   = vertexO;
1167:     break;
1168:   case DM_POLYTOPE_SEGMENT:
1169:     *Nt     = 2;
1170:     *target = segT;
1171:     *size   = segS;
1172:     *cone   = segC;
1173:     *ornt   = segO;
1174:     break;
1175:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1176:     *Nt     = 1;
1177:     *target = tvertT;
1178:     *size   = tvertS;
1179:     *cone   = tvertC;
1180:     *ornt   = tvertO;
1181:     break;
1182:   case DM_POLYTOPE_TRIANGLE:
1183:     *Nt     = 2;
1184:     *target = triT;
1185:     *size   = triS;
1186:     *cone   = triC;
1187:     *ornt   = triO;
1188:     break;
1189:   case DM_POLYTOPE_QUADRILATERAL:
1190:     *Nt     = 3;
1191:     *target = quadT;
1192:     *size   = quadS;
1193:     *cone   = quadC;
1194:     *ornt   = quadO;
1195:     break;
1196:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1197:     *Nt     = 2;
1198:     *target = tsegT;
1199:     *size   = tsegS;
1200:     *cone   = tsegC;
1201:     *ornt   = tsegO;
1202:     break;
1203:   case DM_POLYTOPE_TETRAHEDRON:
1204:     *Nt     = 3;
1205:     *target = tetT;
1206:     *size   = tetS;
1207:     *cone   = tetC;
1208:     *ornt   = tetO;
1209:     break;
1210:   case DM_POLYTOPE_HEXAHEDRON:
1211:     *Nt     = 4;
1212:     *target = hexT;
1213:     *size   = hexS;
1214:     *cone   = hexC;
1215:     *ornt   = hexO;
1216:     break;
1217:   case DM_POLYTOPE_TRI_PRISM:
1218:     *Nt     = 4;
1219:     *target = tripT;
1220:     *size   = tripS;
1221:     *cone   = tripC;
1222:     *ornt   = tripO;
1223:     break;
1224:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
1225:     *Nt     = 2;
1226:     *target = ttriT;
1227:     *size   = ttriS;
1228:     *cone   = ttriC;
1229:     *ornt   = ttriO;
1230:     break;
1231:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1232:     *Nt     = 3;
1233:     *target = tquadT;
1234:     *size   = tquadS;
1235:     *cone   = tquadC;
1236:     *ornt   = tquadO;
1237:     break;
1238:   case DM_POLYTOPE_PYRAMID:
1239:     *Nt     = 5;
1240:     *target = tpyrT;
1241:     *size   = tpyrS;
1242:     *cone   = tpyrC;
1243:     *ornt   = tpyrO;
1244:     break;
1245:   default:
1246:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1247:   }
1248:   PetscFunctionReturn(PETSC_SUCCESS);
1249: }

1251: static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr)
1252: {
1253:   PetscFunctionBegin;
1254:   tr->ops->view                  = DMPlexTransformView_Regular;
1255:   tr->ops->setup                 = DMPlexTransformSetUp_Regular;
1256:   tr->ops->destroy               = DMPlexTransformDestroy_Regular;
1257:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Internal;
1258:   tr->ops->celltransform         = DMPlexTransformCellRefine_Regular;
1259:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular;
1260:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
1261:   PetscFunctionReturn(PETSC_SUCCESS);
1262: }

1264: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr)
1265: {
1266:   DMPlexRefine_Regular *f;

1268:   PetscFunctionBegin;
1270:   PetscCall(PetscNew(&f));
1271:   tr->data = f;

1273:   PetscCall(DMPlexTransformInitialize_Regular(tr));
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }
```