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