Actual source code: plexrefine.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscsf.h>
5: const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "DMPlexCellRefinerTypes", "DM_REFINER_", 0};
7: /*
8: Note that j and invj are non-square:
9: v0 + j x_face = x_cell
10: invj (x_cell - v0) = x_face
11: */
12: static PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
13: {
14: /*
15: 2
16: |\
17: | \
18: | \
19: | \
20: | \
21: | \
22: | \
23: 2 1
24: | \
25: | \
26: | \
27: 0---0-------1
28: v0[Nf][dc]: 3 x 2
29: J[Nf][df][dc]: 3 x 1 x 2
30: invJ[Nf][dc][df]: 3 x 2 x 1
31: detJ[Nf]: 3
32: */
33: static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
34: static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
35: static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0};
36: static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0};
37: /*
38: 3---------2---------2
39: | |
40: | |
41: | |
42: 3 1
43: | |
44: | |
45: | |
46: 0---------0---------1
48: v0[Nf][dc]: 4 x 2
49: J[Nf][df][dc]: 4 x 1 x 2
50: invJ[Nf][dc][df]: 4 x 2 x 1
51: detJ[Nf]: 4
52: */
53: static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 -1.0, 0.0};
54: static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
55: static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
56: static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0};
59: switch (ct) {
60: case DM_POLYTOPE_TRIANGLE: if (Nf) *Nf = 3; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; if (detJ) *detJ = tri_detJ; break;
61: case DM_POLYTOPE_QUADRILATERAL: if (Nf) *Nf = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; if (detJ) *detJ = quad_detJ; break;
62: default:
63: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
64: }
65: return(0);
66: }
68: /* Gets the affine map from the original cell to each subcell */
69: static PetscErrorCode DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
70: {
71: /*
72: 2
73: |\
74: | \
75: | \
76: | \
77: | C \
78: | \
79: | \
80: 2---1---1
81: |\ D / \
82: | 2 0 \
83: |A \ / B \
84: 0---0-------1
85: */
86: static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
87: static PetscReal tri_J[] = {0.5, 0.0,
88: 0.0, 0.5,
90: 0.5, 0.0,
91: 0.0, 0.5,
93: 0.5, 0.0,
94: 0.0, 0.5,
96: 0.0, -0.5,
97: 0.5, 0.5};
98: static PetscReal tri_invJ[] = {2.0, 0.0,
99: 0.0, 2.0,
101: 2.0, 0.0,
102: 0.0, 2.0,
104: 2.0, 0.0,
105: 0.0, 2.0,
107: 2.0, 2.0,
108: -2.0, 0.0};
109: /*
110: 3---------2---------2
111: | | |
112: | D 2 C |
113: | | |
114: 3----3----0----1----1
115: | | |
116: | A 0 B |
117: | | |
118: 0---------0---------1
119: */
120: static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
121: static PetscReal quad_J[] = {0.5, 0.0,
122: 0.0, 0.5,
124: 0.5, 0.0,
125: 0.0, 0.5,
127: 0.5, 0.0,
128: 0.0, 0.5,
130: 0.5, 0.0,
131: 0.0, 0.5};
132: static PetscReal quad_invJ[] = {2.0, 0.0,
133: 0.0, 2.0,
135: 2.0, 0.0,
136: 0.0, 2.0,
138: 2.0, 0.0,
139: 0.0, 2.0,
141: 2.0, 0.0,
142: 0.0, 2.0};
143: /*
144: Bottom (viewed from top) Top
145: 1---------2---------2 7---------2---------6
146: | | | | | |
147: | B 2 C | | H 2 G |
148: | | | | | |
149: 3----3----0----1----1 3----3----0----1----1
150: | | | | | |
151: | A 0 D | | E 0 F |
152: | | | | | |
153: 0---------0---------3 4---------0---------5
154: */
155: 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,
156: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
157: static PetscReal hex_J[] = {0.5, 0.0, 0.0,
158: 0.0, 0.5, 0.0,
159: 0.0, 0.0, 0.5,
161: 0.5, 0.0, 0.0,
162: 0.0, 0.5, 0.0,
163: 0.0, 0.0, 0.5,
165: 0.5, 0.0, 0.0,
166: 0.0, 0.5, 0.0,
167: 0.0, 0.0, 0.5,
169: 0.5, 0.0, 0.0,
170: 0.0, 0.5, 0.0,
171: 0.0, 0.0, 0.5,
173: 0.5, 0.0, 0.0,
174: 0.0, 0.5, 0.0,
175: 0.0, 0.0, 0.5,
177: 0.5, 0.0, 0.0,
178: 0.0, 0.5, 0.0,
179: 0.0, 0.0, 0.5,
181: 0.5, 0.0, 0.0,
182: 0.0, 0.5, 0.0,
183: 0.0, 0.0, 0.5,
185: 0.5, 0.0, 0.0,
186: 0.0, 0.5, 0.0,
187: 0.0, 0.0, 0.5};
188: static PetscReal hex_invJ[] = {2.0, 0.0, 0.0,
189: 0.0, 2.0, 0.0,
190: 0.0, 0.0, 2.0,
192: 2.0, 0.0, 0.0,
193: 0.0, 2.0, 0.0,
194: 0.0, 0.0, 2.0,
196: 2.0, 0.0, 0.0,
197: 0.0, 2.0, 0.0,
198: 0.0, 0.0, 2.0,
200: 2.0, 0.0, 0.0,
201: 0.0, 2.0, 0.0,
202: 0.0, 0.0, 2.0,
204: 2.0, 0.0, 0.0,
205: 0.0, 2.0, 0.0,
206: 0.0, 0.0, 2.0,
208: 2.0, 0.0, 0.0,
209: 0.0, 2.0, 0.0,
210: 0.0, 0.0, 2.0,
212: 2.0, 0.0, 0.0,
213: 0.0, 2.0, 0.0,
214: 0.0, 0.0, 2.0,
216: 2.0, 0.0, 0.0,
217: 0.0, 2.0, 0.0,
218: 0.0, 0.0, 2.0};
221: switch (ct) {
222: case DM_POLYTOPE_TRIANGLE: if (Nc) *Nc = 4; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; break;
223: case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break;
224: case DM_POLYTOPE_HEXAHEDRON: if (Nc) *Nc = 8; if (v0) *v0 = hex_v0; if (J) *J = hex_J; if (invJ) *invJ = hex_invJ; break;
225: default:
226: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
227: }
228: return(0);
229: }
231: /* Should this be here or in the DualSpace somehow? */
232: PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside)
233: {
234: PetscReal sum = 0.0;
235: PetscInt d;
238: *inside = PETSC_TRUE;
239: switch (ct) {
240: case DM_POLYTOPE_TRIANGLE:
241: case DM_POLYTOPE_TETRAHEDRON:
242: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) {
243: if (point[d] < -1.0) {*inside = PETSC_FALSE; break;}
244: sum += point[d];
245: }
246: if (sum > PETSC_SMALL) {*inside = PETSC_FALSE; break;}
247: break;
248: case DM_POLYTOPE_QUADRILATERAL:
249: case DM_POLYTOPE_HEXAHEDRON:
250: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d)
251: if (PetscAbsReal(point[d]) > 1.+PETSC_SMALL) {*inside = PETSC_FALSE; break;}
252: break;
253: default:
254: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
255: }
256: return(0);
257: }
259: /* Regular Refinment of Hybrid Meshes
261: We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
262: to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
263: transformations, such as changing from one type of cell to another, as simplex to hex.
265: To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
266: and the types of the new cells.
268: We need the group multiplication table for group actions from the dihedral group for each cell type.
270: We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
271: we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
272: (face, orient) pairs for each subcell.
273: */
275: /*
276: Input Parameters:
277: + ct - The type of the input cell
278: . co - The orientation of this cell in it parent
279: - cp - The requested cone point of this cell assuming orientation 0
281: Output Parameters:
282: . cpnew - The new cone point, taking inout account the orientation co
283: */
284: PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew)
285: {
286: const PetscInt csize = DMPolytopeTypeGetConeSize(ct);
289: if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;}
290: else {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;}
291: return(0);
292: }
294: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
295: {
296: static PetscReal seg_v[] = {-1.0, 0.0, 1.0};
297: 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};
298: 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};
299: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
300: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
301: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0};
302: static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
303: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0,
304: -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0,
305: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0,
306: -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
307: -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0,
308: -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0,
309: -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
310: -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0};
313: switch (ct) {
314: case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break;
315: case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break;
316: case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break;
317: case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break;
318: case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break;
319: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
320: }
321: return(0);
322: }
324: static PetscErrorCode DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
325: {
326: 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, -1.0/3.0, -1.0/3.0};
327: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
328: -1.0, 0.0, -1.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
329: -1.0, -1.0, 0.0, -1.0/3.0, -1.0, -1.0/3.0, 0.0, -1.0, 0.0,
330: -1.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0,
331: -1.0, -1.0, 1.0, -0.5, -0.5, -0.5};
332: PetscErrorCode ierr;
335: switch (ct) {
336: case DM_POLYTOPE_SEGMENT:
337: case DM_POLYTOPE_QUADRILATERAL:
338: case DM_POLYTOPE_HEXAHEDRON:
339: DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);break;
340: case DM_POLYTOPE_TRIANGLE: *Nv = 7; *subcellV = tri_v; break;
341: case DM_POLYTOPE_TETRAHEDRON: *Nv = 15; *subcellV = tet_v; break;
342: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
343: }
344: return(0);
345: }
347: /*
348: DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type
350: Input Parameters:
351: + cr - The DMPlexCellRefiner object
352: - ct - The cell type
354: Output Parameters:
355: + Nv - The number of refined vertices in the closure of the reference cell of given type
356: - subcellV - The coordinates of these vertices in the reference cell
358: Level: developer
360: .seealso: DMPlexCellRefinerGetSubcellVertices()
361: */
362: static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
363: {
367: (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);
368: return(0);
369: }
371: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
372: {
373: static PetscInt seg_v[] = {0, 1, 1, 2};
374: static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5};
375: static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3};
376: static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9,
377: 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6};
378: 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,
379: 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};
382: if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
383: switch (ct) {
384: case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break;
385: case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break;
386: case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
387: case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break;
388: case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break;
389: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
390: }
391: return(0);
392: }
394: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
395: {
396: static PetscInt tri_v[] = {0, 3, 6, 5, 3, 1, 4, 6, 5, 6, 4, 2};
397: static PetscInt tet_v[] = {0, 3, 4, 1, 7, 8, 14, 10, 6, 12, 11, 5, 3, 4, 14, 10, 2, 5, 11, 9, 1, 8, 14, 4, 13, 12 , 10, 7, 9, 8, 14, 11};
398: PetscErrorCode ierr;
401: switch (ct) {
402: case DM_POLYTOPE_SEGMENT:
403: case DM_POLYTOPE_QUADRILATERAL:
404: case DM_POLYTOPE_HEXAHEDRON:
405: DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);break;
406: case DM_POLYTOPE_TRIANGLE:
407: if (rct != DM_POLYTOPE_QUADRILATERAL) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
408: *Nv = 4; *subcellV = &tri_v[r*(*Nv)]; break;
409: case DM_POLYTOPE_TETRAHEDRON:
410: if (rct != DM_POLYTOPE_HEXAHEDRON) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
411: *Nv = 8; *subcellV = &tet_v[r*(*Nv)]; break;
412: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
413: }
414: return(0);
415: }
417: /*
418: DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type
420: Input Parameters:
421: + cr - The DMPlexCellRefiner object
422: . ct - The cell type
423: . rct - The type of subcell
424: - r - The subcell index
426: Output Parameters:
427: + Nv - The number of refined vertices in the subcell
428: - subcellV - The indices of these vertices in the set of vertices returned by DMPlexCellRefinerGetCellVertices()
430: Level: developer
432: .seealso: DMPlexCellRefinerGetCellVertices()
433: */
434: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
435: {
439: (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);
440: return(0);
441: }
443: /*
444: Input Parameters:
445: + cr - The DMPlexCellRefiner
446: . pct - The cell type of the parent, from whom the new cell is being produced
447: . po - The orientation of the parent cell in its enclosing parent
448: . ct - The type being produced
449: . r - The replica number requested for the produced cell type
450: - o - The relative orientation of the replica
452: Output Parameters:
453: + rnew - The replica number, given the orientation of of the parent
454: - onew - The replica orientation, given the orientation of the parent
455: */
456: static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
457: {
461: (*cr->ops->mapsubcells)(cr, pct, po, ct, r, o, rnew, onew);
462: return(0);
463: }
465: static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
466: {
467: /* We shift any input orientation in order to make it non-negative
468: The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o)
469: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
470: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
471: */
472: PetscInt tri_seg_o[] = {-2, 0,
473: -2, 0,
474: -2, 0,
475: 0, -2,
476: 0, -2,
477: 0, -2};
478: PetscInt tri_seg_r[] = {1, 0, 2,
479: 0, 2, 1,
480: 2, 1, 0,
481: 0, 1, 2,
482: 1, 2, 0,
483: 2, 0, 1};
484: PetscInt tri_tri_o[] = {0, 1, 2, -3, -2, -1,
485: 2, 0, 1, -2, -1, -3,
486: 1, 2, 0, -1, -3, -2,
487: -3, -2, -1, 0, 1, 2,
488: -1, -3, -2, 1, 2, 0,
489: -2, -1, -3, 2, 0, 1};
490: /* orientation if the replica is the center triangle */
491: PetscInt tri_tri_o_c[] = {2, 0, 1, -2, -1, -3,
492: 1, 2, 0, -1, -3, -2,
493: 0, 1, 2, -3, -2, -1,
494: -3, -2, -1, 0, 1, 2,
495: -1, -3, -2, 1, 2, 0,
496: -2, -1, -3, 2, 0, 1};
497: PetscInt tri_tri_r[] = {0, 2, 1, 3,
498: 2, 1, 0, 3,
499: 1, 0, 2, 3,
500: 0, 1, 2, 3,
501: 1, 2, 0, 3,
502: 2, 0, 1, 3};
503: PetscInt quad_seg_r[] = {3, 2, 1, 0,
504: 2, 1, 0, 3,
505: 1, 0, 3, 2,
506: 0, 3, 2, 1,
507: 0, 1, 2, 3,
508: 1, 2, 3, 0,
509: 2, 3, 0, 1,
510: 3, 0, 1, 2};
511: PetscInt quad_quad_o[] = { 0, 1, 2, 3, -4, -3, -2, -1,
512: 4, 0, 1, 2, -3, -2, -1, -4,
513: 3, 4, 0, 1, -2, -1, -4, -3,
514: 2, 3, 4, 0, -1, -4, -3, -2,
515: -4, -3, -2, -1, 0, 1, 2, 3,
516: -1, -4, -3, -2, 1, 2, 3, 0,
517: -2, -1, -4, -3, 2, 3, 0, 1,
518: -3, -2, -1, -4, 3, 0, 1, 2};
519: PetscInt quad_quad_r[] = {0, 3, 2, 1,
520: 3, 2, 1, 0,
521: 2, 1, 0, 3,
522: 1, 0, 3, 2,
523: 0, 1, 2, 3,
524: 1, 2, 3, 0,
525: 2, 3, 0, 1,
526: 3, 0, 1, 2};
527: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
528: 1, 0, -1, -2,
529: -2, -1, 0, 1,
530: -1, -2, 1, 0};
531: PetscInt tquad_tquad_r[] = {1, 0,
532: 1, 0,
533: 0, 1,
534: 0, 1};
537: /* The default is no transformation */
538: *rnew = r;
539: *onew = o;
540: switch (pct) {
541: case DM_POLYTOPE_SEGMENT:
542: if (ct == DM_POLYTOPE_SEGMENT) {
543: if (po == 0 || po == -1) {*rnew = r; *onew = o;}
544: else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;}
545: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po);
546: }
547: break;
548: case DM_POLYTOPE_TRIANGLE:
549: switch (ct) {
550: case DM_POLYTOPE_SEGMENT:
551: if (o == -1) o = 0;
552: if (o == -2) o = 1;
553: *onew = tri_seg_o[(po+3)*2+o];
554: *rnew = tri_seg_r[(po+3)*3+r];
555: break;
556: case DM_POLYTOPE_TRIANGLE:
557: *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3];
558: *rnew = tri_tri_r[(po+3)*4+r];
559: break;
560: default: break;
561: }
562: break;
563: case DM_POLYTOPE_QUADRILATERAL:
564: switch (ct) {
565: case DM_POLYTOPE_SEGMENT:
566: *onew = o;
567: *rnew = quad_seg_r[(po+4)*4+r];
568: break;
569: case DM_POLYTOPE_QUADRILATERAL:
570: *onew = quad_quad_o[(po+4)*8+o+4];
571: *rnew = quad_quad_r[(po+4)*4+r];
572: break;
573: default: break;
574: }
575: break;
576: case DM_POLYTOPE_SEG_PRISM_TENSOR:
577: switch (ct) {
578: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
579: case DM_POLYTOPE_SEG_PRISM_TENSOR:
580: *onew = tquad_tquad_o[(po+2)*4+o+2];
581: *rnew = tquad_tquad_r[(po+2)*2+r];
582: break;
583: default: break;
584: }
585: break;
586: default: break;
587: }
588: return(0);
589: }
591: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
592: {
594: /* We shift any input orientation in order to make it non-negative
595: The orientation array o[po][o] gives the orientation the new replica rnew has to have in order to reproduce the face sequence from (r, o)
596: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
597: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
598: */
599: PetscInt tri_seg_o[] = {0, -2,
600: 0, -2,
601: 0, -2,
602: 0, -2,
603: 0, -2,
604: 0, -2};
605: PetscInt tri_seg_r[] = {2, 1, 0,
606: 1, 0, 2,
607: 0, 2, 1,
608: 0, 1, 2,
609: 1, 2, 0,
610: 2, 0, 1};
611: PetscInt tri_tri_o[] = {0, 1, 2, 3, -4, -3, -2, -1,
612: 3, 0, 1, 2, -3, -2, -1, -4,
613: 1, 2, 3, 0, -1, -4, -3, -2,
614: -4, -3, -2, -1, 0, 1, 2, 3,
615: -1, -4, -3, -2, 1, 2, 3, 0,
616: -3, -2, -1, -4, 3, 0, 1, 2};
617: PetscInt tri_tri_r[] = {0, 2, 1,
618: 2, 1, 0,
619: 1, 0, 2,
620: 0, 1, 2,
621: 1, 2, 0,
622: 2, 0, 1};
623: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
624: 1, 0, -1, -2,
625: -2, -1, 0, 1,
626: -1, -2, 1, 0};
627: PetscInt tquad_tquad_r[] = {1, 0,
628: 1, 0,
629: 0, 1,
630: 0, 1};
631: PetscInt tquad_quad_o[] = {-2, -1, -4, -3, 2, 3, 0, 1,
632: 1, 2, 3, 0, -1, -4, -3, -2,
633: -4, -3, -2, -1, 0, 1, 2, 3,
634: 1, 0, 3, 2, -3, -4, -1, -2};
637: *rnew = r;
638: *onew = o;
639: switch (pct) {
640: case DM_POLYTOPE_TRIANGLE:
641: switch (ct) {
642: case DM_POLYTOPE_SEGMENT:
643: if (o == -1) o = 0;
644: if (o == -2) o = 1;
645: *onew = tri_seg_o[(po+3)*2+o];
646: *rnew = tri_seg_r[(po+3)*3+r];
647: break;
648: case DM_POLYTOPE_QUADRILATERAL:
649: *onew = tri_tri_o[(po+3)*8+o+4];
650: /* TODO See sheet, for fixing po == 1 and 2 */
651: if (po == 2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
652: if (po == 2 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
653: if (po == 1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
654: if (po == 1 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
655: if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
656: if (po == -1 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
657: if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
658: if (po == -2 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
659: *rnew = tri_tri_r[(po+3)*3+r];
660: break;
661: default: break;
662: }
663: break;
664: case DM_POLYTOPE_SEG_PRISM_TENSOR:
665: switch (ct) {
666: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
667: case DM_POLYTOPE_SEG_PRISM_TENSOR:
668: *onew = tquad_tquad_o[(po+2)*4+o+2];
669: *rnew = tquad_tquad_r[(po+2)*2+r];
670: break;
671: case DM_POLYTOPE_QUADRILATERAL:
672: *onew = tquad_quad_o[(po+2)*8+o+4];
673: *rnew = tquad_tquad_r[(po+2)*2+r];
674: break;
675: default: break;
676: }
677: break;
678: default:
679: DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
680: }
681: return(0);
682: }
684: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
685: {
686: return DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
687: }
689: /*@
690: DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type
692: Input Parameter:
693: . source - The cell type for a source point
695: Output Parameter:
696: + Nt - The number of cell types generated by refinement
697: . target - The cell types generated
698: . size - The number of subcells of each type, ordered by dimension
699: . cone - A list of the faces for each subcell of the same type as source
700: - ornt - A list of the face orientations for each subcell of the same type as source
702: Note: The cone array gives the cone of each subcell listed by the first three outputs. For the each cone point, we
703: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
704: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
705: $ the number of cones to be taken, or 0 for the current cell
706: $ the cell cone point number at each level from which it is subdivided
707: $ the replica number r of the subdivision.
708: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
709: $ Nt = 2
710: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
711: $ size = {1, 2}
712: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
713: $ ornt = { 0, 0, 0, 0}
715: Level: developer
717: .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform()
718: @*/
719: PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
720: {
724: (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);
725: return(0);
726: }
728: static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
729: {
730: /* All vertices remain in the refined mesh */
731: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
732: static PetscInt vertexS[] = {1};
733: static PetscInt vertexC[] = {0};
734: static PetscInt vertexO[] = {0};
735: /* Split all edges with a new vertex, making two new 2 edges
736: 0--0--0--1--1
737: */
738: static DMPolytopeType edgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
739: static PetscInt edgeS[] = {1, 2};
740: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
741: static PetscInt edgeO[] = { 0, 0, 0, 0};
742: /* Do not split tensor edges */
743: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
744: static PetscInt tedgeS[] = {1};
745: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
746: static PetscInt tedgeO[] = { 0, 0};
747: /* Add 3 edges inside every triangle, making 4 new triangles.
748: 2
749: |\
750: | \
751: | \
752: 0 1
753: | C \
754: | \
755: | \
756: 2---1---1
757: |\ D / \
758: 1 2 0 0
759: |A \ / B \
760: 0-0-0---1---1
761: */
762: static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
763: static PetscInt triS[] = {3, 4};
764: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
765: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
766: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
767: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
768: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
769: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
770: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
771: static PetscInt triO[] = {0, 0,
772: 0, 0,
773: 0, 0,
774: 0, -2, 0,
775: 0, 0, -2,
776: -2, 0, 0,
777: 0, 0, 0};
778: /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
779: 3----1----2----0----2
780: | | |
781: 0 D 2 C 1
782: | | |
783: 3----3----0----1----1
784: | | |
785: 1 A 0 B 0
786: | | |
787: 0----0----0----1----1
788: */
789: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
790: static PetscInt quadS[] = {1, 4, 4};
791: static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
792: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
793: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
794: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
795: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1,
796: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
797: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2,
798: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
799: static PetscInt quadO[] = {0, 0,
800: 0, 0,
801: 0, 0,
802: 0, 0,
803: 0, 0, -2, 0,
804: 0, 0, 0, -2,
805: -2, 0, 0, 0,
806: 0, -2, 0, 0};
807: /* Add 1 edge inside every tensor quad, making 2 new tensor quads
808: 2----2----1----3----3
809: | | |
810: | | |
811: | | |
812: 4 A 6 B 5
813: | | |
814: | | |
815: | | |
816: 0----0----0----1----1
817: */
818: static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
819: static PetscInt tquadS[] = {1, 2};
820: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
821: 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,
822: 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};
823: static PetscInt tquadO[] = {0, 0,
824: 0, 0, 0, 0,
825: 0, 0, 0, 0};
826: /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
827: 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
828: 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]
829: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
830: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
831: The first four tets just cut off the corners, using the replica notation for new vertices,
832: [v0, (e0, 0), (e2, 0), (e3, 0)]
833: [(e0, 0), v1, (e1, 0), (e4, 0)]
834: [(e2, 0), (e1, 0), v2, (e5, 0)]
835: [(e3, 0), (e4, 0), (e5, 0), v3 ]
836: The next four tets match a vertex to the newly created faces from cutting off those first tets.
837: [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
838: [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
839: [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
840: [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
841: 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
842: [(e2, 0), (e0, 0), (e3, 0)]
843: [(e0, 0), (e1, 0), (e4, 0)]
844: [(e2, 0), (e5, 0), (e1, 0)]
845: [(e3, 0), (e4, 0), (e5, 0)]
846: The next four, from the second group of tets, are
847: [(e2, 0), (e0, 0), (e5, 0)]
848: [(e4, 0), (e0, 0), (e5, 0)]
849: [(e0, 0), (e1, 0), (e5, 0)]
850: [(e5, 0), (e3, 0), (e0, 0)]
851: I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
852: */
853: static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
854: static PetscInt tetS[] = {1, 8, 8};
855: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0,
856: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2,
857: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1,
858: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1,
859: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1,
860: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
861: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1,
862: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0,
863: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0,
864: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0,
865: DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1,
866: DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0,
867: DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2,
868: DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7,
869: DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6,
870: DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3,
871: DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
872: static PetscInt tetO[] = {0, 0,
873: 0, 0, 0,
874: 0, 0, 0,
875: 0, 0, 0,
876: 0, 0, 0,
877: 0, 0, -2,
878: 0, 0, -2,
879: 0, -2, -2,
880: 0, -2, 0,
881: 0, 0, 0, 0,
882: 0, 0, 0, 0,
883: 0, 0, 0, 0,
884: 0, 0, 0, 0,
885: -3, 0, 0, -2,
886: -2, 1, 0, 0,
887: -2, -2, -1, 2,
888: -2, 0, -2, 1};
889: /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
890: 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
891: [v0, v1, v2, v3] f0 bottom
892: [v4, v5, v6, v7] f1 top
893: [v0, v3, v5, v4] f2 front
894: [v2, v1, v7, v6] f3 back
895: [v3, v2, v6, v5] f4 right
896: [v0, v4, v7, v1] f5 left
897: The eight hexes are divided into four on the bottom, and four on the top,
898: [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)]
899: [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)]
900: [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)]
901: [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)]
902: [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)]
903: [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)]
904: [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)]
905: [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7]
906: 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,
907: [(e9, 0), (f2, 0), (c0, 0), (f5, 0)]
908: [(f5, 0), (c0, 0), (f3, 0), (e10, 0)]
909: [(c0, 0), (f4, 0), (e11, 0), (f3, 0)]
910: [(f2, 0), (e8, 0), (f4, 0), (c0, 0)]
911: and on the x-z plane,
912: [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
913: [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
914: [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
915: [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
916: and on the y-z plane,
917: [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
918: [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
919: [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
920: [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
921: */
922: static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
923: static PetscInt hexS[] = {1, 6, 12, 8};
924: static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
925: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
926: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
927: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
928: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
929: DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0,
930: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0,
931: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2,
932: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3,
933: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2,
934: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0,
935: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1,
936: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2,
937: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
938: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3,
939: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
940: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3,
941: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
942: 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,
943: 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,
944: 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,
945: 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,
946: 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,
947: 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,
948: 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,
949: 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};
950: static PetscInt hexO[] = {0, 0,
951: 0, 0,
952: 0, 0,
953: 0, 0,
954: 0, 0,
955: 0, 0,
956: 0, 0, -2, -2,
957: 0, -2, -2, 0,
958: -2, -2, 0, 0,
959: -2, 0, 0, -2,
960: -2, 0, 0, -2,
961: -2, -2, 0, 0,
962: 0, -2, -2, 0,
963: 0, 0, -2, -2,
964: 0, 0, -2, -2,
965: -2, 0, 0, -2,
966: -2, -2, 0, 0,
967: 0, -2, -2, 0,
968: 0, 0, 0, 0, -4, 0,
969: 0, 0, -1, 0, -4, 0,
970: 0, 0, -1, 0, 0, 0,
971: 0, 0, 0, 0, 0, 0,
972: -4, 0, 0, 0, -4, 0,
973: -4, 0, 0, 0, 0, 0,
974: -4, 0, -1, 0, 0, 0,
975: -4, 0, -1, 0, -4, 0};
976: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
977: static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
978: static PetscInt tripS[] = {3, 4, 6, 8};
979: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0,
980: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0,
981: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0,
982: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1,
983: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0,
984: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
985: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
986: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
987: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
988: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
989: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
990: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
991: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
992: 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,
993: 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,
994: 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,
995: 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,
996: 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,
997: 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,
998: 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,
999: 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};
1000: static PetscInt tripO[] = {0, 0,
1001: 0, 0,
1002: 0, 0,
1003: 0, -2, -2,
1004: -2, 0, -2,
1005: -2, -2, 0,
1006: 0, 0, 0,
1007: -2, 0, -2, -2,
1008: -2, 0, -2, -2,
1009: -2, 0, -2, -2,
1010: 0, -2, -2, 0,
1011: 0, -2, -2, 0,
1012: 0, -2, -2, 0,
1013: 0, 0, 0, -1, 0,
1014: 0, 0, 0, 0, -1,
1015: 0, 0, -1, 0, 0,
1016: 2, 0, 0, 0, 0,
1017: -3, 0, 0, -1, 0,
1018: -3, 0, 0, 0, -1,
1019: -3, 0, -1, 0, 0,
1020: -3, 0, 0, 0, 0};
1021: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
1022: 2
1023: |\
1024: | \
1025: | \
1026: 0---1
1028: 2
1030: 0 1
1032: 2
1033: |\
1034: | \
1035: | \
1036: 0---1
1037: */
1038: static DMPolytopeType ttripT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
1039: static PetscInt ttripS[] = {3, 4};
1040: static PetscInt ttripC[] = {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,
1041: 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,
1042: 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,
1043: 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,
1044: 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,
1045: 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,
1046: 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};
1047: static PetscInt ttripO[] = {0, 0, 0, 0,
1048: 0, 0, 0, 0,
1049: 0, 0, 0, 0,
1050: 0, 0, 0, -1, 0,
1051: 0, 0, 0, 0, -1,
1052: 0, 0, -1, 0, 0,
1053: 0, 0, 0, 0, 0};
1054: /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
1055: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1056: static PetscInt tquadpS[] = {1, 4, 4};
1057: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1058: 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,
1059: 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,
1060: 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,
1061: 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,
1062: 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,
1063: 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,
1064: 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,
1065: 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};
1066: static PetscInt tquadpO[] = {0, 0,
1067: 0, 0, 0, 0,
1068: 0, 0, 0, 0,
1069: 0, 0, 0, 0,
1070: 0, 0, 0, 0,
1071: 0, 0, 0, 0, -1, 0,
1072: 0, 0, 0, 0, 0, -1,
1073: 0, 0, -1, 0, 0, 0,
1074: 0, 0, 0, -1, 0, 0};
1077: switch (source) {
1078: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
1079: case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
1080: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1081: case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1082: case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
1083: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1084: case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1085: case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
1086: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1087: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1088: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1089: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1090: }
1091: return(0);
1092: }
1094: static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1095: {
1097: /* Change tensor edges to segments */
1098: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT};
1099: static PetscInt tedgeS[] = {1};
1100: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1101: static PetscInt tedgeO[] = { 0, 0};
1102: /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals.
1103: 2
1104: |\
1105: | \
1106: | \
1107: | \
1108: 0 1
1109: | \
1110: | \
1111: 2 1
1112: |\ / \
1113: | 2 1 \
1114: | \ / \
1115: 1 | 0
1116: | 0 \
1117: | | \
1118: | | \
1119: 0-0-0-----1-----1
1120: */
1121: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1122: static PetscInt triS[] = {1, 3, 3};
1123: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1124: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1125: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1126: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1127: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1128: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1129: static PetscInt triO[] = {0, 0,
1130: 0, 0,
1131: 0, 0,
1132: 0, 0, -2, 0,
1133: 0, 0, 0, -2,
1134: 0, -2, 0, 0};
1135: /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals
1136: 2----2----1----3----3
1137: | | |
1138: | | |
1139: | | |
1140: 4 A 6 B 5
1141: | | |
1142: | | |
1143: | | |
1144: 0----0----0----1----1
1145: */
1146: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1147: static PetscInt tquadS[] = {1, 2};
1148: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1149: /* TODO Fix these */
1150: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1151: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0};
1152: static PetscInt tquadO[] = {0, 0,
1153: 0, 0, -2, -2,
1154: 0, 0, -2, -2};
1155: /* Add 6 triangles inside every cell, making 4 new hexs
1156: TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it
1157: 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
1158: 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]
1159: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1160: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1161: We make a new hex in each corner
1162: [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1163: [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)]
1164: [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)]
1165: [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)]
1166: We create a new face for each edge
1167: [(e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1168: [(f0, 0), (e0, 0), (f1, 0), (c0, 0)]
1169: [(e2, 0), (f0, 0), (c0, 0), (f2, 0)]
1170: [(f3, 0), (e4, 0), (f1, 0), (c0, 0)]
1171: [(e1, 0), (f3, 0), (c0, 0), (f0, 0)]
1172: [(e5, 0), (f3, 0), (c0, 0), (f2, 0)]
1173: I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another.
1174: */
1175: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1176: static PetscInt tetS[] = {1, 4, 6, 4};
1177: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1178: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1179: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1180: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1181: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0,
1182: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1183: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1184: DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3,
1185: DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1186: DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1187: DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
1188: DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2,
1189: DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 0, 2,
1190: DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2};
1191: static PetscInt tetO[] = {0, 0,
1192: 0, 0,
1193: 0, 0,
1194: 0, 0,
1195: 0, 0, -2, -2,
1196: -2, 0, 0, -2,
1197: 0, 0, -2, -2,
1198: -2, 0, 0, -2,
1199: 0, 0, -2, -2,
1200: 0, 0, -2, -2,
1201: 0, 0, 0, 0, 0, 0,
1202: 1, -1, 1, 0, 0, 3,
1203: 0, -4, 1, -1, 0, 3,
1204: 1, -4, 3, -2, -4, 3};
1205: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1206: static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1207: static PetscInt tripS[] = {1, 5, 9, 6};
1208: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1209: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1210: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1211: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1212: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1213: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1214: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2,
1215: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1216: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1217: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1218: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1219: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1220: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1221: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1222: DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 0, 3, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1,
1223: DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 3,
1224: DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 5, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 4, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0,
1225: DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 3, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 0, 6, DM_POLYTOPE_QUADRILATERAL, 1, 4, 2,
1226: DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 3, 3, DM_POLYTOPE_QUADRILATERAL, 0, 6,
1227: DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 8, DM_POLYTOPE_QUADRILATERAL, 1, 3, 2, DM_POLYTOPE_QUADRILATERAL, 0, 7, DM_POLYTOPE_QUADRILATERAL, 1, 4, 3};
1228: static PetscInt tripO[] = {0, 0,
1229: 0, 0,
1230: 0, 0,
1231: 0, 0,
1232: 0, 0,
1233: 0, 0, -2, -2,
1234: -2, 0, 0, -2,
1235: 0, -2, -2, 0,
1236: 0, 0, -2, -2,
1237: 0, 0, -2, -2,
1238: 0, 0, -2, -2,
1239: 0, -2, -2, 0,
1240: 0, -2, -2, 0,
1241: 0, -2, -2, 0,
1242: 0, 0, 0, -1, 0, 1,
1243: 0, 0, 0, 0, 0, -4,
1244: 0, 0, 0, 0, -1, 1,
1245: -4, 0, 0, -1, 0, 1,
1246: -4, 0, 0, 0, 0, -4,
1247: -4, 0, 0, 0, -1, 1};
1248: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms.
1249: 2
1250: |\
1251: | \
1252: | \
1253: 0---1
1255: 2
1257: 0 1
1259: 2
1260: |\
1261: | \
1262: | \
1263: 0---1
1264: */
1265: static DMPolytopeType ttripT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1266: static PetscInt ttripS[] = {1, 3, 3};
1267: static PetscInt ttripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1268: 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,
1269: 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,
1270: 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,
1271: 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, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 1,
1272: 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,
1273: DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 2, DM_POLYTOPE_SEG_PRISM_TENSOR, 0, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 1, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
1274: static PetscInt ttripO[] = {0, 0,
1275: 0, 0, 0, 0,
1276: 0, 0, 0, 0,
1277: 0, 0, 0, 0,
1278: 0, 0, 0, 0, -1, 0,
1279: 0, 0, 0, 0, 0, -1,
1280: 0, 0, 0, -1, 0, 0};
1281: /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms.
1282: 2
1283: |\
1284: | \
1285: | \
1286: 0---1
1288: 2
1290: 0 1
1292: 2
1293: |\
1294: | \
1295: | \
1296: 0---1
1297: */
1298: static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1299: static PetscInt ctripS[] = {1, 3, 3};
1300: static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1301: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1302: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1303: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1304: DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 1,
1305: DM_POLYTOPE_QUADRILATERAL, 1, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 1, 1, DM_POLYTOPE_QUADRILATERAL, 1, 2, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 0, 0,
1306: DM_POLYTOPE_QUADRILATERAL, 1, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 1, 2, DM_POLYTOPE_QUADRILATERAL, 0, 2, DM_POLYTOPE_QUADRILATERAL, 1, 3, 1, DM_POLYTOPE_QUADRILATERAL, 0, 1, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
1307: static PetscInt ctripO[] = {0, 0,
1308: 0, 0, -2, -2,
1309: 0, 0, -2, -2,
1310: 0, 0, -2, -2,
1311: -4, 0, 0, -1, 0, 1,
1312: -4, 0, 0, 0, 0, -4,
1313: -4, 0, 0, 0, -1, 1};
1314: /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */
1315: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1316: static PetscInt tquadpS[] = {1, 4, 4};
1317: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1318: 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,
1319: 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,
1320: 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,
1321: 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,
1322: 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,
1323: 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,
1324: 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,
1325: 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};
1326: static PetscInt tquadpO[] = {0, 0,
1327: 0, 0, 0, 0,
1328: 0, 0, 0, 0,
1329: 0, 0, 0, 0,
1330: 0, 0, 0, 0,
1331: 0, 0, 0, 0, -1, 0,
1332: 0, 0, 0, 0, 0, -1,
1333: 0, 0, -1, 0, 0, 0,
1334: 0, 0, 0, -1, 0, 0};
1335: PetscBool convertTensor = PETSC_TRUE;
1338: if (convertTensor) {
1339: switch (source) {
1340: case DM_POLYTOPE_POINT:
1341: case DM_POLYTOPE_SEGMENT:
1342: case DM_POLYTOPE_QUADRILATERAL:
1343: case DM_POLYTOPE_HEXAHEDRON:
1344: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1345: break;
1346: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1347: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1348: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ctripT; *size = ctripS; *cone = ctripC; *ornt = ctripO; break;
1349: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1350: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1351: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1352: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1353: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1354: }
1355: } else {
1356: switch (source) {
1357: case DM_POLYTOPE_POINT:
1358: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1359: case DM_POLYTOPE_SEGMENT:
1360: case DM_POLYTOPE_QUADRILATERAL:
1361: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1362: case DM_POLYTOPE_HEXAHEDRON:
1363: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1364: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1365: break;
1366: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1367: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1368: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1369: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1370: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1371: }
1372: }
1373: return(0);
1374: }
1376: static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1377: {
1381: switch (source) {
1382: case DM_POLYTOPE_POINT:
1383: case DM_POLYTOPE_SEGMENT:
1384: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1385: case DM_POLYTOPE_TRIANGLE:
1386: case DM_POLYTOPE_TETRAHEDRON:
1387: case DM_POLYTOPE_TRI_PRISM:
1388: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1389: case DM_POLYTOPE_QUADRILATERAL:
1390: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1391: case DM_POLYTOPE_HEXAHEDRON:
1392: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1393: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1394: break;
1395: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1396: }
1397: return(0);
1398: }
1400: static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
1401: {
1402: PetscInt c, cN, *off;
1406: PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
1407: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
1408: const DMPolytopeType ct = (DMPolytopeType) c;
1409: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
1410: const DMPolytopeType ctNew = (DMPolytopeType) cN;
1411: DMPolytopeType *rct;
1412: PetscInt *rsize, *cone, *ornt;
1413: PetscInt Nct, n, i;
1415: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;}
1416: off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
1417: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
1418: const DMPolytopeType ict = (DMPolytopeType) ctOrder[i];
1419: const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];
1421: DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);
1422: if (ict == ct) {
1423: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
1424: if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
1425: break;
1426: }
1427: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
1428: }
1429: }
1430: }
1431: *offset = off;
1432: return(0);
1433: }
1435: static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[])
1436: {
1437: const PetscInt ctSize = DM_NUM_POLYTOPES+1;
1441: if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()");
1442: PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);
1443: PetscArraycpy(cr->ctStart, ctStart, ctSize);
1444: PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);
1445: return(0);
1446: }
1448: /* Construct cell type order since we must loop over cell types in depth order */
1449: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
1450: {
1451: PetscInt *ctO, *ctOInv;
1452: PetscInt c, d, off = 0;
1456: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
1457: for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
1458: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
1459: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
1460: ctO[off++] = c;
1461: }
1462: }
1463: if (DMPolytopeTypeGetDim(ctCell) != 0) {
1464: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
1465: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
1466: ctO[off++] = c;
1467: }
1468: }
1469: for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
1470: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
1471: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
1472: ctO[off++] = c;
1473: }
1474: }
1475: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
1476: if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
1477: ctO[off++] = c;
1478: }
1479: if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
1480: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
1481: ctOInv[ctO[c]] = c;
1482: }
1483: *ctOrder = ctO;
1484: *ctOrderInv = ctOInv;
1485: return(0);
1486: }
1488: PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)
1489: {
1490: DM dm = cr->dm;
1491: DMPolytopeType ctCell;
1492: PetscInt pStart, pEnd, p, c;
1497: if (cr->setupcalled) return(0);
1498: DMPlexGetChart(dm, &pStart, &pEnd);
1499: if (pEnd > pStart) {DMPlexGetCellType(dm, 0, &ctCell);}
1500: else {
1501: PetscInt dim;
1502: DMGetDimension(dm, &dim);
1503: switch (dim) {
1504: case 0: ctCell = DM_POLYTOPE_POINT;break;
1505: case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
1506: case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
1507: case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
1508: default: ctCell = DM_POLYTOPE_TETRAHEDRON;
1509: }
1510: }
1511: DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);
1512: /* Construct sizes and offsets for each cell type */
1513: if (!cr->ctStart) {
1514: PetscInt *ctS, *ctSN, *ctC, *ctCN;
1516: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
1517: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
1518: for (p = pStart; p < pEnd; ++p) {
1519: DMPolytopeType ct;
1520: DMPolytopeType *rct;
1521: PetscInt *rsize, *cone, *ornt;
1522: PetscInt Nct, n;
1524: DMPlexGetCellType(dm, p, &ct);
1525: if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
1526: ++ctC[ct];
1527: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
1528: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
1529: }
1530: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
1531: const PetscInt ct = cr->ctOrder[c];
1532: const PetscInt ctn = cr->ctOrder[c+1];
1534: ctS[ctn] = ctS[ct] + ctC[ct];
1535: ctSN[ctn] = ctSN[ct] + ctCN[ct];
1536: }
1537: PetscFree2(ctC, ctCN);
1538: cr->ctStart = ctS;
1539: cr->ctStartNew = ctSN;
1540: }
1541: CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);
1542: cr->setupcalled = PETSC_TRUE;
1543: return(0);
1544: }
1546: static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v)
1547: {
1551: PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);
1552: return(0);
1553: }
1555: /*
1556: DMPlexCellRefinerView - Views a DMPlexCellRefiner object
1558: Collective on cr
1560: Input Parameters:
1561: + cr - The DMPlexCellRefiner object
1562: - viewer - The PetscViewer object
1564: Level: beginner
1566: .seealso: DMPlexCellRefinerCreate()
1567: */
1568: static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer)
1569: {
1570: PetscBool iascii;
1576: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);}
1577: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1578: PetscViewerASCIIPushTab(viewer);
1579: if (iascii) {DMPlexCellRefinerView_Ascii(cr, viewer);}
1580: PetscViewerASCIIPopTab(viewer);
1581: return(0);
1582: }
1584: PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr)
1585: {
1586: PetscInt c;
1590: if (!*cr) return(0);
1592: PetscObjectDereference((PetscObject) (*cr)->dm);
1593: PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);
1594: PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);
1595: PetscFree((*cr)->offset);
1596: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
1597: PetscFEDestroy(&(*cr)->coordFE[c]);
1598: PetscFEGeomDestroy(&(*cr)->refGeom[c]);
1599: }
1600: PetscFree2((*cr)->coordFE, (*cr)->refGeom);
1601: PetscHeaderDestroy(cr);
1602: return(0);
1603: }
1605: PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr)
1606: {
1607: DMPlexCellRefiner tmp;
1608: PetscErrorCode ierr;
1612: *cr = NULL;
1613: PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);
1614: tmp->setupcalled = PETSC_FALSE;
1616: tmp->dm = dm;
1617: PetscObjectReference((PetscObject) dm);
1618: DMPlexGetCellRefinerType(dm, &tmp->type);
1619: switch (tmp->type) {
1620: case DM_REFINER_REGULAR:
1621: tmp->ops->refine = DMPlexCellRefinerRefine_Regular;
1622: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_Regular;
1623: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_Regular;
1624: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_Regular;
1625: tmp->ops->getaffinetransforms = DMPlexCellRefinerGetAffineTransforms_Regular;
1626: tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular;
1627: break;
1628: case DM_REFINER_TO_BOX:
1629: tmp->ops->refine = DMPlexCellRefinerRefine_ToBox;
1630: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToBox;
1631: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_ToBox;
1632: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox;
1633: break;
1634: case DM_REFINER_TO_SIMPLEX:
1635: tmp->ops->refine = DMPlexCellRefinerRefine_ToSimplex;
1636: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex;
1637: break;
1638: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]);
1639: }
1640: PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);
1641: *cr = tmp;
1642: return(0);
1643: }
1645: /*@
1646: DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell
1648: Input Parameters:
1649: + cr - The DMPlexCellRefiner object
1650: - ct - The cell type
1652: Output Parameters:
1653: + Nc - The number of subcells produced from this cell type
1654: . v0 - The translation of the first vertex for each subcell
1655: . J - The Jacobian for each subcell (map from reference cell to subcell)
1656: - invJ - The inverse Jacobian for each subcell
1658: Level: developer
1660: .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create()
1661: @*/
1662: PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
1663: {
1667: if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner");
1668: (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);
1669: return(0);
1670: }
1672: /*@
1673: DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell
1675: Input Parameters:
1676: + cr - The DMPlexCellRefiner object
1677: - ct - The cell type
1679: Output Parameters:
1680: + Nf - The number of faces for this cell type
1681: . v0 - The translation of the first vertex for each face
1682: . J - The Jacobian for each face (map from original cell to subcell)
1683: . invJ - The inverse Jacobian for each face
1684: - detJ - The determinant of the Jacobian for each face
1686: Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
1688: Level: developer
1690: .seealso: DMPlexCellRefinerGetAffineTransforms(), Create()
1691: @*/
1692: PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
1693: {
1697: if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner");
1698: (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);
1699: return(0);
1700: }
1702: /* Numbering regularly refined meshes
1704: We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute
1705: the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions
1706: about the order of different cell types.
1708: However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous
1709: numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we
1710: will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as
1711: the cell type enumeration within a given depth stratum.
1713: Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we
1714: start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then
1715: offset by the old cell type number and replica number for the insertion.
1716: */
1717: PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
1718: {
1719: DMPolytopeType *rct;
1720: PetscInt *rsize, *cone, *ornt;
1721: PetscInt Nct, n;
1722: PetscInt off = cr->offset[ct*DM_NUM_POLYTOPES+ctNew];
1723: PetscInt ctS = cr->ctStart[ct], ctE = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]];
1724: PetscInt ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]];
1725: PetscInt newp = ctSN;
1726: PetscErrorCode ierr;
1729: if ((p < ctS) || (p >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D is not a %s [%D, %D)", p, DMPolytopeTypes[ct], ctS, ctE);
1730: if (off < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s does not produce type %s", DMPolytopeTypes[ct], DMPolytopeTypes[ctNew]);
1732: newp += off;
1733: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
1734: for (n = 0; n < Nct; ++n) {
1735: if (rct[n] == ctNew) {
1736: if (rsize[n] && r >= rsize[n])
1737: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
1738: newp += (p - ctS) * rsize[n] + r;
1739: break;
1740: }
1741: }
1743: if ((newp < ctSN) || (newp >= ctEN)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %D is not a %s [%D, %D)", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
1744: *pNew = newp;
1745: return(0);
1746: }
1748: static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm)
1749: {
1750: DM dm = cr->dm;
1751: PetscInt pStart, pEnd, p, pNew;
1752: PetscErrorCode ierr;
1755: /* Must create the celltype label here so that we do not automatically try to compute the types */
1756: DMCreateLabel(rdm, "celltype");
1757: DMPlexGetChart(dm, &pStart, &pEnd);
1758: for (p = pStart; p < pEnd; ++p) {
1759: DMPolytopeType ct;
1760: DMPolytopeType *rct;
1761: PetscInt *rsize, *rcone, *rornt;
1762: PetscInt Nct, n, r;
1764: DMPlexGetCellType(dm, p, &ct);
1765: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
1766: for (n = 0; n < Nct; ++n) {
1767: for (r = 0; r < rsize[n]; ++r) {
1768: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
1769: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
1770: DMPlexSetCellType(rdm, pNew, rct[n]);
1771: }
1772: }
1773: }
1774: {
1775: DMLabel ctLabel;
1776: DM_Plex *plex = (DM_Plex *) rdm->data;
1778: DMPlexGetCellTypeLabel(rdm, &ctLabel);
1779: PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
1780: }
1781: return(0);
1782: }
1784: static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm)
1785: {
1786: DM dm = cr->dm;
1787: DMPolytopeType ct;
1788: PetscInt *coneNew, *orntNew;
1789: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1793: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1794: PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
1795: DMPlexGetChart(dm, &pStart, &pEnd);
1796: for (p = pStart; p < pEnd; ++p) {
1797: const PetscInt *cone, *ornt;
1798: PetscInt coff, ooff, c;
1799: DMPolytopeType *rct;
1800: PetscInt *rsize, *rcone, *rornt;
1801: PetscInt Nct, n, r;
1803: DMPlexGetCellType(dm, p, &ct);
1804: DMPlexGetCone(dm, p, &cone);
1805: DMPlexGetConeOrientation(dm, p, &ornt);
1806: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
1807: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1808: const DMPolytopeType ctNew = rct[n];
1809: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1811: for (r = 0; r < rsize[n]; ++r) {
1812: /* pNew is a subcell produced by subdividing p */
1813: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
1814: for (c = 0; c < csizeNew; ++c) {
1815: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1816: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1817: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1818: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1819: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1820: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1821: const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
1822: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1823: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1824: PetscInt lc;
1826: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1827: for (lc = 0; lc < fn; ++lc) {
1828: const PetscInt *ppornt;
1829: PetscInt pcp;
1831: DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);
1832: ppp = pp;
1833: pp = pcone[pcp];
1834: DMPlexGetCellType(dm, pp, &pct);
1835: DMPlexGetCone(dm, pp, &pcone);
1836: DMPlexGetConeOrientation(dm, ppp, &ppornt);
1837: po = ppornt[pcp];
1838: }
1839: pr = rcone[coff++];
1840: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1841: DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);
1842: DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);
1843: orntNew[c] = fo;
1844: }
1845: DMPlexSetCone(rdm, pNew, coneNew);
1846: DMPlexSetConeOrientation(rdm, pNew, orntNew);
1847: }
1848: }
1849: }
1850: PetscFree2(coneNew, orntNew);
1851: DMPlexSymmetrize(rdm);
1852: DMPlexStratify(rdm);
1853: return(0);
1854: }
1856: static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe)
1857: {
1861: if (!cr->coordFE[ct]) {
1862: PetscInt dim, cdim;
1863: PetscBool isSimplex;
1865: switch (ct) {
1866: case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break;
1867: case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break;
1868: case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
1869: case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break;
1870: case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break;
1871: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
1872: }
1873: DMGetCoordinateDim(cr->dm, &cdim);
1874: PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);
1875: {
1876: PetscDualSpace dsp;
1877: PetscQuadrature quad;
1878: DM K;
1879: PetscFEGeom *cg;
1880: PetscReal *Xq, *xq, *wq;
1881: PetscInt Nq, q;
1883: DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);
1884: PetscMalloc1(Nq*cdim, &xq);
1885: for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q];
1886: PetscMalloc1(Nq, &wq);
1887: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
1888: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
1889: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
1890: PetscFESetQuadrature(cr->coordFE[ct], quad);
1892: PetscFEGetDualSpace(cr->coordFE[ct], &dsp);
1893: PetscDualSpaceGetDM(dsp, &K);
1894: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);
1895: cg = cr->refGeom[ct];
1896: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
1897: PetscQuadratureDestroy(&quad);
1898: }
1899: }
1900: *fe = cr->coordFE[ct];
1901: return(0);
1902: }
1904: /*
1905: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
1907: Not collective
1909: Input Parameters:
1910: + cr - The DMPlexCellRefiner
1911: . ct - The type of the parent cell
1912: . rct - The type of the produced cell
1913: . r - The index of the produced cell
1914: - x - The localized coordinates for the parent cell
1916: Output Parameter:
1917: . xr - The localized coordinates for the produced cell
1919: Level: developer
1921: .seealso: DMPlexCellRefinerSetCoordinates()
1922: */
1923: static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1924: {
1925: PetscFE fe = NULL;
1926: PetscInt cdim, Nv, v, *subcellV;
1930: DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);
1931: DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);
1932: PetscFEGetNumComponents(fe, &cdim);
1933: for (v = 0; v < Nv; ++v) {
1934: PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);
1935: }
1936: return(0);
1937: }
1939: static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm)
1940: {
1941: DM dm = cr->dm, cdm;
1942: PetscSection coordSection, coordSectionNew;
1943: Vec coordsLocal, coordsLocalNew;
1944: const PetscScalar *coords;
1945: PetscScalar *coordsNew;
1946: const DMBoundaryType *bd;
1947: const PetscReal *maxCell, *L;
1948: PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1949: PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
1950: PetscErrorCode ierr;
1953: DMGetCoordinateDM(dm, &cdm);
1954: DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
1955: /* Determine if we need to localize coordinates when generating them */
1956: if (isperiodic) {
1957: localizeVertices = PETSC_TRUE;
1958: if (!maxCell) {
1959: PetscBool localized;
1960: DMGetCoordinatesLocalized(dm, &localized);
1961: if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
1962: localizeCells = PETSC_TRUE;
1963: }
1964: }
1966: DMGetCoordinateSection(dm, &coordSection);
1967: PetscSectionGetFieldComponents(coordSection, 0, &dE);
1968: if (maxCell) {
1969: PetscReal maxCellNew[3];
1971: for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
1972: DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
1973: } else {
1974: DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
1975: }
1976: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
1977: PetscSectionSetNumFields(coordSectionNew, 1);
1978: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1979: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1980: if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0, vEndNew);}
1981: else {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}
1983: /* Localization should be inherited */
1984: /* Stefano calculates parent cells for each new cell for localization */
1985: /* Localized cells need coordinates of closure */
1986: for (v = vStartNew; v < vEndNew; ++v) {
1987: PetscSectionSetDof(coordSectionNew, v, dE);
1988: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1989: }
1990: if (localizeCells) {
1991: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1992: for (c = cStart; c < cEnd; ++c) {
1993: PetscInt dof;
1995: PetscSectionGetDof(coordSection, c, &dof);
1996: if (dof) {
1997: DMPolytopeType ct;
1998: DMPolytopeType *rct;
1999: PetscInt *rsize, *rcone, *rornt;
2000: PetscInt dim, cNew, Nct, n, r;
2002: DMPlexGetCellType(dm, c, &ct);
2003: dim = DMPolytopeTypeGetDim(ct);
2004: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2005: /* This allows for different cell types */
2006: for (n = 0; n < Nct; ++n) {
2007: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2008: for (r = 0; r < rsize[n]; ++r) {
2009: PetscInt *closure = NULL;
2010: PetscInt clSize, cl, Nv = 0;
2012: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);
2013: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2014: for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
2015: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2016: PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
2017: PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
2018: }
2019: }
2020: }
2021: }
2022: }
2023: PetscSectionSetUp(coordSectionNew);
2024: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
2025: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
2026: {
2027: VecType vtype;
2028: PetscInt coordSizeNew, bs;
2029: const char *name;
2031: DMGetCoordinatesLocal(dm, &coordsLocal);
2032: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
2033: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
2034: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
2035: PetscObjectGetName((PetscObject) coordsLocal, &name);
2036: PetscObjectSetName((PetscObject) coordsLocalNew, name);
2037: VecGetBlockSize(coordsLocal, &bs);
2038: VecSetBlockSize(coordsLocalNew, bs);
2039: VecGetType(coordsLocal, &vtype);
2040: VecSetType(coordsLocalNew, vtype);
2041: }
2042: VecGetArrayRead(coordsLocal, &coords);
2043: VecGetArray(coordsLocalNew, &coordsNew);
2044: PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
2045: DMPlexGetChart(dm, &pStart, &pEnd);
2046: /* First set coordinates for vertices*/
2047: for (p = pStart; p < pEnd; ++p) {
2048: DMPolytopeType ct;
2049: DMPolytopeType *rct;
2050: PetscInt *rsize, *rcone, *rornt;
2051: PetscInt Nct, n, r;
2052: PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;
2054: DMPlexGetCellType(dm, p, &ct);
2055: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2056: for (n = 0; n < Nct; ++n) {
2057: if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
2058: }
2059: if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2060: PetscInt dof;
2061: PetscSectionGetDof(coordSection, p, &dof);
2062: if (dof) isLocalized = PETSC_TRUE;
2063: }
2064: if (hasVertex) {
2065: PetscScalar *pcoords = NULL;
2066: PetscScalar vcoords[3] = {0., 0., 0.};
2067: PetscInt Nc, Nv, v, d;
2069: DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
2070: if (ct == DM_POLYTOPE_POINT) {
2071: for (d = 0; d < dE; ++d) vcoords[d] = pcoords[d];
2072: } else {
2073: if (localizeVertices) {
2074: PetscScalar anchor[3];
2076: for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
2077: if (!isLocalized) {
2078: Nv = Nc/dE;
2079: for (v = 0; v < Nv; ++v) {DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);}
2080: } else {
2081: Nv = Nc/(2*dE);
2082: for (v = Nv; v < Nv*2; ++v) {DMLocalizeAddCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], vcoords);}
2083: }
2084: } else {
2085: Nv = Nc/dE;
2086: for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) vcoords[d] += pcoords[v*dE+d];
2087: }
2088: for (d = 0; d < dE; ++d) vcoords[d] /= Nv;
2089: }
2090: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
2091: for (n = 0; n < Nct; ++n) {
2092: if (rct[n] != DM_POLYTOPE_POINT) continue;
2093: for (r = 0; r < rsize[n]; ++r) {
2094: PetscInt vNew, off;
2096: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);
2097: PetscSectionGetOffset(coordSectionNew, vNew, &off);
2098: for (d = 0; d < dE; ++d) coordsNew[off+d] = vcoords[d];
2099: }
2100: }
2101: }
2102: }
2103: /* Then set coordinates for cells by localizing */
2104: for (p = pStart; p < pEnd; ++p) {
2105: DMPolytopeType ct;
2106: DMPolytopeType *rct;
2107: PetscInt *rsize, *rcone, *rornt;
2108: PetscInt Nct, n, r;
2109: PetscBool isLocalized = PETSC_FALSE;
2111: DMPlexGetCellType(dm, p, &ct);
2112: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2113: if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2114: PetscInt dof;
2115: PetscSectionGetDof(coordSection, p, &dof);
2116: if (dof) isLocalized = PETSC_TRUE;
2117: }
2118: if (isLocalized) {
2119: const PetscScalar *pcoords;
2121: DMPlexPointLocalRead(cdm, p, coords, &pcoords);
2122: for (n = 0; n < Nct; ++n) {
2123: const PetscInt Nr = rsize[n];
2125: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2126: for (r = 0; r < Nr; ++r) {
2127: PetscInt pNew, offNew;
2129: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2130: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2131: cell to the ones it produces. */
2132: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2133: PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
2134: DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
2135: }
2136: }
2137: }
2138: }
2139: VecRestoreArrayRead(coordsLocal, &coords);
2140: VecRestoreArray(coordsLocalNew, &coordsNew);
2141: DMSetCoordinatesLocal(rdm, coordsLocalNew);
2142: /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
2143: VecDestroy(&coordsLocalNew);
2144: PetscSectionDestroy(&coordSectionNew);
2145: if (!localizeCells) {DMLocalizeCoordinates(rdm);}
2146: return(0);
2147: }
2149: /*@
2150: DMPlexCreateProcessSF - Create an SF which just has process connectivity
2152: Collective on dm
2154: Input Parameters:
2155: + dm - The DM
2156: - sfPoint - The PetscSF which encodes point connectivity
2158: Output Parameters:
2159: + processRanks - A list of process neighbors, or NULL
2160: - sfProcess - An SF encoding the process connectivity, or NULL
2162: Level: developer
2164: .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
2165: @*/
2166: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
2167: {
2168: PetscInt numRoots, numLeaves, l;
2169: const PetscInt *localPoints;
2170: const PetscSFNode *remotePoints;
2171: PetscInt *localPointsNew;
2172: PetscSFNode *remotePointsNew;
2173: PetscInt *ranks, *ranksNew;
2174: PetscMPIInt size;
2175: PetscErrorCode ierr;
2182: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
2183: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2184: PetscMalloc1(numLeaves, &ranks);
2185: for (l = 0; l < numLeaves; ++l) {
2186: ranks[l] = remotePoints[l].rank;
2187: }
2188: PetscSortRemoveDupsInt(&numLeaves, ranks);
2189: PetscMalloc1(numLeaves, &ranksNew);
2190: PetscMalloc1(numLeaves, &localPointsNew);
2191: PetscMalloc1(numLeaves, &remotePointsNew);
2192: for (l = 0; l < numLeaves; ++l) {
2193: ranksNew[l] = ranks[l];
2194: localPointsNew[l] = l;
2195: remotePointsNew[l].index = 0;
2196: remotePointsNew[l].rank = ranksNew[l];
2197: }
2198: PetscFree(ranks);
2199: if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);}
2200: else {PetscFree(ranksNew);}
2201: if (sfProcess) {
2202: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
2203: PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
2204: PetscSFSetFromOptions(*sfProcess);
2205: PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
2206: }
2207: return(0);
2208: }
2210: static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm)
2211: {
2212: DM dm = cr->dm;
2213: DMPlexCellRefiner *crRem;
2214: PetscSF sf, sfNew, sfProcess;
2215: IS processRanks;
2216: MPI_Datatype ctType;
2217: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
2218: const PetscInt *localPoints, *neighbors;
2219: const PetscSFNode *remotePoints;
2220: PetscInt *localPointsNew;
2221: PetscSFNode *remotePointsNew;
2222: PetscInt *ctStartRem, *ctStartNewRem;
2223: PetscInt ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem;
2224: PetscErrorCode ierr;
2227: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
2228: DMGetPointSF(dm, &sf);
2229: DMGetPointSF(rdm, &sfNew);
2230: /* Calculate size of new SF */
2231: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
2232: if (numRoots < 0) return(0);
2233: for (l = 0; l < numLeaves; ++l) {
2234: const PetscInt p = localPoints[l];
2235: DMPolytopeType ct;
2236: DMPolytopeType *rct;
2237: PetscInt *rsize, *rcone, *rornt;
2238: PetscInt Nct, n;
2240: DMPlexGetCellType(dm, p, &ct);
2241: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2242: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
2243: }
2244: /* Communicate ctStart and cStartNew for each remote rank */
2245: DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);
2246: ISGetLocalSize(processRanks, &numNeighbors);
2247: PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);
2248: MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);
2249: MPI_Type_commit(&ctType);
2250: PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);
2251: PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);
2252: PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
2253: PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
2254: MPI_Type_free(&ctType);
2255: PetscSFDestroy(&sfProcess);
2256: PetscMalloc1(numNeighbors, &crRem);
2257: for (n = 0; n < numNeighbors; ++n) {
2258: DMPlexCellRefinerCreate(dm, &crRem[n]);
2259: DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]);
2260: DMPlexCellRefinerSetUp(crRem[n]);
2261: }
2262: PetscFree2(ctStartRem, ctStartNewRem);
2263: /* Calculate new point SF */
2264: PetscMalloc1(numLeavesNew, &localPointsNew);
2265: PetscMalloc1(numLeavesNew, &remotePointsNew);
2266: ISGetIndices(processRanks, &neighbors);
2267: for (l = 0, m = 0; l < numLeaves; ++l) {
2268: PetscInt p = localPoints[l];
2269: PetscInt pRem = remotePoints[l].index;
2270: PetscMPIInt rankRem = remotePoints[l].rank;
2271: DMPolytopeType ct;
2272: DMPolytopeType *rct;
2273: PetscInt *rsize, *rcone, *rornt;
2274: PetscInt neighbor, Nct, n, r;
2276: PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);
2277: if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem);
2278: DMPlexGetCellType(dm, p, &ct);
2279: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2280: for (n = 0; n < Nct; ++n) {
2281: for (r = 0; r < rsize[n]; ++r) {
2282: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2283: DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);
2284: localPointsNew[m] = pNew;
2285: remotePointsNew[m].index = pNewRem;
2286: remotePointsNew[m].rank = rankRem;
2287: ++m;
2288: }
2289: }
2290: }
2291: for (n = 0; n < numNeighbors; ++n) {DMPlexCellRefinerDestroy(&crRem[n]);}
2292: PetscFree(crRem);
2293: if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew);
2294: ISRestoreIndices(processRanks, &neighbors);
2295: ISDestroy(&processRanks);
2296: {
2297: PetscSFNode *rp, *rtmp;
2298: PetscInt *lp, *idx, *ltmp, i;
2300: /* SF needs sorted leaves to correct calculate Gather */
2301: PetscMalloc1(numLeavesNew, &idx);
2302: PetscMalloc1(numLeavesNew, &lp);
2303: PetscMalloc1(numLeavesNew, &rp);
2304: for (i = 0; i < numLeavesNew; ++i) {
2305: if ((localPointsNew[i] < pStartNew) || (localPointsNew[i] >= pEndNew)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %D (%D) not in [%D, %D)", localPointsNew[i], i, pStartNew, pEndNew);
2306: idx[i] = i;
2307: }
2308: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
2309: for (i = 0; i < numLeavesNew; ++i) {
2310: lp[i] = localPointsNew[idx[i]];
2311: rp[i] = remotePointsNew[idx[i]];
2312: }
2313: ltmp = localPointsNew;
2314: localPointsNew = lp;
2315: rtmp = remotePointsNew;
2316: remotePointsNew = rp;
2317: PetscFree(idx);
2318: PetscFree(ltmp);
2319: PetscFree(rtmp);
2320: }
2321: PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
2322: return(0);
2323: }
2325: static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm)
2326: {
2327: DM dm = cr->dm;
2328: PetscInt numLabels, l;
2329: PetscInt pNew;
2333: DMGetNumLabels(dm, &numLabels);
2334: for (l = 0; l < numLabels; ++l) {
2335: DMLabel label, labelNew;
2336: const char *lname;
2337: PetscBool isDepth, isCellType;
2338: IS valueIS;
2339: const PetscInt *values;
2340: PetscInt defVal;
2341: PetscInt numValues, val;
2343: DMGetLabelName(dm, l, &lname);
2344: PetscStrcmp(lname, "depth", &isDepth);
2345: if (isDepth) continue;
2346: PetscStrcmp(lname, "celltype", &isCellType);
2347: if (isCellType) continue;
2348: DMCreateLabel(rdm, lname);
2349: DMGetLabel(dm, lname, &label);
2350: DMGetLabel(rdm, lname, &labelNew);
2351: DMLabelGetDefaultValue(label, &defVal);
2352: DMLabelSetDefaultValue(labelNew, defVal);
2353: DMLabelGetValueIS(label, &valueIS);
2354: ISGetLocalSize(valueIS, &numValues);
2355: ISGetIndices(valueIS, &values);
2356: for (val = 0; val < numValues; ++val) {
2357: IS pointIS;
2358: const PetscInt *points;
2359: PetscInt numPoints, p;
2361: /* Ensure refined label is created with same number of strata as
2362: * original (even if no entries here). */
2363: DMLabelAddStratum(labelNew, values[val]);
2364: DMLabelGetStratumIS(label, values[val], &pointIS);
2365: ISGetLocalSize(pointIS, &numPoints);
2366: ISGetIndices(pointIS, &points);
2367: for (p = 0; p < numPoints; ++p) {
2368: const PetscInt point = points[p];
2369: DMPolytopeType ct;
2370: DMPolytopeType *rct;
2371: PetscInt *rsize, *rcone, *rornt;
2372: PetscInt Nct, n, r;
2374: DMPlexGetCellType(dm, point, &ct);
2375: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2376: for (n = 0; n < Nct; ++n) {
2377: for (r = 0; r < rsize[n]; ++r) {
2378: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);
2379: DMLabelSetValue(labelNew, pNew, values[val]);
2380: }
2381: }
2382: }
2383: ISRestoreIndices(pointIS, &points);
2384: ISDestroy(&pointIS);
2385: }
2386: ISRestoreIndices(valueIS, &values);
2387: ISDestroy(&valueIS);
2388: }
2389: return(0);
2390: }
2392: /* This will only work for interpolated meshes */
2393: PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined)
2394: {
2395: DM rdm;
2396: PetscInt dim, embedDim, depth;
2397: PetscErrorCode ierr;
2401: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
2402: DMSetType(rdm, DMPLEX);
2403: DMGetDimension(dm, &dim);
2404: DMSetDimension(rdm, dim);
2405: DMGetCoordinateDim(dm, &embedDim);
2406: DMSetCoordinateDim(rdm, embedDim);
2407: /* Calculate number of new points of each depth */
2408: DMPlexGetDepth(dm, &depth);
2409: if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement");
2410: /* Step 1: Set chart */
2411: DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);
2412: /* Step 2: Set cone/support sizes (automatically stratifies) */
2413: DMPlexCellRefinerSetConeSizes(cr, rdm);
2414: /* Step 3: Setup refined DM */
2415: DMSetUp(rdm);
2416: /* Step 4: Set cones and supports (automatically symmetrizes) */
2417: DMPlexCellRefinerSetCones(cr, rdm);
2418: /* Step 5: Create pointSF */
2419: DMPlexCellRefinerCreateSF(cr, rdm);
2420: /* Step 6: Create labels */
2421: DMPlexCellRefinerCreateLabels(cr, rdm);
2422: /* Step 7: Set coordinates */
2423: DMPlexCellRefinerSetCoordinates(cr, rdm);
2424: *dmRefined = rdm;
2425: return(0);
2426: }
2428: /*@
2429: DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
2431: Input Parameter:
2432: . dm - The coarse DM
2434: Output Parameter:
2435: . fpointIS - The IS of all the fine points which exist in the original coarse mesh
2437: Level: developer
2439: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexCreateSubpointIS()
2440: @*/
2441: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
2442: {
2443: DMPlexCellRefiner cr;
2444: PetscInt *fpoints;
2445: PetscInt pStart, pEnd, p, vStart, vEnd, v;
2446: PetscErrorCode ierr;
2449: DMPlexGetChart(dm, &pStart, &pEnd);
2450: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2451: DMPlexCellRefinerCreate(dm, &cr);
2452: DMPlexCellRefinerSetUp(cr);
2453: PetscMalloc1(pEnd-pStart, &fpoints);
2454: for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
2455: for (v = vStart; v < vEnd; ++v) {
2456: PetscInt vNew = -1; /* silent overzelous may be used uninitialized */
2458: DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
2459: fpoints[v-pStart] = vNew;
2460: }
2461: DMPlexCellRefinerDestroy(&cr);
2462: ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
2463: return(0);
2464: }
2466: /*@
2467: DMPlexSetRefinementUniform - Set the flag for uniform refinement
2469: Input Parameters:
2470: + dm - The DM
2471: - refinementUniform - The flag for uniform refinement
2473: Level: developer
2475: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
2476: @*/
2477: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
2478: {
2479: DM_Plex *mesh = (DM_Plex*) dm->data;
2483: mesh->refinementUniform = refinementUniform;
2484: return(0);
2485: }
2487: /*@
2488: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
2490: Input Parameter:
2491: . dm - The DM
2493: Output Parameter:
2494: . refinementUniform - The flag for uniform refinement
2496: Level: developer
2498: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
2499: @*/
2500: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
2501: {
2502: DM_Plex *mesh = (DM_Plex*) dm->data;
2507: *refinementUniform = mesh->refinementUniform;
2508: return(0);
2509: }
2511: /*@
2512: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
2514: Input Parameters:
2515: + dm - The DM
2516: - refinementLimit - The maximum cell volume in the refined mesh
2518: Level: developer
2520: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
2521: @*/
2522: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
2523: {
2524: DM_Plex *mesh = (DM_Plex*) dm->data;
2528: mesh->refinementLimit = refinementLimit;
2529: return(0);
2530: }
2532: /*@
2533: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
2535: Input Parameter:
2536: . dm - The DM
2538: Output Parameter:
2539: . refinementLimit - The maximum cell volume in the refined mesh
2541: Level: developer
2543: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
2544: @*/
2545: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
2546: {
2547: DM_Plex *mesh = (DM_Plex*) dm->data;
2552: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
2553: *refinementLimit = mesh->refinementLimit;
2554: return(0);
2555: }
2557: /*@
2558: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
2560: Input Parameters:
2561: + dm - The DM
2562: - refinementFunc - Function giving the maximum cell volume in the refined mesh
2564: Note: The calling sequence is refinementFunc(coords, limit)
2565: $ coords - Coordinates of the current point, usually a cell centroid
2566: $ limit - The maximum cell volume for a cell containing this point
2568: Level: developer
2570: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
2571: @*/
2572: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
2573: {
2574: DM_Plex *mesh = (DM_Plex*) dm->data;
2578: mesh->refinementFunc = refinementFunc;
2579: return(0);
2580: }
2582: /*@
2583: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
2585: Input Parameter:
2586: . dm - The DM
2588: Output Parameter:
2589: . refinementFunc - Function giving the maximum cell volume in the refined mesh
2591: Note: The calling sequence is refinementFunc(coords, limit)
2592: $ coords - Coordinates of the current point, usually a cell centroid
2593: $ limit - The maximum cell volume for a cell containing this point
2595: Level: developer
2597: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
2598: @*/
2599: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
2600: {
2601: DM_Plex *mesh = (DM_Plex*) dm->data;
2606: *refinementFunc = mesh->refinementFunc;
2607: return(0);
2608: }
2610: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
2611: {
2612: PetscBool isUniform;
2613: DMPlexCellRefiner cr;
2614: PetscErrorCode ierr;
2617: DMPlexGetRefinementUniform(dm, &isUniform);
2618: DMViewFromOptions(dm, NULL, "-initref_dm_view");
2619: if (isUniform) {
2620: PetscBool localized;
2622: DMPlexCellRefinerCreate(dm, &cr);
2623: DMPlexCellRefinerSetUp(cr);
2624: DMGetCoordinatesLocalized(dm, &localized);
2625: DMPlexRefineUniform(dm, cr, dmRefined);
2626: DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);
2627: DMCopyBoundary(dm, *dmRefined);
2628: if (localized) {DMLocalizeCoordinates(*dmRefined);}
2629: DMPlexCellRefinerDestroy(&cr);
2630: } else {
2631: DMPlexRefine_Internal(dm, NULL, dmRefined);
2632: }
2633: return(0);
2634: }
2636: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
2637: {
2638: DM cdm = dm;
2639: PetscInt r;
2640: PetscBool isUniform, localized;
2644: DMPlexGetRefinementUniform(dm, &isUniform);
2645: DMGetCoordinatesLocalized(dm, &localized);
2646: if (isUniform) {
2647: for (r = 0; r < nlevels; ++r) {
2648: DMPlexCellRefiner cr;
2650: DMPlexCellRefinerCreate(cdm, &cr);
2651: DMPlexCellRefinerSetUp(cr);
2652: DMPlexRefineUniform(cdm, cr, &dmRefined[r]);
2653: DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);
2654: DMSetRefineLevel(dmRefined[r], cdm->levelup+1);
2655: DMCopyBoundary(cdm, dmRefined[r]);
2656: if (localized) {DMLocalizeCoordinates(dmRefined[r]);}
2657: DMSetCoarseDM(dmRefined[r], cdm);
2658: DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
2659: cdm = dmRefined[r];
2660: DMPlexCellRefinerDestroy(&cr);
2661: }
2662: } else {
2663: for (r = 0; r < nlevels; ++r) {
2664: DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);
2665: DMCopyBoundary(cdm, dmRefined[r]);
2666: if (localized) {DMLocalizeCoordinates(dmRefined[r]);}
2667: DMSetCoarseDM(dmRefined[r], cdm);
2668: cdm = dmRefined[r];
2669: }
2670: }
2671: return(0);
2672: }