Actual source code: plexrefregular.c
1: #include <petsc/private/dmplextransformimpl.h>
3: /* Regular Refinment 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};
88: switch (ct) {
89: case DM_POLYTOPE_TRIANGLE: if (Nf) *Nf = 3; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; if (detJ) *detJ = tri_detJ; break;
90: case DM_POLYTOPE_QUADRILATERAL: if (Nf) *Nf = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; if (detJ) *detJ = quad_detJ; break;
91: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
92: }
93: return(0);
94: }
96: /*@
97: DMPlexRefineRegularGetAffineTransforms - Gets the affine map from the reference cell to each subcell
99: Input Parameters:
100: + cr - The DMPlexCellRefiner object
101: - ct - The cell type
103: Output Parameters:
104: + Nc - The number of subcells produced from this cell type
105: . v0 - The translation of the first vertex for each subcell
106: . J - The Jacobian for each subcell (map from reference cell to subcell)
107: - invJ - The inverse Jacobian for each subcell
109: Level: developer
111: .seealso: DMPlexRefineRegularGetAffineFaceTransforms(), DMPLEXREFINEREGULAR
112: @*/
113: PetscErrorCode DMPlexRefineRegularGetAffineTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
114: {
115: /*
116: 2
117: |\
118: | \
119: | \
120: | \
121: | C \
122: | \
123: | \
124: 2---1---1
125: |\ D / \
126: | 2 0 \
127: |A \ / B \
128: 0---0-------1
129: */
130: static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
131: static PetscReal tri_J[] = {0.5, 0.0,
132: 0.0, 0.5,
134: 0.5, 0.0,
135: 0.0, 0.5,
137: 0.5, 0.0,
138: 0.0, 0.5,
140: 0.0, -0.5,
141: 0.5, 0.5};
142: static PetscReal tri_invJ[] = {2.0, 0.0,
143: 0.0, 2.0,
145: 2.0, 0.0,
146: 0.0, 2.0,
148: 2.0, 0.0,
149: 0.0, 2.0,
151: 2.0, 2.0,
152: -2.0, 0.0};
153: /*
154: 3---------2---------2
155: | | |
156: | D 2 C |
157: | | |
158: 3----3----0----1----1
159: | | |
160: | A 0 B |
161: | | |
162: 0---------0---------1
163: */
164: static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
165: static PetscReal quad_J[] = {0.5, 0.0,
166: 0.0, 0.5,
168: 0.5, 0.0,
169: 0.0, 0.5,
171: 0.5, 0.0,
172: 0.0, 0.5,
174: 0.5, 0.0,
175: 0.0, 0.5};
176: static PetscReal quad_invJ[] = {2.0, 0.0,
177: 0.0, 2.0,
179: 2.0, 0.0,
180: 0.0, 2.0,
182: 2.0, 0.0,
183: 0.0, 2.0,
185: 2.0, 0.0,
186: 0.0, 2.0};
187: /*
188: Bottom (viewed from top) Top
189: 1---------2---------2 7---------2---------6
190: | | | | | |
191: | B 2 C | | H 2 G |
192: | | | | | |
193: 3----3----0----1----1 3----3----0----1----1
194: | | | | | |
195: | A 0 D | | E 0 F |
196: | | | | | |
197: 0---------0---------3 4---------0---------5
198: */
199: 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,
200: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
201: static PetscReal hex_J[] = {0.5, 0.0, 0.0,
202: 0.0, 0.5, 0.0,
203: 0.0, 0.0, 0.5,
205: 0.5, 0.0, 0.0,
206: 0.0, 0.5, 0.0,
207: 0.0, 0.0, 0.5,
209: 0.5, 0.0, 0.0,
210: 0.0, 0.5, 0.0,
211: 0.0, 0.0, 0.5,
213: 0.5, 0.0, 0.0,
214: 0.0, 0.5, 0.0,
215: 0.0, 0.0, 0.5,
217: 0.5, 0.0, 0.0,
218: 0.0, 0.5, 0.0,
219: 0.0, 0.0, 0.5,
221: 0.5, 0.0, 0.0,
222: 0.0, 0.5, 0.0,
223: 0.0, 0.0, 0.5,
225: 0.5, 0.0, 0.0,
226: 0.0, 0.5, 0.0,
227: 0.0, 0.0, 0.5,
229: 0.5, 0.0, 0.0,
230: 0.0, 0.5, 0.0,
231: 0.0, 0.0, 0.5};
232: static PetscReal hex_invJ[] = {2.0, 0.0, 0.0,
233: 0.0, 2.0, 0.0,
234: 0.0, 0.0, 2.0,
236: 2.0, 0.0, 0.0,
237: 0.0, 2.0, 0.0,
238: 0.0, 0.0, 2.0,
240: 2.0, 0.0, 0.0,
241: 0.0, 2.0, 0.0,
242: 0.0, 0.0, 2.0,
244: 2.0, 0.0, 0.0,
245: 0.0, 2.0, 0.0,
246: 0.0, 0.0, 2.0,
248: 2.0, 0.0, 0.0,
249: 0.0, 2.0, 0.0,
250: 0.0, 0.0, 2.0,
252: 2.0, 0.0, 0.0,
253: 0.0, 2.0, 0.0,
254: 0.0, 0.0, 2.0,
256: 2.0, 0.0, 0.0,
257: 0.0, 2.0, 0.0,
258: 0.0, 0.0, 2.0,
260: 2.0, 0.0, 0.0,
261: 0.0, 2.0, 0.0,
262: 0.0, 0.0, 2.0};
265: switch (ct) {
266: case DM_POLYTOPE_TRIANGLE: if (Nc) *Nc = 4; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; break;
267: case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break;
268: case DM_POLYTOPE_HEXAHEDRON: if (Nc) *Nc = 8; if (v0) *v0 = hex_v0; if (J) *J = hex_J; if (invJ) *invJ = hex_invJ; break;
269: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
270: }
271: return(0);
272: }
274: #if 0
275: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
276: {
277: static PetscReal seg_v[] = {-1.0, 0.0, 1.0};
278: 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};
279: 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};
280: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
281: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
282: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0};
283: static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
284: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0,
285: -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0,
286: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0,
287: -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
288: -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0,
289: -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0,
290: -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
291: -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0};
294: switch (ct) {
295: case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break;
296: case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break;
297: case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break;
298: case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break;
299: case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break;
300: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
301: }
302: return(0);
303: }
305: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
306: {
307: static PetscInt seg_v[] = {0, 1, 1, 2};
308: static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5};
309: static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3};
310: static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9,
311: 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6};
312: 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,
313: 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};
316: if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
317: switch (ct) {
318: case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break;
319: case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break;
320: case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
321: case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break;
322: case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break;
323: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
324: }
325: return(0);
326: }
327: #endif
329: static PetscErrorCode DMPlexTransformView_Regular(DMPlexTransform tr, PetscViewer viewer)
330: {
331: PetscBool isascii;
337: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
338: if (isascii) {
339: const char *name;
341: PetscObjectGetName((PetscObject) tr, &name);
342: PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : "");
343: } else {
344: SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
345: }
346: return(0);
347: }
349: static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr)
350: {
352: return(0);
353: }
355: static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr)
356: {
357: DMPlexRefine_Regular *f = (DMPlexRefine_Regular *) tr->data;
358: PetscErrorCode ierr;
361: PetscFree(f);
362: return(0);
363: }
365: PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
366: {
367: static PetscInt seg_seg[] = {1, -1, 0, -1,
368: 0, 0, 1, 0};
369: static PetscInt tri_seg[] = {2, -1, 1, -1, 0, -1,
370: 1, -1, 0, -1, 2, -1,
371: 0, -1, 2, -1, 1, -1,
372: 0, 0, 1, 0, 2, 0,
373: 1, 0, 2, 0, 0, 0,
374: 2, 0, 0, 0, 1, 0};
375: static PetscInt tri_tri[] = {1, -3, 0, -3, 2, -3, 3, -2,
376: 0, -2, 2, -2, 1, -2, 3, -1,
377: 2, -1, 1, -1, 0, -1, 3, -3,
378: 0, 0, 1, 0, 2, 0, 3, 0,
379: 1, 1, 2, 1, 0, 1, 3, 1,
380: 2, 2, 0, 2, 1, 2, 3, 2};
381: static PetscInt quad_seg[] = {1, 0, 0, 0, 3, 0, 2, 0,
382: 0, 0, 3, 0, 2, 0, 1, 0,
383: 3, 0, 2, 0, 1, 0, 0, 0,
384: 2, 0, 1, 0, 0, 0, 3, 0,
385: 0, 0, 1, 0, 2, 0, 3, 0,
386: 1, 0, 2, 0, 3, 0, 0, 0,
387: 2, 0, 3, 0, 0, 0, 1, 0,
388: 3, 0, 0, 0, 1, 0, 2, 0};
389: static PetscInt quad_quad[] = {2, -4, 1, -4, 0, -4, 3, -4,
390: 1, -3, 0, -3, 3, -3, 2, -3,
391: 0, -2, 3, -2, 2, -2, 1, -2,
392: 3, -1, 2, -1, 1, -1, 0, -1,
393: 0, 0, 1, 0, 2, 0, 3, 0,
394: 1, 1, 2, 1, 3, 1, 0, 1,
395: 2, 2, 3, 2, 0, 2, 1, 2,
396: 3, 3, 0, 3, 1, 3, 2, 3};
397: static PetscInt tseg_seg[] = {0, -1,
398: 0, 0,
399: 0, 0,
400: 0, -1};
401: static PetscInt tseg_tseg[] = {1, -2, 0, -2,
402: 1, -1, 0, -1,
403: 0, 0, 1, 0,
404: 0, 1, 1, 1};
405: static PetscInt tet_seg[] = {0, -1,
406: 0, 0,
407: 0, 0,
408: 0, -1,
409: 0, 0,
410: 0, 0,
411: 0, 0,
412: 0, 0,
413: 0, 0,
414: 0, 0,
415: 0, 0,
416: 0, 0,
417: 0, 0,
418: 0, 0,
419: 0, 0,
420: 0, 0,
421: 0, 0,
422: 0, 0,
423: 0, -1,
424: 0, 0,
425: 0, 0,
426: 0, -1,
427: 0, 0,
428: 0, 0};
429: static PetscInt tet_tri[] = {2, -1, 3, -1, 1, -3, 0, -2, 6, 1, 7, -3, 5, 2, 4, -3,
430: 3, -1, 1, -1, 2, -3, 0, -1, 4, 0, 5, 0, 6, 0, 7, 0,
431: 1, -1, 2, -1, 3, -3, 0, -3, 4, 0, 5, 0, 6, 0, 7, 0,
432: 3, -2, 2, -3, 0, -1, 1, -1, 7, -3, 6, 1, 4, 2, 5, -3,
433: 2, -3, 0, -2, 3, -1, 1, -3, 4, 0, 5, 0, 6, 0, 7, 0,
434: 0, -2, 3, -2, 2, -2, 1, -2, 4, 0, 5, 0, 6, 0, 7, 0,
435: 0, -1, 1, -2, 3, -2, 2, -2, 7, 1, 6, -3, 5, -3, 4, 2,
436: 1, -2, 3, -3, 0, -3, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0,
437: 3, -3, 0, -1, 1, -1, 2, -3, 4, 0, 5, 0, 6, 0, 7, 0,
438: 1, -3, 0, -3, 2, -1, 3, -3, 6, -3, 7, 1, 4, -3, 5, 2,
439: 0, -3, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0,
440: 2, -2, 1, -3, 0, -2, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0,
441: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0,
442: 1, 0, 2, 2, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0,
443: 2, 2, 0, 0, 1, 1, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0,
444: 1, 2, 0, 1, 3, 1, 2, 2, 5, 0, 4, 0, 7, -1, 6, -1,
445: 0, 1, 3, 0, 1, 0, 2, 0, 4, 0, 5, 0, 6, 0, 7, 0,
446: 3, 0, 1, 2, 0, 2, 2, 1, 4, 0, 5, 0, 6, 0, 7, 0,
447: 2, 0, 3, 2, 0, 0, 1, 1, 4, -2, 5, -2, 7, 0, 6, 0,
448: 3, 2, 0, 2, 2, 1, 1, 2, 4, 0, 5, 0, 6, 0, 7, 0,
449: 0, 2, 2, 0, 3, 0, 1, 0, 4, 0, 5, 0, 6, 0, 7, 0,
450: 3, 1, 2, 1, 1, 2, 0, 2, 5, -2, 4, -2, 6, -1, 7, -1,
451: 2, 1, 1, 1, 3, 2, 0, 0, 4, 0, 5, 0, 6, 0, 7, 0,
452: 1, 1, 3, 1, 2, 2, 0, 1, 4, 0, 5, 0, 6, 0, 7, 0};
453: static PetscInt tet_tet[] = {2, -12, 3, -12, 1, -12, 0, -12, 6, -9, 7, -9, 5, -12, 4, -12,
454: 3, -11, 1, -11, 2, -11, 0, -11, 4, 0, 5, 0, 6, 0, 7, 0,
455: 1, -10, 2, -10, 3, -10, 0, -10, 4, 0, 5, 0, 6, 0, 7, 0,
456: 3, -9, 2, -9, 0, -9, 1, -9, 7, -9, 6, -9, 4, -12, 5, -12,
457: 2, -8, 0, -8, 3, -8, 1, -8, 4, 0, 5, 0, 6, 0, 7, 0,
458: 0, -7, 3, -7, 2, -7, 1, -7, 4, 0, 5, 0, 6, 0, 7, 0,
459: 0, -6, 1, -6, 3, -6, 2, -6, 4, -3, 5, -3, 7, -6, 6, -6,
460: 1, -5, 3, -5, 0, -5, 2, -5, 4, 0, 5, 0, 6, 0, 7, 0,
461: 3, -4, 0, -4, 1, -4, 2, -4, 4, 0, 5, 0, 6, 0, 7, 0,
462: 1, -3, 0, -3, 2, -3, 3, -3, 5, -3, 4, -3, 6, -6, 7, -6,
463: 0, -2, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0,
464: 2, -1, 1, -1, 0, -1, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0,
465: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0,
466: 1, 1, 2, 1, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0,
467: 2, 2, 0, 2, 1, 2, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0,
468: 1, 3, 0, 3, 3, 3, 2, 3, 5, 0, 4, 0, 7, 0, 6, 0,
469: 0, 4, 3, 4, 1, 4, 2, 4, 4, 0, 5, 0, 6, 0, 7, 0,
470: 3, 5, 1, 5, 0, 5, 2, 5, 4, 0, 5, 0, 6, 0, 7, 0,
471: 2, 6, 3, 6, 0, 6, 1, 6, 6, 6, 7, 6, 4, 6, 5, 6,
472: 3, 7, 0, 7, 2, 7, 1, 7, 4, 0, 5, 0, 6, 0, 7, 0,
473: 0, 8, 2, 8, 3, 8, 1, 8, 4, 0, 5, 0, 6, 0, 7, 0,
474: 3, 9, 2, 9, 1, 9, 0, 9, 7, 6, 6, 6, 5, 6, 4, 6,
475: 2, 10, 1, 10, 3, 10, 0, 10, 4, 0, 5, 0, 6, 0, 7, 0,
476: 1, 11, 3, 11, 2, 11, 0, 11, 4, 0, 5, 0, 6, 0, 7, 0};
477: static PetscInt hex_seg[] = {2, 0, 3, 0, 4, 0, 5, 0, 1, 0, 0, 0,
478: 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0,
479: 5, 0, 4, 0, 1, 0, 0, 0, 3, 0, 2, 0,
480: 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0,
481: 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0,
482: 4, 0, 5, 0, 1, 0, 0, 0, 2, 0, 3, 0,
483: 2, 0, 3, 0, 5, 0, 4, 0, 0, 0, 1, 0,
484: 5, 0, 4, 0, 0, 0, 1, 0, 2, 0, 3, 0,
485: 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0,
486: 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0,
487: 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0,
488: 2, 0, 3, 0, 0, 0, 1, 0, 4, 0, 5, 0,
489: 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0,
490: 1, 0, 0, 0, 5, 0, 4, 0, 2, 0, 3, 0,
491: 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0,
492: 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0,
493: 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0,
494: 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0,
495: 1, 0, 0, 0, 3, 0, 2, 0, 5, 0, 4, 0,
496: 2, 0, 3, 0, 1, 0, 0, 0, 5, 0, 4, 0,
497: 0, 0, 1, 0, 4, 0, 5, 0, 2, 0, 3, 0,
498: 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0,
499: 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0,
500: 0, 0, 1, 0, 2, 0, 3, 0, 5, 0, 4, 0,
501: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0,
502: 0, 0, 1, 0, 5, 0, 4, 0, 2, 0, 3, 0,
503: 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0,
504: 0, 0, 1, 0, 4, 0, 5, 0, 3, 0, 2, 0,
505: 2, 0, 3, 0, 1, 0, 0, 0, 4, 0, 5, 0,
506: 1, 0, 0, 0, 3, 0, 2, 0, 4, 0, 5, 0,
507: 3, 0, 2, 0, 0, 0, 1, 0, 4, 0, 5, 0,
508: 4, 0, 5, 0, 2, 0, 3, 0, 1, 0, 0, 0,
509: 1, 0, 0, 0, 2, 0, 3, 0, 5, 0, 4, 0,
510: 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0,
511: 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0,
512: 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0,
513: 2, 0, 3, 0, 0, 0, 1, 0, 5, 0, 4, 0,
514: 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0,
515: 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0,
516: 4, 0, 5, 0, 3, 0, 2, 0, 0, 0, 1, 0,
517: 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0,
518: 2, 0, 3, 0, 5, 0, 4, 0, 1, 0, 0, 0,
519: 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0,
520: 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0,
521: 3, 0, 2, 0, 4, 0, 5, 0, 1, 0, 0, 0,
522: 5, 0, 4, 0, 1, 0, 0, 0, 2, 0, 3, 0,
523: 4, 0, 5, 0, 0, 0, 1, 0, 2, 0, 3, 0,
524: 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0};
525: 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,
526: 8, -2, 9, -2, 10, -2, 11, -2, 3, -4, 0, -4, 1, -4, 2, -4, 7, 0, 4, 0, 5, 0, 6, 0,
527: 9, 1, 8, 1, 11, 1, 10, 1, 0, 3, 3, 3, 2, 3, 1, 3, 5, 2, 6, 2, 7, 2, 4, 2,
528: 6, 3, 5, 3, 4, 3, 7, 3, 10, -1, 9, -1, 8, -1, 11, -1, 2, -4, 3, -4, 0, -4, 1, -4,
529: 4, 1, 7, 1, 6, 1, 5, 1, 11, 2, 8, 2, 9, 2, 10, 2, 1, 3, 0, 3, 3, 3, 2, 3,
530: 10, -4, 11, -4, 8, -4, 9, -4, 2, 1, 1, 1, 0, 1, 3, 1, 6, -1, 5, -1, 4, -1, 7, -1,
531: 5, -4, 6, -4, 7, -4, 4, -4, 9, 0, 10, 0, 11, 0, 8, 0, 0, -2, 1, -2, 2, -2, 3, -2,
532: 11, 3, 10, 3, 9, 3, 8, 3, 1, -2, 2, -2, 3, -2, 0, -2, 4, -3, 7, -3, 6, -3, 5, -3,
533: 11, -1, 8, -1, 9, -1, 10, -1, 7, 3, 4, 3, 5, 3, 6, 3, 2, 2, 1, 2, 0, 2, 3, 2,
534: 10, 2, 9, 2, 8, 2, 11, 2, 5, 1, 6, 1, 7, 1, 4, 1, 1, -1, 2, -1, 3, -1, 0, -1,
535: 5, 2, 4, 2, 7, 2, 6, 2, 1, 2, 0, 2, 3, 2, 2, 2, 10, -4, 9, -4, 8, -4, 11, -4,
536: 4, -3, 5, -3, 6, -3, 7, -3, 0, -3, 1, -3, 2, -3, 3, -3, 8, -2, 11, -2, 10, -2, 9, -2,
537: 3, 1, 0, 1, 1, 1, 2, 1, 9, -4, 8, -4, 11, -4, 10, -4, 6, 3, 7, 3, 4, 3, 5, 3,
538: 1, 3, 2, 3, 3, 3, 0, 3, 10, 1, 11, 1, 8, 1, 9, 1, 5, -4, 4, -4, 7, -4, 6, -4,
539: 8, 0, 11, 0, 10, 0, 9, 0, 4, -4, 7, -4, 6, -4, 5, -4, 0, 0, 3, 0, 2, 0, 1, 0,
540: 0, 0, 1, 0, 2, 0, 3, 0, 5, -1, 4, -1, 7, -1, 6, -1, 9, -3, 8, -3, 11, -3, 10, -3,
541: 9, -3, 10, -3, 11, -3, 8, -3, 6, -2, 5, -2, 4, -2, 7, -2, 3, -3, 0, -3, 1, -3, 2, -3,
542: 7, 0, 6, 0, 5, 0, 4, 0, 2, -1, 3, -1, 0, -1, 1, -1, 11, 3, 8, 3, 9, 3, 10, 3,
543: 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 10, 2, 11, 2, 8, 2, 9, 2,
544: 6, -1, 7, -1, 4, -1, 5, -1, 3, 0, 2, 0, 1, 0, 0, 0, 9, 1, 10, 1, 11, 1, 8, 1,
545: 2, -4, 1, -4, 0, -4, 3, -4, 11, -2, 10, -2, 9, -2, 8, -2, 7, -2, 6, -2, 5, -2, 4, -2,
546: 1, -1, 0, -1, 3, -1, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 11, -1, 10, -1, 9, -1, 8, -1,
547: 0, -2, 3, -2, 2, -2, 1, -2, 8, 3, 9, 3, 10, 3, 11, 3, 4, 1, 5, 1, 6, 1, 7, 1,
548: 3, -3, 2, -3, 1, -3, 0, -3, 7, -3, 6, -3, 5, -3, 4, -3, 8, 0, 9, 0, 10, 0, 11, 0,
549: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9, 0, 10, 0, 11, 0,
550: 1, 3, 2, 3, 3, 3, 0, 3, 11, -2, 10, -2, 9, -2, 8, -2, 4, 1, 5, 1, 6, 1, 7, 1,
551: 2, 2, 3, 2, 0, 2, 1, 2, 7, -3, 6, -3, 5, -3, 4, -3, 11, -1, 10, -1, 9, -1, 8, -1,
552: 3, 1, 0, 1, 1, 1, 2, 1, 8, 3, 9, 3, 10, 3, 11, 3, 7, -2, 6, -2, 5, -2, 4, -2,
553: 5, 2, 4, 2, 7, 2, 6, 2, 0, -3, 1, -3, 2, -3, 3, -3, 9, 1, 10, 1, 11, 1, 8, 1,
554: 1, -1, 0, -1, 3, -1, 2, -1, 5, -1, 4, -1, 7, -1, 6, -1, 10, 2, 11, 2, 8, 2, 9, 2,
555: 4, -3, 5, -3, 6, -3, 7, -3, 1, 2, 0, 2, 3, 2, 2, 2, 11, 3, 8, 3, 9, 3, 10, 3,
556: 8, 0, 11, 0, 10, 0, 9, 0, 7, 3, 4, 3, 5, 3, 6, 3, 3, -3, 0, -3, 1, -3, 2, -3,
557: 3, -3, 2, -3, 1, -3, 0, -3, 6, 2, 7, 2, 4, 2, 5, 2, 9, -3, 8, -3, 11, -3, 10, -3,
558: 9, -3, 10, -3, 11, -3, 8, -3, 5, 1, 6, 1, 7, 1, 4, 1, 0, 0, 3, 0, 2, 0, 1, 0,
559: 0, -2, 3, -2, 2, -2, 1, -2, 9, -4, 8, -4, 11, -4, 10, -4, 5, -4, 4, -4, 7, -4, 6, -4,
560: 2, -4, 1, -4, 0, -4, 3, -4, 10, 1, 11, 1, 8, 1, 9, 1, 6, 3, 7, 3, 4, 3, 5, 3,
561: 7, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 8, -2, 11, -2, 10, -2, 9, -2,
562: 6, -1, 7, -1, 4, -1, 5, -1, 2, -1, 3, -1, 0, -1, 1, -1, 10, -4, 9, -4, 8, -4, 11, -4,
563: 11, -1, 8, -1, 9, -1, 10, -1, 4, -4, 7, -4, 6, -4, 5, -4, 1, -1, 2, -1, 3, -1, 0, -1,
564: 10, 2, 9, 2, 8, 2, 11, 2, 6, -2, 5, -2, 4, -2, 7, -2, 2, 2, 1, 2, 0, 2, 3, 2,
565: 8, -2, 9, -2, 10, -2, 11, -2, 0, 3, 3, 3, 2, 3, 1, 3, 4, -3, 7, -3, 6, -3, 5, -3,
566: 4, 1, 7, 1, 6, 1, 5, 1, 8, -3, 11, -3, 10, -3, 9, -3, 0, -2, 1, -2, 2, -2, 3, -2,
567: 9, 1, 8, 1, 11, 1, 10, 1, 3, -4, 0, -4, 1, -4, 2, -4, 6, -1, 5, -1, 4, -1, 7, -1,
568: 5, -4, 6, -4, 7, -4, 4, -4, 10, -1, 9, -1, 8, -1, 11, -1, 1, 3, 0, 3, 3, 3, 2, 3,
569: 7, -2, 4, -2, 5, -2, 6, -2, 11, 2, 8, 2, 9, 2, 10, 2, 2, -4, 3, -4, 0, -4, 1, -4,
570: 10, -4, 11, -4, 8, -4, 9, -4, 1, -2, 2, -2, 3, -2, 0, -2, 5, 2, 6, 2, 7, 2, 4, 2,
571: 11, 3, 10, 3, 9, 3, 8, 3, 2, 1, 1, 1, 0, 1, 3, 1, 7, 0, 4, 0, 5, 0, 6, 0,
572: 6, 3, 5, 3, 4, 3, 7, 3, 9, 0, 10, 0, 11, 0, 8, 0, 3, 1, 2, 1, 1, 1, 0, 1};
573: static PetscInt hex_hex[] = {3, -24, 0, -24, 4, -24, 5, -24, 2, -24, 6, -24, 7, -24, 1, -24,
574: 3, -23, 5, -23, 6, -23, 2, -23, 0, -23, 1, -23, 7, -23, 4, -23,
575: 4, -22, 0, -22, 1, -22, 7, -22, 5, -22, 6, -22, 2, -22, 3, -22,
576: 6, -21, 7, -21, 1, -21, 2, -21, 5, -21, 3, -21, 0, -21, 4, -21,
577: 1, -20, 2, -20, 6, -20, 7, -20, 0, -20, 4, -20, 5, -20, 3, -20,
578: 6, -19, 2, -19, 3, -19, 5, -19, 7, -19, 4, -19, 0, -19, 1, -19,
579: 4, -18, 5, -18, 3, -18, 0, -18, 7, -18, 1, -18, 2, -18, 6, -18,
580: 1, -17, 7, -17, 4, -17, 0, -17, 2, -17, 3, -17, 5, -17, 6, -17,
581: 2, -16, 3, -16, 5, -16, 6, -16, 1, -16, 7, -16, 4, -16, 0, -16,
582: 7, -15, 4, -15, 0, -15, 1, -15, 6, -15, 2, -15, 3, -15, 5, -15,
583: 7, -14, 1, -14, 2, -14, 6, -14, 4, -14, 5, -14, 3, -14, 0, -14,
584: 0, -13, 4, -13, 5, -13, 3, -13, 1, -13, 2, -13, 6, -13, 7, -13,
585: 5, -12, 4, -12, 7, -12, 6, -12, 3, -12, 2, -12, 1, -12, 0, -12,
586: 7, -11, 6, -11, 5, -11, 4, -11, 1, -11, 0, -11, 3, -11, 2, -11,
587: 0, -10, 1, -10, 7, -10, 4, -10, 3, -10, 5, -10, 6, -10, 2, -10,
588: 4, -9, 7, -9, 6, -9, 5, -9, 0, -9, 3, -9, 2, -9, 1, -9,
589: 5, -8, 6, -8, 2, -8, 3, -8, 4, -8, 0, -8, 1, -8, 7, -8,
590: 2, -7, 6, -7, 7, -7, 1, -7, 3, -7, 0, -7, 4, -7, 5, -7,
591: 6, -6, 5, -6, 4, -6, 7, -6, 2, -6, 1, -6, 0, -6, 3, -6,
592: 5, -5, 3, -5, 0, -5, 4, -5, 6, -5, 7, -5, 1, -5, 2, -5,
593: 2, -4, 1, -4, 0, -4, 3, -4, 6, -4, 5, -4, 4, -4, 7, -4,
594: 1, -3, 0, -3, 3, -3, 2, -3, 7, -3, 6, -3, 5, -3, 4, -3,
595: 0, -2, 3, -2, 2, -2, 1, -2, 4, -2, 7, -2, 6, -2, 5, -2,
596: 3, -1, 2, -1, 1, -1, 0, -1, 5, -1, 4, -1, 7, -1, 6, -1,
597: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0,
598: 1, 1, 2, 1, 3, 1, 0, 1, 7, 1, 4, 1, 5, 1, 6, 1,
599: 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2,
600: 3, 3, 0, 3, 1, 3, 2, 3, 5, 3, 6, 3, 7, 3, 4, 3,
601: 4, 4, 0, 4, 3, 4, 5, 4, 7, 4, 6, 4, 2, 4, 1, 4,
602: 7, 5, 4, 5, 5, 5, 6, 5, 1, 5, 2, 5, 3, 5, 0, 5,
603: 1, 6, 7, 6, 6, 6, 2, 6, 0, 6, 3, 6, 5, 6, 4, 6,
604: 3, 7, 2, 7, 6, 7, 5, 7, 0, 7, 4, 7, 7, 7, 1, 7,
605: 5, 8, 6, 8, 7, 8, 4, 8, 3, 8, 0, 8, 1, 8, 2, 8,
606: 4, 9, 7, 9, 1, 9, 0, 9, 5, 9, 3, 9, 2, 9, 6, 9,
607: 4, 10, 5, 10, 6, 10, 7, 10, 0, 10, 1, 10, 2, 10, 3, 10,
608: 6, 11, 7, 11, 4, 11, 5, 11, 2, 11, 3, 11, 0, 11, 1, 11,
609: 3, 12, 5, 12, 4, 12, 0, 12, 2, 12, 1, 12, 7, 12, 6, 12,
610: 6, 13, 2, 13, 1, 13, 7, 13, 5, 13, 4, 13, 0, 13, 3, 13,
611: 1, 14, 0, 14, 4, 14, 7, 14, 2, 14, 6, 14, 5, 14, 3, 14,
612: 6, 15, 5, 15, 3, 15, 2, 15, 7, 15, 1, 15, 0, 15, 4, 15,
613: 0, 16, 4, 16, 7, 16, 1, 16, 3, 16, 2, 16, 6, 16, 5, 16,
614: 0, 17, 3, 17, 5, 17, 4, 17, 1, 17, 7, 17, 6, 17, 2, 17,
615: 5, 18, 3, 18, 2, 18, 6, 18, 4, 18, 7, 18, 1, 18, 0, 18,
616: 7, 19, 6, 19, 2, 19, 1, 19, 4, 19, 0, 19, 3, 19, 5, 19,
617: 2, 20, 1, 20, 7, 20, 6, 20, 3, 20, 5, 20, 4, 20, 0, 20,
618: 7, 21, 1, 21, 0, 21, 4, 21, 6, 21, 5, 21, 3, 21, 2, 21,
619: 2, 22, 6, 22, 5, 22, 3, 22, 1, 22, 0, 22, 4, 22, 7, 22,
620: 5, 23, 4, 23, 0, 23, 3, 23, 6, 23, 2, 23, 1, 23, 7, 23};
621: static PetscInt trip_seg[] = {1, 0, 2, 0, 0, 0,
622: 2, 0, 0, 0, 1, 0,
623: 0, 0, 1, 0, 2, 0,
624: 0, -1, 2, -1, 1, -1,
625: 1, -1, 0, -1, 2, -1,
626: 2, -1, 1, -1, 0, -1,
627: 0, 0, 1, 0, 2, 0,
628: 2, 0, 0, 0, 1, 0,
629: 1, 0, 2, 0, 0, 0,
630: 2, -1, 1, -1, 0, -1,
631: 1, -1, 0, -1, 2, -1,
632: 0, -1, 2, -1, 1, -1};
633: static PetscInt trip_tri[] = {1, 1, 2, 1, 0, 1, 3, 1,
634: 2, 2, 0, 2, 1, 2, 3, 2,
635: 0, 0, 1, 0, 2, 0, 3, 0,
636: 2, -1, 1, -1, 0, -1, 3, -3,
637: 0, -2, 2, -2, 1, -2, 3, -1,
638: 1, -3, 0, -3, 2, -3, 3, -2,
639: 0, 0, 1, 0, 2, 0, 3, 0,
640: 2, 2, 0, 2, 1, 2, 3, 2,
641: 1, 1, 2, 1, 0, 1, 3, 1,
642: 1, -3, 0, -3, 2, -3, 3, -2,
643: 0, -2, 2, -2, 1, -2, 3, -1,
644: 2, -1, 1, -1, 0, -1, 3, -3};
645: static PetscInt trip_quad[] = {4, -1, 5, -1, 3, -1, 1, -1, 2, -1, 0, -1,
646: 5, -1, 3, -1, 4, -1, 2, -1, 0, -1, 1, -1,
647: 3, -1, 4, -1, 5, -1, 0, -1, 1, -1, 2, -1,
648: 0, -3, 2, -3, 1, -3, 3, -3, 5, -3, 4, -3,
649: 1, -3, 0, -3, 2, -3, 4, -3, 3, -3, 5, -3,
650: 2, -3, 1, -3, 0, -3, 5, -3, 4, -3, 3, -3,
651: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0,
652: 2, 0, 0, 0, 1, 0, 5, 0, 3, 0, 4, 0,
653: 1, 0, 2, 0, 0, 0, 4, 0, 5, 0, 3, 0,
654: 5, 2, 4, 2, 3, 2, 2, 2, 1, 2, 0, 2,
655: 4, 2, 3, 2, 5, 2, 1, 2, 0, 2, 2, 2,
656: 3, 2, 5, 2, 4, 2, 0, 2, 2, 2, 1, 2};
657: static PetscInt trip_trip[] = {5, -6, 6, -6, 4, -6, 7, -6, 1, -6, 2, -6, 0, -6, 3, -6,
658: 6, -5, 4, -5, 5, -5, 7, -5, 2, -5, 0, -5, 1, -5, 3, -5,
659: 4, -4, 5, -4, 6, -4, 7, -4, 0, -4, 1, -4, 2, -4, 3, -4,
660: 2, -3, 1, -3, 0, -3, 3, -1, 6, -3, 5, -3, 4, -3, 7, -1,
661: 0, -2, 2, -2, 1, -2, 3, -3, 4, -2, 6, -2, 5, -2, 7, -3,
662: 1, -1, 0, -1, 2, -1, 3, -2, 5, -1, 4, -1, 6, -1, 7, -2,
663: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0,
664: 2, 1, 0, 1, 1, 1, 3, 1, 6, 1, 4, 1, 5, 1, 7, 1,
665: 1, 2, 2, 2, 0, 2, 3, 2, 5, 2, 6, 2, 4, 2, 7, 2,
666: 5, 3, 4, 3, 6, 3, 7, 4, 1, 3, 0, 3, 2, 3, 3, 4,
667: 4, 4, 6, 4, 5, 4, 7, 5, 0, 4, 2, 4, 1, 4, 3, 5,
668: 6, 5, 5, 5, 4, 5, 7, 3, 2, 5, 1, 5, 0, 5, 3, 3};
669: static PetscInt ttri_tseg[] = {2, -2, 1, -2, 0, -2,
670: 1, -2, 0, -2, 2, -2,
671: 0, -2, 2, -2, 1, -2,
672: 2, -1, 1, -1, 0, -1,
673: 1, -1, 0, -1, 2, -1,
674: 0, -1, 2, -1, 1, -1,
675: 0, 0, 1, 0, 2, 0,
676: 1, 0, 2, 0, 0, 0,
677: 2, 0, 0, 0, 1, 0,
678: 0, 1, 1, 1, 2, 1,
679: 1, 1, 2, 1, 0, 1,
680: 2, 1, 0, 1, 1, 1};
681: static PetscInt ttri_ttri[] = {1, -6, 0, -6, 2, -6, 3, -5,
682: 0, -5, 2, -5, 1, -5, 3, -4,
683: 2, -4, 1, -4, 0, -4, 3, -6,
684: 1, -3, 0, -3, 2, -3, 3, -2,
685: 0, -2, 2, -2, 1, -2, 3, -1,
686: 2, -1, 1, -1, 0, -1, 3, -3,
687: 0, 0, 1, 0, 2, 0, 3, 0,
688: 1, 1, 2, 1, 0, 1, 3, 1,
689: 2, 2, 0, 2, 1, 2, 3, 2,
690: 0, 3, 1, 3, 2, 3, 3, 3,
691: 1, 4, 2, 4, 0, 4, 3, 4,
692: 2, 5, 0, 5, 1, 5, 3, 5};
693: static PetscInt tquad_tvert[] = {0, -1,
694: 0, -1,
695: 0, -1,
696: 0, -1,
697: 0, 0,
698: 0, 0,
699: 0, 0,
700: 0, 0,
701: 0, 0,
702: 0, 0,
703: 0, 0,
704: 0, 0,
705: 0, -1,
706: 0, -1,
707: 0, -1,
708: 0, -1};
709: static PetscInt tquad_tseg[] = {1, 1, 0, 1, 3, 1, 2, 1,
710: 0, 1, 3, 1, 2, 1, 1, 1,
711: 3, 1, 2, 1, 1, 1, 0, 1,
712: 2, 1, 1, 1, 0, 1, 3, 1,
713: 1, 0, 0, 0, 3, 0, 2, 0,
714: 0, 0, 3, 0, 2, 0, 1, 0,
715: 3, 0, 2, 0, 1, 0, 0, 0,
716: 2, 0, 1, 0, 0, 0, 3, 0,
717: 0, 0, 1, 0, 2, 0, 3, 0,
718: 1, 0, 2, 0, 3, 0, 0, 0,
719: 2, 0, 3, 0, 0, 0, 1, 0,
720: 3, 0, 0, 0, 1, 0, 2, 0,
721: 0, 1, 1, 1, 2, 1, 3, 1,
722: 1, 1, 2, 1, 3, 1, 0, 1,
723: 2, 1, 3, 1, 0, 1, 1, 1,
724: 3, 1, 0, 1, 1, 1, 2, 1};
725: static PetscInt tquad_tquad[] = {2, -8, 1, -8, 0, -8, 3, -8,
726: 1, -7, 0, -7, 3, -7, 2, -7,
727: 0, -6, 3, -6, 2, -6, 1, -6,
728: 3, -5, 2, -5, 1, -5, 0, -5,
729: 2, -4, 1, -4, 0, -4, 3, -4,
730: 1, -3, 0, -3, 3, -3, 2, -3,
731: 0, -2, 3, -2, 2, -2, 1, -2,
732: 3, -1, 2, -1, 1, -1, 0, -1,
733: 0, 0, 1, 0, 2, 0, 3, 0,
734: 1, 1, 2, 1, 3, 1, 0, 1,
735: 2, 2, 3, 2, 0, 2, 1, 2,
736: 3, 3, 0, 3, 1, 3, 2, 3,
737: 0, 4, 1, 4, 2, 4, 3, 4,
738: 1, 5, 2, 5, 3, 5, 0, 5,
739: 2, 6, 3, 6, 0, 6, 1, 6,
740: 3, 7, 0, 7, 1, 7, 2, 7};
743: *rnew = r; *onew = o;
744: if (!so) return(0);
745: switch (sct) {
746: case DM_POLYTOPE_POINT: break;
747: case DM_POLYTOPE_POINT_PRISM_TENSOR:
748: *onew = so < 0 ? -(o+1) : o;
749: break;
750: case DM_POLYTOPE_SEGMENT:
751: switch (tct) {
752: case DM_POLYTOPE_POINT: break;
753: case DM_POLYTOPE_SEGMENT:
754: *rnew = seg_seg[(so+1)*4 + r*2];
755: *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so+1)*4 + r*2 + 1]);
756: break;
757: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
758: }
759: break;
760: case DM_POLYTOPE_TRIANGLE:
761: switch (tct) {
762: case DM_POLYTOPE_SEGMENT:
763: *rnew = tri_seg[(so+3)*6 + r*2];
764: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so+3)*6 + r*2 + 1]);
765: break;
766: case DM_POLYTOPE_TRIANGLE:
767: *rnew = tri_tri[(so+3)*8 + r*2];
768: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so+3)*8 + r*2 + 1]);
769: break;
770: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
771: }
772: break;
773: case DM_POLYTOPE_QUADRILATERAL:
774: switch (tct) {
775: case DM_POLYTOPE_POINT: break;
776: case DM_POLYTOPE_SEGMENT:
777: *rnew = quad_seg[(so+4)*8 + r*2];
778: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so+4)*8 + r*2 + 1]);
779: break;
780: case DM_POLYTOPE_QUADRILATERAL:
781: *rnew = quad_quad[(so+4)*8 + r*2];
782: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so+4)*8 + r*2 + 1]);
783: break;
784: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
785: }
786: break;
787: case DM_POLYTOPE_SEG_PRISM_TENSOR:
788: switch (tct) {
789: case DM_POLYTOPE_POINT_PRISM_TENSOR:
790: *rnew = tseg_seg[(so+2)*2 + r*2];
791: *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so+2)*2 + r*2 + 1]);
792: break;
793: case DM_POLYTOPE_SEG_PRISM_TENSOR:
794: *rnew = tseg_tseg[(so+2)*4 + r*2];
795: *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so+2)*4 + r*2 + 1]);
796: break;
797: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
798: }
799: break;
800: case DM_POLYTOPE_TETRAHEDRON:
801: switch (tct) {
802: case DM_POLYTOPE_SEGMENT:
803: *rnew = tet_seg[(so+12)*2 + r*2];
804: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so+12)*2 + r*2 + 1]);
805: break;
806: case DM_POLYTOPE_TRIANGLE:
807: *rnew = tet_tri[(so+12)*16 + r*2];
808: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so+12)*16 + r*2 + 1]);
809: break;
810: case DM_POLYTOPE_TETRAHEDRON:
811: *rnew = tet_tet[(so+12)*16 + r*2];
812: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so+12)*16 + r*2 + 1]);
813: break;
814: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
815: }
816: break;
817: case DM_POLYTOPE_HEXAHEDRON:
818: switch (tct) {
819: case DM_POLYTOPE_POINT: break;
820: case DM_POLYTOPE_SEGMENT:
821: *rnew = hex_seg[(so+24)*12 + r*2];
822: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so+24)*12 + r*2 + 1]);
823: break;
824: case DM_POLYTOPE_QUADRILATERAL:
825: *rnew = hex_quad[(so+24)*24 + r*2];
826: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so+24)*24 + r*2 + 1]);
827: break;
828: case DM_POLYTOPE_HEXAHEDRON:
829: *rnew = hex_hex[(so+24)*16 + r*2];
830: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so+24)*16 + r*2 + 1]);
831: break;
832: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
833: }
834: break;
835: case DM_POLYTOPE_TRI_PRISM:
836: switch (tct) {
837: case DM_POLYTOPE_SEGMENT:
838: *rnew = trip_seg[(so+6)*6 + r*2];
839: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so+6)*6 + r*2 + 1]);
840: break;
841: case DM_POLYTOPE_TRIANGLE:
842: *rnew = trip_tri[(so+6)*8 + r*2];
843: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so+6)*8 + r*2 + 1]);
844: break;
845: case DM_POLYTOPE_QUADRILATERAL:
846: *rnew = trip_quad[(so+6)*12 + r*2];
847: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so+6)*12 + r*2 + 1]);
848: break;
849: case DM_POLYTOPE_TRI_PRISM:
850: *rnew = trip_trip[(so+6)*16 + r*2];
851: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so+6)*16 + r*2 + 1]);
852: break;
853: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
854: }
855: break;
856: case DM_POLYTOPE_TRI_PRISM_TENSOR:
857: switch (tct) {
858: case DM_POLYTOPE_SEG_PRISM_TENSOR:
859: *rnew = ttri_tseg[(so+6)*6 + r*2];
860: *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so+6)*6 + r*2 + 1]);
861: break;
862: case DM_POLYTOPE_TRI_PRISM_TENSOR:
863: *rnew = ttri_ttri[(so+6)*8 + r*2];
864: *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so+6)*8 + r*2 + 1]);
865: break;
866: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
867: }
868: break;
869: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
870: switch (tct) {
871: case DM_POLYTOPE_POINT_PRISM_TENSOR:
872: *rnew = tquad_tvert[(so+8)*2 + r*2];
873: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so+8)*2 + r*2 + 1]);
874: break;
875: case DM_POLYTOPE_SEG_PRISM_TENSOR:
876: *rnew = tquad_tseg[(so+8)*8 + r*2];
877: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so+8)*8 + r*2 + 1]);
878: break;
879: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
880: *rnew = tquad_tquad[(so+8)*8 + r*2];
881: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so+8)*8 + r*2 + 1]);
882: break;
883: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
884: }
885: break;
886: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
887: }
888: return(0);
889: }
891: PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
892: {
893: /* All vertices remain in the refined mesh */
894: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
895: static PetscInt vertexS[] = {1};
896: static PetscInt vertexC[] = {0};
897: static PetscInt vertexO[] = {0};
898: /* Split all edges with a new vertex, making two new 2 edges
899: 0--0--0--1--1
900: */
901: static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
902: static PetscInt segS[] = {1, 2};
903: 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};
904: static PetscInt segO[] = { 0, 0, 0, 0};
905: /* Do not split tensor edges */
906: static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
907: static PetscInt tvertS[] = {1};
908: static PetscInt tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
909: static PetscInt tvertO[] = { 0, 0};
910: /* Add 3 edges inside every triangle, making 4 new triangles.
911: 2
912: |\
913: | \
914: | \
915: 0 1
916: | C \
917: | \
918: | \
919: 2---1---1
920: |\ D / \
921: 1 2 0 0
922: |A \ / B \
923: 0-0-0---1---1
924: */
925: static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
926: static PetscInt triS[] = {3, 4};
927: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
928: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
929: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
930: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
931: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
932: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
933: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
934: static PetscInt triO[] = {0, 0,
935: 0, 0,
936: 0, 0,
937: 0, -1, 0,
938: 0, 0, -1,
939: -1, 0, 0,
940: 0, 0, 0};
941: /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
942: 3----1----2----0----2
943: | | |
944: 0 D 2 C 1
945: | | |
946: 3----3----0----1----1
947: | | |
948: 1 A 0 B 0
949: | | |
950: 0----0----0----1----1
951: */
952: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
953: static PetscInt quadS[] = {1, 4, 4};
954: static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
955: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
956: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
957: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
958: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1,
959: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
960: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2,
961: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
962: static PetscInt quadO[] = {0, 0,
963: 0, 0,
964: 0, 0,
965: 0, 0,
966: 0, 0, -1, 0,
967: 0, 0, 0, -1,
968: -1, 0, 0, 0,
969: 0, -1, 0, 0};
970: /* Add 1 edge inside every tensor quad, making 2 new tensor quads
971: 2----2----1----3----3
972: | | |
973: | | |
974: | | |
975: 4 A 6 B 5
976: | | |
977: | | |
978: | | |
979: 0----0----0----1----1
980: */
981: static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
982: static PetscInt tsegS[] = {1, 2};
983: static PetscInt tsegC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
984: 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,
985: 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};
986: static PetscInt tsegO[] = {0, 0,
987: 0, 0, 0, 0,
988: 0, 0, 0, 0};
989: /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
990: 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
991: 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]
992: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
993: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
994: The first four tets just cut off the corners, using the replica notation for new vertices,
995: [v0, (e0, 0), (e2, 0), (e3, 0)]
996: [(e0, 0), v1, (e1, 0), (e4, 0)]
997: [(e2, 0), (e1, 0), v2, (e5, 0)]
998: [(e3, 0), (e4, 0), (e5, 0), v3 ]
999: The next four tets match a vertex to the newly created faces from cutting off those first tets.
1000: [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
1001: [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
1002: [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
1003: [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
1004: 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
1005: [(e2, 0), (e0, 0), (e3, 0)]
1006: [(e0, 0), (e1, 0), (e4, 0)]
1007: [(e2, 0), (e5, 0), (e1, 0)]
1008: [(e3, 0), (e4, 0), (e5, 0)]
1009: The next four, from the second group of tets, are
1010: [(e2, 0), (e0, 0), (e5, 0)]
1011: [(e4, 0), (e0, 0), (e5, 0)]
1012: [(e0, 0), (e1, 0), (e5, 0)]
1013: [(e5, 0), (e3, 0), (e0, 0)]
1014: I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
1015: */
1016: static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
1017: static PetscInt tetS[] = {1, 8, 8};
1018: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0,
1019: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1020: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1,
1021: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1022: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1023: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1024: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1,
1025: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1026: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1027: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0,
1028: DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1,
1029: DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0,
1030: DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2,
1031: DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7,
1032: DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6,
1033: DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3,
1034: DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
1035: static PetscInt tetO[] = {0, 0,
1036: 0, 0, 0,
1037: 0, 0, 0,
1038: 0, 0, 0,
1039: 0, 0, 0,
1040: 0, 0, -1,
1041: 0, 0, -1,
1042: 0, -1, -1,
1043: 0, -1, 0,
1044: 0, 0, 0, 0,
1045: 0, 0, 0, 0,
1046: 0, 0, 0, 0,
1047: 0, 0, 0, 0,
1048: -2, 0, 0, -1,
1049: -1, 1, 0, 0,
1050: -1, -1, -3, 2,
1051: -1, 0, -1, 1};
1052: /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
1053: 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
1054: [v0, v1, v2, v3] f0 bottom
1055: [v4, v5, v6, v7] f1 top
1056: [v0, v3, v5, v4] f2 front
1057: [v2, v1, v7, v6] f3 back
1058: [v3, v2, v6, v5] f4 right
1059: [v0, v4, v7, v1] f5 left
1060: The eight hexes are divided into four on the bottom, and four on the top,
1061: [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1062: [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1063: [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1064: [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1065: [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)]
1066: [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)]
1067: [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)]
1068: [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7]
1069: 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,
1070: [(e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1071: [(f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1072: [(c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1073: [(f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1074: and on the x-z plane,
1075: [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
1076: [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
1077: [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
1078: [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
1079: and on the y-z plane,
1080: [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
1081: [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
1082: [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
1083: [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
1084: */
1085: static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1086: static PetscInt hexS[] = {1, 6, 12, 8};
1087: static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1088: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1089: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1090: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1091: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1092: DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0,
1093: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0,
1094: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2,
1095: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3,
1096: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2,
1097: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0,
1098: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1,
1099: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1100: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1101: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3,
1102: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1103: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1104: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1105: 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,
1106: 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,
1107: 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,
1108: 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,
1109: 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,
1110: 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,
1111: 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,
1112: 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};
1113: static PetscInt hexO[] = {0, 0,
1114: 0, 0,
1115: 0, 0,
1116: 0, 0,
1117: 0, 0,
1118: 0, 0,
1119: 0, 0, -1, -1,
1120: 0, -1, -1, 0,
1121: -1, -1, 0, 0,
1122: -1, 0, 0, -1,
1123: -1, 0, 0, -1,
1124: -1, -1, 0, 0,
1125: 0, -1, -1, 0,
1126: 0, 0, -1, -1,
1127: 0, 0, -1, -1,
1128: -1, 0, 0, -1,
1129: -1, -1, 0, 0,
1130: 0, -1, -1, 0,
1131: 0, 0, 0, 0, -2, 0,
1132: 0, 0, -3, 0, -2, 0,
1133: 0, 0, -3, 0, 0, 0,
1134: 0, 0, 0, 0, 0, 0,
1135: -2, 0, 0, 0, -2, 0,
1136: -2, 0, 0, 0, 0, 0,
1137: -2, 0, -3, 0, 0, 0,
1138: -2, 0, -3, 0, -2, 0};
1139: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1140: static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
1141: static PetscInt tripS[] = {3, 4, 6, 8};
1142: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0,
1143: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0,
1144: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0,
1145: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1146: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0,
1147: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1148: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1149: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1150: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1151: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1152: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1153: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1154: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1155: 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,
1156: 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,
1157: 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,
1158: 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,
1159: 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,
1160: 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,
1161: 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,
1162: 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};
1163: static PetscInt tripO[] = {0, 0,
1164: 0, 0,
1165: 0, 0,
1166: 0, -1, -1,
1167: -1, 0, -1,
1168: -1, -1, 0,
1169: 0, 0, 0,
1170: -1, 0, -1, -1,
1171: -1, 0, -1, -1,
1172: -1, 0, -1, -1,
1173: 0, -1, -1, 0,
1174: 0, -1, -1, 0,
1175: 0, -1, -1, 0,
1176: 0, 0, 0, -3, 0,
1177: 0, 0, 0, 0, -3,
1178: 0, 0, -3, 0, 0,
1179: 2, 0, 0, 0, 0,
1180: -2, 0, 0, -3, 0,
1181: -2, 0, 0, 0, -3,
1182: -2, 0, -3, 0, 0,
1183: -2, 0, 0, 0, 0};
1184: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
1185: 2
1186: |\
1187: | \
1188: | \
1189: 0---1
1191: 2
1193: 0 1
1195: 2
1196: |\
1197: | \
1198: | \
1199: 0---1
1200: */
1201: static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
1202: static PetscInt ttriS[] = {3, 4};
1203: 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,
1204: 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,
1205: 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,
1206: 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,
1207: 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,
1208: 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,
1209: 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};
1210: static PetscInt ttriO[] = {0, 0, 0, 0,
1211: 0, 0, 0, 0,
1212: 0, 0, 0, 0,
1213: 0, 0, 0, -1, 0,
1214: 0, 0, 0, 0, -1,
1215: 0, 0, -1, 0, 0,
1216: 0, 0, 0, 0, 0};
1217: /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
1218: static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1219: static PetscInt tquadS[] = {1, 4, 4};
1220: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1221: 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,
1222: 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,
1223: 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,
1224: 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,
1225: 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,
1226: 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,
1227: 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,
1228: 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};
1229: static PetscInt tquadO[] = {0, 0,
1230: 0, 0, 0, 0,
1231: 0, 0, 0, 0,
1232: 0, 0, 0, 0,
1233: 0, 0, 0, 0,
1234: 0, 0, 0, 0, -1, 0,
1235: 0, 0, 0, 0, 0, -1,
1236: 0, 0, -1, 0, 0, 0,
1237: 0, 0, 0, -1, 0, 0};
1238: /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */
1239: static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID};
1240: static PetscInt tpyrS[] = {4, 12, 1, 4, 6};
1241: static PetscInt tpyrC[] = {DM_POLYTOPE_POINT, 2, 1, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0,
1242: DM_POLYTOPE_POINT, 2, 2, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0,
1243: DM_POLYTOPE_POINT, 2, 3, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0,
1244: DM_POLYTOPE_POINT, 2, 4, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0,
1245: /* These four triangle face out of the bottom pyramid */
1246: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0,
1247: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1,
1248: DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1249: DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1250: /* These eight triangles face out of the corner pyramids */
1251: DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 2,
1252: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1253: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1254: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1255: DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
1256: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1,
1257: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 2,
1258: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 3,
1259: /* This quad faces out of the bottom pyramid */
1260: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1261: /* The bottom face of each tet is on the triangular face */
1262: DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_TRIANGLE, 0, 8, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 0,
1263: DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 9, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 1,
1264: DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 10, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2,
1265: DM_POLYTOPE_TRIANGLE, 1, 4, 3, DM_POLYTOPE_TRIANGLE, 0, 11, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3,
1266: /* The front face of all pyramids is toward the front */
1267: DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 11, DM_POLYTOPE_TRIANGLE, 1, 4, 1,
1268: DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 8,
1269: DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 9, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 6,
1270: DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 10, DM_POLYTOPE_TRIANGLE, 1, 3, 1, DM_POLYTOPE_TRIANGLE, 1, 4, 0,
1271: DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, DM_POLYTOPE_TRIANGLE, 1, 4, 2,
1272: DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 1,
1273: };
1274: static PetscInt tpyrO[] = {0, 0,
1275: 0, 0,
1276: 0, 0,
1277: 0, 0,
1278: 0, 0, -1,
1279: 0, 0, -1,
1280: 0, 0, -1,
1281: 0, 0, -1,
1282: 0, -1, 0,
1283: 0, -1, 0,
1284: 0, -1, 0,
1285: 0, -1, 0,
1286: -1, 0, 0,
1287: -1, 0, 0,
1288: -1, 0, 0,
1289: -1, 0, 0,
1290: -1, -1, -1, -1,
1291: 0, -3, -2, -3,
1292: 0, -3, -2, -3,
1293: 0, -3, -2, -3,
1294: 0, -3, -2, -3,
1295: 0, 0, 0, 0, 0,
1296: 0, 0, 0, 0, 0,
1297: 0, 0, 0, 0, 0,
1298: 0, 0, 0, 0, 0,
1299: -2, 0, 0, 0, 0,
1300: 1, 0, 0, 0, 0};
1303: if (rt) *rt = 0;
1304: switch (source) {
1305: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
1306: case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = segT; *size = segS; *cone = segC; *ornt = segO; break;
1307: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tvertT; *size = tvertS; *cone = tvertC; *ornt = tvertO; break;
1308: case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1309: case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
1310: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tsegT; *size = tsegS; *cone = tsegC; *ornt = tsegO; break;
1311: case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1312: case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
1313: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1314: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttriT; *size = ttriS; *cone = ttriC; *ornt = ttriO; break;
1315: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1316: case DM_POLYTOPE_PYRAMID: *Nt = 5; *target = tpyrT; *size = tpyrS; *cone = tpyrC; *ornt = tpyrO; break;
1317: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1318: }
1319: return(0);
1320: }
1322: static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr)
1323: {
1325: tr->ops->view = DMPlexTransformView_Regular;
1326: tr->ops->setup = DMPlexTransformSetUp_Regular;
1327: tr->ops->destroy = DMPlexTransformDestroy_Regular;
1328: tr->ops->celltransform = DMPlexTransformCellRefine_Regular;
1329: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular;
1330: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
1331: return(0);
1332: }
1334: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr)
1335: {
1336: DMPlexRefine_Regular *f;
1337: PetscErrorCode ierr;
1341: PetscNewLog(tr, &f);
1342: tr->data = f;
1344: DMPlexTransformInitialize_Regular(tr);
1345: return(0);
1346: }