Actual source code: plexrefregular.c
1: #include <petsc/private/dmplextransformimpl.h>
3: /* Regular Refinement of Hybrid Meshes
5: We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
6: to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
7: transformations, such as changing from one type of cell to another, as simplex to hex.
9: To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
10: and the types of the new cells.
12: We need the group multiplication table for group actions from the dihedral group for each cell type.
14: We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
15: we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
16: (face, orient) pairs for each subcell.
17: */
19: /*@
20: DMPlexRefineRegularGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell
22: Input Parameters:
23: + cr - The DMPlexCellRefiner object
24: - ct - The cell type
26: Output Parameters:
27: + Nf - The number of faces for this cell type
28: . v0 - The translation of the first vertex for each face
29: . J - The Jacobian for each face (map from original cell to subcell)
30: . invJ - The inverse Jacobian for each face
31: - detJ - The determinant of the Jacobian for each face
33: Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
34: $ v0 + j x_face = x_cell
35: $ invj (x_cell - v0) = x_face
37: Level: developer
39: .seealso: `DMPlexCellRefinerGetAffineTransforms()`, `Create()`
40: @*/
41: PetscErrorCode DMPlexRefineRegularGetAffineFaceTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
42: {
43: /*
44: 2
45: |\
46: | \
47: | \
48: | \
49: | \
50: | \
51: | \
52: 2 1
53: | \
54: | \
55: | \
56: 0---0-------1
57: v0[Nf][dc]: 3 x 2
58: J[Nf][df][dc]: 3 x 1 x 2
59: invJ[Nf][dc][df]: 3 x 2 x 1
60: detJ[Nf]: 3
61: */
62: static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
63: static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
64: static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0};
65: static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0};
66: /*
67: 3---------2---------2
68: | |
69: | |
70: | |
71: 3 1
72: | |
73: | |
74: | |
75: 0---------0---------1
77: v0[Nf][dc]: 4 x 2
78: J[Nf][df][dc]: 4 x 1 x 2
79: invJ[Nf][dc][df]: 4 x 2 x 1
80: detJ[Nf]: 4
81: */
82: static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 - 1.0, 0.0};
83: static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
84: static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
85: static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0};
87: switch (ct) {
88: case DM_POLYTOPE_TRIANGLE:
89: if (Nf) *Nf = 3;
90: if (v0) *v0 = tri_v0;
91: if (J) *J = tri_J;
92: if (invJ) *invJ = tri_invJ;
93: if (detJ) *detJ = tri_detJ;
94: break;
95: case DM_POLYTOPE_QUADRILATERAL:
96: if (Nf) *Nf = 4;
97: if (v0) *v0 = quad_v0;
98: if (J) *J = quad_J;
99: if (invJ) *invJ = quad_invJ;
100: if (detJ) *detJ = quad_detJ;
101: break;
102: default:
103: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
104: }
105: return 0;
106: }
108: /*@
109: DMPlexRefineRegularGetAffineTransforms - Gets the affine map from the reference cell to each subcell
111: Input Parameters:
112: + cr - The DMPlexCellRefiner object
113: - ct - The cell type
115: Output Parameters:
116: + Nc - The number of subcells produced from this cell type
117: . v0 - The translation of the first vertex for each subcell
118: . J - The Jacobian for each subcell (map from reference cell to subcell)
119: - invJ - The inverse Jacobian for each subcell
121: Level: developer
123: .seealso: `DMPlexRefineRegularGetAffineFaceTransforms()`, `DMPLEXREFINEREGULAR`
124: @*/
125: PetscErrorCode DMPlexRefineRegularGetAffineTransforms(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
126: {
127: /*
128: 2
129: |\
130: | \
131: | \
132: | \
133: | C \
134: | \
135: | \
136: 2---1---1
137: |\ D / \
138: | 2 0 \
139: |A \ / B \
140: 0---0-------1
141: */
142: static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
143: static PetscReal tri_J[] = {0.5, 0.0, 0.0, 0.5,
145: 0.5, 0.0, 0.0, 0.5,
147: 0.5, 0.0, 0.0, 0.5,
149: 0.0, -0.5, 0.5, 0.5};
150: static PetscReal tri_invJ[] = {2.0, 0.0, 0.0, 2.0,
152: 2.0, 0.0, 0.0, 2.0,
154: 2.0, 0.0, 0.0, 2.0,
156: 2.0, 2.0, -2.0, 0.0};
157: /*
158: 3---------2---------2
159: | | |
160: | D 2 C |
161: | | |
162: 3----3----0----1----1
163: | | |
164: | A 0 B |
165: | | |
166: 0---------0---------1
167: */
168: static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
169: static PetscReal quad_J[] = {0.5, 0.0, 0.0, 0.5,
171: 0.5, 0.0, 0.0, 0.5,
173: 0.5, 0.0, 0.0, 0.5,
175: 0.5, 0.0, 0.0, 0.5};
176: static PetscReal quad_invJ[] = {2.0, 0.0, 0.0, 2.0,
178: 2.0, 0.0, 0.0, 2.0,
180: 2.0, 0.0, 0.0, 2.0,
182: 2.0, 0.0, 0.0, 2.0};
183: /*
184: Bottom (viewed from top) Top
185: 1---------2---------2 7---------2---------6
186: | | | | | |
187: | B 2 C | | H 2 G |
188: | | | | | |
189: 3----3----0----1----1 3----3----0----1----1
190: | | | | | |
191: | A 0 D | | E 0 F |
192: | | | | | |
193: 0---------0---------3 4---------0---------5
194: */
195: static PetscReal hex_v0[] = {-1.0, -1.0, -1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
196: static PetscReal hex_J[] = {0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
198: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
200: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
202: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
204: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
206: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
208: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5,
210: 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
211: static PetscReal hex_invJ[] = {2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
213: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
215: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
217: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
219: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
221: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
223: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0,
225: 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0};
227: switch (ct) {
228: case DM_POLYTOPE_TRIANGLE:
229: if (Nc) *Nc = 4;
230: if (v0) *v0 = tri_v0;
231: if (J) *J = tri_J;
232: if (invJ) *invJ = tri_invJ;
233: break;
234: case DM_POLYTOPE_QUADRILATERAL:
235: if (Nc) *Nc = 4;
236: if (v0) *v0 = quad_v0;
237: if (J) *J = quad_J;
238: if (invJ) *invJ = quad_invJ;
239: break;
240: case DM_POLYTOPE_HEXAHEDRON:
241: if (Nc) *Nc = 8;
242: if (v0) *v0 = hex_v0;
243: if (J) *J = hex_J;
244: if (invJ) *invJ = hex_invJ;
245: break;
246: default:
247: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
248: }
249: return 0;
250: }
252: #if 0
253: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
254: {
255: static PetscReal seg_v[] = {-1.0, 0.0, 1.0};
256: static PetscReal tri_v[] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
257: static PetscReal quad_v[] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0};
258: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
259: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
260: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0};
261: static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
262: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0,
263: -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0,
264: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0,
265: -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
266: -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0,
267: -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0,
268: -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
269: -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0};
271: switch (ct) {
272: case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break;
273: case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break;
274: case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break;
275: case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break;
276: case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break;
277: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
278: }
279: return 0;
280: }
282: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
283: {
284: static PetscInt seg_v[] = {0, 1, 1, 2};
285: static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5};
286: static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3};
287: static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9,
288: 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6};
289: static PetscInt hex_v[] = {0, 3, 4, 1, 9, 10, 13, 12, 3, 6, 7, 4, 12, 13, 16, 15, 4, 7, 8, 5, 13, 14, 17, 16, 1, 4 , 5 , 2, 10, 11, 14, 13,
290: 9, 12, 13, 10, 18, 19, 22, 21, 10, 13, 14, 11, 19, 20, 23, 22, 13, 16, 17, 14, 22, 23, 26, 25, 12, 15, 16, 13, 21, 22, 25, 24};
293: switch (ct) {
294: case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break;
295: case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break;
296: case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
297: case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break;
298: case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break;
299: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
300: }
301: return 0;
302: }
303: #endif
305: static PetscErrorCode DMPlexTransformView_Regular(DMPlexTransform tr, PetscViewer viewer)
306: {
307: PetscBool isascii;
311: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
312: if (isascii) {
313: const char *name;
315: PetscObjectGetName((PetscObject)tr, &name);
316: PetscViewerASCIIPrintf(viewer, "Regular refinement %s\n", name ? name : "");
317: } else {
318: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
319: }
320: return 0;
321: }
323: static PetscErrorCode DMPlexTransformSetUp_Regular(DMPlexTransform tr)
324: {
325: return 0;
326: }
328: static PetscErrorCode DMPlexTransformDestroy_Regular(DMPlexTransform tr)
329: {
330: DMPlexRefine_Regular *f = (DMPlexRefine_Regular *)tr->data;
332: PetscFree(f);
333: return 0;
334: }
336: PetscErrorCode DMPlexTransformGetSubcellOrientation_Regular(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
337: {
338: static PetscInt seg_seg[] = {1, -1, 0, -1, 0, 0, 1, 0};
339: static PetscInt tri_seg[] = {2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0};
340: static PetscInt tri_tri[] = {1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2};
341: static PetscInt quad_seg[] = {1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0, 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 1, 0, 2, 0};
342: static PetscInt quad_quad[] = {2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0, -3, 3, -3, 2, -3, 0, -2, 3, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 3, 1, 0, 1, 2, 2, 3, 2, 0, 2, 1, 2, 3, 3, 0, 3, 1, 3, 2, 3};
343: static PetscInt tseg_seg[] = {0, -1, 0, 0, 0, 0, 0, -1};
344: static PetscInt tseg_tseg[] = {1, -2, 0, -2, 1, -1, 0, -1, 0, 0, 1, 0, 0, 1, 1, 1};
345: static PetscInt tet_seg[] = {0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0};
346: static PetscInt tet_tri[] = {2, -1, 3, -1, 1, -3, 0, -2, 6, 1, 7, -3, 5, 2, 4, -3, 3, -1, 1, -1, 2, -3, 0, -1, 4, 0, 5, 0, 6, 0, 7, 0, 1, -1, 2, -1, 3, -3, 0, -3, 4, 0, 5, 0, 6, 0, 7, 0, 3, -2, 2, -3, 0, -1, 1, -1, 7, -3, 6, 1, 4, 2, 5, -3,
347: 2, -3, 0, -2, 3, -1, 1, -3, 4, 0, 5, 0, 6, 0, 7, 0, 0, -2, 3, -2, 2, -2, 1, -2, 4, 0, 5, 0, 6, 0, 7, 0, 0, -1, 1, -2, 3, -2, 2, -2, 7, 1, 6, -3, 5, -3, 4, 2, 1, -2, 3, -3, 0, -3, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0,
348: 3, -3, 0, -1, 1, -1, 2, -3, 4, 0, 5, 0, 6, 0, 7, 0, 1, -3, 0, -3, 2, -1, 3, -3, 6, -3, 7, 1, 4, -3, 5, 2, 0, -3, 2, -2, 1, -2, 3, -2, 4, 0, 5, 0, 6, 0, 7, 0, 2, -2, 1, -3, 0, -2, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0,
349: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 0, 2, 2, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0, 2, 2, 0, 0, 1, 1, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 1, 2, 0, 1, 3, 1, 2, 2, 5, 0, 4, 0, 7, -1, 6, -1,
350: 0, 1, 3, 0, 1, 0, 2, 0, 4, 0, 5, 0, 6, 0, 7, 0, 3, 0, 1, 2, 0, 2, 2, 1, 4, 0, 5, 0, 6, 0, 7, 0, 2, 0, 3, 2, 0, 0, 1, 1, 4, -2, 5, -2, 7, 0, 6, 0, 3, 2, 0, 2, 2, 1, 1, 2, 4, 0, 5, 0, 6, 0, 7, 0,
351: 0, 2, 2, 0, 3, 0, 1, 0, 4, 0, 5, 0, 6, 0, 7, 0, 3, 1, 2, 1, 1, 2, 0, 2, 5, -2, 4, -2, 6, -1, 7, -1, 2, 1, 1, 1, 3, 2, 0, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 3, 1, 2, 2, 0, 1, 4, 0, 5, 0, 6, 0, 7, 0};
352: static PetscInt tet_tet[] = {2, -12, 3, -12, 1, -12, 0, -12, 6, -9, 7, -9, 5, -12, 4, -12, 3, -11, 1, -11, 2, -11, 0, -11, 4, 0, 5, 0, 6, 0, 7, 0, 1, -10, 2, -10, 3, -10, 0, -10, 4, 0, 5, 0, 6, 0, 7, 0, 3, -9, 2, -9, 0, -9, 1, -9,
353: 7, -9, 6, -9, 4, -12, 5, -12, 2, -8, 0, -8, 3, -8, 1, -8, 4, 0, 5, 0, 6, 0, 7, 0, 0, -7, 3, -7, 2, -7, 1, -7, 4, 0, 5, 0, 6, 0, 7, 0, 0, -6, 1, -6, 3, -6, 2, -6, 4, -3, 5, -3, 7, -6, 6, -6,
354: 1, -5, 3, -5, 0, -5, 2, -5, 4, 0, 5, 0, 6, 0, 7, 0, 3, -4, 0, -4, 1, -4, 2, -4, 4, 0, 5, 0, 6, 0, 7, 0, 1, -3, 0, -3, 2, -3, 3, -3, 5, -3, 4, -3, 6, -6, 7, -6, 0, -2, 2, -2, 1, -2, 3, -2,
355: 4, 0, 5, 0, 6, 0, 7, 0, 2, -1, 1, -1, 0, -1, 3, -1, 4, 0, 5, 0, 6, 0, 7, 0, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 2, 1, 0, 1, 3, 1, 4, 0, 5, 0, 6, 0, 7, 0,
356: 2, 2, 0, 2, 1, 2, 3, 2, 4, 0, 5, 0, 6, 0, 7, 0, 1, 3, 0, 3, 3, 3, 2, 3, 5, 0, 4, 0, 7, 0, 6, 0, 0, 4, 3, 4, 1, 4, 2, 4, 4, 0, 5, 0, 6, 0, 7, 0, 3, 5, 1, 5, 0, 5, 2, 5,
357: 4, 0, 5, 0, 6, 0, 7, 0, 2, 6, 3, 6, 0, 6, 1, 6, 6, 6, 7, 6, 4, 6, 5, 6, 3, 7, 0, 7, 2, 7, 1, 7, 4, 0, 5, 0, 6, 0, 7, 0, 0, 8, 2, 8, 3, 8, 1, 8, 4, 0, 5, 0, 6, 0, 7, 0,
358: 3, 9, 2, 9, 1, 9, 0, 9, 7, 6, 6, 6, 5, 6, 4, 6, 2, 10, 1, 10, 3, 10, 0, 10, 4, 0, 5, 0, 6, 0, 7, 0, 1, 11, 3, 11, 2, 11, 0, 11, 4, 0, 5, 0, 6, 0, 7, 0};
359: static PetscInt hex_seg[] = {2, 0, 3, 0, 4, 0, 5, 0, 1, 0, 0, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 2, 0, 3, 0,
360: 2, 0, 3, 0, 5, 0, 4, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 4, 0, 5, 0,
361: 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 2, 0, 3, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0,
362: 1, 0, 0, 0, 3, 0, 2, 0, 5, 0, 4, 0, 2, 0, 3, 0, 1, 0, 0, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 0, 0, 1, 0, 2, 0, 3, 0, 5, 0, 4, 0,
363: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 4, 0, 5, 0, 3, 0, 2, 0, 2, 0, 3, 0, 1, 0, 0, 0, 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 4, 0, 5, 0,
364: 3, 0, 2, 0, 0, 0, 1, 0, 4, 0, 5, 0, 4, 0, 5, 0, 2, 0, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0, 5, 0, 4, 0, 5, 0, 4, 0, 2, 0, 3, 0, 0, 0, 1, 0, 1, 0, 0, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 2, 0, 3, 0,
365: 2, 0, 3, 0, 0, 0, 1, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 5, 0, 4, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 4, 0, 5, 0, 3, 0, 2, 0, 0, 0, 1, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 2, 0, 3, 0, 5, 0, 4, 0, 1, 0, 0, 0,
366: 4, 0, 5, 0, 1, 0, 0, 0, 3, 0, 2, 0, 3, 0, 2, 0, 5, 0, 4, 0, 0, 0, 1, 0, 3, 0, 2, 0, 4, 0, 5, 0, 1, 0, 0, 0, 5, 0, 4, 0, 1, 0, 0, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0, 2, 0, 3, 0, 2, 0, 3, 0, 4, 0, 5, 0, 0, 0, 1, 0};
367: static PetscInt hex_quad[] = {7, -2, 4, -2, 5, -2, 6, -2, 8, -3, 11, -3, 10, -3, 9, -3, 3, 1, 2, 1, 1, 1, 0, 1, 8, -2, 9, -2, 10, -2, 11, -2, 3, -4, 0, -4, 1, -4, 2, -4, 7, 0, 4, 0, 5, 0, 6, 0, 9, 1, 8, 1, 11,
368: 1, 10, 1, 0, 3, 3, 3, 2, 3, 1, 3, 5, 2, 6, 2, 7, 2, 4, 2, 6, 3, 5, 3, 4, 3, 7, 3, 10, -1, 9, -1, 8, -1, 11, -1, 2, -4, 3, -4, 0, -4, 1, -4, 4, 1, 7, 1, 6, 1, 5, 1, 11, 2,
369: 8, 2, 9, 2, 10, 2, 1, 3, 0, 3, 3, 3, 2, 3, 10, -4, 11, -4, 8, -4, 9, -4, 2, 1, 1, 1, 0, 1, 3, 1, 6, -1, 5, -1, 4, -1, 7, -1, 5, -4, 6, -4, 7, -4, 4, -4, 9, 0, 10, 0, 11, 0, 8,
370: 0, 0, -2, 1, -2, 2, -2, 3, -2, 11, 3, 10, 3, 9, 3, 8, 3, 1, -2, 2, -2, 3, -2, 0, -2, 4, -3, 7, -3, 6, -3, 5, -3, 11, -1, 8, -1, 9, -1, 10, -1, 7, 3, 4, 3, 5, 3, 6, 3, 2, 2, 1, 2,
371: 0, 2, 3, 2, 10, 2, 9, 2, 8, 2, 11, 2, 5, 1, 6, 1, 7, 1, 4, 1, 1, -1, 2, -1, 3, -1, 0, -1, 5, 2, 4, 2, 7, 2, 6, 2, 1, 2, 0, 2, 3, 2, 2, 2, 10, -4, 9, -4, 8, -4, 11, -4, 4,
372: -3, 5, -3, 6, -3, 7, -3, 0, -3, 1, -3, 2, -3, 3, -3, 8, -2, 11, -2, 10, -2, 9, -2, 3, 1, 0, 1, 1, 1, 2, 1, 9, -4, 8, -4, 11, -4, 10, -4, 6, 3, 7, 3, 4, 3, 5, 3, 1, 3, 2, 3, 3, 3,
373: 0, 3, 10, 1, 11, 1, 8, 1, 9, 1, 5, -4, 4, -4, 7, -4, 6, -4, 8, 0, 11, 0, 10, 0, 9, 0, 4, -4, 7, -4, 6, -4, 5, -4, 0, 0, 3, 0, 2, 0, 1, 0, 0, 0, 1, 0, 2, 0, 3, 0, 5, -1, 4,
374: -1, 7, -1, 6, -1, 9, -3, 8, -3, 11, -3, 10, -3, 9, -3, 10, -3, 11, -3, 8, -3, 6, -2, 5, -2, 4, -2, 7, -2, 3, -3, 0, -3, 1, -3, 2, -3, 7, 0, 6, 0, 5, 0, 4, 0, 2, -1, 3, -1, 0, -1, 1, -1,
375: 11, 3, 8, 3, 9, 3, 10, 3, 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 10, 2, 11, 2, 8, 2, 9, 2, 6, -1, 7, -1, 4, -1, 5, -1, 3, 0, 2, 0, 1, 0, 0, 0, 9, 1, 10, 1, 11,
376: 1, 8, 1, 2, -4, 1, -4, 0, -4, 3, -4, 11, -2, 10, -2, 9, -2, 8, -2, 7, -2, 6, -2, 5, -2, 4, -2, 1, -1, 0, -1, 3, -1, 2, -1, 4, 0, 5, 0, 6, 0, 7, 0, 11, -1, 10, -1, 9, -1, 8, -1, 0, -2,
377: 3, -2, 2, -2, 1, -2, 8, 3, 9, 3, 10, 3, 11, 3, 4, 1, 5, 1, 6, 1, 7, 1, 3, -3, 2, -3, 1, -3, 0, -3, 7, -3, 6, -3, 5, -3, 4, -3, 8, 0, 9, 0, 10, 0, 11, 0, 0, 0, 1, 0, 2, 0, 3,
378: 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9, 0, 10, 0, 11, 0, 1, 3, 2, 3, 3, 3, 0, 3, 11, -2, 10, -2, 9, -2, 8, -2, 4, 1, 5, 1, 6, 1, 7, 1, 2, 2, 3, 2, 0, 2, 1, 2, 7, -3, 6, -3,
379: 5, -3, 4, -3, 11, -1, 10, -1, 9, -1, 8, -1, 3, 1, 0, 1, 1, 1, 2, 1, 8, 3, 9, 3, 10, 3, 11, 3, 7, -2, 6, -2, 5, -2, 4, -2, 5, 2, 4, 2, 7, 2, 6, 2, 0, -3, 1, -3, 2, -3, 3, -3, 9,
380: 1, 10, 1, 11, 1, 8, 1, 1, -1, 0, -1, 3, -1, 2, -1, 5, -1, 4, -1, 7, -1, 6, -1, 10, 2, 11, 2, 8, 2, 9, 2, 4, -3, 5, -3, 6, -3, 7, -3, 1, 2, 0, 2, 3, 2, 2, 2, 11, 3, 8, 3, 9, 3,
381: 10, 3, 8, 0, 11, 0, 10, 0, 9, 0, 7, 3, 4, 3, 5, 3, 6, 3, 3, -3, 0, -3, 1, -3, 2, -3, 3, -3, 2, -3, 1, -3, 0, -3, 6, 2, 7, 2, 4, 2, 5, 2, 9, -3, 8, -3, 11, -3, 10, -3, 9, -3, 10,
382: -3, 11, -3, 8, -3, 5, 1, 6, 1, 7, 1, 4, 1, 0, 0, 3, 0, 2, 0, 1, 0, 0, -2, 3, -2, 2, -2, 1, -2, 9, -4, 8, -4, 11, -4, 10, -4, 5, -4, 4, -4, 7, -4, 6, -4, 2, -4, 1, -4, 0, -4, 3, -4,
383: 10, 1, 11, 1, 8, 1, 9, 1, 6, 3, 7, 3, 4, 3, 5, 3, 7, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 1, 0, 0, 0, 8, -2, 11, -2, 10, -2, 9, -2, 6, -1, 7, -1, 4, -1, 5, -1, 2, -1, 3, -1, 0,
384: -1, 1, -1, 10, -4, 9, -4, 8, -4, 11, -4, 11, -1, 8, -1, 9, -1, 10, -1, 4, -4, 7, -4, 6, -4, 5, -4, 1, -1, 2, -1, 3, -1, 0, -1, 10, 2, 9, 2, 8, 2, 11, 2, 6, -2, 5, -2, 4, -2, 7, -2, 2, 2,
385: 1, 2, 0, 2, 3, 2, 8, -2, 9, -2, 10, -2, 11, -2, 0, 3, 3, 3, 2, 3, 1, 3, 4, -3, 7, -3, 6, -3, 5, -3, 4, 1, 7, 1, 6, 1, 5, 1, 8, -3, 11, -3, 10, -3, 9, -3, 0, -2, 1, -2, 2, -2, 3,
386: -2, 9, 1, 8, 1, 11, 1, 10, 1, 3, -4, 0, -4, 1, -4, 2, -4, 6, -1, 5, -1, 4, -1, 7, -1, 5, -4, 6, -4, 7, -4, 4, -4, 10, -1, 9, -1, 8, -1, 11, -1, 1, 3, 0, 3, 3, 3, 2, 3, 7, -2, 4, -2,
387: 5, -2, 6, -2, 11, 2, 8, 2, 9, 2, 10, 2, 2, -4, 3, -4, 0, -4, 1, -4, 10, -4, 11, -4, 8, -4, 9, -4, 1, -2, 2, -2, 3, -2, 0, -2, 5, 2, 6, 2, 7, 2, 4, 2, 11, 3, 10, 3, 9, 3, 8, 3, 2,
388: 1, 1, 1, 0, 1, 3, 1, 7, 0, 4, 0, 5, 0, 6, 0, 6, 3, 5, 3, 4, 3, 7, 3, 9, 0, 10, 0, 11, 0, 8, 0, 3, 1, 2, 1, 1, 1, 0, 1};
389: static PetscInt hex_hex[] = {3, -24, 0, -24, 4, -24, 5, -24, 2, -24, 6, -24, 7, -24, 1, -24, 3, -23, 5, -23, 6, -23, 2, -23, 0, -23, 1, -23, 7, -23, 4, -23, 4, -22, 0, -22, 1, -22, 7, -22, 5, -22, 6, -22, 2, -22, 3, -22, 6, -21, 7, -21,
390: 1, -21, 2, -21, 5, -21, 3, -21, 0, -21, 4, -21, 1, -20, 2, -20, 6, -20, 7, -20, 0, -20, 4, -20, 5, -20, 3, -20, 6, -19, 2, -19, 3, -19, 5, -19, 7, -19, 4, -19, 0, -19, 1, -19, 4, -18, 5, -18, 3, -18, 0, -18,
391: 7, -18, 1, -18, 2, -18, 6, -18, 1, -17, 7, -17, 4, -17, 0, -17, 2, -17, 3, -17, 5, -17, 6, -17, 2, -16, 3, -16, 5, -16, 6, -16, 1, -16, 7, -16, 4, -16, 0, -16, 7, -15, 4, -15, 0, -15, 1, -15, 6, -15, 2, -15,
392: 3, -15, 5, -15, 7, -14, 1, -14, 2, -14, 6, -14, 4, -14, 5, -14, 3, -14, 0, -14, 0, -13, 4, -13, 5, -13, 3, -13, 1, -13, 2, -13, 6, -13, 7, -13, 5, -12, 4, -12, 7, -12, 6, -12, 3, -12, 2, -12, 1, -12, 0, -12,
393: 7, -11, 6, -11, 5, -11, 4, -11, 1, -11, 0, -11, 3, -11, 2, -11, 0, -10, 1, -10, 7, -10, 4, -10, 3, -10, 5, -10, 6, -10, 2, -10, 4, -9, 7, -9, 6, -9, 5, -9, 0, -9, 3, -9, 2, -9, 1, -9, 5, -8, 6, -8,
394: 2, -8, 3, -8, 4, -8, 0, -8, 1, -8, 7, -8, 2, -7, 6, -7, 7, -7, 1, -7, 3, -7, 0, -7, 4, -7, 5, -7, 6, -6, 5, -6, 4, -6, 7, -6, 2, -6, 1, -6, 0, -6, 3, -6, 5, -5, 3, -5, 0, -5, 4, -5,
395: 6, -5, 7, -5, 1, -5, 2, -5, 2, -4, 1, -4, 0, -4, 3, -4, 6, -4, 5, -4, 4, -4, 7, -4, 1, -3, 0, -3, 3, -3, 2, -3, 7, -3, 6, -3, 5, -3, 4, -3, 0, -2, 3, -2, 2, -2, 1, -2, 4, -2, 7, -2,
396: 6, -2, 5, -2, 3, -1, 2, -1, 1, -1, 0, -1, 5, -1, 4, -1, 7, -1, 6, -1, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 1, 1, 2, 1, 3, 1, 0, 1, 7, 1, 4, 1, 5, 1, 6, 1,
397: 2, 2, 3, 2, 0, 2, 1, 2, 6, 2, 7, 2, 4, 2, 5, 2, 3, 3, 0, 3, 1, 3, 2, 3, 5, 3, 6, 3, 7, 3, 4, 3, 4, 4, 0, 4, 3, 4, 5, 4, 7, 4, 6, 4, 2, 4, 1, 4, 7, 5, 4, 5,
398: 5, 5, 6, 5, 1, 5, 2, 5, 3, 5, 0, 5, 1, 6, 7, 6, 6, 6, 2, 6, 0, 6, 3, 6, 5, 6, 4, 6, 3, 7, 2, 7, 6, 7, 5, 7, 0, 7, 4, 7, 7, 7, 1, 7, 5, 8, 6, 8, 7, 8, 4, 8,
399: 3, 8, 0, 8, 1, 8, 2, 8, 4, 9, 7, 9, 1, 9, 0, 9, 5, 9, 3, 9, 2, 9, 6, 9, 4, 10, 5, 10, 6, 10, 7, 10, 0, 10, 1, 10, 2, 10, 3, 10, 6, 11, 7, 11, 4, 11, 5, 11, 2, 11, 3, 11,
400: 0, 11, 1, 11, 3, 12, 5, 12, 4, 12, 0, 12, 2, 12, 1, 12, 7, 12, 6, 12, 6, 13, 2, 13, 1, 13, 7, 13, 5, 13, 4, 13, 0, 13, 3, 13, 1, 14, 0, 14, 4, 14, 7, 14, 2, 14, 6, 14, 5, 14, 3, 14,
401: 6, 15, 5, 15, 3, 15, 2, 15, 7, 15, 1, 15, 0, 15, 4, 15, 0, 16, 4, 16, 7, 16, 1, 16, 3, 16, 2, 16, 6, 16, 5, 16, 0, 17, 3, 17, 5, 17, 4, 17, 1, 17, 7, 17, 6, 17, 2, 17, 5, 18, 3, 18,
402: 2, 18, 6, 18, 4, 18, 7, 18, 1, 18, 0, 18, 7, 19, 6, 19, 2, 19, 1, 19, 4, 19, 0, 19, 3, 19, 5, 19, 2, 20, 1, 20, 7, 20, 6, 20, 3, 20, 5, 20, 4, 20, 0, 20, 7, 21, 1, 21, 0, 21, 4, 21,
403: 6, 21, 5, 21, 3, 21, 2, 21, 2, 22, 6, 22, 5, 22, 3, 22, 1, 22, 0, 22, 4, 22, 7, 22, 5, 23, 4, 23, 0, 23, 3, 23, 6, 23, 2, 23, 1, 23, 7, 23};
404: static PetscInt trip_seg[] = {1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 0, -1, 2, -1, 1, -1, 1, -1, 0, -1, 2, -1, 2, -1, 1, -1, 0, -1,
405: 0, 0, 1, 0, 2, 0, 2, 0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1};
406: static PetscInt trip_tri[] = {1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2, 0, 0, 1, 0, 2, 0, 3, 0, 2, -1, 1, -1, 0, -1, 3, -3, 0, -2, 2, -2, 1, -2, 3, -1, 1, -3, 0, -3, 2, -3, 3, -2,
407: 0, 0, 1, 0, 2, 0, 3, 0, 2, 2, 0, 2, 1, 2, 3, 2, 1, 1, 2, 1, 0, 1, 3, 1, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3};
408: static PetscInt trip_quad[] = {4, -1, 5, -1, 3, -1, 1, -1, 2, -1, 0, -1, 5, -1, 3, -1, 4, -1, 2, -1, 0, -1, 1, -1, 3, -1, 4, -1, 5, -1, 0, -1, 1, -1, 2, -1, 0, -3, 2, -3, 1, -3, 3, -3, 5, -3, 4, -3,
409: 1, -3, 0, -3, 2, -3, 4, -3, 3, -3, 5, -3, 2, -3, 1, -3, 0, -3, 5, -3, 4, -3, 3, -3, 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 2, 0, 0, 0, 1, 0, 5, 0, 3, 0, 4, 0,
410: 1, 0, 2, 0, 0, 0, 4, 0, 5, 0, 3, 0, 5, 2, 4, 2, 3, 2, 2, 2, 1, 2, 0, 2, 4, 2, 3, 2, 5, 2, 1, 2, 0, 2, 2, 2, 3, 2, 5, 2, 4, 2, 0, 2, 2, 2, 1, 2};
411: static PetscInt trip_trip[] = {5, -6, 6, -6, 4, -6, 7, -6, 1, -6, 2, -6, 0, -6, 3, -6, 6, -5, 4, -5, 5, -5, 7, -5, 2, -5, 0, -5, 1, -5, 3, -5, 4, -4, 5, -4, 6, -4, 7, -4, 0, -4, 1, -4, 2, -4, 3, -4,
412: 2, -3, 1, -3, 0, -3, 3, -1, 6, -3, 5, -3, 4, -3, 7, -1, 0, -2, 2, -2, 1, -2, 3, -3, 4, -2, 6, -2, 5, -2, 7, -3, 1, -1, 0, -1, 2, -1, 3, -2, 5, -1, 4, -1, 6, -1, 7, -2,
413: 0, 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 2, 1, 0, 1, 1, 1, 3, 1, 6, 1, 4, 1, 5, 1, 7, 1, 1, 2, 2, 2, 0, 2, 3, 2, 5, 2, 6, 2, 4, 2, 7, 2,
414: 5, 3, 4, 3, 6, 3, 7, 4, 1, 3, 0, 3, 2, 3, 3, 4, 4, 4, 6, 4, 5, 4, 7, 5, 0, 4, 2, 4, 1, 4, 3, 5, 6, 5, 5, 5, 4, 5, 7, 3, 2, 5, 1, 5, 0, 5, 3, 3};
415: static PetscInt ttri_tseg[] = {2, -2, 1, -2, 0, -2, 1, -2, 0, -2, 2, -2, 0, -2, 2, -2, 1, -2, 2, -1, 1, -1, 0, -1, 1, -1, 0, -1, 2, -1, 0, -1, 2, -1, 1, -1,
416: 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 1, 1, 1, 2, 1, 1, 1, 2, 1, 0, 1, 2, 1, 0, 1, 1, 1};
417: static PetscInt ttri_ttri[] = {1, -6, 0, -6, 2, -6, 3, -5, 0, -5, 2, -5, 1, -5, 3, -4, 2, -4, 1, -4, 0, -4, 3, -6, 1, -3, 0, -3, 2, -3, 3, -2, 0, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 3, -3,
418: 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 0, 1, 3, 1, 2, 2, 0, 2, 1, 2, 3, 2, 0, 3, 1, 3, 2, 3, 3, 3, 1, 4, 2, 4, 0, 4, 3, 4, 2, 5, 0, 5, 1, 5, 3, 5};
419: static PetscInt tquad_tvert[] = {0, -1, 0, -1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, -1, 0, -1};
420: static PetscInt tquad_tseg[] = {1, 1, 0, 1, 3, 1, 2, 1, 0, 1, 3, 1, 2, 1, 1, 1, 3, 1, 2, 1, 1, 1, 0, 1, 2, 1, 1, 1, 0, 1, 3, 1, 1, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 0, 2, 0, 1, 0, 3, 0, 2, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 0, 3, 0,
421: 0, 0, 1, 0, 2, 0, 3, 0, 1, 0, 2, 0, 3, 0, 0, 0, 2, 0, 3, 0, 0, 0, 1, 0, 3, 0, 0, 0, 1, 0, 2, 0, 0, 1, 1, 1, 2, 1, 3, 1, 1, 1, 2, 1, 3, 1, 0, 1, 2, 1, 3, 1, 0, 1, 1, 1, 3, 1, 0, 1, 1, 1, 2, 1};
422: static PetscInt tquad_tquad[] = {2, -8, 1, -8, 0, -8, 3, -8, 1, -7, 0, -7, 3, -7, 2, -7, 0, -6, 3, -6, 2, -6, 1, -6, 3, -5, 2, -5, 1, -5, 0, -5, 2, -4, 1, -4, 0, -4, 3, -4, 1, -3, 0,
423: -3, 3, -3, 2, -3, 0, -2, 3, -2, 2, -2, 1, -2, 3, -1, 2, -1, 1, -1, 0, -1, 0, 0, 1, 0, 2, 0, 3, 0, 1, 1, 2, 1, 3, 1, 0, 1, 2, 2, 3, 2, 0, 2,
424: 1, 2, 3, 3, 0, 3, 1, 3, 2, 3, 0, 4, 1, 4, 2, 4, 3, 4, 1, 5, 2, 5, 3, 5, 0, 5, 2, 6, 3, 6, 0, 6, 1, 6, 3, 7, 0, 7, 1, 7, 2, 7};
427: *rnew = r;
428: *onew = o;
429: if (!so) return 0;
430: switch (sct) {
431: case DM_POLYTOPE_POINT:
432: break;
433: case DM_POLYTOPE_POINT_PRISM_TENSOR:
434: *onew = so < 0 ? -(o + 1) : o;
435: break;
436: case DM_POLYTOPE_SEGMENT:
437: switch (tct) {
438: case DM_POLYTOPE_POINT:
439: break;
440: case DM_POLYTOPE_SEGMENT:
441: *rnew = seg_seg[(so + 1) * 4 + r * 2];
442: *onew = DMPolytopeTypeComposeOrientation(tct, o, seg_seg[(so + 1) * 4 + r * 2 + 1]);
443: break;
444: default:
445: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
446: }
447: break;
448: case DM_POLYTOPE_TRIANGLE:
449: switch (tct) {
450: case DM_POLYTOPE_SEGMENT:
451: *rnew = tri_seg[(so + 3) * 6 + r * 2];
452: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_seg[(so + 3) * 6 + r * 2 + 1]);
453: break;
454: case DM_POLYTOPE_TRIANGLE:
455: *rnew = tri_tri[(so + 3) * 8 + r * 2];
456: *onew = DMPolytopeTypeComposeOrientation(tct, o, tri_tri[(so + 3) * 8 + r * 2 + 1]);
457: break;
458: default:
459: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
460: }
461: break;
462: case DM_POLYTOPE_QUADRILATERAL:
463: switch (tct) {
464: case DM_POLYTOPE_POINT:
465: break;
466: case DM_POLYTOPE_SEGMENT:
467: *rnew = quad_seg[(so + 4) * 8 + r * 2];
468: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_seg[(so + 4) * 8 + r * 2 + 1]);
469: break;
470: case DM_POLYTOPE_QUADRILATERAL:
471: *rnew = quad_quad[(so + 4) * 8 + r * 2];
472: *onew = DMPolytopeTypeComposeOrientation(tct, o, quad_quad[(so + 4) * 8 + r * 2 + 1]);
473: break;
474: default:
475: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
476: }
477: break;
478: case DM_POLYTOPE_SEG_PRISM_TENSOR:
479: switch (tct) {
480: case DM_POLYTOPE_POINT_PRISM_TENSOR:
481: *rnew = tseg_seg[(so + 2) * 2 + r * 2];
482: *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_seg[(so + 2) * 2 + r * 2 + 1]);
483: break;
484: case DM_POLYTOPE_SEG_PRISM_TENSOR:
485: *rnew = tseg_tseg[(so + 2) * 4 + r * 2];
486: *onew = DMPolytopeTypeComposeOrientation(tct, o, tseg_tseg[(so + 2) * 4 + r * 2 + 1]);
487: break;
488: default:
489: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
490: }
491: break;
492: case DM_POLYTOPE_TETRAHEDRON:
493: switch (tct) {
494: case DM_POLYTOPE_SEGMENT:
495: *rnew = tet_seg[(so + 12) * 2 + r * 2];
496: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_seg[(so + 12) * 2 + r * 2 + 1]);
497: break;
498: case DM_POLYTOPE_TRIANGLE:
499: *rnew = tet_tri[(so + 12) * 16 + r * 2];
500: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tri[(so + 12) * 16 + r * 2 + 1]);
501: break;
502: case DM_POLYTOPE_TETRAHEDRON:
503: *rnew = tet_tet[(so + 12) * 16 + r * 2];
504: *onew = DMPolytopeTypeComposeOrientation(tct, o, tet_tet[(so + 12) * 16 + r * 2 + 1]);
505: break;
506: default:
507: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
508: }
509: break;
510: case DM_POLYTOPE_HEXAHEDRON:
511: switch (tct) {
512: case DM_POLYTOPE_POINT:
513: break;
514: case DM_POLYTOPE_SEGMENT:
515: *rnew = hex_seg[(so + 24) * 12 + r * 2];
516: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_seg[(so + 24) * 12 + r * 2 + 1]);
517: break;
518: case DM_POLYTOPE_QUADRILATERAL:
519: *rnew = hex_quad[(so + 24) * 24 + r * 2];
520: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_quad[(so + 24) * 24 + r * 2 + 1]);
521: break;
522: case DM_POLYTOPE_HEXAHEDRON:
523: *rnew = hex_hex[(so + 24) * 16 + r * 2];
524: *onew = DMPolytopeTypeComposeOrientation(tct, o, hex_hex[(so + 24) * 16 + r * 2 + 1]);
525: break;
526: default:
527: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
528: }
529: break;
530: case DM_POLYTOPE_TRI_PRISM:
531: switch (tct) {
532: case DM_POLYTOPE_SEGMENT:
533: *rnew = trip_seg[(so + 6) * 6 + r * 2];
534: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_seg[(so + 6) * 6 + r * 2 + 1]);
535: break;
536: case DM_POLYTOPE_TRIANGLE:
537: *rnew = trip_tri[(so + 6) * 8 + r * 2];
538: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_tri[(so + 6) * 8 + r * 2 + 1]);
539: break;
540: case DM_POLYTOPE_QUADRILATERAL:
541: *rnew = trip_quad[(so + 6) * 12 + r * 2];
542: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_quad[(so + 6) * 12 + r * 2 + 1]);
543: break;
544: case DM_POLYTOPE_TRI_PRISM:
545: *rnew = trip_trip[(so + 6) * 16 + r * 2];
546: *onew = DMPolytopeTypeComposeOrientation(tct, o, trip_trip[(so + 6) * 16 + r * 2 + 1]);
547: break;
548: default:
549: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
550: }
551: break;
552: case DM_POLYTOPE_TRI_PRISM_TENSOR:
553: switch (tct) {
554: case DM_POLYTOPE_SEG_PRISM_TENSOR:
555: *rnew = ttri_tseg[(so + 6) * 6 + r * 2];
556: *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_tseg[(so + 6) * 6 + r * 2 + 1]);
557: break;
558: case DM_POLYTOPE_TRI_PRISM_TENSOR:
559: *rnew = ttri_ttri[(so + 6) * 8 + r * 2];
560: *onew = DMPolytopeTypeComposeOrientation(tct, o, ttri_ttri[(so + 6) * 8 + r * 2 + 1]);
561: break;
562: default:
563: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
564: }
565: break;
566: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
567: switch (tct) {
568: case DM_POLYTOPE_POINT_PRISM_TENSOR:
569: *rnew = tquad_tvert[(so + 8) * 2 + r * 2];
570: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tvert[(so + 8) * 2 + r * 2 + 1]);
571: break;
572: case DM_POLYTOPE_SEG_PRISM_TENSOR:
573: *rnew = tquad_tseg[(so + 8) * 8 + r * 2];
574: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tseg[(so + 8) * 8 + r * 2 + 1]);
575: break;
576: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
577: *rnew = tquad_tquad[(so + 8) * 8 + r * 2];
578: *onew = DMPolytopeTypeComposeOrientation(tct, o, tquad_tquad[(so + 8) * 8 + r * 2 + 1]);
579: break;
580: default:
581: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
582: }
583: break;
584: default:
585: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
586: }
587: return 0;
588: }
590: PetscErrorCode DMPlexTransformCellRefine_Regular(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
591: {
592: /* All vertices remain in the refined mesh */
593: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
594: static PetscInt vertexS[] = {1};
595: static PetscInt vertexC[] = {0};
596: static PetscInt vertexO[] = {0};
597: /* Split all edges with a new vertex, making two new 2 edges
598: 0--0--0--1--1
599: */
600: static DMPolytopeType segT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
601: static PetscInt segS[] = {1, 2};
602: static PetscInt segC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
603: static PetscInt segO[] = {0, 0, 0, 0};
604: /* Do not split tensor edges */
605: static DMPolytopeType tvertT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
606: static PetscInt tvertS[] = {1};
607: static PetscInt tvertC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
608: static PetscInt tvertO[] = {0, 0};
609: /* Add 3 edges inside every triangle, making 4 new triangles.
610: 2
611: |\
612: | \
613: | \
614: 0 1
615: | C \
616: | \
617: | \
618: 2---1---1
619: |\ D / \
620: 1 2 0 0
621: |A \ / B \
622: 0-0-0---1---1
623: */
624: static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
625: static PetscInt triS[] = {3, 4};
626: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
627: static PetscInt triO[] = {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0};
628: /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
629: 3----1----2----0----2
630: | | |
631: 0 D 2 C 1
632: | | |
633: 3----3----0----1----1
634: | | |
635: 1 A 0 B 0
636: | | |
637: 0----0----0----1----1
638: */
639: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
640: static PetscInt quadS[] = {1, 4, 4};
641: static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
642: static PetscInt quadO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, -1, 0, 0};
643: /* Add 1 edge inside every tensor quad, making 2 new tensor quads
644: 2----2----1----3----3
645: | | |
646: | | |
647: | | |
648: 4 A 6 B 5
649: | | |
650: | | |
651: | | |
652: 0----0----0----1----1
653: */
654: static DMPolytopeType tsegT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
655: static PetscInt tsegS[] = {1, 2};
656: static PetscInt tsegC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
657: static PetscInt tsegO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
658: /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
659: The vertices of our reference tet are [(-1, -1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, 1)], which we call [v0, v1, v2, v3]. The first
660: three edges are [v0, v1], [v1, v2], [v2, v0] called e0, e1, and e2, and then three edges to the top point [v0, v3], [v1, v3], [v2, v3]
661: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
662: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
663: The first four tets just cut off the corners, using the replica notation for new vertices,
664: [v0, (e0, 0), (e2, 0), (e3, 0)]
665: [(e0, 0), v1, (e1, 0), (e4, 0)]
666: [(e2, 0), (e1, 0), v2, (e5, 0)]
667: [(e3, 0), (e4, 0), (e5, 0), v3 ]
668: The next four tets match a vertex to the newly created faces from cutting off those first tets.
669: [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
670: [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
671: [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
672: [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
673: We can see that a new edge is introduced in the cell [(e0, 0), (e5, 0)] which we call (-1, 0). The first four faces created are
674: [(e2, 0), (e0, 0), (e3, 0)]
675: [(e0, 0), (e1, 0), (e4, 0)]
676: [(e2, 0), (e5, 0), (e1, 0)]
677: [(e3, 0), (e4, 0), (e5, 0)]
678: The next four, from the second group of tets, are
679: [(e2, 0), (e0, 0), (e5, 0)]
680: [(e4, 0), (e0, 0), (e5, 0)]
681: [(e0, 0), (e1, 0), (e5, 0)]
682: [(e5, 0), (e3, 0), (e0, 0)]
683: I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
684: */
685: static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
686: static PetscInt tetS[] = {1, 8, 8};
687: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
688: static PetscInt tetO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -1, -1, 1, 0, 0, -1, -1, -3, 2, -1, 0, -1, 1};
689: /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
690: The vertices of our reference hex are (-1, -1, -1), (-1, 1, -1), (1, 1, -1), (1, -1, -1), (-1, -1, 1), (1, -1, 1), (1, 1, 1), (-1, 1, 1) which we call [v0, v1, v2, v3, v4, v5, v6, v7]. The fours edges around the bottom [v0, v1], [v1, v2], [v2, v3], [v3, v0] are [e0, e1, e2, e3], and likewise around the top [v4, v5], [v5, v6], [v6, v7], [v7, v4] are [e4, e5, e6, e7]. Finally [v0, v4], [v1, v7], [v2, v6], [v3, v5] are [e9, e10, e11, e8]. The faces of a hex, given in DMPlexGetRawFaces_Internal(), oriented with outward normals, are
691: [v0, v1, v2, v3] f0 bottom
692: [v4, v5, v6, v7] f1 top
693: [v0, v3, v5, v4] f2 front
694: [v2, v1, v7, v6] f3 back
695: [v3, v2, v6, v5] f4 right
696: [v0, v4, v7, v1] f5 left
697: The eight hexes are divided into four on the bottom, and four on the top,
698: [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)]
699: [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)]
700: [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)]
701: [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)]
702: [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)]
703: [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)]
704: [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)]
705: [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7]
706: The 6 internal edges will go from the faces to the central vertex. The 12 internal faces can be divided into groups of 4 by the plane on which they sit. First the faces on the x-y plane are,
707: [(e9, 0), (f2, 0), (c0, 0), (f5, 0)]
708: [(f5, 0), (c0, 0), (f3, 0), (e10, 0)]
709: [(c0, 0), (f4, 0), (e11, 0), (f3, 0)]
710: [(f2, 0), (e8, 0), (f4, 0), (c0, 0)]
711: and on the x-z plane,
712: [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
713: [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
714: [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
715: [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
716: and on the y-z plane,
717: [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
718: [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
719: [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
720: [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
721: */
722: static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
723: static PetscInt hexS[] = {1, 6, 12, 8};
724: static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 5, 3, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_QUADRILATERAL, 0, 11, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 1, 5, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_QUADRILATERAL, 0, 9, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 10, DM_POLYTOPE_QUADRILATERAL, 1, 5, 2};
725: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, -1, -1, -1, 0, 0, 0, -1, -1, 0, 0, 0, -1, -1, 0, 0, -1, -1, -1, 0, 0, -1, -1, -1,
726: 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, -3, 0, -2, 0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 0, -2, 0, -3, 0, 0, 0, -2, 0, -3, 0, -2, 0};
727: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
728: static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
729: static PetscInt tripS[] = {3, 4, 6, 8};
730: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5};
731: static PetscInt tripO[] = {0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, -1, -1, 0, 0, 0, 0, -1, 0, -1, -1, -1, 0, -1, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, -1, -1, 0, 0, -1, -1,
732: 0, 0, 0, 0, -3, 0, 0, 0, 0, 0, -3, 0, 0, -3, 0, 0, 2, 0, 0, 0, 0, -2, 0, 0, -3, 0, -2, 0, 0, 0, -3, -2, 0, -3, 0, 0, -2, 0, 0, 0, 0};
733: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
734: 2
735: |\
736: | \
737: | \
738: 0---1
740: 2
742: 0 1
744: 2
745: |\
746: | \
747: | \
748: 0---1
749: */
750: static DMPolytopeType ttriT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
751: static PetscInt ttriS[] = {3, 4};
752: static PetscInt ttriC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_TRIANGLE, 1, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2};
753: static PetscInt ttriO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0};
754: /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
755: static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
756: static PetscInt tquadS[] = {1, 4, 4};
757: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 5, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 1, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 3, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
758: static PetscInt tquadO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0};
759: /* Add 4 edges, 12 triangles, 1 quad, 4 tetrahedra, and 6 pyramids inside every pyramid. */
760: static DMPolytopeType tpyrT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TETRAHEDRON, DM_POLYTOPE_PYRAMID};
761: static PetscInt tpyrS[] = {4, 12, 1, 4, 6};
762: static PetscInt tpyrC[] = {
763: DM_POLYTOPE_POINT,
764: 2,
765: 1,
766: 1,
767: 0,
768: DM_POLYTOPE_POINT,
769: 1,
770: 0,
771: 0,
772: DM_POLYTOPE_POINT,
773: 2,
774: 2,
775: 1,
776: 0,
777: DM_POLYTOPE_POINT,
778: 1,
779: 0,
780: 0,
781: DM_POLYTOPE_POINT,
782: 2,
783: 3,
784: 1,
785: 0,
786: DM_POLYTOPE_POINT,
787: 1,
788: 0,
789: 0,
790: DM_POLYTOPE_POINT,
791: 2,
792: 4,
793: 1,
794: 0,
795: DM_POLYTOPE_POINT,
796: 1,
797: 0,
798: 0,
799: /* These four triangle face out of the bottom pyramid */
800: DM_POLYTOPE_SEGMENT,
801: 1,
802: 1,
803: 1,
804: DM_POLYTOPE_SEGMENT,
805: 0,
806: 3,
807: DM_POLYTOPE_SEGMENT,
808: 0,
809: 0,
810: DM_POLYTOPE_SEGMENT,
811: 1,
812: 2,
813: 1,
814: DM_POLYTOPE_SEGMENT,
815: 0,
816: 0,
817: DM_POLYTOPE_SEGMENT,
818: 0,
819: 1,
820: DM_POLYTOPE_SEGMENT,
821: 1,
822: 3,
823: 1,
824: DM_POLYTOPE_SEGMENT,
825: 0,
826: 1,
827: DM_POLYTOPE_SEGMENT,
828: 0,
829: 2,
830: DM_POLYTOPE_SEGMENT,
831: 1,
832: 4,
833: 1,
834: DM_POLYTOPE_SEGMENT,
835: 0,
836: 2,
837: DM_POLYTOPE_SEGMENT,
838: 0,
839: 3,
840: /* These eight triangles face out of the corner pyramids */
841: DM_POLYTOPE_SEGMENT,
842: 1,
843: 0,
844: 3,
845: DM_POLYTOPE_SEGMENT,
846: 0,
847: 3,
848: DM_POLYTOPE_SEGMENT,
849: 1,
850: 1,
851: 2,
852: DM_POLYTOPE_SEGMENT,
853: 1,
854: 0,
855: 2,
856: DM_POLYTOPE_SEGMENT,
857: 0,
858: 0,
859: DM_POLYTOPE_SEGMENT,
860: 1,
861: 2,
862: 2,
863: DM_POLYTOPE_SEGMENT,
864: 1,
865: 0,
866: 1,
867: DM_POLYTOPE_SEGMENT,
868: 0,
869: 1,
870: DM_POLYTOPE_SEGMENT,
871: 1,
872: 3,
873: 2,
874: DM_POLYTOPE_SEGMENT,
875: 1,
876: 0,
877: 0,
878: DM_POLYTOPE_SEGMENT,
879: 0,
880: 2,
881: DM_POLYTOPE_SEGMENT,
882: 1,
883: 4,
884: 2,
885: DM_POLYTOPE_SEGMENT,
886: 1,
887: 0,
888: 3,
889: DM_POLYTOPE_SEGMENT,
890: 1,
891: 1,
892: 0,
893: DM_POLYTOPE_SEGMENT,
894: 0,
895: 0,
896: DM_POLYTOPE_SEGMENT,
897: 1,
898: 0,
899: 2,
900: DM_POLYTOPE_SEGMENT,
901: 1,
902: 2,
903: 0,
904: DM_POLYTOPE_SEGMENT,
905: 0,
906: 1,
907: DM_POLYTOPE_SEGMENT,
908: 1,
909: 0,
910: 1,
911: DM_POLYTOPE_SEGMENT,
912: 1,
913: 3,
914: 0,
915: DM_POLYTOPE_SEGMENT,
916: 0,
917: 2,
918: DM_POLYTOPE_SEGMENT,
919: 1,
920: 0,
921: 0,
922: DM_POLYTOPE_SEGMENT,
923: 1,
924: 4,
925: 0,
926: DM_POLYTOPE_SEGMENT,
927: 0,
928: 3,
929: /* This quad faces out of the bottom pyramid */
930: DM_POLYTOPE_SEGMENT,
931: 1,
932: 1,
933: 1,
934: DM_POLYTOPE_SEGMENT,
935: 1,
936: 2,
937: 1,
938: DM_POLYTOPE_SEGMENT,
939: 1,
940: 3,
941: 1,
942: DM_POLYTOPE_SEGMENT,
943: 1,
944: 4,
945: 1,
946: /* The bottom face of each tet is on the triangular face */
947: DM_POLYTOPE_TRIANGLE,
948: 1,
949: 1,
950: 3,
951: DM_POLYTOPE_TRIANGLE,
952: 0,
953: 8,
954: DM_POLYTOPE_TRIANGLE,
955: 0,
956: 4,
957: DM_POLYTOPE_TRIANGLE,
958: 0,
959: 0,
960: DM_POLYTOPE_TRIANGLE,
961: 1,
962: 2,
963: 3,
964: DM_POLYTOPE_TRIANGLE,
965: 0,
966: 9,
967: DM_POLYTOPE_TRIANGLE,
968: 0,
969: 5,
970: DM_POLYTOPE_TRIANGLE,
971: 0,
972: 1,
973: DM_POLYTOPE_TRIANGLE,
974: 1,
975: 3,
976: 3,
977: DM_POLYTOPE_TRIANGLE,
978: 0,
979: 10,
980: DM_POLYTOPE_TRIANGLE,
981: 0,
982: 6,
983: DM_POLYTOPE_TRIANGLE,
984: 0,
985: 2,
986: DM_POLYTOPE_TRIANGLE,
987: 1,
988: 4,
989: 3,
990: DM_POLYTOPE_TRIANGLE,
991: 0,
992: 11,
993: DM_POLYTOPE_TRIANGLE,
994: 0,
995: 7,
996: DM_POLYTOPE_TRIANGLE,
997: 0,
998: 3,
999: /* The front face of all pyramids is toward the front */
1000: DM_POLYTOPE_QUADRILATERAL,
1001: 1,
1002: 0,
1003: 0,
1004: DM_POLYTOPE_TRIANGLE,
1005: 1,
1006: 1,
1007: 0,
1008: DM_POLYTOPE_TRIANGLE,
1009: 0,
1010: 4,
1011: DM_POLYTOPE_TRIANGLE,
1012: 0,
1013: 11,
1014: DM_POLYTOPE_TRIANGLE,
1015: 1,
1016: 4,
1017: 1,
1018: DM_POLYTOPE_QUADRILATERAL,
1019: 1,
1020: 0,
1021: 3,
1022: DM_POLYTOPE_TRIANGLE,
1023: 1,
1024: 1,
1025: 1,
1026: DM_POLYTOPE_TRIANGLE,
1027: 1,
1028: 2,
1029: 0,
1030: DM_POLYTOPE_TRIANGLE,
1031: 0,
1032: 5,
1033: DM_POLYTOPE_TRIANGLE,
1034: 0,
1035: 8,
1036: DM_POLYTOPE_QUADRILATERAL,
1037: 1,
1038: 0,
1039: 2,
1040: DM_POLYTOPE_TRIANGLE,
1041: 0,
1042: 9,
1043: DM_POLYTOPE_TRIANGLE,
1044: 1,
1045: 2,
1046: 1,
1047: DM_POLYTOPE_TRIANGLE,
1048: 1,
1049: 3,
1050: 0,
1051: DM_POLYTOPE_TRIANGLE,
1052: 0,
1053: 6,
1054: DM_POLYTOPE_QUADRILATERAL,
1055: 1,
1056: 0,
1057: 1,
1058: DM_POLYTOPE_TRIANGLE,
1059: 0,
1060: 7,
1061: DM_POLYTOPE_TRIANGLE,
1062: 0,
1063: 10,
1064: DM_POLYTOPE_TRIANGLE,
1065: 1,
1066: 3,
1067: 1,
1068: DM_POLYTOPE_TRIANGLE,
1069: 1,
1070: 4,
1071: 0,
1072: DM_POLYTOPE_QUADRILATERAL,
1073: 0,
1074: 0,
1075: DM_POLYTOPE_TRIANGLE,
1076: 1,
1077: 1,
1078: 2,
1079: DM_POLYTOPE_TRIANGLE,
1080: 1,
1081: 2,
1082: 2,
1083: DM_POLYTOPE_TRIANGLE,
1084: 1,
1085: 3,
1086: 2,
1087: DM_POLYTOPE_TRIANGLE,
1088: 1,
1089: 4,
1090: 2,
1091: DM_POLYTOPE_QUADRILATERAL,
1092: 0,
1093: 0,
1094: DM_POLYTOPE_TRIANGLE,
1095: 0,
1096: 0,
1097: DM_POLYTOPE_TRIANGLE,
1098: 0,
1099: 3,
1100: DM_POLYTOPE_TRIANGLE,
1101: 0,
1102: 2,
1103: DM_POLYTOPE_TRIANGLE,
1104: 0,
1105: 1,
1106: };
1107: static PetscInt tpyrO[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, -1, -1,
1108: -1, 0, -3, -2, -3, 0, -3, -2, -3, 0, -3, -2, -3, 0, -3, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 0};
1110: if (rt) *rt = 0;
1111: switch (source) {
1112: case DM_POLYTOPE_POINT:
1113: *Nt = 1;
1114: *target = vertexT;
1115: *size = vertexS;
1116: *cone = vertexC;
1117: *ornt = vertexO;
1118: break;
1119: case DM_POLYTOPE_SEGMENT:
1120: *Nt = 2;
1121: *target = segT;
1122: *size = segS;
1123: *cone = segC;
1124: *ornt = segO;
1125: break;
1126: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1127: *Nt = 1;
1128: *target = tvertT;
1129: *size = tvertS;
1130: *cone = tvertC;
1131: *ornt = tvertO;
1132: break;
1133: case DM_POLYTOPE_TRIANGLE:
1134: *Nt = 2;
1135: *target = triT;
1136: *size = triS;
1137: *cone = triC;
1138: *ornt = triO;
1139: break;
1140: case DM_POLYTOPE_QUADRILATERAL:
1141: *Nt = 3;
1142: *target = quadT;
1143: *size = quadS;
1144: *cone = quadC;
1145: *ornt = quadO;
1146: break;
1147: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1148: *Nt = 2;
1149: *target = tsegT;
1150: *size = tsegS;
1151: *cone = tsegC;
1152: *ornt = tsegO;
1153: break;
1154: case DM_POLYTOPE_TETRAHEDRON:
1155: *Nt = 3;
1156: *target = tetT;
1157: *size = tetS;
1158: *cone = tetC;
1159: *ornt = tetO;
1160: break;
1161: case DM_POLYTOPE_HEXAHEDRON:
1162: *Nt = 4;
1163: *target = hexT;
1164: *size = hexS;
1165: *cone = hexC;
1166: *ornt = hexO;
1167: break;
1168: case DM_POLYTOPE_TRI_PRISM:
1169: *Nt = 4;
1170: *target = tripT;
1171: *size = tripS;
1172: *cone = tripC;
1173: *ornt = tripO;
1174: break;
1175: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1176: *Nt = 2;
1177: *target = ttriT;
1178: *size = ttriS;
1179: *cone = ttriC;
1180: *ornt = ttriO;
1181: break;
1182: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1183: *Nt = 3;
1184: *target = tquadT;
1185: *size = tquadS;
1186: *cone = tquadC;
1187: *ornt = tquadO;
1188: break;
1189: case DM_POLYTOPE_PYRAMID:
1190: *Nt = 5;
1191: *target = tpyrT;
1192: *size = tpyrS;
1193: *cone = tpyrC;
1194: *ornt = tpyrO;
1195: break;
1196: default:
1197: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1198: }
1199: return 0;
1200: }
1202: static PetscErrorCode DMPlexTransformInitialize_Regular(DMPlexTransform tr)
1203: {
1204: tr->ops->view = DMPlexTransformView_Regular;
1205: tr->ops->setup = DMPlexTransformSetUp_Regular;
1206: tr->ops->destroy = DMPlexTransformDestroy_Regular;
1207: tr->ops->celltransform = DMPlexTransformCellRefine_Regular;
1208: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Regular;
1209: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
1210: return 0;
1211: }
1213: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform tr)
1214: {
1215: DMPlexRefine_Regular *f;
1218: PetscNew(&f);
1219: tr->data = f;
1221: DMPlexTransformInitialize_Regular(tr);
1222: return 0;
1223: }