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