Actual source code: plexrefine.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscsf.h>
5: PetscBool SBRcite = PETSC_FALSE;
6: const char SBRCitation[] = "@article{PlazaCarey2000,\n"
7: " title = {Local refinement of simplicial grids based on the skeleton},\n"
8: " journal = {Applied Numerical Mathematics},\n"
9: " author = {A. Plaza and Graham F. Carey},\n"
10: " volume = {32},\n"
11: " number = {3},\n"
12: " pages = {195--218},\n"
13: " doi = {10.1016/S0168-9274(99)00022-7},\n"
14: " year = {2000}\n}\n";
16: const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "Alfeld2D", "Alfeld3D", "PowellSabin", "BoundaryLayer", "SBR", "DMPlexCellRefinerTypes", "DM_REFINER_", NULL};
18: /*
19: Note that j and invj are non-square:
20: v0 + j x_face = x_cell
21: invj (x_cell - v0) = x_face
22: */
23: static PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
24: {
25: /*
26: 2
27: |\
28: | \
29: | \
30: | \
31: | \
32: | \
33: | \
34: 2 1
35: | \
36: | \
37: | \
38: 0---0-------1
39: v0[Nf][dc]: 3 x 2
40: J[Nf][df][dc]: 3 x 1 x 2
41: invJ[Nf][dc][df]: 3 x 2 x 1
42: detJ[Nf]: 3
43: */
44: static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
45: static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
46: static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0};
47: static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0};
48: /*
49: 3---------2---------2
50: | |
51: | |
52: | |
53: 3 1
54: | |
55: | |
56: | |
57: 0---------0---------1
59: v0[Nf][dc]: 4 x 2
60: J[Nf][df][dc]: 4 x 1 x 2
61: invJ[Nf][dc][df]: 4 x 2 x 1
62: detJ[Nf]: 4
63: */
64: static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 -1.0, 0.0};
65: static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
66: static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
67: static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0};
70: switch (ct) {
71: 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;
72: 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;
73: default:
74: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
75: }
76: return(0);
77: }
79: /* Gets the affine map from the original cell to each subcell */
80: static PetscErrorCode DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
81: {
82: /*
83: 2
84: |\
85: | \
86: | \
87: | \
88: | C \
89: | \
90: | \
91: 2---1---1
92: |\ D / \
93: | 2 0 \
94: |A \ / B \
95: 0---0-------1
96: */
97: static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
98: static PetscReal tri_J[] = {0.5, 0.0,
99: 0.0, 0.5,
101: 0.5, 0.0,
102: 0.0, 0.5,
104: 0.5, 0.0,
105: 0.0, 0.5,
107: 0.0, -0.5,
108: 0.5, 0.5};
109: static PetscReal tri_invJ[] = {2.0, 0.0,
110: 0.0, 2.0,
112: 2.0, 0.0,
113: 0.0, 2.0,
115: 2.0, 0.0,
116: 0.0, 2.0,
118: 2.0, 2.0,
119: -2.0, 0.0};
120: /*
121: 3---------2---------2
122: | | |
123: | D 2 C |
124: | | |
125: 3----3----0----1----1
126: | | |
127: | A 0 B |
128: | | |
129: 0---------0---------1
130: */
131: static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
132: static PetscReal quad_J[] = {0.5, 0.0,
133: 0.0, 0.5,
135: 0.5, 0.0,
136: 0.0, 0.5,
138: 0.5, 0.0,
139: 0.0, 0.5,
141: 0.5, 0.0,
142: 0.0, 0.5};
143: static PetscReal quad_invJ[] = {2.0, 0.0,
144: 0.0, 2.0,
146: 2.0, 0.0,
147: 0.0, 2.0,
149: 2.0, 0.0,
150: 0.0, 2.0,
152: 2.0, 0.0,
153: 0.0, 2.0};
154: /*
155: Bottom (viewed from top) Top
156: 1---------2---------2 7---------2---------6
157: | | | | | |
158: | B 2 C | | H 2 G |
159: | | | | | |
160: 3----3----0----1----1 3----3----0----1----1
161: | | | | | |
162: | A 0 D | | E 0 F |
163: | | | | | |
164: 0---------0---------3 4---------0---------5
165: */
166: 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,
167: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
168: static PetscReal hex_J[] = {0.5, 0.0, 0.0,
169: 0.0, 0.5, 0.0,
170: 0.0, 0.0, 0.5,
172: 0.5, 0.0, 0.0,
173: 0.0, 0.5, 0.0,
174: 0.0, 0.0, 0.5,
176: 0.5, 0.0, 0.0,
177: 0.0, 0.5, 0.0,
178: 0.0, 0.0, 0.5,
180: 0.5, 0.0, 0.0,
181: 0.0, 0.5, 0.0,
182: 0.0, 0.0, 0.5,
184: 0.5, 0.0, 0.0,
185: 0.0, 0.5, 0.0,
186: 0.0, 0.0, 0.5,
188: 0.5, 0.0, 0.0,
189: 0.0, 0.5, 0.0,
190: 0.0, 0.0, 0.5,
192: 0.5, 0.0, 0.0,
193: 0.0, 0.5, 0.0,
194: 0.0, 0.0, 0.5,
196: 0.5, 0.0, 0.0,
197: 0.0, 0.5, 0.0,
198: 0.0, 0.0, 0.5};
199: static PetscReal hex_invJ[] = {2.0, 0.0, 0.0,
200: 0.0, 2.0, 0.0,
201: 0.0, 0.0, 2.0,
203: 2.0, 0.0, 0.0,
204: 0.0, 2.0, 0.0,
205: 0.0, 0.0, 2.0,
207: 2.0, 0.0, 0.0,
208: 0.0, 2.0, 0.0,
209: 0.0, 0.0, 2.0,
211: 2.0, 0.0, 0.0,
212: 0.0, 2.0, 0.0,
213: 0.0, 0.0, 2.0,
215: 2.0, 0.0, 0.0,
216: 0.0, 2.0, 0.0,
217: 0.0, 0.0, 2.0,
219: 2.0, 0.0, 0.0,
220: 0.0, 2.0, 0.0,
221: 0.0, 0.0, 2.0,
223: 2.0, 0.0, 0.0,
224: 0.0, 2.0, 0.0,
225: 0.0, 0.0, 2.0,
227: 2.0, 0.0, 0.0,
228: 0.0, 2.0, 0.0,
229: 0.0, 0.0, 2.0};
232: switch (ct) {
233: case DM_POLYTOPE_TRIANGLE: if (Nc) *Nc = 4; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; break;
234: case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break;
235: case DM_POLYTOPE_HEXAHEDRON: if (Nc) *Nc = 8; if (v0) *v0 = hex_v0; if (J) *J = hex_J; if (invJ) *invJ = hex_invJ; break;
236: default:
237: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
238: }
239: return(0);
240: }
242: /* Should this be here or in the DualSpace somehow? */
243: PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside)
244: {
245: PetscReal sum = 0.0;
246: PetscInt d;
249: *inside = PETSC_TRUE;
250: switch (ct) {
251: case DM_POLYTOPE_TRIANGLE:
252: case DM_POLYTOPE_TETRAHEDRON:
253: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) {
254: if (point[d] < -1.0) {*inside = PETSC_FALSE; break;}
255: sum += point[d];
256: }
257: if (sum > PETSC_SMALL) {*inside = PETSC_FALSE; break;}
258: break;
259: case DM_POLYTOPE_QUADRILATERAL:
260: case DM_POLYTOPE_HEXAHEDRON:
261: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d)
262: if (PetscAbsReal(point[d]) > 1.+PETSC_SMALL) {*inside = PETSC_FALSE; break;}
263: break;
264: default:
265: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
266: }
267: return(0);
268: }
270: /* Regular Refinment of Hybrid Meshes
272: We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
273: to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
274: transformations, such as changing from one type of cell to another, as simplex to hex.
276: To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
277: and the types of the new cells.
279: We need the group multiplication table for group actions from the dihedral group for each cell type.
281: We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
282: we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
283: (face, orient) pairs for each subcell.
284: */
286: /*
287: Input Parameters:
288: + ct - The type of the input cell
289: . co - The orientation of this cell in its parent
290: - cp - The requested cone point of this cell assuming orientation 0
292: Output Parameters:
293: . cpnew - The new cone point, taking into account the orientation co
294: */
295: PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew)
296: {
297: const PetscInt csize = DMPolytopeTypeGetConeSize(ct);
300: if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;}
301: else {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;}
302: return(0);
303: }
305: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
306: {
307: static PetscReal seg_v[] = {-1.0, 0.0, 1.0};
308: 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};
309: 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};
310: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
311: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
312: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0};
313: static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
314: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0,
315: -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0,
316: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0,
317: -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
318: -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0,
319: -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0,
320: -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
321: -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0};
324: switch (ct) {
325: case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break;
326: case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break;
327: case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break;
328: case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break;
329: case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break;
330: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
331: }
332: return(0);
333: }
335: static PetscErrorCode DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
336: {
337: 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};
338: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
339: -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,
340: -1.0, -1.0, 0.0, -1.0/3.0, -1.0, -1.0/3.0, 0.0, -1.0, 0.0,
341: -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,
342: -1.0, -1.0, 1.0, -0.5, -0.5, -0.5};
343: PetscErrorCode ierr;
346: switch (ct) {
347: case DM_POLYTOPE_SEGMENT:
348: case DM_POLYTOPE_QUADRILATERAL:
349: case DM_POLYTOPE_HEXAHEDRON:
350: DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);break;
351: case DM_POLYTOPE_TRIANGLE: *Nv = 7; *subcellV = tri_v; break;
352: case DM_POLYTOPE_TETRAHEDRON: *Nv = 15; *subcellV = tet_v; break;
353: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
354: }
355: return(0);
356: }
358: /*
359: DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type
361: Input Parameters:
362: + cr - The DMPlexCellRefiner object
363: - ct - The cell type
365: Output Parameters:
366: + Nv - The number of refined vertices in the closure of the reference cell of given type
367: - subcellV - The coordinates of these vertices in the reference cell
369: Level: developer
371: .seealso: DMPlexCellRefinerGetSubcellVertices()
372: */
373: static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
374: {
378: if (!cr->ops->getcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
379: (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);
380: return(0);
381: }
383: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
384: {
385: static PetscInt seg_v[] = {0, 1, 1, 2};
386: static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5};
387: static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3};
388: static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9,
389: 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6};
390: 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,
391: 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};
394: if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
395: switch (ct) {
396: case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break;
397: case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break;
398: case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
399: case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break;
400: case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break;
401: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
402: }
403: return(0);
404: }
406: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
407: {
408: static PetscInt tri_v[] = {0, 3, 6, 5, 3, 1, 4, 6, 5, 6, 4, 2};
409: 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};
410: PetscErrorCode ierr;
413: switch (ct) {
414: case DM_POLYTOPE_SEGMENT:
415: case DM_POLYTOPE_QUADRILATERAL:
416: case DM_POLYTOPE_HEXAHEDRON:
417: DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);break;
418: case DM_POLYTOPE_TRIANGLE:
419: if (rct != DM_POLYTOPE_QUADRILATERAL) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
420: *Nv = 4; *subcellV = &tri_v[r*(*Nv)]; break;
421: case DM_POLYTOPE_TETRAHEDRON:
422: if (rct != DM_POLYTOPE_HEXAHEDRON) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
423: *Nv = 8; *subcellV = &tet_v[r*(*Nv)]; break;
424: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
425: }
426: return(0);
427: }
429: /*
430: DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type
432: Input Parameters:
433: + cr - The DMPlexCellRefiner object
434: . ct - The cell type
435: . rct - The type of subcell
436: - r - The subcell index
438: Output Parameters:
439: + Nv - The number of refined vertices in the subcell
440: - subcellV - The indices of these vertices in the set of vertices returned by DMPlexCellRefinerGetCellVertices()
442: Level: developer
444: .seealso: DMPlexCellRefinerGetCellVertices()
445: */
446: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
447: {
451: if (!cr->ops->getsubcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
452: (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);
453: return(0);
454: }
456: /*
457: Input Parameters:
458: + cr - The DMPlexCellRefiner
459: . pct - The cell type of the parent, from whom the new cell is being produced
460: . ct - The type being produced
461: . r - The replica number requested for the produced cell type
462: . Nv - Number of vertices in the closure of the parent cell
463: . dE - Spatial dimension
464: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
466: Output Parameters:
467: . out - The coordinates of the new vertices
468: */
469: static PetscErrorCode DMPlexCellRefinerMapCoordinates(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
470: {
474: if (!cr->ops->mapcoords) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
475: (*cr->ops->mapcoords)(cr, pct, ct, r, Nv, dE, in, out);
476: return(0);
477: }
479: /* Computes new vertex as the barycenter */
480: static PetscErrorCode DMPlexCellRefinerMapCoordinates_Barycenter(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
481: {
482: PetscInt v,d;
485: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
486: for (d = 0; d < dE; ++d) out[d] = 0.0;
487: for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
488: for (d = 0; d < dE; ++d) out[d] /= Nv;
489: return(0);
490: }
492: /*
493: Input Parameters:
494: + cr - The DMPlexCellRefiner
495: . pct - The cell type of the parent, from whom the new cell is being produced
496: . pp - The parent cell
497: . po - The orientation of the parent cell in its enclosing parent
498: . ct - The type being produced
499: . r - The replica number requested for the produced cell type
500: - o - The relative orientation of the replica
502: Output Parameters:
503: + rnew - The replica number, given the orientation of the parent
504: - onew - The replica orientation, given the orientation of the parent
505: */
506: static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
507: {
511: if (!cr->ops->mapsubcells) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
512: (*cr->ops->mapsubcells)(cr, pct, pp, po, ct, r, o, rnew, onew);
513: return(0);
514: }
516: /*
517: This is the group multiplication table for the dihedral group of the cell.
518: */
519: static PetscErrorCode ComposeOrientation_Private(PetscInt n, PetscInt o1, PetscInt o2, PetscInt *o)
520: {
522: if (!n) {*o = 0;}
523: else if (o1 >= 0 && o2 >= 0) {*o = ( o1 + o2) % n;}
524: else if (o1 < 0 && o2 < 0) {*o = (-o1 - o2) % n;}
525: else if (o1 < 0) {*o = -((-(o1+1) + o2) % n + 1);}
526: else if (o2 < 0) {*o = -((-(o2+1) + o1) % n + 1);}
527: return(0);
528: }
530: static PetscErrorCode DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
531: {
535: *rnew = r;
536: ComposeOrientation_Private(DMPolytopeTypeGetConeSize(ct), po, o, onew);
537: return(0);
538: }
540: static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
541: {
542: /* We shift any input orientation in order to make it non-negative
543: 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)
544: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
545: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
546: */
547: PetscInt tri_seg_o[] = {-2, 0,
548: -2, 0,
549: -2, 0,
550: 0, -2,
551: 0, -2,
552: 0, -2};
553: PetscInt tri_seg_r[] = {1, 0, 2,
554: 0, 2, 1,
555: 2, 1, 0,
556: 0, 1, 2,
557: 1, 2, 0,
558: 2, 0, 1};
559: PetscInt tri_tri_o[] = {0, 1, 2, -3, -2, -1,
560: 2, 0, 1, -2, -1, -3,
561: 1, 2, 0, -1, -3, -2,
562: -3, -2, -1, 0, 1, 2,
563: -1, -3, -2, 1, 2, 0,
564: -2, -1, -3, 2, 0, 1};
565: /* orientation if the replica is the center triangle */
566: PetscInt tri_tri_o_c[] = {2, 0, 1, -2, -1, -3,
567: 1, 2, 0, -1, -3, -2,
568: 0, 1, 2, -3, -2, -1,
569: -3, -2, -1, 0, 1, 2,
570: -1, -3, -2, 1, 2, 0,
571: -2, -1, -3, 2, 0, 1};
572: PetscInt tri_tri_r[] = {0, 2, 1, 3,
573: 2, 1, 0, 3,
574: 1, 0, 2, 3,
575: 0, 1, 2, 3,
576: 1, 2, 0, 3,
577: 2, 0, 1, 3};
578: PetscInt quad_seg_r[] = {3, 2, 1, 0,
579: 2, 1, 0, 3,
580: 1, 0, 3, 2,
581: 0, 3, 2, 1,
582: 0, 1, 2, 3,
583: 1, 2, 3, 0,
584: 2, 3, 0, 1,
585: 3, 0, 1, 2};
586: PetscInt quad_quad_o[] = { 0, 1, 2, 3, -4, -3, -2, -1,
587: 4, 0, 1, 2, -3, -2, -1, -4,
588: 3, 4, 0, 1, -2, -1, -4, -3,
589: 2, 3, 4, 0, -1, -4, -3, -2,
590: -4, -3, -2, -1, 0, 1, 2, 3,
591: -1, -4, -3, -2, 1, 2, 3, 0,
592: -2, -1, -4, -3, 2, 3, 0, 1,
593: -3, -2, -1, -4, 3, 0, 1, 2};
594: PetscInt quad_quad_r[] = {0, 3, 2, 1,
595: 3, 2, 1, 0,
596: 2, 1, 0, 3,
597: 1, 0, 3, 2,
598: 0, 1, 2, 3,
599: 1, 2, 3, 0,
600: 2, 3, 0, 1,
601: 3, 0, 1, 2};
602: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
603: 1, 0, -1, -2,
604: -2, -1, 0, 1,
605: -1, -2, 1, 0};
606: PetscInt tquad_tquad_r[] = {1, 0,
607: 1, 0,
608: 0, 1,
609: 0, 1};
612: /* The default is no transformation */
613: *rnew = r;
614: *onew = o;
615: switch (pct) {
616: case DM_POLYTOPE_SEGMENT:
617: if (ct == DM_POLYTOPE_SEGMENT) {
618: if (po == 0 || po == -1) {*rnew = r; *onew = o;}
619: else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;}
620: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po);
621: }
622: break;
623: case DM_POLYTOPE_TRIANGLE:
624: switch (ct) {
625: case DM_POLYTOPE_SEGMENT:
626: if (o == -1) o = 0;
627: if (o == -2) o = 1;
628: *onew = tri_seg_o[(po+3)*2+o];
629: *rnew = tri_seg_r[(po+3)*3+r];
630: break;
631: case DM_POLYTOPE_TRIANGLE:
632: *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3];
633: *rnew = tri_tri_r[(po+3)*4+r];
634: break;
635: default: break;
636: }
637: break;
638: case DM_POLYTOPE_QUADRILATERAL:
639: switch (ct) {
640: case DM_POLYTOPE_SEGMENT:
641: *onew = o;
642: *rnew = quad_seg_r[(po+4)*4+r];
643: break;
644: case DM_POLYTOPE_QUADRILATERAL:
645: *onew = quad_quad_o[(po+4)*8+o+4];
646: *rnew = quad_quad_r[(po+4)*4+r];
647: break;
648: default: break;
649: }
650: break;
651: case DM_POLYTOPE_SEG_PRISM_TENSOR:
652: switch (ct) {
653: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
654: case DM_POLYTOPE_SEG_PRISM_TENSOR:
655: *onew = tquad_tquad_o[(po+2)*4+o+2];
656: *rnew = tquad_tquad_r[(po+2)*2+r];
657: break;
658: default: break;
659: }
660: break;
661: default: break;
662: }
663: return(0);
664: }
666: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
667: {
669: /* We shift any input orientation in order to make it non-negative
670: 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)
671: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
672: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
673: */
674: PetscInt tri_seg_o[] = {0, -2,
675: 0, -2,
676: 0, -2,
677: 0, -2,
678: 0, -2,
679: 0, -2};
680: PetscInt tri_seg_r[] = {2, 1, 0,
681: 1, 0, 2,
682: 0, 2, 1,
683: 0, 1, 2,
684: 1, 2, 0,
685: 2, 0, 1};
686: PetscInt tri_tri_o[] = {0, 1, 2, 3, -4, -3, -2, -1,
687: 3, 0, 1, 2, -3, -2, -1, -4,
688: 1, 2, 3, 0, -1, -4, -3, -2,
689: -4, -3, -2, -1, 0, 1, 2, 3,
690: -1, -4, -3, -2, 1, 2, 3, 0,
691: -3, -2, -1, -4, 3, 0, 1, 2};
692: PetscInt tri_tri_r[] = {0, 2, 1,
693: 2, 1, 0,
694: 1, 0, 2,
695: 0, 1, 2,
696: 1, 2, 0,
697: 2, 0, 1};
698: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
699: 1, 0, -1, -2,
700: -2, -1, 0, 1,
701: -1, -2, 1, 0};
702: PetscInt tquad_tquad_r[] = {1, 0,
703: 1, 0,
704: 0, 1,
705: 0, 1};
706: PetscInt tquad_quad_o[] = {-2, -1, -4, -3, 2, 3, 0, 1,
707: 1, 2, 3, 0, -1, -4, -3, -2,
708: -4, -3, -2, -1, 0, 1, 2, 3,
709: 1, 0, 3, 2, -3, -4, -1, -2};
712: *rnew = r;
713: *onew = o;
714: switch (pct) {
715: case DM_POLYTOPE_TRIANGLE:
716: switch (ct) {
717: case DM_POLYTOPE_SEGMENT:
718: if (o == -1) o = 0;
719: if (o == -2) o = 1;
720: *onew = tri_seg_o[(po+3)*2+o];
721: *rnew = tri_seg_r[(po+3)*3+r];
722: break;
723: case DM_POLYTOPE_QUADRILATERAL:
724: *onew = tri_tri_o[(po+3)*8+o+4];
725: /* TODO See sheet, for fixing po == 1 and 2 */
726: if (po == 2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
727: if (po == 2 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
728: if (po == 1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
729: if (po == 1 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
730: if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
731: if (po == -1 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
732: if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
733: if (po == -2 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
734: *rnew = tri_tri_r[(po+3)*3+r];
735: break;
736: default: break;
737: }
738: break;
739: case DM_POLYTOPE_SEG_PRISM_TENSOR:
740: switch (ct) {
741: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
742: case DM_POLYTOPE_SEG_PRISM_TENSOR:
743: *onew = tquad_tquad_o[(po+2)*4+o+2];
744: *rnew = tquad_tquad_r[(po+2)*2+r];
745: break;
746: case DM_POLYTOPE_QUADRILATERAL:
747: *onew = tquad_quad_o[(po+2)*8+o+4];
748: *rnew = tquad_tquad_r[(po+2)*2+r];
749: break;
750: default: break;
751: }
752: break;
753: default:
754: DMPlexCellRefinerMapSubcells_Regular(cr, pct, pp, po, ct, r, o, rnew, onew);
755: }
756: return(0);
757: }
759: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
760: {
761: return DMPlexCellRefinerMapSubcells_Regular(cr, pct, pp, po, ct, r, o, rnew, onew);
762: }
764: /*@
765: DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type
767: Input Parameters:
768: + source - The cell type for a source point
769: - p - The source point, or PETSC_DETERMINE if the refine is homogeneous
771: Output Parameters:
772: + rt - The refine type for this cell
773: . Nt - The number of cell types generated by refinement
774: . target - The cell types generated
775: . size - The number of subcells of each type, ordered by dimension
776: . cone - A list of the faces for each subcell of the same type as source
777: - ornt - A list of the face orientations for each subcell of the same type as source
779: Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
780: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
781: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
782: $ the number of cones to be taken, or 0 for the current cell
783: $ the cell cone point number at each level from which it is subdivided
784: $ the replica number r of the subdivision.
785: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
786: $ Nt = 2
787: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
788: $ size = {1, 2}
789: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
790: $ ornt = { 0, 0, 0, 0}
792: Level: developer
794: .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform()
795: @*/
796: PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
797: {
801: if (!cr->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
802: (*cr->ops->refine)(cr, source, p, rt, Nt, target, size, cone, ornt);
803: return(0);
804: }
806: static PetscErrorCode DMPlexCellRefinerRefine_None(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
807: {
808: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
809: static PetscInt vertexS[] = {1};
810: static PetscInt vertexC[] = {0};
811: static PetscInt vertexO[] = {0};
812: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
813: static PetscInt edgeS[] = {1};
814: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
815: static PetscInt edgeO[] = {0, 0};
816: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
817: static PetscInt tedgeS[] = {1};
818: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
819: static PetscInt tedgeO[] = {0, 0};
820: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
821: static PetscInt triS[] = {1};
822: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
823: static PetscInt triO[] = {0, 0, 0};
824: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
825: static PetscInt quadS[] = {1};
826: static PetscInt quadC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
827: static PetscInt quadO[] = {0, 0, 0, 0};
828: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
829: static PetscInt tquadS[] = {1};
830: static PetscInt tquadC[] = {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};
831: static PetscInt tquadO[] = {0, 0, 0, 0};
832: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
833: static PetscInt tetS[] = {1};
834: static PetscInt tetC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0};
835: static PetscInt tetO[] = {0, 0, 0, 0};
836: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
837: static PetscInt hexS[] = {1};
838: static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
839: DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
840: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
841: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
842: static PetscInt tripS[] = {1};
843: static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
844: DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
845: static PetscInt tripO[] = {0, 0, 0, 0, 0};
846: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
847: static PetscInt ttripS[] = {1};
848: static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
849: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
850: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
851: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
852: static PetscInt tquadpS[] = {1};
853: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
854: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
855: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
856: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
857: static PetscInt pyrS[] = {1};
858: static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
859: DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
860: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
863: if (rt) *rt = 0;
864: switch (source) {
865: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
866: case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
867: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
868: case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
869: case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
870: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
871: case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
872: case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
873: case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
874: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
875: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
876: case DM_POLYTOPE_PYRAMID: *Nt = 1; *target = pyrT; *size = pyrS; *cone = pyrC; *ornt = pyrO; break;
877: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
878: }
879: return(0);
880: }
882: static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
883: {
884: /* All vertices remain in the refined mesh */
885: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
886: static PetscInt vertexS[] = {1};
887: static PetscInt vertexC[] = {0};
888: static PetscInt vertexO[] = {0};
889: /* Split all edges with a new vertex, making two new 2 edges
890: 0--0--0--1--1
891: */
892: static DMPolytopeType edgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
893: static PetscInt edgeS[] = {1, 2};
894: 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};
895: static PetscInt edgeO[] = { 0, 0, 0, 0};
896: /* Do not split tensor edges */
897: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
898: static PetscInt tedgeS[] = {1};
899: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
900: static PetscInt tedgeO[] = { 0, 0};
901: /* Add 3 edges inside every triangle, making 4 new triangles.
902: 2
903: |\
904: | \
905: | \
906: 0 1
907: | C \
908: | \
909: | \
910: 2---1---1
911: |\ D / \
912: 1 2 0 0
913: |A \ / B \
914: 0-0-0---1---1
915: */
916: static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
917: static PetscInt triS[] = {3, 4};
918: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
919: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
920: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
921: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
922: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
923: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
924: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
925: static PetscInt triO[] = {0, 0,
926: 0, 0,
927: 0, 0,
928: 0, -2, 0,
929: 0, 0, -2,
930: -2, 0, 0,
931: 0, 0, 0};
932: /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
933: 3----1----2----0----2
934: | | |
935: 0 D 2 C 1
936: | | |
937: 3----3----0----1----1
938: | | |
939: 1 A 0 B 0
940: | | |
941: 0----0----0----1----1
942: */
943: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
944: static PetscInt quadS[] = {1, 4, 4};
945: static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
946: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
947: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
948: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
949: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1,
950: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
951: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2,
952: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
953: static PetscInt quadO[] = {0, 0,
954: 0, 0,
955: 0, 0,
956: 0, 0,
957: 0, 0, -2, 0,
958: 0, 0, 0, -2,
959: -2, 0, 0, 0,
960: 0, -2, 0, 0};
961: /* Add 1 edge inside every tensor quad, making 2 new tensor quads
962: 2----2----1----3----3
963: | | |
964: | | |
965: | | |
966: 4 A 6 B 5
967: | | |
968: | | |
969: | | |
970: 0----0----0----1----1
971: */
972: static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
973: static PetscInt tquadS[] = {1, 2};
974: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
975: 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,
976: 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};
977: static PetscInt tquadO[] = {0, 0,
978: 0, 0, 0, 0,
979: 0, 0, 0, 0};
980: /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
981: 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
982: 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]
983: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
984: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
985: The first four tets just cut off the corners, using the replica notation for new vertices,
986: [v0, (e0, 0), (e2, 0), (e3, 0)]
987: [(e0, 0), v1, (e1, 0), (e4, 0)]
988: [(e2, 0), (e1, 0), v2, (e5, 0)]
989: [(e3, 0), (e4, 0), (e5, 0), v3 ]
990: The next four tets match a vertex to the newly created faces from cutting off those first tets.
991: [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
992: [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
993: [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
994: [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
995: 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
996: [(e2, 0), (e0, 0), (e3, 0)]
997: [(e0, 0), (e1, 0), (e4, 0)]
998: [(e2, 0), (e5, 0), (e1, 0)]
999: [(e3, 0), (e4, 0), (e5, 0)]
1000: The next four, from the second group of tets, are
1001: [(e2, 0), (e0, 0), (e5, 0)]
1002: [(e4, 0), (e0, 0), (e5, 0)]
1003: [(e0, 0), (e1, 0), (e5, 0)]
1004: [(e5, 0), (e3, 0), (e0, 0)]
1005: I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
1006: */
1007: static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
1008: static PetscInt tetS[] = {1, 8, 8};
1009: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0,
1010: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1011: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1,
1012: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1013: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1014: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1015: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1,
1016: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1017: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1018: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0,
1019: DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1,
1020: DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0,
1021: DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2,
1022: DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7,
1023: DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6,
1024: DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3,
1025: DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
1026: static PetscInt tetO[] = {0, 0,
1027: 0, 0, 0,
1028: 0, 0, 0,
1029: 0, 0, 0,
1030: 0, 0, 0,
1031: 0, 0, -2,
1032: 0, 0, -2,
1033: 0, -2, -2,
1034: 0, -2, 0,
1035: 0, 0, 0, 0,
1036: 0, 0, 0, 0,
1037: 0, 0, 0, 0,
1038: 0, 0, 0, 0,
1039: -3, 0, 0, -2,
1040: -2, 1, 0, 0,
1041: -2, -2, -1, 2,
1042: -2, 0, -2, 1};
1043: /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
1044: 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
1045: [v0, v1, v2, v3] f0 bottom
1046: [v4, v5, v6, v7] f1 top
1047: [v0, v3, v5, v4] f2 front
1048: [v2, v1, v7, v6] f3 back
1049: [v3, v2, v6, v5] f4 right
1050: [v0, v4, v7, v1] f5 left
1051: The eight hexes are divided into four on the bottom, and four on the top,
1052: [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1053: [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1054: [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1055: [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1056: [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)]
1057: [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)]
1058: [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)]
1059: [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7]
1060: 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,
1061: [(e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1062: [(f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1063: [(c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1064: [(f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1065: and on the x-z plane,
1066: [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
1067: [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
1068: [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
1069: [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
1070: and on the y-z plane,
1071: [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
1072: [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
1073: [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
1074: [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
1075: */
1076: static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1077: static PetscInt hexS[] = {1, 6, 12, 8};
1078: static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1079: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1080: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1081: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1082: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1083: DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0,
1084: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0,
1085: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2,
1086: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3,
1087: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2,
1088: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0,
1089: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1,
1090: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1091: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1092: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3,
1093: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1094: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1095: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1096: 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,
1097: 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,
1098: 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,
1099: 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,
1100: 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,
1101: 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,
1102: 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,
1103: 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};
1104: static PetscInt hexO[] = {0, 0,
1105: 0, 0,
1106: 0, 0,
1107: 0, 0,
1108: 0, 0,
1109: 0, 0,
1110: 0, 0, -2, -2,
1111: 0, -2, -2, 0,
1112: -2, -2, 0, 0,
1113: -2, 0, 0, -2,
1114: -2, 0, 0, -2,
1115: -2, -2, 0, 0,
1116: 0, -2, -2, 0,
1117: 0, 0, -2, -2,
1118: 0, 0, -2, -2,
1119: -2, 0, 0, -2,
1120: -2, -2, 0, 0,
1121: 0, -2, -2, 0,
1122: 0, 0, 0, 0, -4, 0,
1123: 0, 0, -1, 0, -4, 0,
1124: 0, 0, -1, 0, 0, 0,
1125: 0, 0, 0, 0, 0, 0,
1126: -4, 0, 0, 0, -4, 0,
1127: -4, 0, 0, 0, 0, 0,
1128: -4, 0, -1, 0, 0, 0,
1129: -4, 0, -1, 0, -4, 0};
1130: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1131: static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
1132: static PetscInt tripS[] = {3, 4, 6, 8};
1133: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0,
1134: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0,
1135: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0,
1136: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1137: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0,
1138: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1139: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1140: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1141: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1142: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1143: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1144: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1145: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1146: 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,
1147: 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,
1148: 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,
1149: 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,
1150: 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,
1151: 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,
1152: 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,
1153: 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};
1154: static PetscInt tripO[] = {0, 0,
1155: 0, 0,
1156: 0, 0,
1157: 0, -2, -2,
1158: -2, 0, -2,
1159: -2, -2, 0,
1160: 0, 0, 0,
1161: -2, 0, -2, -2,
1162: -2, 0, -2, -2,
1163: -2, 0, -2, -2,
1164: 0, -2, -2, 0,
1165: 0, -2, -2, 0,
1166: 0, -2, -2, 0,
1167: 0, 0, 0, -1, 0,
1168: 0, 0, 0, 0, -1,
1169: 0, 0, -1, 0, 0,
1170: 2, 0, 0, 0, 0,
1171: -3, 0, 0, -1, 0,
1172: -3, 0, 0, 0, -1,
1173: -3, 0, -1, 0, 0,
1174: -3, 0, 0, 0, 0};
1175: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
1176: 2
1177: |\
1178: | \
1179: | \
1180: 0---1
1182: 2
1184: 0 1
1186: 2
1187: |\
1188: | \
1189: | \
1190: 0---1
1191: */
1192: static DMPolytopeType ttripT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
1193: static PetscInt ttripS[] = {3, 4};
1194: 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,
1195: 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,
1196: 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,
1197: 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,
1198: 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,
1199: 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,
1200: 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};
1201: static PetscInt ttripO[] = {0, 0, 0, 0,
1202: 0, 0, 0, 0,
1203: 0, 0, 0, 0,
1204: 0, 0, 0, -1, 0,
1205: 0, 0, 0, 0, -1,
1206: 0, 0, -1, 0, 0,
1207: 0, 0, 0, 0, 0};
1208: /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
1209: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1210: static PetscInt tquadpS[] = {1, 4, 4};
1211: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1212: 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,
1213: 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,
1214: 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,
1215: 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,
1216: 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,
1217: 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,
1218: 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,
1219: 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};
1220: static PetscInt tquadpO[] = {0, 0,
1221: 0, 0, 0, 0,
1222: 0, 0, 0, 0,
1223: 0, 0, 0, 0,
1224: 0, 0, 0, 0,
1225: 0, 0, 0, 0, -1, 0,
1226: 0, 0, 0, 0, 0, -1,
1227: 0, 0, -1, 0, 0, 0,
1228: 0, 0, 0, -1, 0, 0};
1232: if (rt) *rt = 0;
1233: switch (source) {
1234: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
1235: case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
1236: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1237: case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1238: case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
1239: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1240: case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1241: case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
1242: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1243: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1244: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1245: /* TODO Fix pyramids: For now, we just ignore them */
1246: case DM_POLYTOPE_PYRAMID:
1247: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1248: break;
1249: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1250: }
1251: return(0);
1252: }
1254: static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1255: {
1257: /* Change tensor edges to segments */
1258: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT};
1259: static PetscInt tedgeS[] = {1};
1260: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1261: static PetscInt tedgeO[] = { 0, 0};
1262: /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals.
1263: 2
1264: |\
1265: | \
1266: | \
1267: | \
1268: 0 1
1269: | \
1270: | \
1271: 2 1
1272: |\ / \
1273: | 2 1 \
1274: | \ / \
1275: 1 | 0
1276: | 0 \
1277: | | \
1278: | | \
1279: 0-0-0-----1-----1
1280: */
1281: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1282: static PetscInt triS[] = {1, 3, 3};
1283: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1284: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1285: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1286: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1287: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1288: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1289: static PetscInt triO[] = {0, 0,
1290: 0, 0,
1291: 0, 0,
1292: 0, 0, -2, 0,
1293: 0, 0, 0, -2,
1294: 0, -2, 0, 0};
1295: /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals
1296: 2----2----1----3----3
1297: | | |
1298: | | |
1299: | | |
1300: 4 A 6 B 5
1301: | | |
1302: | | |
1303: | | |
1304: 0----0----0----1----1
1305: */
1306: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1307: static PetscInt tquadS[] = {1, 2};
1308: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1309: /* TODO Fix these */
1310: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1311: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0};
1312: static PetscInt tquadO[] = {0, 0,
1313: 0, 0, -2, -2,
1314: 0, 0, -2, -2};
1315: /* Add 6 triangles inside every cell, making 4 new hexs
1316: TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it
1317: 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
1318: 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]
1319: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1320: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1321: We make a new hex in each corner
1322: [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1323: [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)]
1324: [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)]
1325: [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)]
1326: We create a new face for each edge
1327: [(e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1328: [(f0, 0), (e0, 0), (f1, 0), (c0, 0)]
1329: [(e2, 0), (f0, 0), (c0, 0), (f2, 0)]
1330: [(f3, 0), (e4, 0), (f1, 0), (c0, 0)]
1331: [(e1, 0), (f3, 0), (c0, 0), (f0, 0)]
1332: [(e5, 0), (f3, 0), (c0, 0), (f2, 0)]
1333: I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another.
1334: */
1335: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1336: static PetscInt tetS[] = {1, 4, 6, 4};
1337: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1338: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1339: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1340: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1341: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0,
1342: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1343: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1344: DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3,
1345: DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1346: DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1347: 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,
1348: 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,
1349: 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,
1350: 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};
1351: static PetscInt tetO[] = {0, 0,
1352: 0, 0,
1353: 0, 0,
1354: 0, 0,
1355: 0, 0, -2, -2,
1356: -2, 0, 0, -2,
1357: 0, 0, -2, -2,
1358: -2, 0, 0, -2,
1359: 0, 0, -2, -2,
1360: 0, 0, -2, -2,
1361: 0, 0, 0, 0, 0, 0,
1362: 1, -1, 1, 0, 0, 3,
1363: 0, -4, 1, -1, 0, 3,
1364: 1, -4, 3, -2, -4, 3};
1365: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1366: static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1367: static PetscInt tripS[] = {1, 5, 9, 6};
1368: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1369: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1370: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1371: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1372: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1373: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1374: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2,
1375: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1376: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1377: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1378: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1379: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1380: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1381: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1382: 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,
1383: 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,
1384: 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,
1385: 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,
1386: 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,
1387: 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};
1388: static PetscInt tripO[] = {0, 0,
1389: 0, 0,
1390: 0, 0,
1391: 0, 0,
1392: 0, 0,
1393: 0, 0, -2, -2,
1394: -2, 0, 0, -2,
1395: 0, -2, -2, 0,
1396: 0, 0, -2, -2,
1397: 0, 0, -2, -2,
1398: 0, 0, -2, -2,
1399: 0, -2, -2, 0,
1400: 0, -2, -2, 0,
1401: 0, -2, -2, 0,
1402: 0, 0, 0, -1, 0, 1,
1403: 0, 0, 0, 0, 0, -4,
1404: 0, 0, 0, 0, -1, 1,
1405: -4, 0, 0, -1, 0, 1,
1406: -4, 0, 0, 0, 0, -4,
1407: -4, 0, 0, 0, -1, 1};
1408: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms.
1409: 2
1410: |\
1411: | \
1412: | \
1413: 0---1
1415: 2
1417: 0 1
1419: 2
1420: |\
1421: | \
1422: | \
1423: 0---1
1424: */
1425: static DMPolytopeType ttripT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1426: static PetscInt ttripS[] = {1, 3, 3};
1427: static PetscInt ttripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1428: 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,
1429: 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,
1430: 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,
1431: 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,
1432: 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,
1433: 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};
1434: static PetscInt ttripO[] = {0, 0,
1435: 0, 0, 0, 0,
1436: 0, 0, 0, 0,
1437: 0, 0, 0, 0,
1438: 0, 0, 0, 0, -1, 0,
1439: 0, 0, 0, 0, 0, -1,
1440: 0, 0, 0, -1, 0, 0};
1441: /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms.
1442: 2
1443: |\
1444: | \
1445: | \
1446: 0---1
1448: 2
1450: 0 1
1452: 2
1453: |\
1454: | \
1455: | \
1456: 0---1
1457: */
1458: static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1459: static PetscInt ctripS[] = {1, 3, 3};
1460: static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1461: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1462: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1463: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1464: 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,
1465: 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,
1466: 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};
1467: static PetscInt ctripO[] = {0, 0,
1468: 0, 0, -2, -2,
1469: 0, 0, -2, -2,
1470: 0, 0, -2, -2,
1471: -4, 0, 0, -1, 0, 1,
1472: -4, 0, 0, 0, 0, -4,
1473: -4, 0, 0, 0, -1, 1};
1474: /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */
1475: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1476: static PetscInt tquadpS[] = {1, 4, 4};
1477: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1478: 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,
1479: 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,
1480: 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,
1481: 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,
1482: 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,
1483: 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,
1484: 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,
1485: 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};
1486: static PetscInt tquadpO[] = {0, 0,
1487: 0, 0, 0, 0,
1488: 0, 0, 0, 0,
1489: 0, 0, 0, 0,
1490: 0, 0, 0, 0,
1491: 0, 0, 0, 0, -1, 0,
1492: 0, 0, 0, 0, 0, -1,
1493: 0, 0, -1, 0, 0, 0,
1494: 0, 0, 0, -1, 0, 0};
1495: PetscBool convertTensor = PETSC_TRUE;
1498: if (rt) *rt = 0;
1499: if (convertTensor) {
1500: switch (source) {
1501: case DM_POLYTOPE_POINT:
1502: case DM_POLYTOPE_SEGMENT:
1503: case DM_POLYTOPE_QUADRILATERAL:
1504: case DM_POLYTOPE_HEXAHEDRON:
1505: DMPlexCellRefinerRefine_Regular(cr, source, p, rt, Nt, target, size, cone, ornt);
1506: break;
1507: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1508: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1509: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ctripT; *size = ctripS; *cone = ctripC; *ornt = ctripO; break;
1510: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1511: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1512: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1513: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1514: /* TODO Fix pyramids: For now, we just ignore them */
1515: case DM_POLYTOPE_PYRAMID:
1516: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1517: break;
1518: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1519: }
1520: } else {
1521: switch (source) {
1522: case DM_POLYTOPE_POINT:
1523: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1524: case DM_POLYTOPE_SEGMENT:
1525: case DM_POLYTOPE_QUADRILATERAL:
1526: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1527: case DM_POLYTOPE_HEXAHEDRON:
1528: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1529: DMPlexCellRefinerRefine_Regular(cr, source, p, rt, Nt, target, size, cone, ornt);
1530: break;
1531: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1532: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1533: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1534: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1535: /* TODO Fix pyramids: For now, we just ignore them */
1536: case DM_POLYTOPE_PYRAMID:
1537: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1538: break;
1539: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1540: }
1541: }
1542: return(0);
1543: }
1545: static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1546: {
1550: if (rt) *rt = 0;
1551: switch (source) {
1552: case DM_POLYTOPE_POINT:
1553: case DM_POLYTOPE_SEGMENT:
1554: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1555: case DM_POLYTOPE_TRIANGLE:
1556: case DM_POLYTOPE_TETRAHEDRON:
1557: case DM_POLYTOPE_TRI_PRISM:
1558: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1559: case DM_POLYTOPE_QUADRILATERAL:
1560: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1561: case DM_POLYTOPE_HEXAHEDRON:
1562: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1563: DMPlexCellRefinerRefine_Regular(cr, source, p, rt, Nt, target, size, cone, ornt);
1564: break;
1565: /* TODO Fix pyramids: For now, we just ignore them */
1566: case DM_POLYTOPE_PYRAMID:
1567: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1568: break;
1569: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1570: }
1571: return(0);
1572: }
1574: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1575: {
1577: /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
1578: 2
1579: |\
1580: |\\
1581: | |\
1582: | \ \
1583: | | \
1584: | \ \
1585: | | \
1586: 2 \ \
1587: | | 1
1588: | 2 \
1589: | | \
1590: | /\ \
1591: | 0 1 |
1592: | / \ |
1593: |/ \|
1594: 0---0----1
1595: */
1596: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1597: static PetscInt triS[] = {1, 3, 3};
1598: static PetscInt triC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1599: DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1600: DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1601: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1602: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1,
1603: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
1604: static PetscInt triO[] = {0, 0,
1605: 0, 0,
1606: 0, 0,
1607: 0, 0, -2,
1608: 0, 0, -2,
1609: 0, 0, -2};
1612: if (rt) *rt = 0;
1613: switch (source) {
1614: case DM_POLYTOPE_POINT:
1615: case DM_POLYTOPE_SEGMENT:
1616: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1617: case DM_POLYTOPE_QUADRILATERAL:
1618: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1619: case DM_POLYTOPE_TETRAHEDRON:
1620: case DM_POLYTOPE_HEXAHEDRON:
1621: case DM_POLYTOPE_TRI_PRISM:
1622: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1623: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1624: case DM_POLYTOPE_PYRAMID:
1625: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1626: break;
1627: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1628: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1629: }
1630: return(0);
1631: }
1633: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1634: {
1636: /* Add 6 triangles inside every cell, making 4 new tets
1637: 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
1638: 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]
1639: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1640: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1641: We make a new tet on each face
1642: [v0, v1, v2, (c0, 0)]
1643: [v0, v3, v1, (c0, 0)]
1644: [v0, v2, v3, (c0, 0)]
1645: [v2, v1, v3, (c0, 0)]
1646: We create a new face for each edge
1647: [v0, (c0, 0), v1 ]
1648: [v0, v2, (c0, 0)]
1649: [v2, v1, (c0, 0)]
1650: [v0, (c0, 0), v3 ]
1651: [v1, v3, (c0, 0)]
1652: [v3, v2, (c0, 0)]
1653: */
1654: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
1655: static PetscInt tetS[] = {1, 4, 6, 4};
1656: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1657: DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1658: DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1659: DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1660: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0,
1661: DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1662: DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1663: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0,
1664: DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1,
1665: DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1666: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2,
1667: DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4,
1668: DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5,
1669: DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
1670: static PetscInt tetO[] = {0, 0,
1671: 0, 0,
1672: 0, 0,
1673: 0, 0,
1674: 0, -2, -2,
1675: -2, 0, -2,
1676: -2, 0, -2,
1677: 0, -2, -2,
1678: -2, 0, -2,
1679: -2, 0, -2,
1680: 0, 0, 0, 0,
1681: 0, 0, -3, 0,
1682: 0, -3, -3, 0,
1683: 0, -3, -1, -1};
1686: if (rt) *rt = 0;
1687: switch (source) {
1688: case DM_POLYTOPE_POINT:
1689: case DM_POLYTOPE_SEGMENT:
1690: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1691: case DM_POLYTOPE_TRIANGLE:
1692: case DM_POLYTOPE_QUADRILATERAL:
1693: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1694: case DM_POLYTOPE_HEXAHEDRON:
1695: case DM_POLYTOPE_TRI_PRISM:
1696: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1697: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1698: case DM_POLYTOPE_PYRAMID:
1699: DMPlexCellRefinerRefine_None(cr, source, p, rt, Nt, target, size, cone, ornt);
1700: break;
1701: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1702: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1703: }
1704: return(0);
1705: }
1707: typedef struct {
1708: PetscInt n;
1709: PetscReal r;
1710: PetscScalar *h;
1711: PetscInt *Nt;
1712: DMPolytopeType **target;
1713: PetscInt **size;
1714: PetscInt **cone;
1715: PetscInt **ornt;
1716: } PlexRefiner_BL;
1718: static PetscErrorCode DMPlexCellRefinerSetUp_BL(DMPlexCellRefiner cr)
1719: {
1720: PlexRefiner_BL *crbl;
1722: PetscInt i,n;
1723: PetscReal r;
1724: PetscInt c1,c2,o1,o2;
1727: PetscNew(&crbl);
1728: cr->data = crbl;
1729: crbl->n = 1; /* 1 split -> 2 new cells */
1730: crbl->r = 1; /* linear progression */
1732: /* TODO: add setfromoptions to the refiners? */
1733: PetscOptionsGetInt(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_splits", &crbl->n, NULL);
1734: if (crbl->n < 1) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Number of splits %D must be positive",crbl->n);
1735: PetscOptionsGetReal(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_progression", &crbl->r, NULL);
1736: n = crbl->n;
1737: r = crbl->r;
1739: /* we only split DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR and DM_POLYTOPE_QUAD_PRISM_TENSOR */
1740: PetscMalloc5(4,&crbl->Nt,4,&crbl->target,4,&crbl->size,4,&crbl->cone,4,&crbl->ornt);
1742: /* progression */
1743: PetscMalloc1(n,&crbl->h);
1744: if (r > 1) {
1745: PetscReal d = (r-1.)/(PetscPowRealInt(r,n+1)-1.);
1747: crbl->h[0] = d;
1748: for (i = 1; i < n; i++) {
1749: d *= r;
1750: crbl->h[i] = crbl->h[i-1] + d;
1751: }
1752: } else { /* linear */
1753: for (i = 0; i < n; i++) crbl->h[i] = (i + 1.)/(n+1); /* linear */
1754: }
1756: /* DM_POLYTOPE_POINT_PRISM_TENSOR produces n points and n+1 tensor segments */
1757: c1 = 14+6*(n-1);
1758: o1 = 2*(n+1);
1759: crbl->Nt[0] = 2;
1761: PetscMalloc4(crbl->Nt[0],&crbl->target[0],crbl->Nt[0],&crbl->size[0],c1,&crbl->cone[0],o1,&crbl->ornt[0]);
1763: crbl->target[0][0] = DM_POLYTOPE_POINT;
1764: crbl->target[0][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1766: crbl->size[0][0] = n;
1767: crbl->size[0][1] = n+1;
1769: /* the tensor segments */
1770: crbl->cone[0][0] = DM_POLYTOPE_POINT;
1771: crbl->cone[0][1] = 1;
1772: crbl->cone[0][2] = 0;
1773: crbl->cone[0][3] = 0;
1774: crbl->cone[0][4] = DM_POLYTOPE_POINT;
1775: crbl->cone[0][5] = 0;
1776: crbl->cone[0][6] = 0;
1777: for (i = 0; i < n-1; i++) {
1778: crbl->cone[0][7+6*i+0] = DM_POLYTOPE_POINT;
1779: crbl->cone[0][7+6*i+1] = 0;
1780: crbl->cone[0][7+6*i+2] = i;
1781: crbl->cone[0][7+6*i+3] = DM_POLYTOPE_POINT;
1782: crbl->cone[0][7+6*i+4] = 0;
1783: crbl->cone[0][7+6*i+5] = i+1;
1784: }
1785: crbl->cone[0][7+6*(n-1)+0] = DM_POLYTOPE_POINT;
1786: crbl->cone[0][7+6*(n-1)+1] = 0;
1787: crbl->cone[0][7+6*(n-1)+2] = n-1;
1788: crbl->cone[0][7+6*(n-1)+3] = DM_POLYTOPE_POINT;
1789: crbl->cone[0][7+6*(n-1)+4] = 1;
1790: crbl->cone[0][7+6*(n-1)+5] = 1;
1791: crbl->cone[0][7+6*(n-1)+6] = 0;
1792: for (i = 0; i < o1; i++) crbl->ornt[0][i] = 0;
1794: /* DM_POLYTOPE_SEG_PRISM_TENSOR produces n segments and n+1 tensor quads */
1795: c1 = 8*n;
1796: c2 = 30+14*(n-1);
1797: o1 = 2*n;
1798: o2 = 4*(n+1);
1799: crbl->Nt[1] = 2;
1801: PetscMalloc4(crbl->Nt[1],&crbl->target[1],crbl->Nt[1],&crbl->size[1],c1+c2,&crbl->cone[1],o1+o2,&crbl->ornt[1]);
1803: crbl->target[1][0] = DM_POLYTOPE_SEGMENT;
1804: crbl->target[1][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1806: crbl->size[1][0] = n;
1807: crbl->size[1][1] = n+1;
1809: /* the segments */
1810: for (i = 0; i < n; i++) {
1811: crbl->cone[1][8*i+0] = DM_POLYTOPE_POINT;
1812: crbl->cone[1][8*i+1] = 1;
1813: crbl->cone[1][8*i+2] = 2;
1814: crbl->cone[1][8*i+3] = i;
1815: crbl->cone[1][8*i+4] = DM_POLYTOPE_POINT;
1816: crbl->cone[1][8*i+5] = 1;
1817: crbl->cone[1][8*i+6] = 3;
1818: crbl->cone[1][8*i+7] = i;
1819: }
1821: /* the tensor quads */
1822: crbl->cone[1][c1+ 0] = DM_POLYTOPE_SEGMENT;
1823: crbl->cone[1][c1+ 1] = 1;
1824: crbl->cone[1][c1+ 2] = 0;
1825: crbl->cone[1][c1+ 3] = 0;
1826: crbl->cone[1][c1+ 4] = DM_POLYTOPE_SEGMENT;
1827: crbl->cone[1][c1+ 5] = 0;
1828: crbl->cone[1][c1+ 6] = 0;
1829: crbl->cone[1][c1+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1830: crbl->cone[1][c1+ 8] = 1;
1831: crbl->cone[1][c1+ 9] = 2;
1832: crbl->cone[1][c1+10] = 0;
1833: crbl->cone[1][c1+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1834: crbl->cone[1][c1+12] = 1;
1835: crbl->cone[1][c1+13] = 3;
1836: crbl->cone[1][c1+14] = 0;
1837: for (i = 0; i < n-1; i++) {
1838: crbl->cone[1][c1+15+14*i+ 0] = DM_POLYTOPE_SEGMENT;
1839: crbl->cone[1][c1+15+14*i+ 1] = 0;
1840: crbl->cone[1][c1+15+14*i+ 2] = i;
1841: crbl->cone[1][c1+15+14*i+ 3] = DM_POLYTOPE_SEGMENT;
1842: crbl->cone[1][c1+15+14*i+ 4] = 0;
1843: crbl->cone[1][c1+15+14*i+ 5] = i+1;
1844: crbl->cone[1][c1+15+14*i+ 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1845: crbl->cone[1][c1+15+14*i+ 7] = 1;
1846: crbl->cone[1][c1+15+14*i+ 8] = 2;
1847: crbl->cone[1][c1+15+14*i+ 9] = i+1;
1848: crbl->cone[1][c1+15+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1849: crbl->cone[1][c1+15+14*i+11] = 1;
1850: crbl->cone[1][c1+15+14*i+12] = 3;
1851: crbl->cone[1][c1+15+14*i+13] = i+1;
1852: }
1853: crbl->cone[1][c1+15+14*(n-1)+ 0] = DM_POLYTOPE_SEGMENT;
1854: crbl->cone[1][c1+15+14*(n-1)+ 1] = 0;
1855: crbl->cone[1][c1+15+14*(n-1)+ 2] = n-1;
1856: crbl->cone[1][c1+15+14*(n-1)+ 3] = DM_POLYTOPE_SEGMENT;
1857: crbl->cone[1][c1+15+14*(n-1)+ 4] = 1;
1858: crbl->cone[1][c1+15+14*(n-1)+ 5] = 1;
1859: crbl->cone[1][c1+15+14*(n-1)+ 6] = 0;
1860: crbl->cone[1][c1+15+14*(n-1)+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1861: crbl->cone[1][c1+15+14*(n-1)+ 8] = 1;
1862: crbl->cone[1][c1+15+14*(n-1)+ 9] = 2;
1863: crbl->cone[1][c1+15+14*(n-1)+10] = n;
1864: crbl->cone[1][c1+15+14*(n-1)+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1865: crbl->cone[1][c1+15+14*(n-1)+12] = 1;
1866: crbl->cone[1][c1+15+14*(n-1)+13] = 3;
1867: crbl->cone[1][c1+15+14*(n-1)+14] = n;
1868: for (i = 0; i < o1+o2; i++) crbl->ornt[1][i] = 0;
1870: /* DM_POLYTOPE_TRI_PRISM_TENSOR produces n triangles and n+1 tensor triangular prisms */
1871: c1 = 12*n;
1872: c2 = 38+18*(n-1);
1873: o1 = 3*n;
1874: o2 = 5*(n+1);
1875: crbl->Nt[2] = 2;
1877: PetscMalloc4(crbl->Nt[2],&crbl->target[2],crbl->Nt[2],&crbl->size[2],c1+c2,&crbl->cone[2],o1+o2,&crbl->ornt[2]);
1879: crbl->target[2][0] = DM_POLYTOPE_TRIANGLE;
1880: crbl->target[2][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;
1882: crbl->size[2][0] = n;
1883: crbl->size[2][1] = n+1;
1885: /* the triangles */
1886: for (i = 0; i < n; i++) {
1887: crbl->cone[2][12*i+ 0] = DM_POLYTOPE_SEGMENT;
1888: crbl->cone[2][12*i+ 1] = 1;
1889: crbl->cone[2][12*i+ 2] = 2;
1890: crbl->cone[2][12*i+ 3] = i;
1891: crbl->cone[2][12*i+ 4] = DM_POLYTOPE_SEGMENT;
1892: crbl->cone[2][12*i+ 5] = 1;
1893: crbl->cone[2][12*i+ 6] = 3;
1894: crbl->cone[2][12*i+ 7] = i;
1895: crbl->cone[2][12*i+ 8] = DM_POLYTOPE_SEGMENT;
1896: crbl->cone[2][12*i+ 9] = 1;
1897: crbl->cone[2][12*i+10] = 4;
1898: crbl->cone[2][12*i+11] = i;
1899: }
1901: /* the triangular prisms */
1902: crbl->cone[2][c1+ 0] = DM_POLYTOPE_TRIANGLE;
1903: crbl->cone[2][c1+ 1] = 1;
1904: crbl->cone[2][c1+ 2] = 0;
1905: crbl->cone[2][c1+ 3] = 0;
1906: crbl->cone[2][c1+ 4] = DM_POLYTOPE_TRIANGLE;
1907: crbl->cone[2][c1+ 5] = 0;
1908: crbl->cone[2][c1+ 6] = 0;
1909: crbl->cone[2][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1910: crbl->cone[2][c1+ 8] = 1;
1911: crbl->cone[2][c1+ 9] = 2;
1912: crbl->cone[2][c1+10] = 0;
1913: crbl->cone[2][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1914: crbl->cone[2][c1+12] = 1;
1915: crbl->cone[2][c1+13] = 3;
1916: crbl->cone[2][c1+14] = 0;
1917: crbl->cone[2][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1918: crbl->cone[2][c1+16] = 1;
1919: crbl->cone[2][c1+17] = 4;
1920: crbl->cone[2][c1+18] = 0;
1921: for (i = 0; i < n-1; i++) {
1922: crbl->cone[2][c1+19+18*i+ 0] = DM_POLYTOPE_TRIANGLE;
1923: crbl->cone[2][c1+19+18*i+ 1] = 0;
1924: crbl->cone[2][c1+19+18*i+ 2] = i;
1925: crbl->cone[2][c1+19+18*i+ 3] = DM_POLYTOPE_TRIANGLE;
1926: crbl->cone[2][c1+19+18*i+ 4] = 0;
1927: crbl->cone[2][c1+19+18*i+ 5] = i+1;
1928: crbl->cone[2][c1+19+18*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1929: crbl->cone[2][c1+19+18*i+ 7] = 1;
1930: crbl->cone[2][c1+19+18*i+ 8] = 2;
1931: crbl->cone[2][c1+19+18*i+ 9] = i+1;
1932: crbl->cone[2][c1+19+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1933: crbl->cone[2][c1+19+18*i+11] = 1;
1934: crbl->cone[2][c1+19+18*i+12] = 3;
1935: crbl->cone[2][c1+19+18*i+13] = i+1;
1936: crbl->cone[2][c1+19+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1937: crbl->cone[2][c1+19+18*i+15] = 1;
1938: crbl->cone[2][c1+19+18*i+16] = 4;
1939: crbl->cone[2][c1+19+18*i+17] = i+1;
1940: }
1941: crbl->cone[2][c1+19+18*(n-1)+ 0] = DM_POLYTOPE_TRIANGLE;
1942: crbl->cone[2][c1+19+18*(n-1)+ 1] = 0;
1943: crbl->cone[2][c1+19+18*(n-1)+ 2] = n-1;
1944: crbl->cone[2][c1+19+18*(n-1)+ 3] = DM_POLYTOPE_TRIANGLE;
1945: crbl->cone[2][c1+19+18*(n-1)+ 4] = 1;
1946: crbl->cone[2][c1+19+18*(n-1)+ 5] = 1;
1947: crbl->cone[2][c1+19+18*(n-1)+ 6] = 0;
1948: crbl->cone[2][c1+19+18*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1949: crbl->cone[2][c1+19+18*(n-1)+ 8] = 1;
1950: crbl->cone[2][c1+19+18*(n-1)+ 9] = 2;
1951: crbl->cone[2][c1+19+18*(n-1)+10] = n;
1952: crbl->cone[2][c1+19+18*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1953: crbl->cone[2][c1+19+18*(n-1)+12] = 1;
1954: crbl->cone[2][c1+19+18*(n-1)+13] = 3;
1955: crbl->cone[2][c1+19+18*(n-1)+14] = n;
1956: crbl->cone[2][c1+19+18*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1957: crbl->cone[2][c1+19+18*(n-1)+16] = 1;
1958: crbl->cone[2][c1+19+18*(n-1)+17] = 4;
1959: crbl->cone[2][c1+19+18*(n-1)+18] = n;
1960: for (i = 0; i < o1+o2; i++) crbl->ornt[2][i] = 0;
1962: /* DM_POLYTOPE_QUAD_PRISM_TENSOR produces n quads and n+1 tensor quad prisms */
1963: c1 = 16*n;
1964: c2 = 46+22*(n-1);
1965: o1 = 4*n;
1966: o2 = 6*(n+1);
1967: crbl->Nt[3] = 2;
1969: PetscMalloc4(crbl->Nt[3],&crbl->target[3],crbl->Nt[3],&crbl->size[3],c1+c2,&crbl->cone[3],o1+o2,&crbl->ornt[3]);
1971: crbl->target[3][0] = DM_POLYTOPE_QUADRILATERAL;
1972: crbl->target[3][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;
1974: crbl->size[3][0] = n;
1975: crbl->size[3][1] = n+1;
1977: /* the quads */
1978: for (i = 0; i < n; i++) {
1979: crbl->cone[3][16*i+ 0] = DM_POLYTOPE_SEGMENT;
1980: crbl->cone[3][16*i+ 1] = 1;
1981: crbl->cone[3][16*i+ 2] = 2;
1982: crbl->cone[3][16*i+ 3] = i;
1983: crbl->cone[3][16*i+ 4] = DM_POLYTOPE_SEGMENT;
1984: crbl->cone[3][16*i+ 5] = 1;
1985: crbl->cone[3][16*i+ 6] = 3;
1986: crbl->cone[3][16*i+ 7] = i;
1987: crbl->cone[3][16*i+ 8] = DM_POLYTOPE_SEGMENT;
1988: crbl->cone[3][16*i+ 9] = 1;
1989: crbl->cone[3][16*i+10] = 4;
1990: crbl->cone[3][16*i+11] = i;
1991: crbl->cone[3][16*i+12] = DM_POLYTOPE_SEGMENT;
1992: crbl->cone[3][16*i+13] = 1;
1993: crbl->cone[3][16*i+14] = 5;
1994: crbl->cone[3][16*i+15] = i;
1995: }
1997: /* the quad prisms */
1998: crbl->cone[3][c1+ 0] = DM_POLYTOPE_QUADRILATERAL;
1999: crbl->cone[3][c1+ 1] = 1;
2000: crbl->cone[3][c1+ 2] = 0;
2001: crbl->cone[3][c1+ 3] = 0;
2002: crbl->cone[3][c1+ 4] = DM_POLYTOPE_QUADRILATERAL;
2003: crbl->cone[3][c1+ 5] = 0;
2004: crbl->cone[3][c1+ 6] = 0;
2005: crbl->cone[3][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2006: crbl->cone[3][c1+ 8] = 1;
2007: crbl->cone[3][c1+ 9] = 2;
2008: crbl->cone[3][c1+10] = 0;
2009: crbl->cone[3][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2010: crbl->cone[3][c1+12] = 1;
2011: crbl->cone[3][c1+13] = 3;
2012: crbl->cone[3][c1+14] = 0;
2013: crbl->cone[3][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2014: crbl->cone[3][c1+16] = 1;
2015: crbl->cone[3][c1+17] = 4;
2016: crbl->cone[3][c1+18] = 0;
2017: crbl->cone[3][c1+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2018: crbl->cone[3][c1+20] = 1;
2019: crbl->cone[3][c1+21] = 5;
2020: crbl->cone[3][c1+22] = 0;
2021: for (i = 0; i < n-1; i++) {
2022: crbl->cone[3][c1+23+22*i+ 0] = DM_POLYTOPE_QUADRILATERAL;
2023: crbl->cone[3][c1+23+22*i+ 1] = 0;
2024: crbl->cone[3][c1+23+22*i+ 2] = i;
2025: crbl->cone[3][c1+23+22*i+ 3] = DM_POLYTOPE_QUADRILATERAL;
2026: crbl->cone[3][c1+23+22*i+ 4] = 0;
2027: crbl->cone[3][c1+23+22*i+ 5] = i+1;
2028: crbl->cone[3][c1+23+22*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2029: crbl->cone[3][c1+23+22*i+ 7] = 1;
2030: crbl->cone[3][c1+23+22*i+ 8] = 2;
2031: crbl->cone[3][c1+23+22*i+ 9] = i+1;
2032: crbl->cone[3][c1+23+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2033: crbl->cone[3][c1+23+22*i+11] = 1;
2034: crbl->cone[3][c1+23+22*i+12] = 3;
2035: crbl->cone[3][c1+23+22*i+13] = i+1;
2036: crbl->cone[3][c1+23+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2037: crbl->cone[3][c1+23+22*i+15] = 1;
2038: crbl->cone[3][c1+23+22*i+16] = 4;
2039: crbl->cone[3][c1+23+22*i+17] = i+1;
2040: crbl->cone[3][c1+23+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2041: crbl->cone[3][c1+23+22*i+19] = 1;
2042: crbl->cone[3][c1+23+22*i+20] = 5;
2043: crbl->cone[3][c1+23+22*i+21] = i+1;
2044: }
2045: crbl->cone[3][c1+23+22*(n-1)+ 0] = DM_POLYTOPE_QUADRILATERAL;
2046: crbl->cone[3][c1+23+22*(n-1)+ 1] = 0;
2047: crbl->cone[3][c1+23+22*(n-1)+ 2] = n-1;
2048: crbl->cone[3][c1+23+22*(n-1)+ 3] = DM_POLYTOPE_QUADRILATERAL;
2049: crbl->cone[3][c1+23+22*(n-1)+ 4] = 1;
2050: crbl->cone[3][c1+23+22*(n-1)+ 5] = 1;
2051: crbl->cone[3][c1+23+22*(n-1)+ 6] = 0;
2052: crbl->cone[3][c1+23+22*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2053: crbl->cone[3][c1+23+22*(n-1)+ 8] = 1;
2054: crbl->cone[3][c1+23+22*(n-1)+ 9] = 2;
2055: crbl->cone[3][c1+23+22*(n-1)+10] = n;
2056: crbl->cone[3][c1+23+22*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2057: crbl->cone[3][c1+23+22*(n-1)+12] = 1;
2058: crbl->cone[3][c1+23+22*(n-1)+13] = 3;
2059: crbl->cone[3][c1+23+22*(n-1)+14] = n;
2060: crbl->cone[3][c1+23+22*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2061: crbl->cone[3][c1+23+22*(n-1)+16] = 1;
2062: crbl->cone[3][c1+23+22*(n-1)+17] = 4;
2063: crbl->cone[3][c1+23+22*(n-1)+18] = n;
2064: crbl->cone[3][c1+23+22*(n-1)+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2065: crbl->cone[3][c1+23+22*(n-1)+20] = 1;
2066: crbl->cone[3][c1+23+22*(n-1)+21] = 5;
2067: crbl->cone[3][c1+23+22*(n-1)+22] = n;
2068: for (i = 0; i < o1+o2; i++) crbl->ornt[3][i] = 0;
2069: return(0);
2070: }
2072: static PetscErrorCode DMPlexCellRefinerDestroy_BL(DMPlexCellRefiner cr)
2073: {
2074: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2078: PetscFree4(crbl->target[0],crbl->size[0],crbl->cone[0],crbl->ornt[0]);
2079: PetscFree4(crbl->target[1],crbl->size[1],crbl->cone[1],crbl->ornt[1]);
2080: PetscFree4(crbl->target[2],crbl->size[2],crbl->cone[2],crbl->ornt[2]);
2081: PetscFree4(crbl->target[3],crbl->size[3],crbl->cone[3],crbl->ornt[3]);
2082: PetscFree5(crbl->Nt,crbl->target,crbl->size,crbl->cone,crbl->ornt);
2083: PetscFree(crbl->h);
2084: PetscFree(cr->data);
2085: return(0);
2086: }
2088: static PetscErrorCode DMPlexCellRefinerRefine_BL(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
2089: {
2090: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2091: PetscErrorCode ierr;
2094: if (rt) *rt = 0;
2095: switch (source) {
2096: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2097: *Nt = crbl->Nt[0];
2098: *target = crbl->target[0];
2099: *size = crbl->size[0];
2100: *cone = crbl->cone[0];
2101: *ornt = crbl->ornt[0];
2102: break;
2103: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2104: *Nt = crbl->Nt[1];
2105: *target = crbl->target[1];
2106: *size = crbl->size[1];
2107: *cone = crbl->cone[1];
2108: *ornt = crbl->ornt[1];
2109: break;
2110: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2111: *Nt = crbl->Nt[2];
2112: *target = crbl->target[2];
2113: *size = crbl->size[2];
2114: *cone = crbl->cone[2];
2115: *ornt = crbl->ornt[2];
2116: break;
2117: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2118: *Nt = crbl->Nt[3];
2119: *target = crbl->target[3];
2120: *size = crbl->size[3];
2121: *cone = crbl->cone[3];
2122: *ornt = crbl->ornt[3];
2123: break;
2124: default:
2125: DMPlexCellRefinerRefine_None(cr,source,p,rt,Nt,target,size,cone,ornt);
2126: }
2127: return(0);
2128: }
2130: static PetscErrorCode DMPlexCellRefinerMapSubcells_BL(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
2131: {
2132: /* We shift any input orientation in order to make it non-negative
2133: 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)
2134: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
2135: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
2136: */
2137: PetscInt tquad_seg_o[] = { 0, 1, -2, -1,
2138: 0, 1, -2, -1,
2139: -2, -1, 0, 1,
2140: -2, -1, 0, 1};
2141: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
2142: 1, 0, -1, -2,
2143: -2, -1, 0, 1,
2144: -1, -2, 1, 0};
2145: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2146: const PetscInt n = crbl->n;
2150: *rnew = r;
2151: *onew = o;
2152: switch (pct) {
2153: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2154: if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
2155: if (po == 0 || po == -1) {*rnew = r; *onew = o;}
2156: else if (po == 1 || po == -2) {*rnew = n - r; *onew = (o == 0 || o == -1) ? -2 : 0;}
2157: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for tensor segment", po);
2158: }
2159: break;
2160: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2161: switch (ct) {
2162: case DM_POLYTOPE_SEGMENT:
2163: *onew = tquad_seg_o[(po+2)*4+o+2];
2164: *rnew = r;
2165: break;
2166: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2167: *onew = tquad_tquad_o[(po+2)*4+o+2];
2168: *rnew = r;
2169: break;
2170: default: break;
2171: }
2172: break;
2173: default:
2174: DMPlexCellRefinerMapSubcells_None(cr, pct, pp, po, ct, r, o, rnew, onew);
2175: }
2176: return(0);
2177: }
2179: static PetscErrorCode DMPlexCellRefinerMapCoordinates_BL(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2180: {
2181: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2182: PetscInt d;
2183: PetscErrorCode ierr;
2186: switch (pct) {
2187: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2188: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
2189: if (Nv != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent vertices %D",Nv);
2190: if (r >= crbl->n || r < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Invalid replica %D, must be in [0,%D)",r,crbl->n);
2191: for (d = 0; d < dE; d++) out[d] = in[d] + crbl->h[r] * (in[d + dE] - in[d]);
2192: break;
2193: default:
2194: DMPlexCellRefinerMapCoordinates_Barycenter(cr,pct,ct,r,Nv,dE,in,out);
2195: }
2196: return(0);
2197: }
2199: typedef struct {
2200: DMLabel splitPoints; /* List of edges to be bisected (1) and cells to be divided (2) */
2201: PetscSection secEdgeLen; /* Section for edge length field */
2202: PetscReal *edgeLen; /* Storage for edge length field */
2203: PetscInt *splitArray; /* Array for communication of split points label */
2204: } PlexRefiner_SBR;
2206: typedef struct _p_PointQueue *PointQueue;
2207: struct _p_PointQueue {
2208: PetscInt size; /* Size of the storage array */
2209: PetscInt *points; /* Array of mesh points */
2210: PetscInt front; /* Index of the front of the queue */
2211: PetscInt back; /* Index of the back of the queue */
2212: PetscInt num; /* Number of enqueued points */
2213: };
2215: static PetscErrorCode PointQueueCreate(PetscInt size, PointQueue *queue)
2216: {
2217: PointQueue q;
2221: if (size < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Queue size %D must be non-negative", size);
2222: PetscCalloc1(1, &q);
2223: q->size = size;
2224: PetscMalloc1(q->size, &q->points);
2225: q->num = 0;
2226: q->front = 0;
2227: q->back = q->size-1;
2228: *queue = q;
2229: return(0);
2230: }
2232: static PetscErrorCode PointQueueDestroy(PointQueue *queue)
2233: {
2234: PointQueue q = *queue;
2238: PetscFree(q->points);
2239: PetscFree(q);
2240: *queue = NULL;
2241: return(0);
2242: }
2244: static PetscErrorCode PointQueueEnsureSize(PointQueue queue)
2245: {
2249: if (queue->num < queue->size) return(0);
2250: queue->size *= 2;
2251: PetscRealloc(queue->size * sizeof(PetscInt), &queue->points);
2252: return(0);
2253: }
2255: static PetscErrorCode PointQueueEnqueue(PointQueue queue, PetscInt p)
2256: {
2260: PointQueueEnsureSize(queue);
2261: queue->back = (queue->back + 1) % queue->size;
2262: queue->points[queue->back] = p;
2263: ++queue->num;
2264: return(0);
2265: }
2267: static PetscErrorCode PointQueueDequeue(PointQueue queue, PetscInt *p)
2268: {
2270: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot dequeue from an empty queue");
2271: *p = queue->points[queue->front];
2272: queue->front = (queue->front + 1) % queue->size;
2273: --queue->num;
2274: return(0);
2275: }
2277: #if 0
2278: static PetscErrorCode PointQueueFront(PointQueue queue, PetscInt *p)
2279: {
2281: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the front of an empty queue");
2282: *p = queue->points[queue->front];
2283: return(0);
2284: }
2286: static PetscErrorCode PointQueueBack(PointQueue queue, PetscInt *p)
2287: {
2289: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the back of an empty queue");
2290: *p = queue->points[queue->back];
2291: return(0);
2292: }
2293: #endif
2295: PETSC_STATIC_INLINE PetscBool PointQueueEmpty(PointQueue queue)
2296: {
2297: if (!queue->num) return PETSC_TRUE;
2298: return PETSC_FALSE;
2299: }
2301: static PetscErrorCode SBRGetEdgeLen_Private(DMPlexCellRefiner cr, PetscInt edge, PetscReal *len)
2302: {
2303: PlexRefiner_SBR *sbr = (PlexRefiner_SBR *) cr->data;
2304: DM dm = cr->dm;
2305: PetscInt off;
2306: PetscErrorCode ierr;
2309: PetscSectionGetOffset(sbr->secEdgeLen, edge, &off);
2310: if (sbr->edgeLen[off] <= 0.0) {
2311: DM cdm;
2312: Vec coordsLocal;
2313: const PetscScalar *coords;
2314: const PetscInt *cone;
2315: PetscScalar *cA, *cB;
2316: PetscInt coneSize, cdim;
2318: DMGetCoordinateDM(dm, &cdm);
2319: DMPlexGetCone(dm, edge, &cone);
2320: DMPlexGetConeSize(dm, edge, &coneSize);
2321: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %D cone size must be 2, not %D", edge, coneSize);
2322: DMGetCoordinateDim(dm, &cdim);
2323: DMGetCoordinatesLocal(dm, &coordsLocal);
2324: VecGetArrayRead(coordsLocal, &coords);
2325: DMPlexPointLocalRead(cdm, cone[0], coords, &cA);
2326: DMPlexPointLocalRead(cdm, cone[1], coords, &cB);
2327: sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
2328: VecRestoreArrayRead(coordsLocal, &coords);
2329: }
2330: *len = sbr->edgeLen[off];
2331: return(0);
2332: }
2334: /* Mark local edges that should be split */
2335: /* TODO This will not work in 3D */
2336: static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexCellRefiner cr, PointQueue queue)
2337: {
2338: PlexRefiner_SBR *sbr = (PlexRefiner_SBR *) cr->data;
2339: DM dm = cr->dm;
2340: PetscErrorCode ierr;
2343: while (!PointQueueEmpty(queue)) {
2344: PetscInt p = -1;
2345: const PetscInt *support;
2346: PetscInt supportSize, s;
2348: PointQueueDequeue(queue, &p);
2349: DMPlexGetSupport(dm, p, &support);
2350: DMPlexGetSupportSize(dm, p, &supportSize);
2351: for (s = 0; s < supportSize; ++s) {
2352: const PetscInt cell = support[s];
2353: const PetscInt *cone;
2354: PetscInt coneSize, c;
2355: PetscInt cval, eval, maxedge;
2356: PetscReal len, maxlen;
2358: DMLabelGetValue(sbr->splitPoints, cell, &cval);
2359: if (cval == 2) continue;
2360: DMPlexGetCone(dm, cell, &cone);
2361: DMPlexGetConeSize(dm, cell, &coneSize);
2362: SBRGetEdgeLen_Private(cr, cone[0], &maxlen);
2363: maxedge = cone[0];
2364: for (c = 1; c < coneSize; ++c) {
2365: SBRGetEdgeLen_Private(cr, cone[c], &len);
2366: if (len > maxlen) {maxlen = len; maxedge = cone[c];}
2367: }
2368: DMLabelGetValue(sbr->splitPoints, maxedge, &eval);
2369: if (eval != 1) {
2370: DMLabelSetValue(sbr->splitPoints, maxedge, 1);
2371: PointQueueEnqueue(queue, maxedge);
2372: }
2373: DMLabelSetValue(sbr->splitPoints, cell, 2);
2374: }
2375: }
2376: return(0);
2377: }
2379: static PetscErrorCode SBRInitializeComm(DMPlexCellRefiner cr, PetscSF pointSF)
2380: {
2381: PlexRefiner_SBR *sbr = (PlexRefiner_SBR *) cr->data;
2382: DMLabel splitPoints = sbr->splitPoints;
2383: PetscInt *splitArray = sbr->splitArray;
2384: const PetscInt *degree;
2385: const PetscInt *points;
2386: PetscInt Nl, l, pStart, pEnd, p, val;
2387: PetscErrorCode ierr;
2390: /* Add in leaves */
2391: PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
2392: for (l = 0; l < Nl; ++l) {
2393: DMLabelGetValue(splitPoints, points[l], &val);
2394: if (val > 0) splitArray[points[l]] = val;
2395: }
2396: /* Add in shared roots */
2397: PetscSFComputeDegreeBegin(pointSF, °ree);
2398: PetscSFComputeDegreeEnd(pointSF, °ree);
2399: DMPlexGetChart(cr->dm, &pStart, &pEnd);
2400: for (p = pStart; p < pEnd; ++p) {
2401: if (degree[p]) {
2402: DMLabelGetValue(splitPoints, p, &val);
2403: if (val > 0) splitArray[p] = val;
2404: }
2405: }
2406: return(0);
2407: }
2409: static PetscErrorCode SBRFinalizeComm(DMPlexCellRefiner cr, PetscSF pointSF, PointQueue queue)
2410: {
2411: PlexRefiner_SBR *sbr = (PlexRefiner_SBR *) cr->data;
2412: DMLabel splitPoints = sbr->splitPoints;
2413: PetscInt *splitArray = sbr->splitArray;
2414: const PetscInt *degree;
2415: const PetscInt *points;
2416: PetscInt Nl, l, pStart, pEnd, p, val;
2417: PetscErrorCode ierr;
2420: /* Read out leaves */
2421: PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
2422: for (l = 0; l < Nl; ++l) {
2423: const PetscInt p = points[l];
2424: const PetscInt cval = splitArray[p];
2426: if (cval) {
2427: DMLabelGetValue(splitPoints, p, &val);
2428: if (val <= 0) {
2429: DMLabelSetValue(splitPoints, p, cval);
2430: PointQueueEnqueue(queue, p);
2431: }
2432: }
2433: }
2434: /* Read out shared roots */
2435: PetscSFComputeDegreeBegin(pointSF, °ree);
2436: PetscSFComputeDegreeEnd(pointSF, °ree);
2437: DMPlexGetChart(cr->dm, &pStart, &pEnd);
2438: for (p = pStart; p < pEnd; ++p) {
2439: if (degree[p]) {
2440: const PetscInt cval = splitArray[p];
2442: if (cval) {
2443: DMLabelGetValue(splitPoints, p, &val);
2444: if (val <= 0) {
2445: DMLabelSetValue(splitPoints, p, cval);
2446: PointQueueEnqueue(queue, p);
2447: }
2448: }
2449: }
2450: }
2451: return(0);
2452: }
2454: static PetscErrorCode DMPlexCellRefinerSetUp_SBR(DMPlexCellRefiner cr)
2455: {
2456: PlexRefiner_SBR *sbr;
2457: PetscSF pointSF;
2458: PointQueue queue = NULL;
2459: IS refineIS;
2460: const PetscInt *refineCells;
2461: PetscMPIInt size;
2462: PetscInt pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
2463: PetscBool empty;
2464: PetscErrorCode ierr;
2467: PetscCitationsRegister(SBRCitation, &SBRcite);
2468: PetscNew(&sbr);
2469: cr->data = sbr;
2470: DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints);
2471: /* Create edge lengths */
2472: DMPlexGetDepthStratum(cr->dm, 1, &eStart, &eEnd);
2473: PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen);
2474: PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd);
2475: for (e = eStart; e < eEnd; ++e) {
2476: PetscSectionSetDof(sbr->secEdgeLen, e, 1);
2477: }
2478: PetscSectionSetUp(sbr->secEdgeLen);
2479: PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize);
2480: PetscCalloc1(edgeLenSize, &sbr->edgeLen);
2481: /* Add edges of cells that are marked for refinement to edge queue */
2482: if (!cr->adaptLabel) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "CellRefiner must have an adaptation label in order to use SBR algorithm");
2483: PointQueueCreate(1024, &queue);
2484: DMLabelGetStratumIS(cr->adaptLabel, DM_ADAPT_REFINE, &refineIS);
2485: DMLabelGetStratumSize(cr->adaptLabel, DM_ADAPT_REFINE, &Nc);
2486: if (refineIS) {ISGetIndices(refineIS, &refineCells);}
2487: for (c = 0; c < Nc; ++c) {
2488: PetscInt *closure = NULL;
2489: PetscInt Ncl, cl;
2491: DMLabelSetValue(sbr->splitPoints, refineCells[c], 2);
2492: DMPlexGetTransitiveClosure(cr->dm, refineCells[c], PETSC_TRUE, &Ncl, &closure);
2493: for (cl = 0; cl < Ncl; cl += 2) {
2494: const PetscInt edge = closure[cl];
2496: if (edge >= eStart && edge < eEnd) {
2497: DMLabelSetValue(sbr->splitPoints, edge, 1);
2498: PointQueueEnqueue(queue, edge);
2499: }
2500: }
2501: DMPlexRestoreTransitiveClosure(cr->dm, refineCells[c], PETSC_TRUE, &Ncl, &closure);
2502: }
2503: if (refineIS) {ISRestoreIndices(refineIS, &refineCells);}
2504: ISDestroy(&refineIS);
2505: /* Setup communication */
2506: MPI_Comm_size(PetscObjectComm((PetscObject) cr->dm), &size);
2507: DMGetPointSF(cr->dm, &pointSF);
2508: if (size > 1) {
2509: PetscInt pStart, pEnd;
2511: DMPlexGetChart(cr->dm, &pStart, &pEnd);
2512: PetscCalloc1(pEnd-pStart, &sbr->splitArray);
2513: }
2514: /* While edge queue is not empty: */
2515: empty = PointQueueEmpty(queue);
2516: MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) cr->dm));
2517: while (!empty) {
2518: SBRSplitLocalEdges_Private(cr, queue);
2519: /* Communicate marked edges
2520: An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
2521: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2523: TODO: We could use in-place communication with a different SF
2524: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2525: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2527: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2528: values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new
2529: edge to the queue.
2530: */
2531: if (size > 1) {
2532: SBRInitializeComm(cr, pointSF);
2533: PetscSFReduceBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
2534: PetscSFReduceEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
2535: PetscSFBcastBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
2536: PetscSFBcastEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
2537: SBRFinalizeComm(cr, pointSF, queue);
2538: }
2539: empty = PointQueueEmpty(queue);
2540: MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) cr->dm));
2541: }
2542: PetscFree(sbr->splitArray);
2543: /* Calculate refineType for each cell */
2544: DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &cr->refineType);
2545: DMPlexGetChart(cr->dm, &pStart, &pEnd);
2546: for (p = pStart; p < pEnd; ++p) {
2547: DMPolytopeType ct;
2548: PetscInt val;
2550: DMPlexGetCellType(cr->dm, p, &ct);
2551: switch (ct) {
2552: case DM_POLYTOPE_POINT:
2553: DMLabelSetValue(cr->refineType, p, 0);break;
2554: case DM_POLYTOPE_SEGMENT:
2555: DMLabelGetValue(sbr->splitPoints, p, &val);
2556: if (val == 1) {DMLabelSetValue(cr->refineType, p, 2);}
2557: else {DMLabelSetValue(cr->refineType, p, 1);}
2558: break;
2559: case DM_POLYTOPE_TRIANGLE:
2560: DMLabelGetValue(sbr->splitPoints, p, &val);
2561: if (val == 2) {
2562: const PetscInt *cone;
2563: PetscReal lens[3];
2564: PetscInt vals[3], i;
2566: DMPlexGetCone(cr->dm, p, &cone);
2567: for (i = 0; i < 3; ++i) {
2568: DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]);
2569: vals[i] = vals[i] < 0 ? 0 : vals[i];
2570: SBRGetEdgeLen_Private(cr, cone[i], &lens[i]);
2571: }
2572: if (vals[0] && vals[1] && vals[2]) {DMLabelSetValue(cr->refineType, p, 4);}
2573: else if (vals[0] && vals[1]) {DMLabelSetValue(cr->refineType, p, lens[0] > lens[1] ? 5 : 6);}
2574: else if (vals[1] && vals[2]) {DMLabelSetValue(cr->refineType, p, lens[1] > lens[2] ? 7 : 8);}
2575: else if (vals[2] && vals[0]) {DMLabelSetValue(cr->refineType, p, lens[2] > lens[0] ? 9: 10);}
2576: else if (vals[0]) {DMLabelSetValue(cr->refineType, p, 11);}
2577: else if (vals[1]) {DMLabelSetValue(cr->refineType, p, 12);}
2578: else if (vals[2]) {DMLabelSetValue(cr->refineType, p, 13);}
2579: else SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D does not fit any refinement type (%D, %D, %D)", p, vals[0], vals[1], vals[2]);
2580: } else {DMLabelSetValue(cr->refineType, p, 3);}
2581: break;
2582: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
2583: }
2584: DMLabelGetValue(sbr->splitPoints, p, &val);
2585: }
2586: /* Cleanup */
2587: PointQueueDestroy(&queue);
2588: return(0);
2589: }
2591: static PetscErrorCode DMPlexCellRefinerDestroy_SBR(DMPlexCellRefiner cr)
2592: {
2593: PlexRefiner_SBR *sbr = (PlexRefiner_SBR *) cr->data;
2594: PetscErrorCode ierr;
2597: PetscFree(sbr->edgeLen);
2598: PetscSectionDestroy(&sbr->secEdgeLen);
2599: DMLabelDestroy(&sbr->splitPoints);
2600: PetscFree(cr->data);
2601: return(0);
2602: }
2604: static PetscErrorCode DMPlexCellRefinerRefine_SBR(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
2605: {
2606: /* Add 1 edge inside this triangle, making 2 new triangles.
2607: 2
2608: |\
2609: | \
2610: | \
2611: | \
2612: | 1
2613: | \
2614: | B \
2615: 2 1
2616: | / \
2617: | ____/ 0
2618: |/ A \
2619: 0-----0-----1
2620: */
2621: static DMPolytopeType triT10[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
2622: static PetscInt triS10[] = {1, 2};
2623: static PetscInt triC10[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
2624: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2625: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0};
2626: static PetscInt triO10[] = {0, 0,
2627: 0, 0, -2,
2628: 0, 0, 0};
2629: /* Add 1 edge inside this triangle, making 2 new triangles.
2630: 2
2631: |\
2632: | \
2633: | \
2634: 0 \
2635: | \
2636: | \
2637: | \
2638: 2 A 1
2639: |\ \
2640: 1 ---\ \
2641: | B ----\\
2642: 0-----0-----1
2643: */
2644: static PetscInt triC11[] = {DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0,
2645: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2646: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0};
2647: static PetscInt triO11[] = {0, 0,
2648: 0, 0, -2,
2649: 0, 0, 0};
2650: /* Add 1 edge inside this triangle, making 2 new triangles.
2651: 2
2652: |\
2653: |\\
2654: || \
2655: | \ \
2656: | | \
2657: | | \
2658: | | \
2659: 2 \ 1
2660: | A | \
2661: | | B \
2662: | \ \
2663: 0-----0-----1
2664: */
2665: static PetscInt triC12[] = {DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 1, 0, 0,
2666: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2667: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0};
2668: static PetscInt triO12[] = {0, 0,
2669: 0, 0, -2,
2670: 0, 0, 0};
2671: /* Add 2 edges inside this triangle, making 3 new triangles.
2672: 2
2673: |\
2674: | \
2675: | \
2676: 0 \
2677: | 1
2678: | \
2679: | B \
2680: 2-------1
2681: | C / \
2682: 1 ____/ 0
2683: |/ A \
2684: 0-----0-----1
2685: */
2686: static DMPolytopeType triT20[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
2687: static PetscInt triS20[] = {2, 3};
2688: static PetscInt triC20[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
2689: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
2690: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2691: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1,
2692: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2693: static PetscInt triO20[] = {0, 0,
2694: 0, 0,
2695: 0, 0, -2,
2696: 0, 0, -2,
2697: 0, 0, 0};
2698: /* Add 1 edge inside this triangle, making 2 new triangles.
2699: 2
2700: |\
2701: | \
2702: | \
2703: 0 \
2704: | 1
2705: | \
2706: | B \
2707: 2 1
2708: | /|\
2709: 1 ____/ / 0
2710: |/ A / C \
2711: 0-----0-----1
2712: */
2713: static PetscInt triC21[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
2714: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 0, 0,
2715: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
2716: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2717: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2718: static PetscInt triO21[] = {0, 0,
2719: 0, 0,
2720: 0, -2, -2,
2721: 0, 0, 0,
2722: 0, 0, 0};
2723: /* Add 2 edges inside this triangle, making 3 new triangles.
2724: 2
2725: |\
2726: | \
2727: | \
2728: 0 \
2729: | \
2730: | \
2731: | \
2732: 2 A 1
2733: |\ \
2734: 1 ---\ \
2735: |B \_C----\\
2736: 0-----0-----1
2737: */
2738: static PetscInt triC22[] = {DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0,
2739: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
2740: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2741: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1,
2742: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2743: static PetscInt triO22[] = {0, 0,
2744: 0, 0,
2745: 0, 0, -2,
2746: 0, 0, -2,
2747: 0, 0, 0};
2748: /* Add 1 edge inside this triangle, making 2 new triangles.
2749: 2
2750: |\
2751: | \
2752: | \
2753: 0 \
2754: | \
2755: | C \
2756: | \
2757: 2-------1
2758: |\ A \
2759: 1 ---\ \
2760: | B ----\\
2761: 0-----0-----1
2762: */
2763: static PetscInt triC23[] = {DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0,
2764: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 1, 0,
2765: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
2766: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2767: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2768: static PetscInt triO23[] = {0, 0,
2769: 0, 0,
2770: 0, -2, -2,
2771: 0, 0, 0,
2772: 0, 0, 0};
2773: /* Add 2 edges inside this triangle, making 3 new triangles.
2774: 2
2775: |\
2776: |\\
2777: || \
2778: | \ \
2779: | | \
2780: | | \
2781: | | \
2782: 2 \ C 1
2783: | A | / \
2784: | | |B \
2785: | \/ \
2786: 0-----0-----1
2787: */
2788: static PetscInt triC24[] = {DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 1, 0, 0,
2789: DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
2790: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2791: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1,
2792: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2793: static PetscInt triO24[] = {0, 0,
2794: 0, 0,
2795: 0, 0, -2,
2796: 0, 0, -2,
2797: 0, 0, 0};
2798: /* Add 1 edge inside this triangle, making 2 new triangles.
2799: 2
2800: |\
2801: |\\
2802: || \
2803: | \ \
2804: | | \
2805: | | \
2806: | | \
2807: 2 A \ 1
2808: |\ | \
2809: | \__| B \
2810: | C \\ \
2811: 0-----0-----1
2812: */
2813: static PetscInt triC25[] = {DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 1, 0, 0,
2814: DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 2, 0,
2815: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
2816: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
2817: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
2818: static PetscInt triO25[] = {0, 0,
2819: 0, 0,
2820: 0, -2, -2,
2821: 0, 0, 0,
2822: 0, 0, 0};
2823: DMLabel rtype = cr->refineType;
2824: PetscInt val;
2829: if (p < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
2830: DMLabelGetValue(rtype, p, &val);
2831: if (rt) *rt = val;
2832: switch (source) {
2833: case DM_POLYTOPE_POINT:
2834: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2835: case DM_POLYTOPE_QUADRILATERAL:
2836: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2837: case DM_POLYTOPE_TETRAHEDRON:
2838: case DM_POLYTOPE_HEXAHEDRON:
2839: case DM_POLYTOPE_TRI_PRISM:
2840: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2841: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2842: case DM_POLYTOPE_PYRAMID:
2843: DMPlexCellRefinerRefine_None(cr, source, p, NULL, Nt, target, size, cone, ornt);
2844: break;
2845: case DM_POLYTOPE_SEGMENT:
2846: if (val == 1) {DMPlexCellRefinerRefine_None(cr, source, p, NULL, Nt, target, size, cone, ornt);}
2847: else {DMPlexCellRefinerRefine_Regular(cr, source, p, NULL, Nt, target, size, cone, ornt);}
2848: break;
2849: case DM_POLYTOPE_TRIANGLE:
2850: switch (val) {
2851: case 12: *Nt = 2; *target = triT10; *size = triS10; *cone = triC10; *ornt = triO10; break;
2852: case 13: *Nt = 2; *target = triT10; *size = triS10; *cone = triC11; *ornt = triO11; break;
2853: case 11: *Nt = 2; *target = triT10; *size = triS10; *cone = triC12; *ornt = triO12; break;
2854: case 5: *Nt = 2; *target = triT20; *size = triS20; *cone = triC24; *ornt = triO24; break;
2855: case 6: *Nt = 2; *target = triT20; *size = triS20; *cone = triC21; *ornt = triO21; break;
2856: case 7: *Nt = 2; *target = triT20; *size = triS20; *cone = triC20; *ornt = triO20; break;
2857: case 8: *Nt = 2; *target = triT20; *size = triS20; *cone = triC23; *ornt = triO23; break;
2858: case 9: *Nt = 2; *target = triT20; *size = triS20; *cone = triC22; *ornt = triO22; break;
2859: case 10: *Nt = 2; *target = triT20; *size = triS20; *cone = triC25; *ornt = triO25; break;
2860: case 4: DMPlexCellRefinerRefine_Regular(cr, source, p, NULL, Nt, target, size, cone, ornt); break;
2861: default: DMPlexCellRefinerRefine_None(cr, source, p, NULL, Nt, target, size, cone, ornt);
2862: }
2863: break;
2864: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
2865: }
2866: return(0);
2867: }
2869: static PetscErrorCode DMPlexCellRefinerMapSubcells_SBR(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt pp, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
2870: {
2871: /* We shift any input orientation in order to make it non-negative
2872: 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)
2873: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
2874: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
2875: */
2876: PetscInt rt;
2877: PetscErrorCode ierr;
2880: DMLabelGetValue(cr->refineType, pp, &rt);
2881: *rnew = r;
2882: *onew = o;
2883: switch (rt) {
2884: case 5:
2885: case 6:
2886: case 7:
2887: case 8:
2888: case 9:
2889: case 10:
2890: switch (ct) {
2891: case DM_POLYTOPE_SEGMENT:
2892: break;
2893: case DM_POLYTOPE_TRIANGLE:
2894: break;
2895: default: break;
2896: }
2897: break;
2898: case 11:
2899: case 12:
2900: case 13:
2901: switch (ct) {
2902: case DM_POLYTOPE_SEGMENT:
2903: break;
2904: case DM_POLYTOPE_TRIANGLE:
2905: *onew = po < 0 ? -(o+1) : o;
2906: *rnew = po < 0 ? (r+1)%2 : r;
2907: break;
2908: default: break;
2909: }
2910: break;
2911: case 2:
2912: case 4:
2913: DMPlexCellRefinerMapSubcells_Regular(cr, pct, pp, po, ct, r, o, rnew, onew);
2914: break;
2915: default: DMPlexCellRefinerMapSubcells_None(cr, pct, pp, po, ct, r, o, rnew, onew);
2916: }
2917: return(0);
2918: }
2920: static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
2921: {
2922: PetscInt c, cN, *off;
2926: if (cr->refineType) {
2927: IS rtIS;
2928: const PetscInt *reftypes;
2929: PetscInt Nrt, r;
2931: DMLabelGetNumValues(cr->refineType, &Nrt);
2932: DMLabelGetValueIS(cr->refineType, &rtIS);
2933: ISGetIndices(rtIS, &reftypes);
2934: PetscCalloc1(Nrt*DM_NUM_POLYTOPES, &off);
2935: for (r = 0; r < Nrt; ++r) {
2936: const PetscInt rt = reftypes[r];
2937: IS rtIS;
2938: const PetscInt *points;
2939: DMPolytopeType ct;
2940: PetscInt p;
2942: DMLabelGetStratumIS(cr->refineType, rt, &rtIS);
2943: ISGetIndices(rtIS, &points);
2944: p = points[0];
2945: ISRestoreIndices(rtIS, &points);
2946: ISDestroy(&rtIS);
2947: DMPlexGetCellType(cr->dm, p, &ct);
2948: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
2949: const DMPolytopeType ctNew = (DMPolytopeType) cN;
2950: DMPolytopeType *rct;
2951: PetscInt *rsize, *cone, *ornt;
2952: PetscInt Nct, n, s;
2954: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[r*DM_NUM_POLYTOPES+ctNew] = -1; break;}
2955: off[r*DM_NUM_POLYTOPES+ctNew] = 0;
2956: for (s = 0; s <= r; ++s) {
2957: const PetscInt st = reftypes[s];
2958: DMPolytopeType sct;
2959: PetscInt q, qrt;
2961: DMLabelGetStratumIS(cr->refineType, st, &rtIS);
2962: ISGetIndices(rtIS, &points);
2963: q = points[0];
2964: ISRestoreIndices(rtIS, &points);
2965: ISDestroy(&rtIS);
2966: DMPlexGetCellType(cr->dm, q, &sct);
2967: DMPlexCellRefinerRefine(cr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
2968: if (st != qrt) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %D of point %D does not match predicted type %D", qrt, q, st);
2969: if (st == rt) {
2970: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
2971: if (n == Nct) off[r*DM_NUM_POLYTOPES+ctNew] = -1;
2972: break;
2973: }
2974: for (n = 0; n < Nct; ++n) {
2975: if (rct[n] == ctNew) {
2976: PetscInt sn;
2978: DMLabelGetStratumSize(cr->refineType, st, &sn);
2979: off[r*DM_NUM_POLYTOPES+ctNew] += sn * rsize[n];
2980: }
2981: }
2982: }
2983: }
2984: }
2985: #if 0
2986: {
2987: PetscInt cols = 8;
2988: PetscInt f, g;
2989: PetscPrintf(PETSC_COMM_SELF, "Offsets\n");
2990: PetscPrintf(PETSC_COMM_SELF, " ");
2991: for (g = 0; g < cols; ++g) {
2992: PetscPrintf(PETSC_COMM_SELF, " % 14s", DMPolytopeTypes[g]);
2993: }
2994: PetscPrintf(PETSC_COMM_SELF, "\n");
2995: for (f = 0; f < Nrt; ++f) {
2996: PetscPrintf(PETSC_COMM_SELF, "%2d |", reftypes[f]);
2997: for (g = 0; g < cols; ++g) {
2998: PetscPrintf(PETSC_COMM_SELF, " % 14d", PetscRealPart(off[f*DM_NUM_POLYTOPES+g]));
2999: }
3000: PetscPrintf(PETSC_COMM_SELF, " |\n");
3001: }
3002: }
3003: #endif
3004: ISRestoreIndices(rtIS, &reftypes);
3005: ISDestroy(&rtIS);
3006: } else {
3007: PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
3008: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
3009: const DMPolytopeType ct = (DMPolytopeType) c;
3010: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
3011: const DMPolytopeType ctNew = (DMPolytopeType) cN;
3012: DMPolytopeType *rct;
3013: PetscInt *rsize, *cone, *ornt;
3014: PetscInt Nct, n, i;
3016: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;}
3017: off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
3018: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
3019: const DMPolytopeType ict = (DMPolytopeType) ctOrder[i];
3020: const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];
3022: DMPlexCellRefinerRefine(cr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt);
3023: if (ict == ct) {
3024: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
3025: if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
3026: break;
3027: }
3028: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
3029: }
3030: }
3031: }
3032: }
3033: *offset = off;
3034: return(0);
3035: }
3037: static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[])
3038: {
3039: const PetscInt ctSize = DM_NUM_POLYTOPES+1;
3043: if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()");
3044: PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);
3045: PetscArraycpy(cr->ctStart, ctStart, ctSize);
3046: PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);
3047: return(0);
3048: }
3050: /* Construct cell type order since we must loop over cell types in depth order */
3051: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
3052: {
3053: PetscInt *ctO, *ctOInv;
3054: PetscInt c, d, off = 0;
3058: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
3059: for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
3060: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
3061: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
3062: ctO[off++] = c;
3063: }
3064: }
3065: if (DMPolytopeTypeGetDim(ctCell) != 0) {
3066: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
3067: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
3068: ctO[off++] = c;
3069: }
3070: }
3071: for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
3072: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
3073: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
3074: ctO[off++] = c;
3075: }
3076: }
3077: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
3078: if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
3079: ctO[off++] = c;
3080: }
3081: if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
3082: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
3083: ctOInv[ctO[c]] = c;
3084: }
3085: *ctOrder = ctO;
3086: *ctOrderInv = ctOInv;
3087: return(0);
3088: }
3090: PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)
3091: {
3092: DM dm = cr->dm;
3093: DMPolytopeType ctCell;
3094: PetscInt pStart, pEnd, p, c;
3099: if (cr->setupcalled) return(0);
3100: if (cr->ops->setup) {
3101: (*cr->ops->setup)(cr);
3102: }
3103: DMPlexGetChart(dm, &pStart, &pEnd);
3104: if (pEnd > pStart) {
3105: DMPlexGetCellType(dm, 0, &ctCell);
3106: } else {
3107: PetscInt dim;
3108: DMGetDimension(dm, &dim);
3109: switch (dim) {
3110: case 0: ctCell = DM_POLYTOPE_POINT;break;
3111: case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
3112: case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
3113: case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
3114: default: ctCell = DM_POLYTOPE_UNKNOWN;
3115: }
3116: }
3117: DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);
3118: /* Construct sizes and offsets for each cell type */
3119: if (!cr->ctStart) {
3120: PetscInt *ctS, *ctSN, *ctC, *ctCN;
3122: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
3123: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
3124: for (p = pStart; p < pEnd; ++p) {
3125: DMPolytopeType ct;
3126: DMPolytopeType *rct;
3127: PetscInt *rsize, *cone, *ornt;
3128: PetscInt Nct, n;
3130: DMPlexGetCellType(dm, p, &ct);
3131: if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
3132: ++ctC[ct];
3133: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
3134: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
3135: }
3136: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
3137: const PetscInt ct = cr->ctOrder[c];
3138: const PetscInt ctn = cr->ctOrder[c+1];
3140: ctS[ctn] = ctS[ct] + ctC[ct];
3141: ctSN[ctn] = ctSN[ct] + ctCN[ct];
3142: }
3143: PetscFree2(ctC, ctCN);
3144: cr->ctStart = ctS;
3145: cr->ctStartNew = ctSN;
3146: }
3147: CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);
3148: cr->setupcalled = PETSC_TRUE;
3149: return(0);
3150: }
3152: static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v)
3153: {
3157: PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);
3158: return(0);
3159: }
3161: /*
3162: DMPlexCellRefinerView - Views a DMPlexCellRefiner object
3164: Collective on cr
3166: Input Parameters:
3167: + cr - The DMPlexCellRefiner object
3168: - viewer - The PetscViewer object
3170: Level: beginner
3172: .seealso: DMPlexCellRefinerCreate()
3173: */
3174: static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer)
3175: {
3176: PetscBool iascii;
3182: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);}
3183: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3184: PetscViewerASCIIPushTab(viewer);
3185: if (iascii) {DMPlexCellRefinerView_Ascii(cr, viewer);}
3186: PetscViewerASCIIPopTab(viewer);
3187: return(0);
3188: }
3190: PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr)
3191: {
3192: PetscInt c;
3196: if (!*cr) return(0);
3198: if ((*cr)->ops->destroy) {
3199: ((*cr)->ops->destroy)(*cr);
3200: }
3201: PetscObjectDereference((PetscObject) (*cr)->dm);
3202: DMLabelDestroy(&(*cr)->refineType);
3203: PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);
3204: PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);
3205: PetscFree((*cr)->offset);
3206: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
3207: PetscFEDestroy(&(*cr)->coordFE[c]);
3208: PetscFEGeomDestroy(&(*cr)->refGeom[c]);
3209: }
3210: PetscFree2((*cr)->coordFE, (*cr)->refGeom);
3211: PetscHeaderDestroy(cr);
3212: return(0);
3213: }
3215: PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr)
3216: {
3217: DMPlexCellRefiner tmp;
3218: PetscErrorCode ierr;
3222: *cr = NULL;
3223: PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);
3224: tmp->setupcalled = PETSC_FALSE;
3226: tmp->dm = dm;
3227: PetscObjectReference((PetscObject) dm);
3228: DMPlexGetCellRefinerType(dm, &tmp->type);
3229: switch (tmp->type) {
3230: case DM_REFINER_REGULAR:
3231: tmp->ops->refine = DMPlexCellRefinerRefine_Regular;
3232: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_Regular;
3233: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_Regular;
3234: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_Regular;
3235: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3236: tmp->ops->getaffinetransforms = DMPlexCellRefinerGetAffineTransforms_Regular;
3237: tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular;
3238: break;
3239: case DM_REFINER_TO_BOX:
3240: tmp->ops->refine = DMPlexCellRefinerRefine_ToBox;
3241: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToBox;
3242: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_ToBox;
3243: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox;
3244: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3245: break;
3246: case DM_REFINER_TO_SIMPLEX:
3247: tmp->ops->refine = DMPlexCellRefinerRefine_ToSimplex;
3248: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex;
3249: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3250: break;
3251: case DM_REFINER_ALFELD2D:
3252: tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld2D;
3253: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
3254: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3255: break;
3256: case DM_REFINER_ALFELD3D:
3257: tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld3D;
3258: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
3259: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3260: break;
3261: case DM_REFINER_BOUNDARYLAYER:
3262: tmp->ops->setup = DMPlexCellRefinerSetUp_BL;
3263: tmp->ops->destroy = DMPlexCellRefinerDestroy_BL;
3264: tmp->ops->refine = DMPlexCellRefinerRefine_BL;
3265: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_BL;
3266: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_BL;
3267: break;
3268: case DM_REFINER_SBR:
3269: tmp->ops->setup = DMPlexCellRefinerSetUp_SBR;
3270: tmp->ops->destroy = DMPlexCellRefinerDestroy_SBR;
3271: tmp->ops->refine = DMPlexCellRefinerRefine_SBR;
3272: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_SBR;
3273: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
3274: break;
3275: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]);
3276: }
3277: PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);
3278: *cr = tmp;
3279: return(0);
3280: }
3282: /*@
3283: DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell
3285: Input Parameters:
3286: + cr - The DMPlexCellRefiner object
3287: - ct - The cell type
3289: Output Parameters:
3290: + Nc - The number of subcells produced from this cell type
3291: . v0 - The translation of the first vertex for each subcell
3292: . J - The Jacobian for each subcell (map from reference cell to subcell)
3293: - invJ - The inverse Jacobian for each subcell
3295: Level: developer
3297: .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create()
3298: @*/
3299: PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
3300: {
3304: if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner");
3305: (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);
3306: return(0);
3307: }
3309: /*@
3310: DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell
3312: Input Parameters:
3313: + cr - The DMPlexCellRefiner object
3314: - ct - The cell type
3316: Output Parameters:
3317: + Nf - The number of faces for this cell type
3318: . v0 - The translation of the first vertex for each face
3319: . J - The Jacobian for each face (map from original cell to subcell)
3320: . invJ - The inverse Jacobian for each face
3321: - detJ - The determinant of the Jacobian for each face
3323: Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
3325: Level: developer
3327: .seealso: DMPlexCellRefinerGetAffineTransforms(), Create()
3328: @*/
3329: PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
3330: {
3334: if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner");
3335: (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);
3336: return(0);
3337: }
3339: /* Numbering regularly refined meshes
3341: We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute
3342: the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions
3343: about the order of different cell types.
3345: However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous
3346: numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we
3347: will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as
3348: the cell type enumeration within a given depth stratum.
3350: Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we
3351: start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then
3352: offset by the old cell type number and replica number for the insertion.
3353: */
3355: static PetscErrorCode DMPlexCellRefinerGetReducedPointNumber(DMPlexCellRefiner cr, PetscInt rt, PetscInt p, PetscInt *rp)
3356: {
3357: IS rtIS;
3358: const PetscInt *points;
3359: PetscInt n;
3360: PetscErrorCode ierr;
3363: /* TODO Move this inside the DMLabel so that I do not have to create the IS */
3364: DMLabelGetStratumIS(cr->refineType, rt, &rtIS);
3365: ISGetLocalSize(rtIS, &n);
3366: ISGetIndices(rtIS, &points);
3367: PetscFindInt(p, n, points, rp);
3368: ISRestoreIndices(rtIS, &points);
3369: ISDestroy(&rtIS);
3370: return(0);
3371: }
3373: /*
3374: DMPlexCellRefinerGetNewPoint - Get the number of a point in the refined mesh based on information from the original mesh.
3376: Not collective
3378: Input Parameters:
3379: + cr - The cell refiner
3380: . ct - The type of the original point which produces the new point
3381: . ctNew - The type of the new point
3382: . p - The original point which produces the new point
3383: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
3385: Output Parameters:
3386: . pNew - The new point number
3388: Level: developer
3390: .seealso: DMPlexCellRefinerRefine()
3391: */
3392: PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
3393: {
3394: DMPolytopeType *rct;
3395: PetscInt *rsize, *cone, *ornt;
3396: PetscInt rt, Nct, n, off, rp;
3397: PetscInt ctS = cr->ctStart[ct], ctE = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]];
3398: PetscInt ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]];
3399: PetscInt newp = ctSN, cind;
3403: 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);
3404: DMPlexCellRefinerRefine(cr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
3405: if (cr->refineType) {
3406: /* TODO Make this a function in DMLabel */
3407: {
3408: IS rtIS;
3409: const PetscInt *reftypes;
3410: PetscInt Nrt;
3412: DMLabelGetNumValues(cr->refineType, &Nrt);
3413: DMLabelGetValueIS(cr->refineType, &rtIS);
3414: ISGetIndices(rtIS, &reftypes);
3415: for (cind = 0; cind < Nrt; ++cind) if (reftypes[cind] == rt) break;
3416: if (cind >= Nrt) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unable to locate refine type %D", rt);
3417: ISRestoreIndices(rtIS, &reftypes);
3418: ISDestroy(&rtIS);
3419: }
3420: DMPlexCellRefinerGetReducedPointNumber(cr, rt, p, &rp);
3421: if (rp < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %D does not have refine type %D", DMPolytopeTypes[ct], p, rt);
3422: } else {
3423: cind = ct;
3424: rp = p - ctS;
3425: }
3426: off = cr->offset[cind*DM_NUM_POLYTOPES + ctNew];
3427: if (off < 0) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%D) of point %D does not produce type %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew]);
3428: newp += off;
3429: for (n = 0; n < Nct; ++n) {
3430: if (rct[n] == ctNew) {
3431: if (rsize[n] && r >= rsize[n])
3432: 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]);
3433: newp += rp * rsize[n] + r;
3434: break;
3435: }
3436: }
3438: 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);
3439: *pNew = newp;
3440: return(0);
3441: }
3443: static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm)
3444: {
3445: DM dm = cr->dm;
3446: PetscInt pStart, pEnd, p, pNew;
3447: PetscErrorCode ierr;
3450: /* Must create the celltype label here so that we do not automatically try to compute the types */
3451: DMCreateLabel(rdm, "celltype");
3452: DMPlexGetChart(dm, &pStart, &pEnd);
3453: for (p = pStart; p < pEnd; ++p) {
3454: DMPolytopeType ct;
3455: DMPolytopeType *rct;
3456: PetscInt *rsize, *rcone, *rornt;
3457: PetscInt Nct, n, r;
3459: DMPlexGetCellType(dm, p, &ct);
3460: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3461: for (n = 0; n < Nct; ++n) {
3462: for (r = 0; r < rsize[n]; ++r) {
3463: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3464: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
3465: DMPlexSetCellType(rdm, pNew, rct[n]);
3466: }
3467: }
3468: }
3469: {
3470: DMLabel ctLabel;
3471: DM_Plex *plex = (DM_Plex *) rdm->data;
3473: DMPlexGetCellTypeLabel(rdm, &ctLabel);
3474: PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
3475: }
3476: return(0);
3477: }
3479: static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm)
3480: {
3481: DM dm = cr->dm;
3482: DMPolytopeType ct;
3483: PetscInt *coneNew, *orntNew;
3484: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
3488: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
3489: PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
3490: DMPlexGetChart(dm, &pStart, &pEnd);
3491: for (p = pStart; p < pEnd; ++p) {
3492: const PetscInt *cone, *ornt;
3493: PetscInt coff, ooff, c;
3494: DMPolytopeType *rct;
3495: PetscInt *rsize, *rcone, *rornt;
3496: PetscInt Nct, n, r;
3497: DMPlexGetCellType(dm, p, &ct);
3498: DMPlexGetCone(dm, p, &cone);
3499: DMPlexGetConeOrientation(dm, p, &ornt);
3500: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3501: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
3502: const DMPolytopeType ctNew = rct[n];
3503: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
3505: for (r = 0; r < rsize[n]; ++r) {
3506: /* pNew is a subcell produced by subdividing p */
3507: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3508: for (c = 0; c < csizeNew; ++c) {
3509: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
3510: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
3511: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
3512: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
3513: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
3514: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
3515: const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
3516: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
3517: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
3518: PetscInt lc;
3520: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
3521: for (lc = 0; lc < fn; ++lc) {
3522: const PetscInt *ppornt;
3523: PetscInt pcp;
3525: DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);
3526: ppp = pp;
3527: pp = pcone[pcp];
3528: DMPlexGetCellType(dm, pp, &pct);
3529: DMPlexGetCone(dm, pp, &pcone);
3530: DMPlexGetConeOrientation(dm, ppp, &ppornt);
3531: if (po < 0 && pct != DM_POLYTOPE_POINT) {
3532: const PetscInt pornt = ppornt[pcp];
3533: const PetscInt pcsize = DMPolytopeTypeGetConeSize(pct);
3534: const PetscInt pcstart = pornt < 0 ? -(pornt+1) : pornt;
3535: const PetscInt rcstart = (pcstart+pcsize-1)%pcsize;
3536: po = pornt < 0 ? -(rcstart+1) : rcstart;
3537: } else {
3538: po = ppornt[pcp];
3539: }
3540: }
3541: pr = rcone[coff++];
3542: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
3543: DMPlexCellRefinerMapSubcells(cr, pct, pp, po, ft, pr, fo, &pr, &fo);
3544: DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);
3545: orntNew[c] = fo;
3546: }
3547: DMPlexSetCone(rdm, pNew, coneNew);
3548: DMPlexSetConeOrientation(rdm, pNew, orntNew);
3549: }
3550: }
3551: }
3552: PetscFree2(coneNew, orntNew);
3553: DMViewFromOptions(rdm, NULL, "-rdm_view");
3554: DMPlexSymmetrize(rdm);
3555: DMPlexStratify(rdm);
3556: return(0);
3557: }
3559: static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe)
3560: {
3564: if (!cr->coordFE[ct]) {
3565: PetscInt dim, cdim;
3566: PetscBool isSimplex;
3568: switch (ct) {
3569: case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break;
3570: case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break;
3571: case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
3572: case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break;
3573: case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break;
3574: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
3575: }
3576: DMGetCoordinateDim(cr->dm, &cdim);
3577: PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);
3578: {
3579: PetscDualSpace dsp;
3580: PetscQuadrature quad;
3581: DM K;
3582: PetscFEGeom *cg;
3583: PetscReal *Xq, *xq, *wq;
3584: PetscInt Nq, q;
3586: DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);
3587: PetscMalloc1(Nq*cdim, &xq);
3588: for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q];
3589: PetscMalloc1(Nq, &wq);
3590: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
3591: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
3592: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
3593: PetscFESetQuadrature(cr->coordFE[ct], quad);
3595: PetscFEGetDualSpace(cr->coordFE[ct], &dsp);
3596: PetscDualSpaceGetDM(dsp, &K);
3597: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);
3598: cg = cr->refGeom[ct];
3599: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
3600: PetscQuadratureDestroy(&quad);
3601: }
3602: }
3603: *fe = cr->coordFE[ct];
3604: return(0);
3605: }
3607: /*
3608: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
3610: Not collective
3612: Input Parameters:
3613: + cr - The DMPlexCellRefiner
3614: . ct - The type of the parent cell
3615: . rct - The type of the produced cell
3616: . r - The index of the produced cell
3617: - x - The localized coordinates for the parent cell
3619: Output Parameter:
3620: . xr - The localized coordinates for the produced cell
3622: Level: developer
3624: .seealso: DMPlexCellRefinerSetCoordinates()
3625: */
3626: static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
3627: {
3628: PetscFE fe = NULL;
3629: PetscInt cdim, Nv, v, *subcellV;
3633: DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);
3634: DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);
3635: PetscFEGetNumComponents(fe, &cdim);
3636: for (v = 0; v < Nv; ++v) {
3637: PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);
3638: }
3639: return(0);
3640: }
3642: static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm)
3643: {
3644: DM dm = cr->dm, cdm;
3645: PetscSection coordSection, coordSectionNew;
3646: Vec coordsLocal, coordsLocalNew;
3647: const PetscScalar *coords;
3648: PetscScalar *coordsNew;
3649: const DMBoundaryType *bd;
3650: const PetscReal *maxCell, *L;
3651: PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
3652: PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
3653: PetscErrorCode ierr;
3656: DMGetCoordinateDM(dm, &cdm);
3657: DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
3658: /* Determine if we need to localize coordinates when generating them */
3659: if (isperiodic) {
3660: localizeVertices = PETSC_TRUE;
3661: if (!maxCell) {
3662: PetscBool localized;
3663: DMGetCoordinatesLocalized(dm, &localized);
3664: if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
3665: localizeCells = PETSC_TRUE;
3666: }
3667: }
3669: DMGetCoordinateSection(dm, &coordSection);
3670: PetscSectionGetFieldComponents(coordSection, 0, &dE);
3671: if (maxCell) {
3672: PetscReal maxCellNew[3];
3674: for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
3675: DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
3676: } else {
3677: DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
3678: }
3679: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
3680: PetscSectionSetNumFields(coordSectionNew, 1);
3681: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
3682: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
3683: if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0, vEndNew);}
3684: else {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}
3686: /* Localization should be inherited */
3687: /* Stefano calculates parent cells for each new cell for localization */
3688: /* Localized cells need coordinates of closure */
3689: for (v = vStartNew; v < vEndNew; ++v) {
3690: PetscSectionSetDof(coordSectionNew, v, dE);
3691: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
3692: }
3693: if (localizeCells) {
3694: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3695: for (c = cStart; c < cEnd; ++c) {
3696: PetscInt dof;
3698: PetscSectionGetDof(coordSection, c, &dof);
3699: if (dof) {
3700: DMPolytopeType ct;
3701: DMPolytopeType *rct;
3702: PetscInt *rsize, *rcone, *rornt;
3703: PetscInt dim, cNew, Nct, n, r;
3705: DMPlexGetCellType(dm, c, &ct);
3706: dim = DMPolytopeTypeGetDim(ct);
3707: DMPlexCellRefinerRefine(cr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3708: /* This allows for different cell types */
3709: for (n = 0; n < Nct; ++n) {
3710: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
3711: for (r = 0; r < rsize[n]; ++r) {
3712: PetscInt *closure = NULL;
3713: PetscInt clSize, cl, Nv = 0;
3715: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);
3716: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
3717: for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
3718: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
3719: PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
3720: PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
3721: }
3722: }
3723: }
3724: }
3725: }
3726: PetscSectionSetUp(coordSectionNew);
3727: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
3728: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
3729: {
3730: VecType vtype;
3731: PetscInt coordSizeNew, bs;
3732: const char *name;
3734: DMGetCoordinatesLocal(dm, &coordsLocal);
3735: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
3736: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
3737: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
3738: PetscObjectGetName((PetscObject) coordsLocal, &name);
3739: PetscObjectSetName((PetscObject) coordsLocalNew, name);
3740: VecGetBlockSize(coordsLocal, &bs);
3741: VecSetBlockSize(coordsLocalNew, bs);
3742: VecGetType(coordsLocal, &vtype);
3743: VecSetType(coordsLocalNew, vtype);
3744: }
3745: VecGetArrayRead(coordsLocal, &coords);
3746: VecGetArray(coordsLocalNew, &coordsNew);
3747: PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
3748: DMPlexGetChart(dm, &pStart, &pEnd);
3749: /* First set coordinates for vertices*/
3750: for (p = pStart; p < pEnd; ++p) {
3751: DMPolytopeType ct;
3752: DMPolytopeType *rct;
3753: PetscInt *rsize, *rcone, *rornt;
3754: PetscInt Nct, n, r;
3755: PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;
3757: DMPlexGetCellType(dm, p, &ct);
3758: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3759: for (n = 0; n < Nct; ++n) {
3760: if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
3761: }
3762: if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
3763: PetscInt dof;
3764: PetscSectionGetDof(coordSection, p, &dof);
3765: if (dof) isLocalized = PETSC_TRUE;
3766: }
3767: if (hasVertex) {
3768: const PetscScalar *icoords = NULL;
3769: PetscScalar *pcoords = NULL;
3770: PetscInt Nc, Nv, v, d;
3772: DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
3774: icoords = pcoords;
3775: Nv = Nc/dE;
3776: if (ct != DM_POLYTOPE_POINT) {
3777: if (localizeVertices) {
3778: PetscScalar anchor[3];
3780: for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
3781: if (!isLocalized) {
3782: for (v = 0; v < Nv; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
3783: } else {
3784: Nv = Nc/(2*dE);
3785: icoords = pcoords + Nv*dE;
3786: for (v = Nv; v < Nv*2; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
3787: }
3788: }
3789: }
3790: for (n = 0; n < Nct; ++n) {
3791: if (rct[n] != DM_POLYTOPE_POINT) continue;
3792: for (r = 0; r < rsize[n]; ++r) {
3793: PetscScalar vcoords[3];
3794: PetscInt vNew, off;
3796: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);
3797: PetscSectionGetOffset(coordSectionNew, vNew, &off);
3798: DMPlexCellRefinerMapCoordinates(cr, ct, rct[n], r, Nv, dE, icoords, vcoords);
3799: DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);
3800: }
3801: }
3802: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
3803: }
3804: }
3805: /* Then set coordinates for cells by localizing */
3806: for (p = pStart; p < pEnd; ++p) {
3807: DMPolytopeType ct;
3808: DMPolytopeType *rct;
3809: PetscInt *rsize, *rcone, *rornt;
3810: PetscInt Nct, n, r;
3811: PetscBool isLocalized = PETSC_FALSE;
3813: DMPlexGetCellType(dm, p, &ct);
3814: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3815: if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
3816: PetscInt dof;
3817: PetscSectionGetDof(coordSection, p, &dof);
3818: if (dof) isLocalized = PETSC_TRUE;
3819: }
3820: if (isLocalized) {
3821: const PetscScalar *pcoords;
3823: DMPlexPointLocalRead(cdm, p, coords, &pcoords);
3824: for (n = 0; n < Nct; ++n) {
3825: const PetscInt Nr = rsize[n];
3827: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
3828: for (r = 0; r < Nr; ++r) {
3829: PetscInt pNew, offNew;
3831: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
3832: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
3833: cell to the ones it produces. */
3834: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3835: PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
3836: DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
3837: }
3838: }
3839: }
3840: }
3841: VecRestoreArrayRead(coordsLocal, &coords);
3842: VecRestoreArray(coordsLocalNew, &coordsNew);
3843: DMSetCoordinatesLocal(rdm, coordsLocalNew);
3844: /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
3845: VecDestroy(&coordsLocalNew);
3846: PetscSectionDestroy(&coordSectionNew);
3847: if (!localizeCells) {DMLocalizeCoordinates(rdm);}
3848: return(0);
3849: }
3851: /*@
3852: DMPlexCreateProcessSF - Create an SF which just has process connectivity
3854: Collective on dm
3856: Input Parameters:
3857: + dm - The DM
3858: - sfPoint - The PetscSF which encodes point connectivity
3860: Output Parameters:
3861: + processRanks - A list of process neighbors, or NULL
3862: - sfProcess - An SF encoding the process connectivity, or NULL
3864: Level: developer
3866: .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
3867: @*/
3868: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
3869: {
3870: PetscInt numRoots, numLeaves, l;
3871: const PetscInt *localPoints;
3872: const PetscSFNode *remotePoints;
3873: PetscInt *localPointsNew;
3874: PetscSFNode *remotePointsNew;
3875: PetscInt *ranks, *ranksNew;
3876: PetscMPIInt size;
3877: PetscErrorCode ierr;
3884: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
3885: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3886: PetscMalloc1(numLeaves, &ranks);
3887: for (l = 0; l < numLeaves; ++l) {
3888: ranks[l] = remotePoints[l].rank;
3889: }
3890: PetscSortRemoveDupsInt(&numLeaves, ranks);
3891: PetscMalloc1(numLeaves, &ranksNew);
3892: PetscMalloc1(numLeaves, &localPointsNew);
3893: PetscMalloc1(numLeaves, &remotePointsNew);
3894: for (l = 0; l < numLeaves; ++l) {
3895: ranksNew[l] = ranks[l];
3896: localPointsNew[l] = l;
3897: remotePointsNew[l].index = 0;
3898: remotePointsNew[l].rank = ranksNew[l];
3899: }
3900: PetscFree(ranks);
3901: if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);}
3902: else {PetscFree(ranksNew);}
3903: if (sfProcess) {
3904: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
3905: PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
3906: PetscSFSetFromOptions(*sfProcess);
3907: PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
3908: }
3909: return(0);
3910: }
3912: static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm)
3913: {
3914: DM dm = cr->dm;
3915: DMPlexCellRefiner *crRem;
3916: PetscSF sf, sfNew, sfProcess;
3917: IS processRanks;
3918: MPI_Datatype ctType;
3919: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
3920: const PetscInt *localPoints, *neighbors;
3921: const PetscSFNode *remotePoints;
3922: PetscInt *localPointsNew;
3923: PetscSFNode *remotePointsNew;
3924: PetscInt *ctStartRem, *ctStartNewRem;
3925: PetscInt ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem;
3926: /* Brute force algorithm */
3927: PetscSF rsf;
3928: PetscSection s;
3929: const PetscInt *rootdegree;
3930: PetscInt *rootPointsNew, *remoteOffsets;
3931: PetscInt numPointsNew, pStart, pEnd, p;
3932: PetscErrorCode ierr;
3935: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
3936: DMGetPointSF(dm, &sf);
3937: DMGetPointSF(rdm, &sfNew);
3938: /* Calculate size of new SF */
3939: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
3940: if (numRoots < 0) return(0);
3941: for (l = 0; l < numLeaves; ++l) {
3942: const PetscInt p = localPoints[l];
3943: DMPolytopeType ct;
3944: DMPolytopeType *rct;
3945: PetscInt *rsize, *rcone, *rornt;
3946: PetscInt Nct, n;
3948: DMPlexGetCellType(dm, p, &ct);
3949: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3950: for (n = 0; n < Nct; ++n) {
3951: numLeavesNew += rsize[n];
3952: }
3953: }
3954: if (cr->refineType) {
3955: DMPlexGetChart(dm, &pStart, &pEnd);
3956: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);
3957: PetscSectionSetChart(s, pStart, pEnd);
3958: for (p = pStart; p < pEnd; ++p) {
3959: DMPolytopeType ct;
3960: DMPolytopeType *rct;
3961: PetscInt *rsize, *rcone, *rornt;
3962: PetscInt Nct, n;
3964: DMPlexGetCellType(dm, p, &ct);
3965: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3966: for (n = 0; n < Nct; ++n) {
3967: PetscSectionAddDof(s, p, rsize[n]);
3968: }
3969: }
3970: PetscSectionSetUp(s);
3971: PetscSectionGetStorageSize(s, &numPointsNew);
3972: PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
3973: PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
3974: PetscFree(remoteOffsets);
3975: PetscSFComputeDegreeBegin(sf, &rootdegree);
3976: PetscSFComputeDegreeEnd(sf, &rootdegree);
3977: PetscMalloc1(numPointsNew, &rootPointsNew);
3978: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
3979: for (p = pStart; p < pEnd; ++p) {
3980: DMPolytopeType ct;
3981: DMPolytopeType *rct;
3982: PetscInt *rsize, *rcone, *rornt;
3983: PetscInt Nct, n, r, off;
3985: if (!rootdegree[p-pStart]) continue;
3986: PetscSectionGetOffset(s, p, &off);
3987: DMPlexGetCellType(dm, p, &ct);
3988: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
3989: for (n = 0, m = 0; n < Nct; ++n) {
3990: for (r = 0; r < rsize[n]; ++r, ++m) {
3991: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3992: rootPointsNew[off+m] = pNew;
3993: }
3994: }
3995: }
3996: PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
3997: PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
3998: PetscSFDestroy(&rsf);
3999: PetscMalloc1(numLeavesNew, &localPointsNew);
4000: PetscMalloc1(numLeavesNew, &remotePointsNew);
4001: for (l = 0, m = 0; l < numLeaves; ++l) {
4002: const PetscInt p = localPoints[l];
4003: DMPolytopeType ct;
4004: DMPolytopeType *rct;
4005: PetscInt *rsize, *rcone, *rornt;
4006: PetscInt Nct, n, r, q, off;
4008: PetscSectionGetOffset(s, p, &off);
4009: DMPlexGetCellType(dm, p, &ct);
4010: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
4011: for (n = 0, q = 0; n < Nct; ++n) {
4012: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
4013: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
4014: localPointsNew[m] = pNew;
4015: remotePointsNew[m].index = rootPointsNew[off+q];
4016: remotePointsNew[m].rank = remotePoints[l].rank;
4017: }
4018: }
4019: }
4020: PetscSectionDestroy(&s);
4021: PetscFree(rootPointsNew);
4022: } else {
4023: /* Communicate ctStart and cStartNew for each remote rank */
4024: DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);
4025: ISGetLocalSize(processRanks, &numNeighbors);
4026: PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);
4027: MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);
4028: MPI_Type_commit(&ctType);
4029: PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem,MPI_REPLACE);
4030: PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem,MPI_REPLACE);
4031: PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem,MPI_REPLACE);
4032: PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem,MPI_REPLACE);
4033: MPI_Type_free(&ctType);
4034: PetscSFDestroy(&sfProcess);
4035: PetscMalloc1(numNeighbors, &crRem);
4036: for (n = 0; n < numNeighbors; ++n) {
4037: DMPlexCellRefinerCreate(dm, &crRem[n]);
4038: DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]);
4039: DMPlexCellRefinerSetUp(crRem[n]);
4040: }
4041: PetscFree2(ctStartRem, ctStartNewRem);
4042: /* Calculate new point SF */
4043: PetscMalloc1(numLeavesNew, &localPointsNew);
4044: PetscMalloc1(numLeavesNew, &remotePointsNew);
4045: ISGetIndices(processRanks, &neighbors);
4046: for (l = 0, m = 0; l < numLeaves; ++l) {
4047: PetscInt p = localPoints[l];
4048: PetscInt pRem = remotePoints[l].index;
4049: PetscMPIInt rankRem = remotePoints[l].rank;
4050: DMPolytopeType ct;
4051: DMPolytopeType *rct;
4052: PetscInt *rsize, *rcone, *rornt;
4053: PetscInt neighbor, Nct, n, r;
4055: PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);
4056: if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem);
4057: DMPlexGetCellType(dm, p, &ct);
4058: DMPlexCellRefinerRefine(cr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
4059: for (n = 0; n < Nct; ++n) {
4060: for (r = 0; r < rsize[n]; ++r) {
4061: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
4062: DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);
4063: localPointsNew[m] = pNew;
4064: remotePointsNew[m].index = pNewRem;
4065: remotePointsNew[m].rank = rankRem;
4066: ++m;
4067: }
4068: }
4069: }
4070: for (n = 0; n < numNeighbors; ++n) {DMPlexCellRefinerDestroy(&crRem[n]);}
4071: PetscFree(crRem);
4072: if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew);
4073: ISRestoreIndices(processRanks, &neighbors);
4074: ISDestroy(&processRanks);
4075: }
4076: {
4077: PetscSFNode *rp, *rtmp;
4078: PetscInt *lp, *idx, *ltmp, i;
4080: /* SF needs sorted leaves to correct calculate Gather */
4081: PetscMalloc1(numLeavesNew, &idx);
4082: PetscMalloc1(numLeavesNew, &lp);
4083: PetscMalloc1(numLeavesNew, &rp);
4084: for (i = 0; i < numLeavesNew; ++i) {
4085: 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);
4086: idx[i] = i;
4087: }
4088: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
4089: for (i = 0; i < numLeavesNew; ++i) {
4090: lp[i] = localPointsNew[idx[i]];
4091: rp[i] = remotePointsNew[idx[i]];
4092: }
4093: ltmp = localPointsNew;
4094: localPointsNew = lp;
4095: rtmp = remotePointsNew;
4096: remotePointsNew = rp;
4097: PetscFree(idx);
4098: PetscFree(ltmp);
4099: PetscFree(rtmp);
4100: }
4101: PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
4102: return(0);
4103: }
4105: static PetscErrorCode RefineLabel_Internal(DMPlexCellRefiner cr, DMLabel label, DMLabel labelNew)
4106: {
4107: DM dm = cr->dm;
4108: IS valueIS;
4109: const PetscInt *values;
4110: PetscInt defVal, Nv, val;
4111: PetscErrorCode ierr;
4114: DMLabelGetDefaultValue(label, &defVal);
4115: DMLabelSetDefaultValue(labelNew, defVal);
4116: DMLabelGetValueIS(label, &valueIS);
4117: ISGetLocalSize(valueIS, &Nv);
4118: ISGetIndices(valueIS, &values);
4119: for (val = 0; val < Nv; ++val) {
4120: IS pointIS;
4121: const PetscInt *points;
4122: PetscInt numPoints, p;
4124: /* Ensure refined label is created with same number of strata as
4125: * original (even if no entries here). */
4126: DMLabelAddStratum(labelNew, values[val]);
4127: DMLabelGetStratumIS(label, values[val], &pointIS);
4128: ISGetLocalSize(pointIS, &numPoints);
4129: ISGetIndices(pointIS, &points);
4130: for (p = 0; p < numPoints; ++p) {
4131: const PetscInt point = points[p];
4132: DMPolytopeType ct;
4133: DMPolytopeType *rct;
4134: PetscInt *rsize, *rcone, *rornt;
4135: PetscInt Nct, n, r, pNew=0;
4137: DMPlexGetCellType(dm, point, &ct);
4138: DMPlexCellRefinerRefine(cr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
4139: for (n = 0; n < Nct; ++n) {
4140: for (r = 0; r < rsize[n]; ++r) {
4141: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);
4142: DMLabelSetValue(labelNew, pNew, values[val]);
4143: }
4144: }
4145: }
4146: ISRestoreIndices(pointIS, &points);
4147: ISDestroy(&pointIS);
4148: }
4149: ISRestoreIndices(valueIS, &values);
4150: ISDestroy(&valueIS);
4151: return(0);
4152: }
4154: static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm)
4155: {
4156: DM dm = cr->dm;
4157: PetscInt numLabels, l;
4161: DMGetNumLabels(dm, &numLabels);
4162: for (l = 0; l < numLabels; ++l) {
4163: DMLabel label, labelNew;
4164: const char *lname;
4165: PetscBool isDepth, isCellType;
4167: DMGetLabelName(dm, l, &lname);
4168: PetscStrcmp(lname, "depth", &isDepth);
4169: if (isDepth) continue;
4170: PetscStrcmp(lname, "celltype", &isCellType);
4171: if (isCellType) continue;
4172: DMCreateLabel(rdm, lname);
4173: DMGetLabel(dm, lname, &label);
4174: DMGetLabel(rdm, lname, &labelNew);
4175: RefineLabel_Internal(cr, label, labelNew);
4176: }
4177: return(0);
4178: }
4180: /* This will only work for interpolated meshes */
4181: PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined)
4182: {
4183: DM rdm;
4184: PetscInt dim, embedDim, depth;
4185: PetscErrorCode ierr;
4189: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
4190: DMSetType(rdm, DMPLEX);
4191: DMGetDimension(dm, &dim);
4192: DMSetDimension(rdm, dim);
4193: DMGetCoordinateDim(dm, &embedDim);
4194: DMSetCoordinateDim(rdm, embedDim);
4195: /* Calculate number of new points of each depth */
4196: DMPlexGetDepth(dm, &depth);
4197: if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement");
4198: /* Step 1: Set chart */
4199: DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);
4200: /* Step 2: Set cone/support sizes (automatically stratifies) */
4201: DMPlexCellRefinerSetConeSizes(cr, rdm);
4202: /* Step 3: Setup refined DM */
4203: DMSetUp(rdm);
4204: /* Step 4: Set cones and supports (automatically symmetrizes) */
4205: DMPlexCellRefinerSetCones(cr, rdm);
4206: /* Step 5: Create pointSF */
4207: DMPlexCellRefinerCreateSF(cr, rdm);
4208: /* Step 6: Create labels */
4209: DMPlexCellRefinerCreateLabels(cr, rdm);
4210: /* Step 7: Set coordinates */
4211: DMPlexCellRefinerSetCoordinates(cr, rdm);
4212: *dmRefined = rdm;
4213: return(0);
4214: }
4216: /*@
4217: DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
4219: Input Parameter:
4220: . dm - The coarse DM
4222: Output Parameter:
4223: . fpointIS - The IS of all the fine points which exist in the original coarse mesh
4225: Level: developer
4227: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
4228: @*/
4229: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
4230: {
4231: DMPlexCellRefiner cr;
4232: PetscInt *fpoints;
4233: PetscInt pStart, pEnd, p, vStart, vEnd, v;
4234: PetscErrorCode ierr;
4237: DMPlexGetChart(dm, &pStart, &pEnd);
4238: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
4239: DMPlexCellRefinerCreate(dm, &cr);
4240: DMPlexCellRefinerSetUp(cr);
4241: PetscMalloc1(pEnd-pStart, &fpoints);
4242: for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
4243: for (v = vStart; v < vEnd; ++v) {
4244: PetscInt vNew = -1; /* silent overzelous may be used uninitialized */
4246: DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
4247: fpoints[v-pStart] = vNew;
4248: }
4249: DMPlexCellRefinerDestroy(&cr);
4250: ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
4251: return(0);
4252: }
4254: /*@
4255: DMPlexSetRefinementUniform - Set the flag for uniform refinement
4257: Input Parameters:
4258: + dm - The DM
4259: - refinementUniform - The flag for uniform refinement
4261: Level: developer
4263: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
4264: @*/
4265: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
4266: {
4267: DM_Plex *mesh = (DM_Plex*) dm->data;
4271: mesh->refinementUniform = refinementUniform;
4272: return(0);
4273: }
4275: /*@
4276: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
4278: Input Parameter:
4279: . dm - The DM
4281: Output Parameter:
4282: . refinementUniform - The flag for uniform refinement
4284: Level: developer
4286: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
4287: @*/
4288: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
4289: {
4290: DM_Plex *mesh = (DM_Plex*) dm->data;
4295: *refinementUniform = mesh->refinementUniform;
4296: return(0);
4297: }
4299: /*@
4300: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
4302: Input Parameters:
4303: + dm - The DM
4304: - refinementLimit - The maximum cell volume in the refined mesh
4306: Level: developer
4308: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
4309: @*/
4310: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
4311: {
4312: DM_Plex *mesh = (DM_Plex*) dm->data;
4316: mesh->refinementLimit = refinementLimit;
4317: return(0);
4318: }
4320: /*@
4321: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
4323: Input Parameter:
4324: . dm - The DM
4326: Output Parameter:
4327: . refinementLimit - The maximum cell volume in the refined mesh
4329: Level: developer
4331: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
4332: @*/
4333: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
4334: {
4335: DM_Plex *mesh = (DM_Plex*) dm->data;
4340: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
4341: *refinementLimit = mesh->refinementLimit;
4342: return(0);
4343: }
4345: /*@
4346: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
4348: Input Parameters:
4349: + dm - The DM
4350: - refinementFunc - Function giving the maximum cell volume in the refined mesh
4352: Note: The calling sequence is refinementFunc(coords, limit)
4353: $ coords - Coordinates of the current point, usually a cell centroid
4354: $ limit - The maximum cell volume for a cell containing this point
4356: Level: developer
4358: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
4359: @*/
4360: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
4361: {
4362: DM_Plex *mesh = (DM_Plex*) dm->data;
4366: mesh->refinementFunc = refinementFunc;
4367: return(0);
4368: }
4370: /*@
4371: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
4373: Input Parameter:
4374: . dm - The DM
4376: Output Parameter:
4377: . refinementFunc - Function giving the maximum cell volume in the refined mesh
4379: Note: The calling sequence is refinementFunc(coords, limit)
4380: $ coords - Coordinates of the current point, usually a cell centroid
4381: $ limit - The maximum cell volume for a cell containing this point
4383: Level: developer
4385: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
4386: @*/
4387: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
4388: {
4389: DM_Plex *mesh = (DM_Plex*) dm->data;
4394: *refinementFunc = mesh->refinementFunc;
4395: return(0);
4396: }
4398: static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm)
4399: {
4400: DM dm = cr->dm;
4401: PetscInt Nf, f, Nds, s;
4405: DMGetNumFields(dm, &Nf);
4406: for (f = 0; f < Nf; ++f) {
4407: DMLabel label, labelNew;
4408: PetscObject obj;
4409: const char *lname;
4411: DMGetField(rdm, f, &label, &obj);
4412: if (!label) continue;
4413: PetscObjectGetName((PetscObject) label, &lname);
4414: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
4415: RefineLabel_Internal(cr, label, labelNew);
4416: DMSetField_Internal(rdm, f, labelNew, obj);
4417: DMLabelDestroy(&labelNew);
4418: }
4419: DMGetNumDS(dm, &Nds);
4420: for (s = 0; s < Nds; ++s) {
4421: DMLabel label, labelNew;
4422: const char *lname;
4424: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
4425: if (!label) continue;
4426: PetscObjectGetName((PetscObject) label, &lname);
4427: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
4428: RefineLabel_Internal(cr, label, labelNew);
4429: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
4430: DMLabelDestroy(&labelNew);
4431: }
4432: return(0);
4433: }
4435: PetscErrorCode DMPlexCellRefinerAdaptLabel(DM dm, DMLabel adaptLabel, DM *dmRefined)
4436: {
4437: DMPlexCellRefiner cr;
4438: DM cdm, rcdm;
4439: PetscErrorCode ierr;
4442: DMPlexCellRefinerCreate(dm, &cr);
4443: cr->adaptLabel = adaptLabel;
4444: DMPlexCellRefinerSetUp(cr);
4445: DMPlexRefineUniform(dm, cr, dmRefined);
4446: DMCopyDisc(dm, *dmRefined);
4447: DMGetCoordinateDM(dm, &cdm);
4448: DMGetCoordinateDM(*dmRefined, &rcdm);
4449: DMCopyDisc(cdm, rcdm);
4450: RefineDiscLabels_Internal(cr, *dmRefined);
4451: DMCopyBoundary(dm, *dmRefined);
4452: DMPlexCellRefinerDestroy(&cr);
4453: return(0);
4454: }
4456: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
4457: {
4458: PetscBool isUniform;
4459: DMPlexCellRefiner cr;
4460: PetscErrorCode ierr;
4463: DMPlexGetRefinementUniform(dm, &isUniform);
4464: DMViewFromOptions(dm, NULL, "-initref_dm_view");
4465: if (isUniform) {
4466: DM cdm, rcdm;
4467: PetscBool localized;
4469: DMPlexCellRefinerCreate(dm, &cr);
4470: DMPlexCellRefinerSetUp(cr);
4471: DMGetCoordinatesLocalized(dm, &localized);
4472: DMPlexRefineUniform(dm, cr, dmRefined);
4473: DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);
4474: DMCopyDisc(dm, *dmRefined);
4475: DMGetCoordinateDM(dm, &cdm);
4476: DMGetCoordinateDM(*dmRefined, &rcdm);
4477: DMCopyDisc(cdm, rcdm);
4478: RefineDiscLabels_Internal(cr, *dmRefined);
4479: DMCopyBoundary(dm, *dmRefined);
4480: DMPlexCellRefinerDestroy(&cr);
4481: } else {
4482: DMPlexRefine_Internal(dm, NULL, dmRefined);
4483: }
4484: if (*dmRefined) {
4485: ((DM_Plex *) (*dmRefined)->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
4486: ((DM_Plex *) (*dmRefined)->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
4487: }
4488: return(0);
4489: }
4491: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
4492: {
4493: DM cdm = dm;
4494: PetscInt r;
4495: PetscBool isUniform, localized;
4499: DMPlexGetRefinementUniform(dm, &isUniform);
4500: DMGetCoordinatesLocalized(dm, &localized);
4501: if (isUniform) {
4502: for (r = 0; r < nlevels; ++r) {
4503: DMPlexCellRefiner cr;
4504: DM codm, rcodm;
4506: DMPlexCellRefinerCreate(cdm, &cr);
4507: DMPlexCellRefinerSetUp(cr);
4508: DMPlexRefineUniform(cdm, cr, &dmRefined[r]);
4509: DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);
4510: DMSetRefineLevel(dmRefined[r], cdm->levelup+1);
4511: DMCopyDisc(cdm, dmRefined[r]);
4512: DMGetCoordinateDM(dm, &codm);
4513: DMGetCoordinateDM(dmRefined[r], &rcodm);
4514: DMCopyDisc(codm, rcodm);
4515: RefineDiscLabels_Internal(cr, dmRefined[r]);
4516: DMCopyBoundary(cdm, dmRefined[r]);
4517: DMSetCoarseDM(dmRefined[r], cdm);
4518: DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
4519: if (dmRefined[r]) {
4520: ((DM_Plex *) (dmRefined[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
4521: ((DM_Plex *) (dmRefined[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
4522: }
4523: cdm = dmRefined[r];
4524: DMPlexCellRefinerDestroy(&cr);
4525: }
4526: } else {
4527: for (r = 0; r < nlevels; ++r) {
4528: DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);
4529: DMCopyDisc(cdm, dmRefined[r]);
4530: DMCopyBoundary(cdm, dmRefined[r]);
4531: if (localized) {DMLocalizeCoordinates(dmRefined[r]);}
4532: DMSetCoarseDM(dmRefined[r], cdm);
4533: if (dmRefined[r]) {
4534: ((DM_Plex *) (dmRefined[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM;
4535: ((DM_Plex *) (dmRefined[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2;
4536: }
4537: cdm = dmRefined[r];
4538: }
4539: }
4540: return(0);
4541: }