Actual source code: plexrefine.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscsf.h>
5: const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "Alfeld2D", "Alfeld3D", "PowellSabin", "BoundaryLayer", "DMPlexCellRefinerTypes", "DM_REFINER_", NULL};
7: /*
8: Note that j and invj are non-square:
9: v0 + j x_face = x_cell
10: invj (x_cell - v0) = x_face
11: */
12: static PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
13: {
14: /*
15: 2
16: |\
17: | \
18: | \
19: | \
20: | \
21: | \
22: | \
23: 2 1
24: | \
25: | \
26: | \
27: 0---0-------1
28: v0[Nf][dc]: 3 x 2
29: J[Nf][df][dc]: 3 x 1 x 2
30: invJ[Nf][dc][df]: 3 x 2 x 1
31: detJ[Nf]: 3
32: */
33: static PetscReal tri_v0[] = {0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
34: static PetscReal tri_J[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
35: static PetscReal tri_invJ[] = {1.0, 0.0, -0.5, 0.5, 0.0, -1.0};
36: static PetscReal tri_detJ[] = {1.0, 1.414213562373095, 1.0};
37: /*
38: 3---------2---------2
39: | |
40: | |
41: | |
42: 3 1
43: | |
44: | |
45: | |
46: 0---------0---------1
48: v0[Nf][dc]: 4 x 2
49: J[Nf][df][dc]: 4 x 1 x 2
50: invJ[Nf][dc][df]: 4 x 2 x 1
51: detJ[Nf]: 4
52: */
53: static PetscReal quad_v0[] = {0.0, -1.0, 1.0, 0.0, 0.0, 1.0 -1.0, 0.0};
54: static PetscReal quad_J[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
55: static PetscReal quad_invJ[] = {1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, -1.0};
56: static PetscReal quad_detJ[] = {1.0, 1.0, 1.0, 1.0};
59: switch (ct) {
60: case DM_POLYTOPE_TRIANGLE: if (Nf) *Nf = 3; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; if (detJ) *detJ = tri_detJ; break;
61: case DM_POLYTOPE_QUADRILATERAL: if (Nf) *Nf = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; if (detJ) *detJ = quad_detJ; break;
62: default:
63: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
64: }
65: return(0);
66: }
68: /* Gets the affine map from the original cell to each subcell */
69: static PetscErrorCode DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
70: {
71: /*
72: 2
73: |\
74: | \
75: | \
76: | \
77: | C \
78: | \
79: | \
80: 2---1---1
81: |\ D / \
82: | 2 0 \
83: |A \ / B \
84: 0---0-------1
85: */
86: static PetscReal tri_v0[] = {-1.0, -1.0, 0.0, -1.0, -1.0, 0.0, 0.0, -1.0};
87: static PetscReal tri_J[] = {0.5, 0.0,
88: 0.0, 0.5,
90: 0.5, 0.0,
91: 0.0, 0.5,
93: 0.5, 0.0,
94: 0.0, 0.5,
96: 0.0, -0.5,
97: 0.5, 0.5};
98: static PetscReal tri_invJ[] = {2.0, 0.0,
99: 0.0, 2.0,
101: 2.0, 0.0,
102: 0.0, 2.0,
104: 2.0, 0.0,
105: 0.0, 2.0,
107: 2.0, 2.0,
108: -2.0, 0.0};
109: /*
110: 3---------2---------2
111: | | |
112: | D 2 C |
113: | | |
114: 3----3----0----1----1
115: | | |
116: | A 0 B |
117: | | |
118: 0---------0---------1
119: */
120: static PetscReal quad_v0[] = {-1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
121: static PetscReal quad_J[] = {0.5, 0.0,
122: 0.0, 0.5,
124: 0.5, 0.0,
125: 0.0, 0.5,
127: 0.5, 0.0,
128: 0.0, 0.5,
130: 0.5, 0.0,
131: 0.0, 0.5};
132: static PetscReal quad_invJ[] = {2.0, 0.0,
133: 0.0, 2.0,
135: 2.0, 0.0,
136: 0.0, 2.0,
138: 2.0, 0.0,
139: 0.0, 2.0,
141: 2.0, 0.0,
142: 0.0, 2.0};
143: /*
144: Bottom (viewed from top) Top
145: 1---------2---------2 7---------2---------6
146: | | | | | |
147: | B 2 C | | H 2 G |
148: | | | | | |
149: 3----3----0----1----1 3----3----0----1----1
150: | | | | | |
151: | A 0 D | | E 0 F |
152: | | | | | |
153: 0---------0---------3 4---------0---------5
154: */
155: static PetscReal hex_v0[] = {-1.0, -1.0, -1.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, -1.0,
156: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
157: static PetscReal hex_J[] = {0.5, 0.0, 0.0,
158: 0.0, 0.5, 0.0,
159: 0.0, 0.0, 0.5,
161: 0.5, 0.0, 0.0,
162: 0.0, 0.5, 0.0,
163: 0.0, 0.0, 0.5,
165: 0.5, 0.0, 0.0,
166: 0.0, 0.5, 0.0,
167: 0.0, 0.0, 0.5,
169: 0.5, 0.0, 0.0,
170: 0.0, 0.5, 0.0,
171: 0.0, 0.0, 0.5,
173: 0.5, 0.0, 0.0,
174: 0.0, 0.5, 0.0,
175: 0.0, 0.0, 0.5,
177: 0.5, 0.0, 0.0,
178: 0.0, 0.5, 0.0,
179: 0.0, 0.0, 0.5,
181: 0.5, 0.0, 0.0,
182: 0.0, 0.5, 0.0,
183: 0.0, 0.0, 0.5,
185: 0.5, 0.0, 0.0,
186: 0.0, 0.5, 0.0,
187: 0.0, 0.0, 0.5};
188: static PetscReal hex_invJ[] = {2.0, 0.0, 0.0,
189: 0.0, 2.0, 0.0,
190: 0.0, 0.0, 2.0,
192: 2.0, 0.0, 0.0,
193: 0.0, 2.0, 0.0,
194: 0.0, 0.0, 2.0,
196: 2.0, 0.0, 0.0,
197: 0.0, 2.0, 0.0,
198: 0.0, 0.0, 2.0,
200: 2.0, 0.0, 0.0,
201: 0.0, 2.0, 0.0,
202: 0.0, 0.0, 2.0,
204: 2.0, 0.0, 0.0,
205: 0.0, 2.0, 0.0,
206: 0.0, 0.0, 2.0,
208: 2.0, 0.0, 0.0,
209: 0.0, 2.0, 0.0,
210: 0.0, 0.0, 2.0,
212: 2.0, 0.0, 0.0,
213: 0.0, 2.0, 0.0,
214: 0.0, 0.0, 2.0,
216: 2.0, 0.0, 0.0,
217: 0.0, 2.0, 0.0,
218: 0.0, 0.0, 2.0};
221: switch (ct) {
222: case DM_POLYTOPE_TRIANGLE: if (Nc) *Nc = 4; if (v0) *v0 = tri_v0; if (J) *J = tri_J; if (invJ) *invJ = tri_invJ; break;
223: case DM_POLYTOPE_QUADRILATERAL: if (Nc) *Nc = 4; if (v0) *v0 = quad_v0; if (J) *J = quad_J; if (invJ) *invJ = quad_invJ; break;
224: case DM_POLYTOPE_HEXAHEDRON: if (Nc) *Nc = 8; if (v0) *v0 = hex_v0; if (J) *J = hex_J; if (invJ) *invJ = hex_invJ; break;
225: default:
226: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
227: }
228: return(0);
229: }
231: /* Should this be here or in the DualSpace somehow? */
232: PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside)
233: {
234: PetscReal sum = 0.0;
235: PetscInt d;
238: *inside = PETSC_TRUE;
239: switch (ct) {
240: case DM_POLYTOPE_TRIANGLE:
241: case DM_POLYTOPE_TETRAHEDRON:
242: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d) {
243: if (point[d] < -1.0) {*inside = PETSC_FALSE; break;}
244: sum += point[d];
245: }
246: if (sum > PETSC_SMALL) {*inside = PETSC_FALSE; break;}
247: break;
248: case DM_POLYTOPE_QUADRILATERAL:
249: case DM_POLYTOPE_HEXAHEDRON:
250: for (d = 0; d < DMPolytopeTypeGetDim(ct); ++d)
251: if (PetscAbsReal(point[d]) > 1.+PETSC_SMALL) {*inside = PETSC_FALSE; break;}
252: break;
253: default:
254: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported polytope type %s", DMPolytopeTypes[ct]);
255: }
256: return(0);
257: }
259: /* Regular Refinment of Hybrid Meshes
261: We would like to express regular refinement as a small set of rules that can be applied on every point of the Plex
262: to automatically generate a refined Plex. In fact, we would like these rules to be general enough to encompass other
263: transformations, such as changing from one type of cell to another, as simplex to hex.
265: To start, we can create a function that takes an original cell type and returns the number of new cells replacing it
266: and the types of the new cells.
268: We need the group multiplication table for group actions from the dihedral group for each cell type.
270: We need an operator which takes in a cell, and produces a new set of cells with new faces and correct orientations. I think
271: we can just write this operator for faces with identity, and then compose the face orientation actions to get the actual
272: (face, orient) pairs for each subcell.
273: */
275: /*
276: Input Parameters:
277: + ct - The type of the input cell
278: . co - The orientation of this cell in its parent
279: - cp - The requested cone point of this cell assuming orientation 0
281: Output Parameters:
282: . cpnew - The new cone point, taking into account the orientation co
283: */
284: PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew)
285: {
286: const PetscInt csize = DMPolytopeTypeGetConeSize(ct);
289: if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;}
290: else {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;}
291: return(0);
292: }
294: static PetscErrorCode DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
295: {
296: static PetscReal seg_v[] = {-1.0, 0.0, 1.0};
297: static PetscReal tri_v[] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0};
298: static PetscReal quad_v[] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0};
299: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
300: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
301: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0};
302: static PetscReal hex_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
303: -1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, -1.0,
304: -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0,
305: -1.0, -1.0, 0.0, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0,
306: -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
307: -1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0,
308: -1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, 1.0,
309: -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
310: -1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0};
313: switch (ct) {
314: case DM_POLYTOPE_SEGMENT: *Nv = 3; *subcellV = seg_v; break;
315: case DM_POLYTOPE_TRIANGLE: *Nv = 6; *subcellV = tri_v; break;
316: case DM_POLYTOPE_QUADRILATERAL: *Nv = 9; *subcellV = quad_v; break;
317: case DM_POLYTOPE_TETRAHEDRON: *Nv = 10; *subcellV = tet_v; break;
318: case DM_POLYTOPE_HEXAHEDRON: *Nv = 27; *subcellV = hex_v; break;
319: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
320: }
321: return(0);
322: }
324: static PetscErrorCode DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
325: {
326: static PetscReal tri_v[] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, -1.0/3.0, -1.0/3.0};
327: static PetscReal tet_v[] = {-1.0, -1.0, -1.0, 0.0, -1.0, -1.0, 1.0, -1.0, -1.0,
328: -1.0, 0.0, -1.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0, -1.0, -1.0, 1.0, -1.0,
329: -1.0, -1.0, 0.0, -1.0/3.0, -1.0, -1.0/3.0, 0.0, -1.0, 0.0,
330: -1.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0/3.0, -1.0, 0.0, 0.0,
331: -1.0, -1.0, 1.0, -0.5, -0.5, -0.5};
332: PetscErrorCode ierr;
335: switch (ct) {
336: case DM_POLYTOPE_SEGMENT:
337: case DM_POLYTOPE_QUADRILATERAL:
338: case DM_POLYTOPE_HEXAHEDRON:
339: DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);break;
340: case DM_POLYTOPE_TRIANGLE: *Nv = 7; *subcellV = tri_v; break;
341: case DM_POLYTOPE_TETRAHEDRON: *Nv = 15; *subcellV = tet_v; break;
342: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
343: }
344: return(0);
345: }
347: /*
348: DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type
350: Input Parameters:
351: + cr - The DMPlexCellRefiner object
352: - ct - The cell type
354: Output Parameters:
355: + Nv - The number of refined vertices in the closure of the reference cell of given type
356: - subcellV - The coordinates of these vertices in the reference cell
358: Level: developer
360: .seealso: DMPlexCellRefinerGetSubcellVertices()
361: */
362: static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
363: {
367: if (!cr->ops->getcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
368: (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);
369: return(0);
370: }
372: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
373: {
374: static PetscInt seg_v[] = {0, 1, 1, 2};
375: static PetscInt tri_v[] = {0, 3, 5, 3, 1, 4, 5, 4, 2, 3, 4, 5};
376: static PetscInt quad_v[] = {0, 4, 8, 7, 4, 1, 5, 8, 8, 5, 2, 6, 7, 8, 6, 3};
377: static PetscInt tet_v[] = {0, 3, 1, 6, 3, 2, 4, 8, 1, 4, 5, 7, 6, 8, 7, 9,
378: 1, 6, 3, 7, 8, 4, 3, 7, 7, 3, 1, 4, 7, 3, 8, 6};
379: 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,
380: 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};
383: if (ct != rct) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
384: switch (ct) {
385: case DM_POLYTOPE_SEGMENT: *Nv = 2; *subcellV = &seg_v[r*(*Nv)]; break;
386: case DM_POLYTOPE_TRIANGLE: *Nv = 3; *subcellV = &tri_v[r*(*Nv)]; break;
387: case DM_POLYTOPE_QUADRILATERAL: *Nv = 4; *subcellV = &quad_v[r*(*Nv)]; break;
388: case DM_POLYTOPE_TETRAHEDRON: *Nv = 4; *subcellV = &tet_v[r*(*Nv)]; break;
389: case DM_POLYTOPE_HEXAHEDRON: *Nv = 8; *subcellV = &hex_v[r*(*Nv)]; break;
390: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
391: }
392: return(0);
393: }
395: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
396: {
397: static PetscInt tri_v[] = {0, 3, 6, 5, 3, 1, 4, 6, 5, 6, 4, 2};
398: 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};
399: PetscErrorCode ierr;
402: switch (ct) {
403: case DM_POLYTOPE_SEGMENT:
404: case DM_POLYTOPE_QUADRILATERAL:
405: case DM_POLYTOPE_HEXAHEDRON:
406: DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);break;
407: case DM_POLYTOPE_TRIANGLE:
408: if (rct != DM_POLYTOPE_QUADRILATERAL) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
409: *Nv = 4; *subcellV = &tri_v[r*(*Nv)]; break;
410: case DM_POLYTOPE_TETRAHEDRON:
411: if (rct != DM_POLYTOPE_HEXAHEDRON) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
412: *Nv = 8; *subcellV = &tet_v[r*(*Nv)]; break;
413: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No subcell vertices for cell type %s", DMPolytopeTypes[ct]);
414: }
415: return(0);
416: }
418: /*
419: DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type
421: Input Parameters:
422: + cr - The DMPlexCellRefiner object
423: . ct - The cell type
424: . rct - The type of subcell
425: - r - The subcell index
427: Output Parameters:
428: + Nv - The number of refined vertices in the subcell
429: - subcellV - The indices of these vertices in the set of vertices returned by DMPlexCellRefinerGetCellVertices()
431: Level: developer
433: .seealso: DMPlexCellRefinerGetCellVertices()
434: */
435: static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
436: {
440: if (!cr->ops->getsubcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
441: (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);
442: return(0);
443: }
445: /*
446: Input Parameters:
447: + cr - The DMPlexCellRefiner
448: . pct - The cell type of the parent, from whom the new cell is being produced
449: . ct - The type being produced
450: . r - The replica number requested for the produced cell type
451: . Nv - Number of vertices in the closure of the parent cell
452: . dE - Spatial dimension
453: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
455: Output Parameters:
456: . out - The coordinates of the new vertices
457: */
458: static PetscErrorCode DMPlexCellRefinerMapCoordinates(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
459: {
463: if (!cr->ops->mapcoords) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
464: (*cr->ops->mapcoords)(cr, pct, ct, r, Nv, dE, in, out);
465: return(0);
466: }
468: /* Computes new vertex as the barycenter */
469: static PetscErrorCode DMPlexCellRefinerMapCoordinates_Barycenter(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
470: {
471: PetscInt v,d;
474: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
475: for (d = 0; d < dE; ++d) out[d] = 0.0;
476: for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
477: for (d = 0; d < dE; ++d) out[d] /= Nv;
478: return(0);
479: }
481: /*
482: Input Parameters:
483: + cr - The DMPlexCellRefiner
484: . pct - The cell type of the parent, from whom the new cell is being produced
485: . po - The orientation of the parent cell in its enclosing parent
486: . ct - The type being produced
487: . r - The replica number requested for the produced cell type
488: - o - The relative orientation of the replica
490: Output Parameters:
491: + rnew - The replica number, given the orientation of the parent
492: - onew - The replica orientation, given the orientation of the parent
493: */
494: static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
495: {
499: if (!cr->ops->mapsubcells) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
500: (*cr->ops->mapsubcells)(cr, pct, po, ct, r, o, rnew, onew);
501: return(0);
502: }
504: /*
505: This is the group multiplication table for the dihedral group of the cell.
506: */
507: static PetscErrorCode ComposeOrientation_Private(PetscInt n, PetscInt o1, PetscInt o2, PetscInt *o)
508: {
510: if (!n) {*o = 0;}
511: else if (o1 >= 0 && o2 >= 0) {*o = ( o1 + o2) % n;}
512: else if (o1 < 0 && o2 < 0) {*o = (-o1 - o2) % n;}
513: else if (o1 < 0) {*o = -((-(o1+1) + o2) % n + 1);}
514: else if (o2 < 0) {*o = -((-(o2+1) + o1) % n + 1);}
515: return(0);
516: }
518: static PetscErrorCode DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
519: {
523: *rnew = r;
524: ComposeOrientation_Private(DMPolytopeTypeGetConeSize(ct), po, o, onew);
525: return(0);
526: }
528: static PetscErrorCode DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
529: {
530: /* We shift any input orientation in order to make it non-negative
531: 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)
532: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
533: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
534: */
535: PetscInt tri_seg_o[] = {-2, 0,
536: -2, 0,
537: -2, 0,
538: 0, -2,
539: 0, -2,
540: 0, -2};
541: PetscInt tri_seg_r[] = {1, 0, 2,
542: 0, 2, 1,
543: 2, 1, 0,
544: 0, 1, 2,
545: 1, 2, 0,
546: 2, 0, 1};
547: PetscInt tri_tri_o[] = {0, 1, 2, -3, -2, -1,
548: 2, 0, 1, -2, -1, -3,
549: 1, 2, 0, -1, -3, -2,
550: -3, -2, -1, 0, 1, 2,
551: -1, -3, -2, 1, 2, 0,
552: -2, -1, -3, 2, 0, 1};
553: /* orientation if the replica is the center triangle */
554: PetscInt tri_tri_o_c[] = {2, 0, 1, -2, -1, -3,
555: 1, 2, 0, -1, -3, -2,
556: 0, 1, 2, -3, -2, -1,
557: -3, -2, -1, 0, 1, 2,
558: -1, -3, -2, 1, 2, 0,
559: -2, -1, -3, 2, 0, 1};
560: PetscInt tri_tri_r[] = {0, 2, 1, 3,
561: 2, 1, 0, 3,
562: 1, 0, 2, 3,
563: 0, 1, 2, 3,
564: 1, 2, 0, 3,
565: 2, 0, 1, 3};
566: PetscInt quad_seg_r[] = {3, 2, 1, 0,
567: 2, 1, 0, 3,
568: 1, 0, 3, 2,
569: 0, 3, 2, 1,
570: 0, 1, 2, 3,
571: 1, 2, 3, 0,
572: 2, 3, 0, 1,
573: 3, 0, 1, 2};
574: PetscInt quad_quad_o[] = { 0, 1, 2, 3, -4, -3, -2, -1,
575: 4, 0, 1, 2, -3, -2, -1, -4,
576: 3, 4, 0, 1, -2, -1, -4, -3,
577: 2, 3, 4, 0, -1, -4, -3, -2,
578: -4, -3, -2, -1, 0, 1, 2, 3,
579: -1, -4, -3, -2, 1, 2, 3, 0,
580: -2, -1, -4, -3, 2, 3, 0, 1,
581: -3, -2, -1, -4, 3, 0, 1, 2};
582: PetscInt quad_quad_r[] = {0, 3, 2, 1,
583: 3, 2, 1, 0,
584: 2, 1, 0, 3,
585: 1, 0, 3, 2,
586: 0, 1, 2, 3,
587: 1, 2, 3, 0,
588: 2, 3, 0, 1,
589: 3, 0, 1, 2};
590: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
591: 1, 0, -1, -2,
592: -2, -1, 0, 1,
593: -1, -2, 1, 0};
594: PetscInt tquad_tquad_r[] = {1, 0,
595: 1, 0,
596: 0, 1,
597: 0, 1};
600: /* The default is no transformation */
601: *rnew = r;
602: *onew = o;
603: switch (pct) {
604: case DM_POLYTOPE_SEGMENT:
605: if (ct == DM_POLYTOPE_SEGMENT) {
606: if (po == 0 || po == -1) {*rnew = r; *onew = o;}
607: else if (po == 1 || po == -2) {*rnew = (r+1)%2; *onew = (o == 0 || o == -1) ? -2 : 0;}
608: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for segment", po);
609: }
610: break;
611: case DM_POLYTOPE_TRIANGLE:
612: switch (ct) {
613: case DM_POLYTOPE_SEGMENT:
614: if (o == -1) o = 0;
615: if (o == -2) o = 1;
616: *onew = tri_seg_o[(po+3)*2+o];
617: *rnew = tri_seg_r[(po+3)*3+r];
618: break;
619: case DM_POLYTOPE_TRIANGLE:
620: *onew = r == 3 && po < 0 ? tri_tri_o_c[((po+3)%3)*6+o+3] : tri_tri_o[(po+3)*6+o+3];
621: *rnew = tri_tri_r[(po+3)*4+r];
622: break;
623: default: break;
624: }
625: break;
626: case DM_POLYTOPE_QUADRILATERAL:
627: switch (ct) {
628: case DM_POLYTOPE_SEGMENT:
629: *onew = o;
630: *rnew = quad_seg_r[(po+4)*4+r];
631: break;
632: case DM_POLYTOPE_QUADRILATERAL:
633: *onew = quad_quad_o[(po+4)*8+o+4];
634: *rnew = quad_quad_r[(po+4)*4+r];
635: break;
636: default: break;
637: }
638: break;
639: case DM_POLYTOPE_SEG_PRISM_TENSOR:
640: switch (ct) {
641: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
642: case DM_POLYTOPE_SEG_PRISM_TENSOR:
643: *onew = tquad_tquad_o[(po+2)*4+o+2];
644: *rnew = tquad_tquad_r[(po+2)*2+r];
645: break;
646: default: break;
647: }
648: break;
649: default: break;
650: }
651: return(0);
652: }
654: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
655: {
657: /* We shift any input orientation in order to make it non-negative
658: 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)
659: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
660: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
661: */
662: PetscInt tri_seg_o[] = {0, -2,
663: 0, -2,
664: 0, -2,
665: 0, -2,
666: 0, -2,
667: 0, -2};
668: PetscInt tri_seg_r[] = {2, 1, 0,
669: 1, 0, 2,
670: 0, 2, 1,
671: 0, 1, 2,
672: 1, 2, 0,
673: 2, 0, 1};
674: PetscInt tri_tri_o[] = {0, 1, 2, 3, -4, -3, -2, -1,
675: 3, 0, 1, 2, -3, -2, -1, -4,
676: 1, 2, 3, 0, -1, -4, -3, -2,
677: -4, -3, -2, -1, 0, 1, 2, 3,
678: -1, -4, -3, -2, 1, 2, 3, 0,
679: -3, -2, -1, -4, 3, 0, 1, 2};
680: PetscInt tri_tri_r[] = {0, 2, 1,
681: 2, 1, 0,
682: 1, 0, 2,
683: 0, 1, 2,
684: 1, 2, 0,
685: 2, 0, 1};
686: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
687: 1, 0, -1, -2,
688: -2, -1, 0, 1,
689: -1, -2, 1, 0};
690: PetscInt tquad_tquad_r[] = {1, 0,
691: 1, 0,
692: 0, 1,
693: 0, 1};
694: PetscInt tquad_quad_o[] = {-2, -1, -4, -3, 2, 3, 0, 1,
695: 1, 2, 3, 0, -1, -4, -3, -2,
696: -4, -3, -2, -1, 0, 1, 2, 3,
697: 1, 0, 3, 2, -3, -4, -1, -2};
700: *rnew = r;
701: *onew = o;
702: switch (pct) {
703: case DM_POLYTOPE_TRIANGLE:
704: switch (ct) {
705: case DM_POLYTOPE_SEGMENT:
706: if (o == -1) o = 0;
707: if (o == -2) o = 1;
708: *onew = tri_seg_o[(po+3)*2+o];
709: *rnew = tri_seg_r[(po+3)*3+r];
710: break;
711: case DM_POLYTOPE_QUADRILATERAL:
712: *onew = tri_tri_o[(po+3)*8+o+4];
713: /* TODO See sheet, for fixing po == 1 and 2 */
714: if (po == 2 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
715: if (po == 2 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
716: if (po == 1 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
717: if (po == 1 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
718: if (po == -1 && r == 2 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+3)%4)+4];
719: if (po == -1 && r == 2 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+5)%4)];
720: if (po == -2 && r == 1 && o >= 0) *onew = tri_tri_o[(po+3)*8+((o+1)%4)+4];
721: if (po == -2 && r == 1 && o < 0) *onew = tri_tri_o[(po+3)*8+((o+7)%4)];
722: *rnew = tri_tri_r[(po+3)*3+r];
723: break;
724: default: break;
725: }
726: break;
727: case DM_POLYTOPE_SEG_PRISM_TENSOR:
728: switch (ct) {
729: /* DM_POLYTOPE_POINT_PRISM_TENSOR does not change */
730: case DM_POLYTOPE_SEG_PRISM_TENSOR:
731: *onew = tquad_tquad_o[(po+2)*4+o+2];
732: *rnew = tquad_tquad_r[(po+2)*2+r];
733: break;
734: case DM_POLYTOPE_QUADRILATERAL:
735: *onew = tquad_quad_o[(po+2)*8+o+4];
736: *rnew = tquad_tquad_r[(po+2)*2+r];
737: break;
738: default: break;
739: }
740: break;
741: default:
742: DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
743: }
744: return(0);
745: }
747: static PetscErrorCode DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
748: {
749: return DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);
750: }
752: /*@
753: DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type
755: Input Parameter:
756: . source - The cell type for a source point
758: Output Parameter:
759: + Nt - The number of cell types generated by refinement
760: . target - The cell types generated
761: . size - The number of subcells of each type, ordered by dimension
762: . cone - A list of the faces for each subcell of the same type as source
763: - ornt - A list of the face orientations for each subcell of the same type as source
765: Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
766: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
767: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
768: $ the number of cones to be taken, or 0 for the current cell
769: $ the cell cone point number at each level from which it is subdivided
770: $ the replica number r of the subdivision.
771: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
772: $ Nt = 2
773: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
774: $ size = {1, 2}
775: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
776: $ ornt = { 0, 0, 0, 0}
778: Level: developer
780: .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform()
781: @*/
782: PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
783: {
787: if (!cr->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
788: (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);
789: return(0);
790: }
792: static PetscErrorCode DMPlexCellRefinerRefine_None(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
793: {
794: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
795: static PetscInt vertexS[] = {1};
796: static PetscInt vertexC[] = {0};
797: static PetscInt vertexO[] = {0};
798: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
799: static PetscInt edgeS[] = {1};
800: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
801: static PetscInt edgeO[] = {0, 0};
802: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
803: static PetscInt tedgeS[] = {1};
804: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
805: static PetscInt tedgeO[] = {0, 0};
806: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
807: static PetscInt triS[] = {1};
808: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
809: static PetscInt triO[] = {0, 0, 0};
810: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
811: static PetscInt quadS[] = {1};
812: 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};
813: static PetscInt quadO[] = {0, 0, 0, 0};
814: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
815: static PetscInt tquadS[] = {1};
816: 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};
817: static PetscInt tquadO[] = {0, 0, 0, 0};
818: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
819: static PetscInt tetS[] = {1};
820: 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};
821: static PetscInt tetO[] = {0, 0, 0, 0};
822: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
823: static PetscInt hexS[] = {1};
824: static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
825: DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
826: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
827: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
828: static PetscInt tripS[] = {1};
829: static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
830: DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
831: static PetscInt tripO[] = {0, 0, 0, 0, 0};
832: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
833: static PetscInt ttripS[] = {1};
834: static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
835: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
836: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
837: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
838: static PetscInt tquadpS[] = {1};
839: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
840: 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};
841: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
844: switch (source) {
845: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
846: case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
847: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
848: case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
849: case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
850: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
851: case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
852: case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
853: case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
854: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
855: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
856: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
857: }
858: return(0);
859: }
861: static PetscErrorCode DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
862: {
863: /* All vertices remain in the refined mesh */
864: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
865: static PetscInt vertexS[] = {1};
866: static PetscInt vertexC[] = {0};
867: static PetscInt vertexO[] = {0};
868: /* Split all edges with a new vertex, making two new 2 edges
869: 0--0--0--1--1
870: */
871: static DMPolytopeType edgeT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT};
872: static PetscInt edgeS[] = {1, 2};
873: 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};
874: static PetscInt edgeO[] = { 0, 0, 0, 0};
875: /* Do not split tensor edges */
876: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
877: static PetscInt tedgeS[] = {1};
878: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
879: static PetscInt tedgeO[] = { 0, 0};
880: /* Add 3 edges inside every triangle, making 4 new triangles.
881: 2
882: |\
883: | \
884: | \
885: 0 1
886: | C \
887: | \
888: | \
889: 2---1---1
890: |\ D / \
891: 1 2 0 0
892: |A \ / B \
893: 0-0-0---1---1
894: */
895: static DMPolytopeType triT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
896: static PetscInt triS[] = {3, 4};
897: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
898: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
899: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 0, 0,
900: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
901: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
902: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
903: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2};
904: static PetscInt triO[] = {0, 0,
905: 0, 0,
906: 0, 0,
907: 0, -2, 0,
908: 0, 0, -2,
909: -2, 0, 0,
910: 0, 0, 0};
911: /* Add a vertex in the center of each quadrilateral, and 4 edges inside, making 4 new quads.
912: 3----1----2----0----2
913: | | |
914: 0 D 2 C 1
915: | | |
916: 3----3----0----1----1
917: | | |
918: 1 A 0 B 0
919: | | |
920: 0----0----0----1----1
921: */
922: static DMPolytopeType quadT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
923: static PetscInt quadS[] = {1, 4, 4};
924: static PetscInt quadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
925: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
926: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
927: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
928: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1,
929: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
930: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2,
931: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0};
932: static PetscInt quadO[] = {0, 0,
933: 0, 0,
934: 0, 0,
935: 0, 0,
936: 0, 0, -2, 0,
937: 0, 0, 0, -2,
938: -2, 0, 0, 0,
939: 0, -2, 0, 0};
940: /* Add 1 edge inside every tensor quad, making 2 new tensor quads
941: 2----2----1----3----3
942: | | |
943: | | |
944: | | |
945: 4 A 6 B 5
946: | | |
947: | | |
948: | | |
949: 0----0----0----1----1
950: */
951: static DMPolytopeType tquadT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR};
952: static PetscInt tquadS[] = {1, 2};
953: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
954: 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,
955: 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};
956: static PetscInt tquadO[] = {0, 0,
957: 0, 0, 0, 0,
958: 0, 0, 0, 0};
959: /* Add 1 edge and 8 triangles inside every cell, making 8 new tets
960: 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
961: 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]
962: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
963: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
964: The first four tets just cut off the corners, using the replica notation for new vertices,
965: [v0, (e0, 0), (e2, 0), (e3, 0)]
966: [(e0, 0), v1, (e1, 0), (e4, 0)]
967: [(e2, 0), (e1, 0), v2, (e5, 0)]
968: [(e3, 0), (e4, 0), (e5, 0), v3 ]
969: The next four tets match a vertex to the newly created faces from cutting off those first tets.
970: [(e2, 0), (e3, 0), (e0, 0), (e5, 0)]
971: [(e4, 0), (e1, 0), (e0, 0), (e5, 0)]
972: [(e5, 0), (e0, 0), (e2, 0), (e1, 0)]
973: [(e5, 0), (e0, 0), (e4, 0), (e3, 0)]
974: 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
975: [(e2, 0), (e0, 0), (e3, 0)]
976: [(e0, 0), (e1, 0), (e4, 0)]
977: [(e2, 0), (e5, 0), (e1, 0)]
978: [(e3, 0), (e4, 0), (e5, 0)]
979: The next four, from the second group of tets, are
980: [(e2, 0), (e0, 0), (e5, 0)]
981: [(e4, 0), (e0, 0), (e5, 0)]
982: [(e0, 0), (e1, 0), (e5, 0)]
983: [(e5, 0), (e3, 0), (e0, 0)]
984: I could write a program to generate these orientations by comparing the faces from GetRawFaces() with my existing table.
985: */
986: static DMPolytopeType tetT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
987: static PetscInt tetS[] = {1, 8, 8};
988: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 2, 2, 1, 0,
989: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2,
990: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1,
991: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 0, 1,
992: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1,
993: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
994: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 1,
995: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 0,
996: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 0,
997: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 0,
998: DM_POLYTOPE_TRIANGLE, 1, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 1, 2, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 1,
999: DM_POLYTOPE_TRIANGLE, 1, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 2, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 0,
1000: DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 1, DM_POLYTOPE_TRIANGLE, 1, 2, 2, DM_POLYTOPE_TRIANGLE, 1, 3, 2,
1001: DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 3, DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 7,
1002: DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 1, 3, 3, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 6,
1003: DM_POLYTOPE_TRIANGLE, 0, 4, DM_POLYTOPE_TRIANGLE, 0, 6, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 1, 0, 3,
1004: DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 7, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 1, 1, 3};
1005: static PetscInt tetO[] = {0, 0,
1006: 0, 0, 0,
1007: 0, 0, 0,
1008: 0, 0, 0,
1009: 0, 0, 0,
1010: 0, 0, -2,
1011: 0, 0, -2,
1012: 0, -2, -2,
1013: 0, -2, 0,
1014: 0, 0, 0, 0,
1015: 0, 0, 0, 0,
1016: 0, 0, 0, 0,
1017: 0, 0, 0, 0,
1018: -3, 0, 0, -2,
1019: -2, 1, 0, 0,
1020: -2, -2, -1, 2,
1021: -2, 0, -2, 1};
1022: /* Add a vertex in the center of each cell, add 6 edges and 12 quads inside every cell, making 8 new hexes
1023: 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
1024: [v0, v1, v2, v3] f0 bottom
1025: [v4, v5, v6, v7] f1 top
1026: [v0, v3, v5, v4] f2 front
1027: [v2, v1, v7, v6] f3 back
1028: [v3, v2, v6, v5] f4 right
1029: [v0, v4, v7, v1] f5 left
1030: The eight hexes are divided into four on the bottom, and four on the top,
1031: [v0, (e0, 0), (f0, 0), (e3, 0), (e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1032: [(e0, 0), v1, (e1, 0), (f0, 0), (f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1033: [(f0, 0), (e1, 0), v2, (e2, 0), (c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1034: [(e3, 0), (f0, 0), (e2, 0), v3, (f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1035: [(e9, 0), (f5, 0), (c0, 0), (f2, 0), v4, (e4, 0), (f1, 0), (e7, 0)]
1036: [(f2, 0), (c0, 0), (f4, 0), (e8, 0), (e4, 0), v5, (e5, 0), (f1, 0)]
1037: [(c0, 0), (f3, 0), (e11, 0), (f4, 0), (f1, 0), (e5, 0), v6, (e6, 0)]
1038: [(f5, 0), (e10, 0), (f3, 0), (c0, 0), (e7, 0), (f1, 0), (e6, 0), v7]
1039: 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,
1040: [(e9, 0), (f2, 0), (c0, 0), (f5, 0)]
1041: [(f5, 0), (c0, 0), (f3, 0), (e10, 0)]
1042: [(c0, 0), (f4, 0), (e11, 0), (f3, 0)]
1043: [(f2, 0), (e8, 0), (f4, 0), (c0, 0)]
1044: and on the x-z plane,
1045: [(f0, 0), (e0, 0), (f5, 0), (c0, 0)]
1046: [(c0, 0), (f5, 0), (e7, 0), (f1, 0)]
1047: [(f4, 0), (c0, 0), (f1, 0), (e5, 0)]
1048: [(e2, 0), (f0, 0), (c0, 0), (f4, 0)]
1049: and on the y-z plane,
1050: [(e3, 0), (f2, 0), (c0, 0), (f0, 0)]
1051: [(f2, 0), (e4, 0), (f1, 0), (c0, 0)]
1052: [(c0, 0), (f1, 0), (e6, 0), (f3, 0)]
1053: [(f0, 0), (c0, 0), (f3, 0), (e1, 0)]
1054: */
1055: static DMPolytopeType hexT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1056: static PetscInt hexS[] = {1, 6, 12, 8};
1057: static PetscInt hexC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1058: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1059: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1060: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1061: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1062: DM_POLYTOPE_POINT, 1, 5, 0, DM_POLYTOPE_POINT, 0, 0,
1063: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 0,
1064: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 5, 2,
1065: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3,
1066: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 2,
1067: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 5, 3, DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 0, 0,
1068: DM_POLYTOPE_SEGMENT, 0, 5, DM_POLYTOPE_SEGMENT, 1, 5, 1, DM_POLYTOPE_SEGMENT, 1, 1, 3, DM_POLYTOPE_SEGMENT, 0, 1,
1069: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1070: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1071: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 3,
1072: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1073: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1074: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1075: 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,
1076: 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,
1077: 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,
1078: 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,
1079: 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,
1080: 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,
1081: 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,
1082: 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};
1083: static PetscInt hexO[] = {0, 0,
1084: 0, 0,
1085: 0, 0,
1086: 0, 0,
1087: 0, 0,
1088: 0, 0,
1089: 0, 0, -2, -2,
1090: 0, -2, -2, 0,
1091: -2, -2, 0, 0,
1092: -2, 0, 0, -2,
1093: -2, 0, 0, -2,
1094: -2, -2, 0, 0,
1095: 0, -2, -2, 0,
1096: 0, 0, -2, -2,
1097: 0, 0, -2, -2,
1098: -2, 0, 0, -2,
1099: -2, -2, 0, 0,
1100: 0, -2, -2, 0,
1101: 0, 0, 0, 0, -4, 0,
1102: 0, 0, -1, 0, -4, 0,
1103: 0, 0, -1, 0, 0, 0,
1104: 0, 0, 0, 0, 0, 0,
1105: -4, 0, 0, 0, -4, 0,
1106: -4, 0, 0, 0, 0, 0,
1107: -4, 0, -1, 0, 0, 0,
1108: -4, 0, -1, 0, -4, 0};
1109: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1110: static DMPolytopeType tripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM};
1111: static PetscInt tripS[] = {3, 4, 6, 8};
1112: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 1, 3, 0,
1113: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 1, 4, 0,
1114: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 1, 2, 0,
1115: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1116: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 0,
1117: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1118: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1119: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1120: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 4, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1121: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1122: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1123: DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 4, 2, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1124: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1125: 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,
1126: 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,
1127: 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,
1128: 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,
1129: 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,
1130: 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,
1131: 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,
1132: 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};
1133: static PetscInt tripO[] = {0, 0,
1134: 0, 0,
1135: 0, 0,
1136: 0, -2, -2,
1137: -2, 0, -2,
1138: -2, -2, 0,
1139: 0, 0, 0,
1140: -2, 0, -2, -2,
1141: -2, 0, -2, -2,
1142: -2, 0, -2, -2,
1143: 0, -2, -2, 0,
1144: 0, -2, -2, 0,
1145: 0, -2, -2, 0,
1146: 0, 0, 0, -1, 0,
1147: 0, 0, 0, 0, -1,
1148: 0, 0, -1, 0, 0,
1149: 2, 0, 0, 0, 0,
1150: -3, 0, 0, -1, 0,
1151: -3, 0, 0, 0, -1,
1152: -3, 0, -1, 0, 0,
1153: -3, 0, 0, 0, 0};
1154: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new prisms.
1155: 2
1156: |\
1157: | \
1158: | \
1159: 0---1
1161: 2
1163: 0 1
1165: 2
1166: |\
1167: | \
1168: | \
1169: 0---1
1170: */
1171: static DMPolytopeType ttripT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR};
1172: static PetscInt ttripS[] = {3, 4};
1173: 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,
1174: 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,
1175: 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,
1176: 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,
1177: 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,
1178: 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,
1179: 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};
1180: static PetscInt ttripO[] = {0, 0, 0, 0,
1181: 0, 0, 0, 0,
1182: 0, 0, 0, 0,
1183: 0, 0, 0, -1, 0,
1184: 0, 0, 0, 0, -1,
1185: 0, 0, -1, 0, 0,
1186: 0, 0, 0, 0, 0};
1187: /* Add 1 edge and 4 tensor quads inside every tensor quad prism, making 4 new prisms. */
1188: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1189: static PetscInt tquadpS[] = {1, 4, 4};
1190: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1191: 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,
1192: 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,
1193: 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,
1194: 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,
1195: 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,
1196: 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,
1197: 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,
1198: 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};
1199: static PetscInt tquadpO[] = {0, 0,
1200: 0, 0, 0, 0,
1201: 0, 0, 0, 0,
1202: 0, 0, 0, 0,
1203: 0, 0, 0, 0,
1204: 0, 0, 0, 0, -1, 0,
1205: 0, 0, 0, 0, 0, -1,
1206: 0, 0, -1, 0, 0, 0,
1207: 0, 0, 0, -1, 0, 0};
1210: switch (source) {
1211: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
1212: case DM_POLYTOPE_SEGMENT: *Nt = 2; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
1213: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1214: case DM_POLYTOPE_TRIANGLE: *Nt = 2; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1215: case DM_POLYTOPE_QUADRILATERAL: *Nt = 3; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
1216: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1217: case DM_POLYTOPE_TETRAHEDRON: *Nt = 3; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1218: case DM_POLYTOPE_HEXAHEDRON: *Nt = 4; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
1219: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1220: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 2; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1221: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1222: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1223: }
1224: return(0);
1225: }
1227: static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1228: {
1230: /* Change tensor edges to segments */
1231: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_SEGMENT};
1232: static PetscInt tedgeS[] = {1};
1233: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
1234: static PetscInt tedgeO[] = { 0, 0};
1235: /* Add 1 vertex, 3 edges inside every triangle, making 3 new quadrilaterals.
1236: 2
1237: |\
1238: | \
1239: | \
1240: | \
1241: 0 1
1242: | \
1243: | \
1244: 2 1
1245: |\ / \
1246: | 2 1 \
1247: | \ / \
1248: 1 | 0
1249: | 0 \
1250: | | \
1251: | | \
1252: 0-0-0-----1-----1
1253: */
1254: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1255: static PetscInt triS[] = {1, 3, 3};
1256: static PetscInt triC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1257: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1258: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1259: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1260: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1261: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0};
1262: static PetscInt triO[] = {0, 0,
1263: 0, 0,
1264: 0, 0,
1265: 0, 0, -2, 0,
1266: 0, 0, 0, -2,
1267: 0, -2, 0, 0};
1268: /* Add 1 edge inside every tensor quad, making 2 new quadrilaterals
1269: 2----2----1----3----3
1270: | | |
1271: | | |
1272: | | |
1273: 4 A 6 B 5
1274: | | |
1275: | | |
1276: | | |
1277: 0----0----0----1----1
1278: */
1279: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL};
1280: static PetscInt tquadS[] = {1, 2};
1281: static PetscInt tquadC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1282: /* TODO Fix these */
1283: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1284: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 0};
1285: static PetscInt tquadO[] = {0, 0,
1286: 0, 0, -2, -2,
1287: 0, 0, -2, -2};
1288: /* Add 6 triangles inside every cell, making 4 new hexs
1289: TODO: Need different SubcellMap(). Need to make a struct with the function pointers in it
1290: 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
1291: 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]
1292: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1293: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1294: We make a new hex in each corner
1295: [v0, (e0, 0), (f0, 0), (e2, 0), (e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1296: [v1, (e4, 0), (f3, 0), (e1, 0), (e0, 0), (f0, 0), (c0, 0), (f1, 0)]
1297: [v2, (e1, 0), (f3, 0), (e5, 0), (e2, 0), (f2, 0), (c0, 0), (f0, 0)]
1298: [v3, (e4, 0), (f1, 0), (e3, 0), (e5, 0), (f2, 0), (c0, 0), (f3, 0)]
1299: We create a new face for each edge
1300: [(e3, 0), (f2, 0), (c0, 0), (f1, 0)]
1301: [(f0, 0), (e0, 0), (f1, 0), (c0, 0)]
1302: [(e2, 0), (f0, 0), (c0, 0), (f2, 0)]
1303: [(f3, 0), (e4, 0), (f1, 0), (c0, 0)]
1304: [(e1, 0), (f3, 0), (c0, 0), (f0, 0)]
1305: [(e5, 0), (f3, 0), (c0, 0), (f2, 0)]
1306: I could write a program to generate these from the first hex by acting with the symmetry group to take one subcell into another.
1307: */
1308: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1309: static PetscInt tetS[] = {1, 4, 6, 4};
1310: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1311: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1312: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1313: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1314: DM_POLYTOPE_SEGMENT, 1, 2, 2, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0,
1315: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1316: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1317: DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 3,
1318: DM_POLYTOPE_SEGMENT, 1, 3, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 0, 1,
1319: DM_POLYTOPE_SEGMENT, 1, 3, 2, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 1,
1320: 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,
1321: 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,
1322: 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,
1323: 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};
1324: static PetscInt tetO[] = {0, 0,
1325: 0, 0,
1326: 0, 0,
1327: 0, 0,
1328: 0, 0, -2, -2,
1329: -2, 0, 0, -2,
1330: 0, 0, -2, -2,
1331: -2, 0, 0, -2,
1332: 0, 0, -2, -2,
1333: 0, 0, -2, -2,
1334: 0, 0, 0, 0, 0, 0,
1335: 1, -1, 1, 0, 0, 3,
1336: 0, -4, 1, -1, 0, 3,
1337: 1, -4, 3, -2, -4, 3};
1338: /* Add 3 quads inside every triangular prism, making 4 new prisms. */
1339: static DMPolytopeType tripT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1340: static PetscInt tripS[] = {1, 5, 9, 6};
1341: static PetscInt tripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1342: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1343: DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_POINT, 0, 0,
1344: DM_POLYTOPE_POINT, 1, 3, 0, DM_POLYTOPE_POINT, 0, 0,
1345: DM_POLYTOPE_POINT, 1, 4, 0, DM_POLYTOPE_POINT, 0, 0,
1346: DM_POLYTOPE_SEGMENT, 1, 2, 3, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 1,
1347: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 1, 3, 3, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 2,
1348: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 1, DM_POLYTOPE_SEGMENT, 1, 4, 3,
1349: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1350: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1351: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1352: DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 2,
1353: DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 2,
1354: DM_POLYTOPE_SEGMENT, 0, 4, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 2,
1355: 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,
1356: 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,
1357: 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,
1358: 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,
1359: 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,
1360: 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};
1361: static PetscInt tripO[] = {0, 0,
1362: 0, 0,
1363: 0, 0,
1364: 0, 0,
1365: 0, 0,
1366: 0, 0, -2, -2,
1367: -2, 0, 0, -2,
1368: 0, -2, -2, 0,
1369: 0, 0, -2, -2,
1370: 0, 0, -2, -2,
1371: 0, 0, -2, -2,
1372: 0, -2, -2, 0,
1373: 0, -2, -2, 0,
1374: 0, -2, -2, 0,
1375: 0, 0, 0, -1, 0, 1,
1376: 0, 0, 0, 0, 0, -4,
1377: 0, 0, 0, 0, -1, 1,
1378: -4, 0, 0, -1, 0, 1,
1379: -4, 0, 0, 0, 0, -4,
1380: -4, 0, 0, 0, -1, 1};
1381: /* Add 3 tensor quads inside every tensor triangular prism, making 4 new tensor triangular prisms.
1382: 2
1383: |\
1384: | \
1385: | \
1386: 0---1
1388: 2
1390: 0 1
1392: 2
1393: |\
1394: | \
1395: | \
1396: 0---1
1397: */
1398: static DMPolytopeType ttripT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1399: static PetscInt ttripS[] = {1, 3, 3};
1400: static PetscInt ttripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1401: 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,
1402: 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,
1403: 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,
1404: 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,
1405: 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,
1406: 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};
1407: static PetscInt ttripO[] = {0, 0,
1408: 0, 0, 0, 0,
1409: 0, 0, 0, 0,
1410: 0, 0, 0, 0,
1411: 0, 0, 0, 0, -1, 0,
1412: 0, 0, 0, 0, 0, -1,
1413: 0, 0, 0, -1, 0, 0};
1414: /* TODO Add 3 quads inside every tensor triangular prism, making 4 new triangular prisms.
1415: 2
1416: |\
1417: | \
1418: | \
1419: 0---1
1421: 2
1423: 0 1
1425: 2
1426: |\
1427: | \
1428: | \
1429: 0---1
1430: */
1431: static DMPolytopeType ctripT[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_HEXAHEDRON};
1432: static PetscInt ctripS[] = {1, 3, 3};
1433: static PetscInt ctripC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1434: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0,
1435: DM_POLYTOPE_SEGMENT, 1, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 3, 0,
1436: DM_POLYTOPE_SEGMENT, 1, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 2, DM_POLYTOPE_SEGMENT, 1, 4, 0,
1437: 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,
1438: 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,
1439: 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};
1440: static PetscInt ctripO[] = {0, 0,
1441: 0, 0, -2, -2,
1442: 0, 0, -2, -2,
1443: 0, 0, -2, -2,
1444: -4, 0, 0, -1, 0, 1,
1445: -4, 0, 0, 0, 0, -4,
1446: -4, 0, 0, 0, -1, 1};
1447: /* Add 1 edge and 4 quads inside every tensor quad prism, making 4 new hexahedra. */
1448: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_QUAD_PRISM_TENSOR};
1449: static PetscInt tquadpS[] = {1, 4, 4};
1450: static PetscInt tquadpC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
1451: 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,
1452: 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,
1453: 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,
1454: 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,
1455: 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,
1456: 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,
1457: 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,
1458: 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};
1459: static PetscInt tquadpO[] = {0, 0,
1460: 0, 0, 0, 0,
1461: 0, 0, 0, 0,
1462: 0, 0, 0, 0,
1463: 0, 0, 0, 0,
1464: 0, 0, 0, 0, -1, 0,
1465: 0, 0, 0, 0, 0, -1,
1466: 0, 0, -1, 0, 0, 0,
1467: 0, 0, 0, -1, 0, 0};
1468: PetscBool convertTensor = PETSC_TRUE;
1471: if (convertTensor) {
1472: switch (source) {
1473: case DM_POLYTOPE_POINT:
1474: case DM_POLYTOPE_SEGMENT:
1475: case DM_POLYTOPE_QUADRILATERAL:
1476: case DM_POLYTOPE_HEXAHEDRON:
1477: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1478: break;
1479: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
1480: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 2; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
1481: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ctripT; *size = ctripS; *cone = ctripC; *ornt = ctripO; break;
1482: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 3; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
1483: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1484: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1485: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1486: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1487: }
1488: } else {
1489: switch (source) {
1490: case DM_POLYTOPE_POINT:
1491: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1492: case DM_POLYTOPE_SEGMENT:
1493: case DM_POLYTOPE_QUADRILATERAL:
1494: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1495: case DM_POLYTOPE_HEXAHEDRON:
1496: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1497: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1498: break;
1499: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1500: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1501: case DM_POLYTOPE_TRI_PRISM: *Nt = 4; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
1502: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 3; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
1503: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1504: }
1505: }
1506: return(0);
1507: }
1509: static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1510: {
1514: switch (source) {
1515: case DM_POLYTOPE_POINT:
1516: case DM_POLYTOPE_SEGMENT:
1517: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1518: case DM_POLYTOPE_TRIANGLE:
1519: case DM_POLYTOPE_TETRAHEDRON:
1520: case DM_POLYTOPE_TRI_PRISM:
1521: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1522: case DM_POLYTOPE_QUADRILATERAL:
1523: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1524: case DM_POLYTOPE_HEXAHEDRON:
1525: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1526: DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);
1527: break;
1528: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1529: }
1530: return(0);
1531: }
1533: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1534: {
1536: /* Add 1 vertex, 3 edges inside every triangle, making 3 new triangles.
1537: 2
1538: |\
1539: |\\
1540: | |\
1541: | \ \
1542: | | \
1543: | \ \
1544: | | \
1545: 2 \ \
1546: | | 1
1547: | 2 \
1548: | | \
1549: | /\ \
1550: | 0 1 |
1551: | / \ |
1552: |/ \|
1553: 0---0----1
1554: */
1555: static DMPolytopeType triT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
1556: static PetscInt triS[] = {1, 3, 3};
1557: static PetscInt triC[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1558: DM_POLYTOPE_POINT, 2, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1559: DM_POLYTOPE_POINT, 2, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1560: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 0,
1561: DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 1,
1562: DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 2};
1563: static PetscInt triO[] = {0, 0,
1564: 0, 0,
1565: 0, 0,
1566: 0, 0, -2,
1567: 0, 0, -2,
1568: 0, 0, -2};
1571: switch (source) {
1572: case DM_POLYTOPE_POINT:
1573: case DM_POLYTOPE_SEGMENT:
1574: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1575: case DM_POLYTOPE_QUADRILATERAL:
1576: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1577: case DM_POLYTOPE_TETRAHEDRON:
1578: case DM_POLYTOPE_HEXAHEDRON:
1579: case DM_POLYTOPE_TRI_PRISM:
1580: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1581: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1582: DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1583: break;
1584: case DM_POLYTOPE_TRIANGLE: *Nt = 3; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
1585: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1586: }
1587: return(0);
1588: }
1590: static PetscErrorCode DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1591: {
1593: /* Add 6 triangles inside every cell, making 4 new tets
1594: 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
1595: 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]
1596: called e3, e4, and e5. The faces of a tet, given in DMPlexGetRawFaces_Internal() are
1597: [v0, v1, v2], [v0, v3, v1], [v0, v2, v3], [v2, v1, v3]
1598: We make a new tet on each face
1599: [v0, v1, v2, (c0, 0)]
1600: [v0, v3, v1, (c0, 0)]
1601: [v0, v2, v3, (c0, 0)]
1602: [v2, v1, v3, (c0, 0)]
1603: We create a new face for each edge
1604: [v0, (c0, 0), v1 ]
1605: [v0, v2, (c0, 0)]
1606: [v2, v1, (c0, 0)]
1607: [v0, (c0, 0), v3 ]
1608: [v1, v3, (c0, 0)]
1609: [v3, v2, (c0, 0)]
1610: */
1611: static DMPolytopeType tetT[] = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE, DM_POLYTOPE_TETRAHEDRON};
1612: static PetscInt tetS[] = {1, 4, 6, 4};
1613: static PetscInt tetC[] = {DM_POLYTOPE_POINT, 3, 0, 0, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1614: DM_POLYTOPE_POINT, 3, 0, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1615: DM_POLYTOPE_POINT, 3, 0, 2, 0, 0, DM_POLYTOPE_POINT, 0, 0,
1616: DM_POLYTOPE_POINT, 3, 1, 0, 1, 0, DM_POLYTOPE_POINT, 0, 0,
1617: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 2, 0, 0, 0,
1618: DM_POLYTOPE_SEGMENT, 2, 0, 2, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 0,
1619: DM_POLYTOPE_SEGMENT, 2, 0, 1, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 0, 2,
1620: DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 2, 1, 0, 0,
1621: DM_POLYTOPE_SEGMENT, 2, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 3, DM_POLYTOPE_SEGMENT, 0, 1,
1622: DM_POLYTOPE_SEGMENT, 2, 2, 1, 0, DM_POLYTOPE_SEGMENT, 0, 2, DM_POLYTOPE_SEGMENT, 0, 3,
1623: DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 2,
1624: DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 0, DM_POLYTOPE_TRIANGLE, 0, 4,
1625: DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 0, 1, DM_POLYTOPE_TRIANGLE, 0, 3, DM_POLYTOPE_TRIANGLE, 0, 5,
1626: DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 0, 2, DM_POLYTOPE_TRIANGLE, 0, 5, DM_POLYTOPE_TRIANGLE, 0, 4};
1627: static PetscInt tetO[] = {0, 0,
1628: 0, 0,
1629: 0, 0,
1630: 0, 0,
1631: 0, -2, -2,
1632: -2, 0, -2,
1633: -2, 0, -2,
1634: 0, -2, -2,
1635: -2, 0, -2,
1636: -2, 0, -2,
1637: 0, 0, 0, 0,
1638: 0, 0, -3, 0,
1639: 0, -3, -3, 0,
1640: 0, -3, -1, -1};
1643: switch (source) {
1644: case DM_POLYTOPE_POINT:
1645: case DM_POLYTOPE_SEGMENT:
1646: case DM_POLYTOPE_POINT_PRISM_TENSOR:
1647: case DM_POLYTOPE_TRIANGLE:
1648: case DM_POLYTOPE_QUADRILATERAL:
1649: case DM_POLYTOPE_SEG_PRISM_TENSOR:
1650: case DM_POLYTOPE_HEXAHEDRON:
1651: case DM_POLYTOPE_TRI_PRISM:
1652: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1653: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1654: DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);
1655: break;
1656: case DM_POLYTOPE_TETRAHEDRON: *Nt = 4; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
1657: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1658: }
1659: return(0);
1660: }
1662: typedef struct {
1663: PetscInt n;
1664: PetscReal r;
1665: PetscScalar *h;
1666: PetscInt *Nt;
1667: DMPolytopeType **target;
1668: PetscInt **size;
1669: PetscInt **cone;
1670: PetscInt **ornt;
1671: } PlexRefiner_BL;
1673: static PetscErrorCode DMPlexCellRefinerSetUp_BL(DMPlexCellRefiner cr)
1674: {
1675: PlexRefiner_BL *crbl;
1677: PetscInt i,n;
1678: PetscReal r;
1679: PetscInt c1,c2,o1,o2;
1682: PetscNew(&crbl);
1683: cr->data = crbl;
1684: crbl->n = 1; /* 1 split -> 2 new cells */
1685: crbl->r = 1; /* linear progression */
1687: /* TODO: add setfromoptions to the refiners? */
1688: PetscOptionsGetInt(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_splits", &crbl->n, NULL);
1689: if (crbl->n < 1) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Number of splits %D must be positive",crbl->n);
1690: PetscOptionsGetReal(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_progression", &crbl->r, NULL);
1691: n = crbl->n;
1692: r = crbl->r;
1694: /* we only split DM_POLYTOPE_POINT_PRISM_TENSOR, DM_POLYTOPE_SEG_PRISM_TENSOR, DM_POLYTOPE_TRI_PRISM_TENSOR and DM_POLYTOPE_QUAD_PRISM_TENSOR */
1695: PetscMalloc5(4,&crbl->Nt,4,&crbl->target,4,&crbl->size,4,&crbl->cone,4,&crbl->ornt);
1697: /* progression */
1698: PetscMalloc1(n,&crbl->h);
1699: if (r > 1) {
1700: PetscReal d = (r-1.)/(PetscPowRealInt(r,n+1)-1.);
1702: crbl->h[0] = d;
1703: for (i = 1; i < n; i++) {
1704: d *= r;
1705: crbl->h[i] = crbl->h[i-1] + d;
1706: }
1707: } else { /* linear */
1708: for (i = 0; i < n; i++) crbl->h[i] = (i + 1.)/(n+1); /* linear */
1709: }
1711: /* DM_POLYTOPE_POINT_PRISM_TENSOR produces n points and n+1 tensor segments */
1712: c1 = 14+6*(n-1);
1713: o1 = 2*(n+1);
1714: crbl->Nt[0] = 2;
1716: PetscMalloc4(crbl->Nt[0],&crbl->target[0],crbl->Nt[0],&crbl->size[0],c1,&crbl->cone[0],o1,&crbl->ornt[0]);
1718: crbl->target[0][0] = DM_POLYTOPE_POINT;
1719: crbl->target[0][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1721: crbl->size[0][0] = n;
1722: crbl->size[0][1] = n+1;
1724: /* the tensor segments */
1725: crbl->cone[0][0] = DM_POLYTOPE_POINT;
1726: crbl->cone[0][1] = 1;
1727: crbl->cone[0][2] = 0;
1728: crbl->cone[0][3] = 0;
1729: crbl->cone[0][4] = DM_POLYTOPE_POINT;
1730: crbl->cone[0][5] = 0;
1731: crbl->cone[0][6] = 0;
1732: for (i = 0; i < n-1; i++) {
1733: crbl->cone[0][7+6*i+0] = DM_POLYTOPE_POINT;
1734: crbl->cone[0][7+6*i+1] = 0;
1735: crbl->cone[0][7+6*i+2] = i;
1736: crbl->cone[0][7+6*i+3] = DM_POLYTOPE_POINT;
1737: crbl->cone[0][7+6*i+4] = 0;
1738: crbl->cone[0][7+6*i+5] = i+1;
1739: }
1740: crbl->cone[0][7+6*(n-1)+0] = DM_POLYTOPE_POINT;
1741: crbl->cone[0][7+6*(n-1)+1] = 0;
1742: crbl->cone[0][7+6*(n-1)+2] = n-1;
1743: crbl->cone[0][7+6*(n-1)+3] = DM_POLYTOPE_POINT;
1744: crbl->cone[0][7+6*(n-1)+4] = 1;
1745: crbl->cone[0][7+6*(n-1)+5] = 1;
1746: crbl->cone[0][7+6*(n-1)+6] = 0;
1747: for (i = 0; i < o1; i++) crbl->ornt[0][i] = 0;
1749: /* DM_POLYTOPE_SEG_PRISM_TENSOR produces n segments and n+1 tensor quads */
1750: c1 = 8*n;
1751: c2 = 30+14*(n-1);
1752: o1 = 2*n;
1753: o2 = 4*(n+1);
1754: crbl->Nt[1] = 2;
1756: PetscMalloc4(crbl->Nt[1],&crbl->target[1],crbl->Nt[1],&crbl->size[1],c1+c2,&crbl->cone[1],o1+o2,&crbl->ornt[1]);
1758: crbl->target[1][0] = DM_POLYTOPE_SEGMENT;
1759: crbl->target[1][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1761: crbl->size[1][0] = n;
1762: crbl->size[1][1] = n+1;
1764: /* the segments */
1765: for (i = 0; i < n; i++) {
1766: crbl->cone[1][8*i+0] = DM_POLYTOPE_POINT;
1767: crbl->cone[1][8*i+1] = 1;
1768: crbl->cone[1][8*i+2] = 2;
1769: crbl->cone[1][8*i+3] = i;
1770: crbl->cone[1][8*i+4] = DM_POLYTOPE_POINT;
1771: crbl->cone[1][8*i+5] = 1;
1772: crbl->cone[1][8*i+6] = 3;
1773: crbl->cone[1][8*i+7] = i;
1774: }
1776: /* the tensor quads */
1777: crbl->cone[1][c1+ 0] = DM_POLYTOPE_SEGMENT;
1778: crbl->cone[1][c1+ 1] = 1;
1779: crbl->cone[1][c1+ 2] = 0;
1780: crbl->cone[1][c1+ 3] = 0;
1781: crbl->cone[1][c1+ 4] = DM_POLYTOPE_SEGMENT;
1782: crbl->cone[1][c1+ 5] = 0;
1783: crbl->cone[1][c1+ 6] = 0;
1784: crbl->cone[1][c1+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1785: crbl->cone[1][c1+ 8] = 1;
1786: crbl->cone[1][c1+ 9] = 2;
1787: crbl->cone[1][c1+10] = 0;
1788: crbl->cone[1][c1+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1789: crbl->cone[1][c1+12] = 1;
1790: crbl->cone[1][c1+13] = 3;
1791: crbl->cone[1][c1+14] = 0;
1792: for (i = 0; i < n-1; i++) {
1793: crbl->cone[1][c1+15+14*i+ 0] = DM_POLYTOPE_SEGMENT;
1794: crbl->cone[1][c1+15+14*i+ 1] = 0;
1795: crbl->cone[1][c1+15+14*i+ 2] = i;
1796: crbl->cone[1][c1+15+14*i+ 3] = DM_POLYTOPE_SEGMENT;
1797: crbl->cone[1][c1+15+14*i+ 4] = 0;
1798: crbl->cone[1][c1+15+14*i+ 5] = i+1;
1799: crbl->cone[1][c1+15+14*i+ 6] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1800: crbl->cone[1][c1+15+14*i+ 7] = 1;
1801: crbl->cone[1][c1+15+14*i+ 8] = 2;
1802: crbl->cone[1][c1+15+14*i+ 9] = i+1;
1803: crbl->cone[1][c1+15+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1804: crbl->cone[1][c1+15+14*i+11] = 1;
1805: crbl->cone[1][c1+15+14*i+12] = 3;
1806: crbl->cone[1][c1+15+14*i+13] = i+1;
1807: }
1808: crbl->cone[1][c1+15+14*(n-1)+ 0] = DM_POLYTOPE_SEGMENT;
1809: crbl->cone[1][c1+15+14*(n-1)+ 1] = 0;
1810: crbl->cone[1][c1+15+14*(n-1)+ 2] = n-1;
1811: crbl->cone[1][c1+15+14*(n-1)+ 3] = DM_POLYTOPE_SEGMENT;
1812: crbl->cone[1][c1+15+14*(n-1)+ 4] = 1;
1813: crbl->cone[1][c1+15+14*(n-1)+ 5] = 1;
1814: crbl->cone[1][c1+15+14*(n-1)+ 6] = 0;
1815: crbl->cone[1][c1+15+14*(n-1)+ 7] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1816: crbl->cone[1][c1+15+14*(n-1)+ 8] = 1;
1817: crbl->cone[1][c1+15+14*(n-1)+ 9] = 2;
1818: crbl->cone[1][c1+15+14*(n-1)+10] = n;
1819: crbl->cone[1][c1+15+14*(n-1)+11] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1820: crbl->cone[1][c1+15+14*(n-1)+12] = 1;
1821: crbl->cone[1][c1+15+14*(n-1)+13] = 3;
1822: crbl->cone[1][c1+15+14*(n-1)+14] = n;
1823: for (i = 0; i < o1+o2; i++) crbl->ornt[1][i] = 0;
1825: /* DM_POLYTOPE_TRI_PRISM_TENSOR produces n triangles and n+1 tensor triangular prisms */
1826: c1 = 12*n;
1827: c2 = 38+18*(n-1);
1828: o1 = 3*n;
1829: o2 = 5*(n+1);
1830: crbl->Nt[2] = 2;
1832: PetscMalloc4(crbl->Nt[2],&crbl->target[2],crbl->Nt[2],&crbl->size[2],c1+c2,&crbl->cone[2],o1+o2,&crbl->ornt[2]);
1834: crbl->target[2][0] = DM_POLYTOPE_TRIANGLE;
1835: crbl->target[2][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;
1837: crbl->size[2][0] = n;
1838: crbl->size[2][1] = n+1;
1840: /* the triangles */
1841: for (i = 0; i < n; i++) {
1842: crbl->cone[2][12*i+ 0] = DM_POLYTOPE_SEGMENT;
1843: crbl->cone[2][12*i+ 1] = 1;
1844: crbl->cone[2][12*i+ 2] = 2;
1845: crbl->cone[2][12*i+ 3] = i;
1846: crbl->cone[2][12*i+ 4] = DM_POLYTOPE_SEGMENT;
1847: crbl->cone[2][12*i+ 5] = 1;
1848: crbl->cone[2][12*i+ 6] = 3;
1849: crbl->cone[2][12*i+ 7] = i;
1850: crbl->cone[2][12*i+ 8] = DM_POLYTOPE_SEGMENT;
1851: crbl->cone[2][12*i+ 9] = 1;
1852: crbl->cone[2][12*i+10] = 4;
1853: crbl->cone[2][12*i+11] = i;
1854: }
1856: /* the triangular prisms */
1857: crbl->cone[2][c1+ 0] = DM_POLYTOPE_TRIANGLE;
1858: crbl->cone[2][c1+ 1] = 1;
1859: crbl->cone[2][c1+ 2] = 0;
1860: crbl->cone[2][c1+ 3] = 0;
1861: crbl->cone[2][c1+ 4] = DM_POLYTOPE_TRIANGLE;
1862: crbl->cone[2][c1+ 5] = 0;
1863: crbl->cone[2][c1+ 6] = 0;
1864: crbl->cone[2][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1865: crbl->cone[2][c1+ 8] = 1;
1866: crbl->cone[2][c1+ 9] = 2;
1867: crbl->cone[2][c1+10] = 0;
1868: crbl->cone[2][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1869: crbl->cone[2][c1+12] = 1;
1870: crbl->cone[2][c1+13] = 3;
1871: crbl->cone[2][c1+14] = 0;
1872: crbl->cone[2][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1873: crbl->cone[2][c1+16] = 1;
1874: crbl->cone[2][c1+17] = 4;
1875: crbl->cone[2][c1+18] = 0;
1876: for (i = 0; i < n-1; i++) {
1877: crbl->cone[2][c1+19+18*i+ 0] = DM_POLYTOPE_TRIANGLE;
1878: crbl->cone[2][c1+19+18*i+ 1] = 0;
1879: crbl->cone[2][c1+19+18*i+ 2] = i;
1880: crbl->cone[2][c1+19+18*i+ 3] = DM_POLYTOPE_TRIANGLE;
1881: crbl->cone[2][c1+19+18*i+ 4] = 0;
1882: crbl->cone[2][c1+19+18*i+ 5] = i+1;
1883: crbl->cone[2][c1+19+18*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1884: crbl->cone[2][c1+19+18*i+ 7] = 1;
1885: crbl->cone[2][c1+19+18*i+ 8] = 2;
1886: crbl->cone[2][c1+19+18*i+ 9] = i+1;
1887: crbl->cone[2][c1+19+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1888: crbl->cone[2][c1+19+18*i+11] = 1;
1889: crbl->cone[2][c1+19+18*i+12] = 3;
1890: crbl->cone[2][c1+19+18*i+13] = i+1;
1891: crbl->cone[2][c1+19+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1892: crbl->cone[2][c1+19+18*i+15] = 1;
1893: crbl->cone[2][c1+19+18*i+16] = 4;
1894: crbl->cone[2][c1+19+18*i+17] = i+1;
1895: }
1896: crbl->cone[2][c1+19+18*(n-1)+ 0] = DM_POLYTOPE_TRIANGLE;
1897: crbl->cone[2][c1+19+18*(n-1)+ 1] = 0;
1898: crbl->cone[2][c1+19+18*(n-1)+ 2] = n-1;
1899: crbl->cone[2][c1+19+18*(n-1)+ 3] = DM_POLYTOPE_TRIANGLE;
1900: crbl->cone[2][c1+19+18*(n-1)+ 4] = 1;
1901: crbl->cone[2][c1+19+18*(n-1)+ 5] = 1;
1902: crbl->cone[2][c1+19+18*(n-1)+ 6] = 0;
1903: crbl->cone[2][c1+19+18*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1904: crbl->cone[2][c1+19+18*(n-1)+ 8] = 1;
1905: crbl->cone[2][c1+19+18*(n-1)+ 9] = 2;
1906: crbl->cone[2][c1+19+18*(n-1)+10] = n;
1907: crbl->cone[2][c1+19+18*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1908: crbl->cone[2][c1+19+18*(n-1)+12] = 1;
1909: crbl->cone[2][c1+19+18*(n-1)+13] = 3;
1910: crbl->cone[2][c1+19+18*(n-1)+14] = n;
1911: crbl->cone[2][c1+19+18*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1912: crbl->cone[2][c1+19+18*(n-1)+16] = 1;
1913: crbl->cone[2][c1+19+18*(n-1)+17] = 4;
1914: crbl->cone[2][c1+19+18*(n-1)+18] = n;
1915: for (i = 0; i < o1+o2; i++) crbl->ornt[2][i] = 0;
1917: /* DM_POLYTOPE_QUAD_PRISM_TENSOR produces n quads and n+1 tensor quad prisms */
1918: c1 = 16*n;
1919: c2 = 46+22*(n-1);
1920: o1 = 4*n;
1921: o2 = 6*(n+1);
1922: crbl->Nt[3] = 2;
1924: PetscMalloc4(crbl->Nt[3],&crbl->target[3],crbl->Nt[3],&crbl->size[3],c1+c2,&crbl->cone[3],o1+o2,&crbl->ornt[3]);
1926: crbl->target[3][0] = DM_POLYTOPE_QUADRILATERAL;
1927: crbl->target[3][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;
1929: crbl->size[3][0] = n;
1930: crbl->size[3][1] = n+1;
1932: /* the quads */
1933: for (i = 0; i < n; i++) {
1934: crbl->cone[3][16*i+ 0] = DM_POLYTOPE_SEGMENT;
1935: crbl->cone[3][16*i+ 1] = 1;
1936: crbl->cone[3][16*i+ 2] = 2;
1937: crbl->cone[3][16*i+ 3] = i;
1938: crbl->cone[3][16*i+ 4] = DM_POLYTOPE_SEGMENT;
1939: crbl->cone[3][16*i+ 5] = 1;
1940: crbl->cone[3][16*i+ 6] = 3;
1941: crbl->cone[3][16*i+ 7] = i;
1942: crbl->cone[3][16*i+ 8] = DM_POLYTOPE_SEGMENT;
1943: crbl->cone[3][16*i+ 9] = 1;
1944: crbl->cone[3][16*i+10] = 4;
1945: crbl->cone[3][16*i+11] = i;
1946: crbl->cone[3][16*i+12] = DM_POLYTOPE_SEGMENT;
1947: crbl->cone[3][16*i+13] = 1;
1948: crbl->cone[3][16*i+14] = 5;
1949: crbl->cone[3][16*i+15] = i;
1950: }
1952: /* the quad prisms */
1953: crbl->cone[3][c1+ 0] = DM_POLYTOPE_QUADRILATERAL;
1954: crbl->cone[3][c1+ 1] = 1;
1955: crbl->cone[3][c1+ 2] = 0;
1956: crbl->cone[3][c1+ 3] = 0;
1957: crbl->cone[3][c1+ 4] = DM_POLYTOPE_QUADRILATERAL;
1958: crbl->cone[3][c1+ 5] = 0;
1959: crbl->cone[3][c1+ 6] = 0;
1960: crbl->cone[3][c1+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1961: crbl->cone[3][c1+ 8] = 1;
1962: crbl->cone[3][c1+ 9] = 2;
1963: crbl->cone[3][c1+10] = 0;
1964: crbl->cone[3][c1+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1965: crbl->cone[3][c1+12] = 1;
1966: crbl->cone[3][c1+13] = 3;
1967: crbl->cone[3][c1+14] = 0;
1968: crbl->cone[3][c1+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1969: crbl->cone[3][c1+16] = 1;
1970: crbl->cone[3][c1+17] = 4;
1971: crbl->cone[3][c1+18] = 0;
1972: crbl->cone[3][c1+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1973: crbl->cone[3][c1+20] = 1;
1974: crbl->cone[3][c1+21] = 5;
1975: crbl->cone[3][c1+22] = 0;
1976: for (i = 0; i < n-1; i++) {
1977: crbl->cone[3][c1+23+22*i+ 0] = DM_POLYTOPE_QUADRILATERAL;
1978: crbl->cone[3][c1+23+22*i+ 1] = 0;
1979: crbl->cone[3][c1+23+22*i+ 2] = i;
1980: crbl->cone[3][c1+23+22*i+ 3] = DM_POLYTOPE_QUADRILATERAL;
1981: crbl->cone[3][c1+23+22*i+ 4] = 0;
1982: crbl->cone[3][c1+23+22*i+ 5] = i+1;
1983: crbl->cone[3][c1+23+22*i+ 6] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1984: crbl->cone[3][c1+23+22*i+ 7] = 1;
1985: crbl->cone[3][c1+23+22*i+ 8] = 2;
1986: crbl->cone[3][c1+23+22*i+ 9] = i+1;
1987: crbl->cone[3][c1+23+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1988: crbl->cone[3][c1+23+22*i+11] = 1;
1989: crbl->cone[3][c1+23+22*i+12] = 3;
1990: crbl->cone[3][c1+23+22*i+13] = i+1;
1991: crbl->cone[3][c1+23+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1992: crbl->cone[3][c1+23+22*i+15] = 1;
1993: crbl->cone[3][c1+23+22*i+16] = 4;
1994: crbl->cone[3][c1+23+22*i+17] = i+1;
1995: crbl->cone[3][c1+23+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1996: crbl->cone[3][c1+23+22*i+19] = 1;
1997: crbl->cone[3][c1+23+22*i+20] = 5;
1998: crbl->cone[3][c1+23+22*i+21] = i+1;
1999: }
2000: crbl->cone[3][c1+23+22*(n-1)+ 0] = DM_POLYTOPE_QUADRILATERAL;
2001: crbl->cone[3][c1+23+22*(n-1)+ 1] = 0;
2002: crbl->cone[3][c1+23+22*(n-1)+ 2] = n-1;
2003: crbl->cone[3][c1+23+22*(n-1)+ 3] = DM_POLYTOPE_QUADRILATERAL;
2004: crbl->cone[3][c1+23+22*(n-1)+ 4] = 1;
2005: crbl->cone[3][c1+23+22*(n-1)+ 5] = 1;
2006: crbl->cone[3][c1+23+22*(n-1)+ 6] = 0;
2007: crbl->cone[3][c1+23+22*(n-1)+ 7] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2008: crbl->cone[3][c1+23+22*(n-1)+ 8] = 1;
2009: crbl->cone[3][c1+23+22*(n-1)+ 9] = 2;
2010: crbl->cone[3][c1+23+22*(n-1)+10] = n;
2011: crbl->cone[3][c1+23+22*(n-1)+11] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2012: crbl->cone[3][c1+23+22*(n-1)+12] = 1;
2013: crbl->cone[3][c1+23+22*(n-1)+13] = 3;
2014: crbl->cone[3][c1+23+22*(n-1)+14] = n;
2015: crbl->cone[3][c1+23+22*(n-1)+15] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2016: crbl->cone[3][c1+23+22*(n-1)+16] = 1;
2017: crbl->cone[3][c1+23+22*(n-1)+17] = 4;
2018: crbl->cone[3][c1+23+22*(n-1)+18] = n;
2019: crbl->cone[3][c1+23+22*(n-1)+19] = DM_POLYTOPE_SEG_PRISM_TENSOR;
2020: crbl->cone[3][c1+23+22*(n-1)+20] = 1;
2021: crbl->cone[3][c1+23+22*(n-1)+21] = 5;
2022: crbl->cone[3][c1+23+22*(n-1)+22] = n;
2023: for (i = 0; i < o1+o2; i++) crbl->ornt[3][i] = 0;
2024: return(0);
2025: }
2027: static PetscErrorCode DMPlexCellRefinerDestroy_BL(DMPlexCellRefiner cr)
2028: {
2029: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2033: PetscFree4(crbl->target[0],crbl->size[0],crbl->cone[0],crbl->ornt[0]);
2034: PetscFree4(crbl->target[1],crbl->size[1],crbl->cone[1],crbl->ornt[1]);
2035: PetscFree4(crbl->target[2],crbl->size[2],crbl->cone[2],crbl->ornt[2]);
2036: PetscFree4(crbl->target[3],crbl->size[3],crbl->cone[3],crbl->ornt[3]);
2037: PetscFree5(crbl->Nt,crbl->target,crbl->size,crbl->cone,crbl->ornt);
2038: PetscFree(crbl->h);
2039: PetscFree(cr->data);
2040: return(0);
2041: }
2043: static PetscErrorCode DMPlexCellRefinerRefine_BL(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
2044: {
2045: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2046: PetscErrorCode ierr;
2049: switch (source) {
2050: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2051: *Nt = crbl->Nt[0];
2052: *target = crbl->target[0];
2053: *size = crbl->size[0];
2054: *cone = crbl->cone[0];
2055: *ornt = crbl->ornt[0];
2056: break;
2057: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2058: *Nt = crbl->Nt[1];
2059: *target = crbl->target[1];
2060: *size = crbl->size[1];
2061: *cone = crbl->cone[1];
2062: *ornt = crbl->ornt[1];
2063: break;
2064: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2065: *Nt = crbl->Nt[2];
2066: *target = crbl->target[2];
2067: *size = crbl->size[2];
2068: *cone = crbl->cone[2];
2069: *ornt = crbl->ornt[2];
2070: break;
2071: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2072: *Nt = crbl->Nt[3];
2073: *target = crbl->target[3];
2074: *size = crbl->size[3];
2075: *cone = crbl->cone[3];
2076: *ornt = crbl->ornt[3];
2077: break;
2078: default:
2079: DMPlexCellRefinerRefine_None(cr,source,Nt,target,size,cone,ornt);
2080: }
2081: return(0);
2082: }
2084: static PetscErrorCode DMPlexCellRefinerMapSubcells_BL(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
2085: {
2086: /* We shift any input orientation in order to make it non-negative
2087: 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)
2088: The replica array r[po][r] gives the new replica number rnew given that the parent point has orientation po
2089: Overall, replica (r, o) in a parent with orientation 0 matches replica (rnew, onew) in a parent with orientation po
2090: */
2091: PetscInt tquad_seg_o[] = { 0, 1, -2, -1,
2092: 0, 1, -2, -1,
2093: -2, -1, 0, 1,
2094: -2, -1, 0, 1};
2095: PetscInt tquad_tquad_o[] = { 0, 1, -2, -1,
2096: 1, 0, -1, -2,
2097: -2, -1, 0, 1,
2098: -1, -2, 1, 0};
2099: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2100: const PetscInt n = crbl->n;
2104: *rnew = r;
2105: *onew = o;
2106: switch (pct) {
2107: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2108: if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR) {
2109: if (po == 0 || po == -1) {*rnew = r; *onew = o;}
2110: else if (po == 1 || po == -2) {*rnew = n - r; *onew = (o == 0 || o == -1) ? -2 : 0;}
2111: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid orientation %D for tensor segment", po);
2112: }
2113: break;
2114: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2115: switch (ct) {
2116: case DM_POLYTOPE_SEGMENT:
2117: *onew = tquad_seg_o[(po+2)*4+o+2];
2118: *rnew = r;
2119: break;
2120: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2121: *onew = tquad_tquad_o[(po+2)*4+o+2];
2122: *rnew = r;
2123: break;
2124: default: break;
2125: }
2126: break;
2127: default:
2128: DMPlexCellRefinerMapSubcells_None(cr, pct, po, ct, r, o, rnew, onew);
2129: }
2130: return(0);
2131: }
2133: static PetscErrorCode DMPlexCellRefinerMapCoordinates_BL(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
2134: {
2135: PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2136: PetscInt d;
2137: PetscErrorCode ierr;
2140: switch (pct) {
2141: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2142: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
2143: if (Nv != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent vertices %D",Nv);
2144: if (r >= crbl->n || r < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Invalid replica %D, must be in [0,%D)",r,crbl->n);
2145: for (d = 0; d < dE; d++) out[d] = in[d] + crbl->h[r] * (in[d + dE] - in[d]);
2146: break;
2147: default:
2148: DMPlexCellRefinerMapCoordinates_Barycenter(cr,pct,ct,r,Nv,dE,in,out);
2149: }
2150: return(0);
2151: }
2153: static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
2154: {
2155: PetscInt c, cN, *off;
2159: PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
2160: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
2161: const DMPolytopeType ct = (DMPolytopeType) c;
2162: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
2163: const DMPolytopeType ctNew = (DMPolytopeType) cN;
2164: DMPolytopeType *rct;
2165: PetscInt *rsize, *cone, *ornt;
2166: PetscInt Nct, n, i;
2168: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; break;}
2169: off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
2170: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
2171: const DMPolytopeType ict = (DMPolytopeType) ctOrder[i];
2172: const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];
2174: DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);
2175: if (ict == ct) {
2176: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
2177: if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
2178: break;
2179: }
2180: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
2181: }
2182: }
2183: }
2184: *offset = off;
2185: return(0);
2186: }
2188: static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[])
2189: {
2190: const PetscInt ctSize = DM_NUM_POLYTOPES+1;
2194: if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()");
2195: PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);
2196: PetscArraycpy(cr->ctStart, ctStart, ctSize);
2197: PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);
2198: return(0);
2199: }
2201: /* Construct cell type order since we must loop over cell types in depth order */
2202: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
2203: {
2204: PetscInt *ctO, *ctOInv;
2205: PetscInt c, d, off = 0;
2209: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
2210: for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
2211: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2212: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
2213: ctO[off++] = c;
2214: }
2215: }
2216: if (DMPolytopeTypeGetDim(ctCell) != 0) {
2217: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2218: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
2219: ctO[off++] = c;
2220: }
2221: }
2222: for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
2223: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2224: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
2225: ctO[off++] = c;
2226: }
2227: }
2228: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2229: if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
2230: ctO[off++] = c;
2231: }
2232: if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
2233: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
2234: ctOInv[ctO[c]] = c;
2235: }
2236: *ctOrder = ctO;
2237: *ctOrderInv = ctOInv;
2238: return(0);
2239: }
2241: PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)
2242: {
2243: DM dm = cr->dm;
2244: DMPolytopeType ctCell;
2245: PetscInt pStart, pEnd, p, c;
2250: if (cr->setupcalled) return(0);
2251: if (cr->ops->setup) {
2252: (*cr->ops->setup)(cr);
2253: }
2254: DMPlexGetChart(dm, &pStart, &pEnd);
2255: if (pEnd > pStart) {
2256: DMPlexGetCellType(dm, 0, &ctCell);
2257: } else {
2258: PetscInt dim;
2259: DMGetDimension(dm, &dim);
2260: switch (dim) {
2261: case 0: ctCell = DM_POLYTOPE_POINT;break;
2262: case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
2263: case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
2264: case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
2265: default: ctCell = DM_POLYTOPE_UNKNOWN;
2266: }
2267: }
2268: DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);
2269: /* Construct sizes and offsets for each cell type */
2270: if (!cr->ctStart) {
2271: PetscInt *ctS, *ctSN, *ctC, *ctCN;
2273: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
2274: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
2275: for (p = pStart; p < pEnd; ++p) {
2276: DMPolytopeType ct;
2277: DMPolytopeType *rct;
2278: PetscInt *rsize, *cone, *ornt;
2279: PetscInt Nct, n;
2281: DMPlexGetCellType(dm, p, &ct);
2282: if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
2283: ++ctC[ct];
2284: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
2285: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
2286: }
2287: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
2288: const PetscInt ct = cr->ctOrder[c];
2289: const PetscInt ctn = cr->ctOrder[c+1];
2291: ctS[ctn] = ctS[ct] + ctC[ct];
2292: ctSN[ctn] = ctSN[ct] + ctCN[ct];
2293: }
2294: PetscFree2(ctC, ctCN);
2295: cr->ctStart = ctS;
2296: cr->ctStartNew = ctSN;
2297: }
2298: CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);
2299: cr->setupcalled = PETSC_TRUE;
2300: return(0);
2301: }
2303: static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v)
2304: {
2308: PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);
2309: return(0);
2310: }
2312: /*
2313: DMPlexCellRefinerView - Views a DMPlexCellRefiner object
2315: Collective on cr
2317: Input Parameters:
2318: + cr - The DMPlexCellRefiner object
2319: - viewer - The PetscViewer object
2321: Level: beginner
2323: .seealso: DMPlexCellRefinerCreate()
2324: */
2325: static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer)
2326: {
2327: PetscBool iascii;
2333: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);}
2334: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2335: PetscViewerASCIIPushTab(viewer);
2336: if (iascii) {DMPlexCellRefinerView_Ascii(cr, viewer);}
2337: PetscViewerASCIIPopTab(viewer);
2338: return(0);
2339: }
2341: PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr)
2342: {
2343: PetscInt c;
2347: if (!*cr) return(0);
2349: if ((*cr)->ops->destroy) {
2350: ((*cr)->ops->destroy)(*cr);
2351: }
2352: PetscObjectDereference((PetscObject) (*cr)->dm);
2353: PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);
2354: PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);
2355: PetscFree((*cr)->offset);
2356: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
2357: PetscFEDestroy(&(*cr)->coordFE[c]);
2358: PetscFEGeomDestroy(&(*cr)->refGeom[c]);
2359: }
2360: PetscFree2((*cr)->coordFE, (*cr)->refGeom);
2361: PetscHeaderDestroy(cr);
2362: return(0);
2363: }
2365: PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr)
2366: {
2367: DMPlexCellRefiner tmp;
2368: PetscErrorCode ierr;
2372: *cr = NULL;
2373: PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);
2374: tmp->setupcalled = PETSC_FALSE;
2376: tmp->dm = dm;
2377: PetscObjectReference((PetscObject) dm);
2378: DMPlexGetCellRefinerType(dm, &tmp->type);
2379: switch (tmp->type) {
2380: case DM_REFINER_REGULAR:
2381: tmp->ops->refine = DMPlexCellRefinerRefine_Regular;
2382: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_Regular;
2383: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_Regular;
2384: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_Regular;
2385: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
2386: tmp->ops->getaffinetransforms = DMPlexCellRefinerGetAffineTransforms_Regular;
2387: tmp->ops->getaffinefacetransforms = DMPlexCellRefinerGetAffineFaceTransforms_Regular;
2388: break;
2389: case DM_REFINER_TO_BOX:
2390: tmp->ops->refine = DMPlexCellRefinerRefine_ToBox;
2391: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToBox;
2392: tmp->ops->getcellvertices = DMPlexCellRefinerGetCellVertices_ToBox;
2393: tmp->ops->getsubcellvertices = DMPlexCellRefinerGetSubcellVertices_ToBox;
2394: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
2395: break;
2396: case DM_REFINER_TO_SIMPLEX:
2397: tmp->ops->refine = DMPlexCellRefinerRefine_ToSimplex;
2398: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_ToSimplex;
2399: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
2400: break;
2401: case DM_REFINER_ALFELD2D:
2402: tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld2D;
2403: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
2404: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
2405: break;
2406: case DM_REFINER_ALFELD3D:
2407: tmp->ops->refine = DMPlexCellRefinerRefine_Alfeld3D;
2408: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_None;
2409: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_Barycenter;
2410: break;
2411: case DM_REFINER_BOUNDARYLAYER:
2412: tmp->ops->setup = DMPlexCellRefinerSetUp_BL;
2413: tmp->ops->destroy = DMPlexCellRefinerDestroy_BL;
2414: tmp->ops->refine = DMPlexCellRefinerRefine_BL;
2415: tmp->ops->mapsubcells = DMPlexCellRefinerMapSubcells_BL;
2416: tmp->ops->mapcoords = DMPlexCellRefinerMapCoordinates_BL;
2417: break;
2418: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid cell refiner type %s", DMPlexCellRefinerTypes[tmp->type]);
2419: }
2420: PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);
2421: *cr = tmp;
2422: return(0);
2423: }
2425: /*@
2426: DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell
2428: Input Parameters:
2429: + cr - The DMPlexCellRefiner object
2430: - ct - The cell type
2432: Output Parameters:
2433: + Nc - The number of subcells produced from this cell type
2434: . v0 - The translation of the first vertex for each subcell
2435: . J - The Jacobian for each subcell (map from reference cell to subcell)
2436: - invJ - The inverse Jacobian for each subcell
2438: Level: developer
2440: .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create()
2441: @*/
2442: PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
2443: {
2447: if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner");
2448: (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);
2449: return(0);
2450: }
2452: /*@
2453: DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell
2455: Input Parameters:
2456: + cr - The DMPlexCellRefiner object
2457: - ct - The cell type
2459: Output Parameters:
2460: + Nf - The number of faces for this cell type
2461: . v0 - The translation of the first vertex for each face
2462: . J - The Jacobian for each face (map from original cell to subcell)
2463: . invJ - The inverse Jacobian for each face
2464: - detJ - The determinant of the Jacobian for each face
2466: Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
2468: Level: developer
2470: .seealso: DMPlexCellRefinerGetAffineTransforms(), Create()
2471: @*/
2472: PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
2473: {
2477: if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner");
2478: (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);
2479: return(0);
2480: }
2482: /* Numbering regularly refined meshes
2484: We want the numbering of the new mesh to respect the same depth stratification as the old mesh. We first compute
2485: the number of new points at each depth. This means that offsets for each depth can be computed, making no assumptions
2486: about the order of different cell types.
2488: However, when we want to order different depth strata, it will be very useful to make assumptions about contiguous
2489: numbering of different cell types, especially if we want to compute new numberings without communication. Therefore, we
2490: will require that cells are numbering contiguously for each cell type, and that those blocks come in the same order as
2491: the cell type enumeration within a given depth stratum.
2493: Thus, at each depth, each cell type will add a certain number of points at that depth. To get the new point number, we
2494: start at the new depth offset, run through all prior cell types incrementing by the total addition from that type, then
2495: offset by the old cell type number and replica number for the insertion.
2496: */
2497: PetscErrorCode DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
2498: {
2499: DMPolytopeType *rct;
2500: PetscInt *rsize, *cone, *ornt;
2501: PetscInt Nct, n;
2502: PetscInt off = cr->offset[ct*DM_NUM_POLYTOPES+ctNew];
2503: PetscInt ctS = cr->ctStart[ct], ctE = cr->ctStart[cr->ctOrder[cr->ctOrderInv[ct]+1]];
2504: PetscInt ctSN = cr->ctStartNew[ctNew], ctEN = cr->ctStartNew[cr->ctOrder[cr->ctOrderInv[ctNew]+1]];
2505: PetscInt newp = ctSN;
2509: 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);
2510: if (off < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s does not produce type %s", DMPolytopeTypes[ct], DMPolytopeTypes[ctNew]);
2512: newp += off;
2513: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);
2514: for (n = 0; n < Nct; ++n) {
2515: if (rct[n] == ctNew) {
2516: if (rsize[n] && r >= rsize[n])
2517: 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]);
2518: newp += (p - ctS) * rsize[n] + r;
2519: break;
2520: }
2521: }
2523: 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);
2524: *pNew = newp;
2525: return(0);
2526: }
2528: static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm)
2529: {
2530: DM dm = cr->dm;
2531: PetscInt pStart, pEnd, p, pNew;
2532: PetscErrorCode ierr;
2535: /* Must create the celltype label here so that we do not automatically try to compute the types */
2536: DMCreateLabel(rdm, "celltype");
2537: DMPlexGetChart(dm, &pStart, &pEnd);
2538: for (p = pStart; p < pEnd; ++p) {
2539: DMPolytopeType ct;
2540: DMPolytopeType *rct;
2541: PetscInt *rsize, *rcone, *rornt;
2542: PetscInt Nct, n, r;
2544: DMPlexGetCellType(dm, p, &ct);
2545: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2546: for (n = 0; n < Nct; ++n) {
2547: for (r = 0; r < rsize[n]; ++r) {
2548: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2549: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
2550: DMPlexSetCellType(rdm, pNew, rct[n]);
2551: }
2552: }
2553: }
2554: {
2555: DMLabel ctLabel;
2556: DM_Plex *plex = (DM_Plex *) rdm->data;
2558: DMPlexGetCellTypeLabel(rdm, &ctLabel);
2559: PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
2560: }
2561: return(0);
2562: }
2564: static PetscErrorCode DMPlexCellRefinerSetCones(DMPlexCellRefiner cr, DM rdm)
2565: {
2566: DM dm = cr->dm;
2567: DMPolytopeType ct;
2568: PetscInt *coneNew, *orntNew;
2569: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
2573: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
2574: PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
2575: DMPlexGetChart(dm, &pStart, &pEnd);
2576: for (p = pStart; p < pEnd; ++p) {
2577: const PetscInt *cone, *ornt;
2578: PetscInt coff, ooff, c;
2579: DMPolytopeType *rct;
2580: PetscInt *rsize, *rcone, *rornt;
2581: PetscInt Nct, n, r;
2582: DMPlexGetCellType(dm, p, &ct);
2583: DMPlexGetCone(dm, p, &cone);
2584: DMPlexGetConeOrientation(dm, p, &ornt);
2585: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2586: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
2587: const DMPolytopeType ctNew = rct[n];
2588: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
2590: for (r = 0; r < rsize[n]; ++r) {
2591: /* pNew is a subcell produced by subdividing p */
2592: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2593: for (c = 0; c < csizeNew; ++c) {
2594: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
2595: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
2596: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
2597: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
2598: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
2599: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
2600: const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
2601: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
2602: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
2603: PetscInt lc;
2605: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
2606: for (lc = 0; lc < fn; ++lc) {
2607: const PetscInt *ppornt;
2608: PetscInt pcp;
2610: DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);
2611: ppp = pp;
2612: pp = pcone[pcp];
2613: DMPlexGetCellType(dm, pp, &pct);
2614: DMPlexGetCone(dm, pp, &pcone);
2615: DMPlexGetConeOrientation(dm, ppp, &ppornt);
2616: if (po < 0 && pct != DM_POLYTOPE_POINT) {
2617: const PetscInt pornt = ppornt[pcp];
2618: const PetscInt pcsize = DMPolytopeTypeGetConeSize(pct);
2619: const PetscInt pcstart = pornt < 0 ? -(pornt+1) : pornt;
2620: const PetscInt rcstart = (pcstart+pcsize-1)%pcsize;
2621: po = pornt < 0 ? -(rcstart+1) : rcstart;
2622: } else {
2623: po = ppornt[pcp];
2624: }
2625: }
2626: pr = rcone[coff++];
2627: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
2628: DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);
2629: DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);
2630: orntNew[c] = fo;
2631: }
2632: DMPlexSetCone(rdm, pNew, coneNew);
2633: DMPlexSetConeOrientation(rdm, pNew, orntNew);
2634: }
2635: }
2636: }
2637: PetscFree2(coneNew, orntNew);
2638: DMPlexSymmetrize(rdm);
2639: DMPlexStratify(rdm);
2640: return(0);
2641: }
2643: static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe)
2644: {
2648: if (!cr->coordFE[ct]) {
2649: PetscInt dim, cdim;
2650: PetscBool isSimplex;
2652: switch (ct) {
2653: case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break;
2654: case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break;
2655: case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
2656: case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break;
2657: case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break;
2658: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
2659: }
2660: DMGetCoordinateDim(cr->dm, &cdim);
2661: PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);
2662: {
2663: PetscDualSpace dsp;
2664: PetscQuadrature quad;
2665: DM K;
2666: PetscFEGeom *cg;
2667: PetscReal *Xq, *xq, *wq;
2668: PetscInt Nq, q;
2670: DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);
2671: PetscMalloc1(Nq*cdim, &xq);
2672: for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q];
2673: PetscMalloc1(Nq, &wq);
2674: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
2675: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
2676: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
2677: PetscFESetQuadrature(cr->coordFE[ct], quad);
2679: PetscFEGetDualSpace(cr->coordFE[ct], &dsp);
2680: PetscDualSpaceGetDM(dsp, &K);
2681: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);
2682: cg = cr->refGeom[ct];
2683: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2684: PetscQuadratureDestroy(&quad);
2685: }
2686: }
2687: *fe = cr->coordFE[ct];
2688: return(0);
2689: }
2691: /*
2692: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
2694: Not collective
2696: Input Parameters:
2697: + cr - The DMPlexCellRefiner
2698: . ct - The type of the parent cell
2699: . rct - The type of the produced cell
2700: . r - The index of the produced cell
2701: - x - The localized coordinates for the parent cell
2703: Output Parameter:
2704: . xr - The localized coordinates for the produced cell
2706: Level: developer
2708: .seealso: DMPlexCellRefinerSetCoordinates()
2709: */
2710: static PetscErrorCode DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
2711: {
2712: PetscFE fe = NULL;
2713: PetscInt cdim, Nv, v, *subcellV;
2717: DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);
2718: DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);
2719: PetscFEGetNumComponents(fe, &cdim);
2720: for (v = 0; v < Nv; ++v) {
2721: PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);
2722: }
2723: return(0);
2724: }
2726: static PetscErrorCode DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr, DM rdm)
2727: {
2728: DM dm = cr->dm, cdm;
2729: PetscSection coordSection, coordSectionNew;
2730: Vec coordsLocal, coordsLocalNew;
2731: const PetscScalar *coords;
2732: PetscScalar *coordsNew;
2733: const DMBoundaryType *bd;
2734: const PetscReal *maxCell, *L;
2735: PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
2736: PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
2737: PetscErrorCode ierr;
2740: DMGetCoordinateDM(dm, &cdm);
2741: DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
2742: /* Determine if we need to localize coordinates when generating them */
2743: if (isperiodic) {
2744: localizeVertices = PETSC_TRUE;
2745: if (!maxCell) {
2746: PetscBool localized;
2747: DMGetCoordinatesLocalized(dm, &localized);
2748: if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
2749: localizeCells = PETSC_TRUE;
2750: }
2751: }
2753: DMGetCoordinateSection(dm, &coordSection);
2754: PetscSectionGetFieldComponents(coordSection, 0, &dE);
2755: if (maxCell) {
2756: PetscReal maxCellNew[3];
2758: for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
2759: DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
2760: } else {
2761: DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
2762: }
2763: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
2764: PetscSectionSetNumFields(coordSectionNew, 1);
2765: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
2766: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
2767: if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0, vEndNew);}
2768: else {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}
2770: /* Localization should be inherited */
2771: /* Stefano calculates parent cells for each new cell for localization */
2772: /* Localized cells need coordinates of closure */
2773: for (v = vStartNew; v < vEndNew; ++v) {
2774: PetscSectionSetDof(coordSectionNew, v, dE);
2775: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
2776: }
2777: if (localizeCells) {
2778: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2779: for (c = cStart; c < cEnd; ++c) {
2780: PetscInt dof;
2782: PetscSectionGetDof(coordSection, c, &dof);
2783: if (dof) {
2784: DMPolytopeType ct;
2785: DMPolytopeType *rct;
2786: PetscInt *rsize, *rcone, *rornt;
2787: PetscInt dim, cNew, Nct, n, r;
2789: DMPlexGetCellType(dm, c, &ct);
2790: dim = DMPolytopeTypeGetDim(ct);
2791: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2792: /* This allows for different cell types */
2793: for (n = 0; n < Nct; ++n) {
2794: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
2795: for (r = 0; r < rsize[n]; ++r) {
2796: PetscInt *closure = NULL;
2797: PetscInt clSize, cl, Nv = 0;
2799: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);
2800: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2801: for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
2802: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
2803: PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
2804: PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
2805: }
2806: }
2807: }
2808: }
2809: }
2810: PetscSectionSetUp(coordSectionNew);
2811: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
2812: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
2813: {
2814: VecType vtype;
2815: PetscInt coordSizeNew, bs;
2816: const char *name;
2818: DMGetCoordinatesLocal(dm, &coordsLocal);
2819: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
2820: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
2821: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
2822: PetscObjectGetName((PetscObject) coordsLocal, &name);
2823: PetscObjectSetName((PetscObject) coordsLocalNew, name);
2824: VecGetBlockSize(coordsLocal, &bs);
2825: VecSetBlockSize(coordsLocalNew, bs);
2826: VecGetType(coordsLocal, &vtype);
2827: VecSetType(coordsLocalNew, vtype);
2828: }
2829: VecGetArrayRead(coordsLocal, &coords);
2830: VecGetArray(coordsLocalNew, &coordsNew);
2831: PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
2832: DMPlexGetChart(dm, &pStart, &pEnd);
2833: /* First set coordinates for vertices*/
2834: for (p = pStart; p < pEnd; ++p) {
2835: DMPolytopeType ct;
2836: DMPolytopeType *rct;
2837: PetscInt *rsize, *rcone, *rornt;
2838: PetscInt Nct, n, r;
2839: PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;
2841: DMPlexGetCellType(dm, p, &ct);
2842: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2843: for (n = 0; n < Nct; ++n) {
2844: if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
2845: }
2846: if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2847: PetscInt dof;
2848: PetscSectionGetDof(coordSection, p, &dof);
2849: if (dof) isLocalized = PETSC_TRUE;
2850: }
2851: if (hasVertex) {
2852: const PetscScalar *icoords = NULL;
2853: PetscScalar *pcoords = NULL;
2854: PetscInt Nc, Nv, v, d;
2856: DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
2858: icoords = pcoords;
2859: Nv = Nc/dE;
2860: if (ct != DM_POLYTOPE_POINT) {
2861: if (localizeVertices) {
2862: PetscScalar anchor[3];
2864: for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
2865: if (!isLocalized) {
2866: for (v = 0; v < Nv; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
2867: } else {
2868: Nv = Nc/(2*dE);
2869: icoords = pcoords + Nv*dE;
2870: for (v = Nv; v < Nv*2; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
2871: }
2872: }
2873: }
2874: for (n = 0; n < Nct; ++n) {
2875: if (rct[n] != DM_POLYTOPE_POINT) continue;
2876: for (r = 0; r < rsize[n]; ++r) {
2877: PetscScalar vcoords[3];
2878: PetscInt vNew, off;
2880: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);
2881: PetscSectionGetOffset(coordSectionNew, vNew, &off);
2882: DMPlexCellRefinerMapCoordinates(cr, ct, rct[n], r, Nv, dE, icoords, vcoords);
2883: DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);
2884: }
2885: }
2886: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
2887: }
2888: }
2889: /* Then set coordinates for cells by localizing */
2890: for (p = pStart; p < pEnd; ++p) {
2891: DMPolytopeType ct;
2892: DMPolytopeType *rct;
2893: PetscInt *rsize, *rcone, *rornt;
2894: PetscInt Nct, n, r;
2895: PetscBool isLocalized = PETSC_FALSE;
2897: DMPlexGetCellType(dm, p, &ct);
2898: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
2899: if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2900: PetscInt dof;
2901: PetscSectionGetDof(coordSection, p, &dof);
2902: if (dof) isLocalized = PETSC_TRUE;
2903: }
2904: if (isLocalized) {
2905: const PetscScalar *pcoords;
2907: DMPlexPointLocalRead(cdm, p, coords, &pcoords);
2908: for (n = 0; n < Nct; ++n) {
2909: const PetscInt Nr = rsize[n];
2911: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2912: for (r = 0; r < Nr; ++r) {
2913: PetscInt pNew, offNew;
2915: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
2916: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
2917: cell to the ones it produces. */
2918: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
2919: PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
2920: DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
2921: }
2922: }
2923: }
2924: }
2925: VecRestoreArrayRead(coordsLocal, &coords);
2926: VecRestoreArray(coordsLocalNew, &coordsNew);
2927: DMSetCoordinatesLocal(rdm, coordsLocalNew);
2928: /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
2929: VecDestroy(&coordsLocalNew);
2930: PetscSectionDestroy(&coordSectionNew);
2931: if (!localizeCells) {DMLocalizeCoordinates(rdm);}
2932: return(0);
2933: }
2935: /*@
2936: DMPlexCreateProcessSF - Create an SF which just has process connectivity
2938: Collective on dm
2940: Input Parameters:
2941: + dm - The DM
2942: - sfPoint - The PetscSF which encodes point connectivity
2944: Output Parameters:
2945: + processRanks - A list of process neighbors, or NULL
2946: - sfProcess - An SF encoding the process connectivity, or NULL
2948: Level: developer
2950: .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
2951: @*/
2952: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
2953: {
2954: PetscInt numRoots, numLeaves, l;
2955: const PetscInt *localPoints;
2956: const PetscSFNode *remotePoints;
2957: PetscInt *localPointsNew;
2958: PetscSFNode *remotePointsNew;
2959: PetscInt *ranks, *ranksNew;
2960: PetscMPIInt size;
2961: PetscErrorCode ierr;
2968: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
2969: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2970: PetscMalloc1(numLeaves, &ranks);
2971: for (l = 0; l < numLeaves; ++l) {
2972: ranks[l] = remotePoints[l].rank;
2973: }
2974: PetscSortRemoveDupsInt(&numLeaves, ranks);
2975: PetscMalloc1(numLeaves, &ranksNew);
2976: PetscMalloc1(numLeaves, &localPointsNew);
2977: PetscMalloc1(numLeaves, &remotePointsNew);
2978: for (l = 0; l < numLeaves; ++l) {
2979: ranksNew[l] = ranks[l];
2980: localPointsNew[l] = l;
2981: remotePointsNew[l].index = 0;
2982: remotePointsNew[l].rank = ranksNew[l];
2983: }
2984: PetscFree(ranks);
2985: if (processRanks) {ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);}
2986: else {PetscFree(ranksNew);}
2987: if (sfProcess) {
2988: PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);
2989: PetscObjectSetName((PetscObject) *sfProcess, "Process SF");
2990: PetscSFSetFromOptions(*sfProcess);
2991: PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
2992: }
2993: return(0);
2994: }
2996: static PetscErrorCode DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr, DM rdm)
2997: {
2998: DM dm = cr->dm;
2999: DMPlexCellRefiner *crRem;
3000: PetscSF sf, sfNew, sfProcess;
3001: IS processRanks;
3002: MPI_Datatype ctType;
3003: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
3004: const PetscInt *localPoints, *neighbors;
3005: const PetscSFNode *remotePoints;
3006: PetscInt *localPointsNew;
3007: PetscSFNode *remotePointsNew;
3008: PetscInt *ctStartRem, *ctStartNewRem;
3009: PetscInt ctSize = DM_NUM_POLYTOPES+1, numNeighbors, n, pStartNew, pEndNew, pNew, pNewRem;
3010: PetscErrorCode ierr;
3013: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
3014: DMGetPointSF(dm, &sf);
3015: DMGetPointSF(rdm, &sfNew);
3016: /* Calculate size of new SF */
3017: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
3018: if (numRoots < 0) return(0);
3019: for (l = 0; l < numLeaves; ++l) {
3020: const PetscInt p = localPoints[l];
3021: DMPolytopeType ct;
3022: DMPolytopeType *rct;
3023: PetscInt *rsize, *rcone, *rornt;
3024: PetscInt Nct, n;
3026: DMPlexGetCellType(dm, p, &ct);
3027: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3028: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
3029: }
3030: /* Communicate ctStart and cStartNew for each remote rank */
3031: DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);
3032: ISGetLocalSize(processRanks, &numNeighbors);
3033: PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);
3034: MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);
3035: MPI_Type_commit(&ctType);
3036: PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);
3037: PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);
3038: PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
3039: PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);
3040: MPI_Type_free(&ctType);
3041: PetscSFDestroy(&sfProcess);
3042: PetscMalloc1(numNeighbors, &crRem);
3043: for (n = 0; n < numNeighbors; ++n) {
3044: DMPlexCellRefinerCreate(dm, &crRem[n]);
3045: DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]);
3046: DMPlexCellRefinerSetUp(crRem[n]);
3047: }
3048: PetscFree2(ctStartRem, ctStartNewRem);
3049: /* Calculate new point SF */
3050: PetscMalloc1(numLeavesNew, &localPointsNew);
3051: PetscMalloc1(numLeavesNew, &remotePointsNew);
3052: ISGetIndices(processRanks, &neighbors);
3053: for (l = 0, m = 0; l < numLeaves; ++l) {
3054: PetscInt p = localPoints[l];
3055: PetscInt pRem = remotePoints[l].index;
3056: PetscMPIInt rankRem = remotePoints[l].rank;
3057: DMPolytopeType ct;
3058: DMPolytopeType *rct;
3059: PetscInt *rsize, *rcone, *rornt;
3060: PetscInt neighbor, Nct, n, r;
3062: PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);
3063: if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem);
3064: DMPlexGetCellType(dm, p, &ct);
3065: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3066: for (n = 0; n < Nct; ++n) {
3067: for (r = 0; r < rsize[n]; ++r) {
3068: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);
3069: DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);
3070: localPointsNew[m] = pNew;
3071: remotePointsNew[m].index = pNewRem;
3072: remotePointsNew[m].rank = rankRem;
3073: ++m;
3074: }
3075: }
3076: }
3077: for (n = 0; n < numNeighbors; ++n) {DMPlexCellRefinerDestroy(&crRem[n]);}
3078: PetscFree(crRem);
3079: if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew);
3080: ISRestoreIndices(processRanks, &neighbors);
3081: ISDestroy(&processRanks);
3082: {
3083: PetscSFNode *rp, *rtmp;
3084: PetscInt *lp, *idx, *ltmp, i;
3086: /* SF needs sorted leaves to correct calculate Gather */
3087: PetscMalloc1(numLeavesNew, &idx);
3088: PetscMalloc1(numLeavesNew, &lp);
3089: PetscMalloc1(numLeavesNew, &rp);
3090: for (i = 0; i < numLeavesNew; ++i) {
3091: 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);
3092: idx[i] = i;
3093: }
3094: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
3095: for (i = 0; i < numLeavesNew; ++i) {
3096: lp[i] = localPointsNew[idx[i]];
3097: rp[i] = remotePointsNew[idx[i]];
3098: }
3099: ltmp = localPointsNew;
3100: localPointsNew = lp;
3101: rtmp = remotePointsNew;
3102: remotePointsNew = rp;
3103: PetscFree(idx);
3104: PetscFree(ltmp);
3105: PetscFree(rtmp);
3106: }
3107: PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
3108: return(0);
3109: }
3111: static PetscErrorCode RefineLabel_Internal(DMPlexCellRefiner cr, DMLabel label, DMLabel labelNew)
3112: {
3113: DM dm = cr->dm;
3114: IS valueIS;
3115: const PetscInt *values;
3116: PetscInt defVal, Nv, val;
3117: PetscErrorCode ierr;
3120: DMLabelGetDefaultValue(label, &defVal);
3121: DMLabelSetDefaultValue(labelNew, defVal);
3122: DMLabelGetValueIS(label, &valueIS);
3123: ISGetLocalSize(valueIS, &Nv);
3124: ISGetIndices(valueIS, &values);
3125: for (val = 0; val < Nv; ++val) {
3126: IS pointIS;
3127: const PetscInt *points;
3128: PetscInt numPoints, p;
3130: /* Ensure refined label is created with same number of strata as
3131: * original (even if no entries here). */
3132: DMLabelAddStratum(labelNew, values[val]);
3133: DMLabelGetStratumIS(label, values[val], &pointIS);
3134: ISGetLocalSize(pointIS, &numPoints);
3135: ISGetIndices(pointIS, &points);
3136: for (p = 0; p < numPoints; ++p) {
3137: const PetscInt point = points[p];
3138: DMPolytopeType ct;
3139: DMPolytopeType *rct;
3140: PetscInt *rsize, *rcone, *rornt;
3141: PetscInt Nct, n, r, pNew;
3143: DMPlexGetCellType(dm, point, &ct);
3144: DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);
3145: for (n = 0; n < Nct; ++n) {
3146: for (r = 0; r < rsize[n]; ++r) {
3147: DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);
3148: DMLabelSetValue(labelNew, pNew, values[val]);
3149: }
3150: }
3151: }
3152: ISRestoreIndices(pointIS, &points);
3153: ISDestroy(&pointIS);
3154: }
3155: ISRestoreIndices(valueIS, &values);
3156: ISDestroy(&valueIS);
3157: return(0);
3158: }
3160: static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm)
3161: {
3162: DM dm = cr->dm;
3163: PetscInt numLabels, l;
3167: DMGetNumLabels(dm, &numLabels);
3168: for (l = 0; l < numLabels; ++l) {
3169: DMLabel label, labelNew;
3170: const char *lname;
3171: PetscBool isDepth, isCellType;
3173: DMGetLabelName(dm, l, &lname);
3174: PetscStrcmp(lname, "depth", &isDepth);
3175: if (isDepth) continue;
3176: PetscStrcmp(lname, "celltype", &isCellType);
3177: if (isCellType) continue;
3178: DMCreateLabel(rdm, lname);
3179: DMGetLabel(dm, lname, &label);
3180: DMGetLabel(rdm, lname, &labelNew);
3181: RefineLabel_Internal(cr, label, labelNew);
3182: }
3183: return(0);
3184: }
3186: /* This will only work for interpolated meshes */
3187: PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined)
3188: {
3189: DM rdm;
3190: PetscInt dim, embedDim, depth;
3191: PetscErrorCode ierr;
3195: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
3196: DMSetType(rdm, DMPLEX);
3197: DMGetDimension(dm, &dim);
3198: DMSetDimension(rdm, dim);
3199: DMGetCoordinateDim(dm, &embedDim);
3200: DMSetCoordinateDim(rdm, embedDim);
3201: /* Calculate number of new points of each depth */
3202: DMPlexGetDepth(dm, &depth);
3203: if (depth >= 0 && dim != depth) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated for regular refinement");
3204: /* Step 1: Set chart */
3205: DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);
3206: /* Step 2: Set cone/support sizes (automatically stratifies) */
3207: DMPlexCellRefinerSetConeSizes(cr, rdm);
3208: /* Step 3: Setup refined DM */
3209: DMSetUp(rdm);
3210: /* Step 4: Set cones and supports (automatically symmetrizes) */
3211: DMPlexCellRefinerSetCones(cr, rdm);
3212: /* Step 5: Create pointSF */
3213: DMPlexCellRefinerCreateSF(cr, rdm);
3214: /* Step 6: Create labels */
3215: DMPlexCellRefinerCreateLabels(cr, rdm);
3216: /* Step 7: Set coordinates */
3217: DMPlexCellRefinerSetCoordinates(cr, rdm);
3218: *dmRefined = rdm;
3219: return(0);
3220: }
3222: /*@
3223: DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
3225: Input Parameter:
3226: . dm - The coarse DM
3228: Output Parameter:
3229: . fpointIS - The IS of all the fine points which exist in the original coarse mesh
3231: Level: developer
3233: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
3234: @*/
3235: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
3236: {
3237: DMPlexCellRefiner cr;
3238: PetscInt *fpoints;
3239: PetscInt pStart, pEnd, p, vStart, vEnd, v;
3240: PetscErrorCode ierr;
3243: DMPlexGetChart(dm, &pStart, &pEnd);
3244: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3245: DMPlexCellRefinerCreate(dm, &cr);
3246: DMPlexCellRefinerSetUp(cr);
3247: PetscMalloc1(pEnd-pStart, &fpoints);
3248: for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1;
3249: for (v = vStart; v < vEnd; ++v) {
3250: PetscInt vNew = -1; /* silent overzelous may be used uninitialized */
3252: DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);
3253: fpoints[v-pStart] = vNew;
3254: }
3255: DMPlexCellRefinerDestroy(&cr);
3256: ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);
3257: return(0);
3258: }
3260: /*@
3261: DMPlexSetRefinementUniform - Set the flag for uniform refinement
3263: Input Parameters:
3264: + dm - The DM
3265: - refinementUniform - The flag for uniform refinement
3267: Level: developer
3269: .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3270: @*/
3271: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
3272: {
3273: DM_Plex *mesh = (DM_Plex*) dm->data;
3277: mesh->refinementUniform = refinementUniform;
3278: return(0);
3279: }
3281: /*@
3282: DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
3284: Input Parameter:
3285: . dm - The DM
3287: Output Parameter:
3288: . refinementUniform - The flag for uniform refinement
3290: Level: developer
3292: .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3293: @*/
3294: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
3295: {
3296: DM_Plex *mesh = (DM_Plex*) dm->data;
3301: *refinementUniform = mesh->refinementUniform;
3302: return(0);
3303: }
3305: /*@
3306: DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
3308: Input Parameters:
3309: + dm - The DM
3310: - refinementLimit - The maximum cell volume in the refined mesh
3312: Level: developer
3314: .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3315: @*/
3316: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
3317: {
3318: DM_Plex *mesh = (DM_Plex*) dm->data;
3322: mesh->refinementLimit = refinementLimit;
3323: return(0);
3324: }
3326: /*@
3327: DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
3329: Input Parameter:
3330: . dm - The DM
3332: Output Parameter:
3333: . refinementLimit - The maximum cell volume in the refined mesh
3335: Level: developer
3337: .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3338: @*/
3339: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
3340: {
3341: DM_Plex *mesh = (DM_Plex*) dm->data;
3346: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
3347: *refinementLimit = mesh->refinementLimit;
3348: return(0);
3349: }
3351: /*@
3352: DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
3354: Input Parameters:
3355: + dm - The DM
3356: - refinementFunc - Function giving the maximum cell volume in the refined mesh
3358: Note: The calling sequence is refinementFunc(coords, limit)
3359: $ coords - Coordinates of the current point, usually a cell centroid
3360: $ limit - The maximum cell volume for a cell containing this point
3362: Level: developer
3364: .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3365: @*/
3366: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
3367: {
3368: DM_Plex *mesh = (DM_Plex*) dm->data;
3372: mesh->refinementFunc = refinementFunc;
3373: return(0);
3374: }
3376: /*@
3377: DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
3379: Input Parameter:
3380: . dm - The DM
3382: Output Parameter:
3383: . refinementFunc - Function giving the maximum cell volume in the refined mesh
3385: Note: The calling sequence is refinementFunc(coords, limit)
3386: $ coords - Coordinates of the current point, usually a cell centroid
3387: $ limit - The maximum cell volume for a cell containing this point
3389: Level: developer
3391: .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3392: @*/
3393: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
3394: {
3395: DM_Plex *mesh = (DM_Plex*) dm->data;
3400: *refinementFunc = mesh->refinementFunc;
3401: return(0);
3402: }
3404: static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm)
3405: {
3406: DM dm = cr->dm;
3407: PetscInt Nf, f, Nds, s;
3411: DMGetNumFields(dm, &Nf);
3412: for (f = 0; f < Nf; ++f) {
3413: DMLabel label, labelNew;
3414: PetscObject obj;
3415: const char *lname;
3417: DMGetField(rdm, f, &label, &obj);
3418: if (!label) continue;
3419: PetscObjectGetName((PetscObject) label, &lname);
3420: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
3421: RefineLabel_Internal(cr, label, labelNew);
3422: DMSetField_Internal(rdm, f, labelNew, obj);
3423: DMLabelDestroy(&labelNew);
3424: }
3425: DMGetNumDS(dm, &Nds);
3426: for (s = 0; s < Nds; ++s) {
3427: DMLabel label, labelNew;
3428: const char *lname;
3430: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
3431: if (!label) continue;
3432: PetscObjectGetName((PetscObject) label, &lname);
3433: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
3434: RefineLabel_Internal(cr, label, labelNew);
3435: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
3436: DMLabelDestroy(&labelNew);
3437: }
3438: return(0);
3439: }
3441: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
3442: {
3443: PetscBool isUniform;
3444: DMPlexCellRefiner cr;
3445: PetscErrorCode ierr;
3448: DMPlexGetRefinementUniform(dm, &isUniform);
3449: DMViewFromOptions(dm, NULL, "-initref_dm_view");
3450: if (isUniform) {
3451: DM cdm, rcdm;
3452: PetscBool localized;
3454: DMPlexCellRefinerCreate(dm, &cr);
3455: DMPlexCellRefinerSetUp(cr);
3456: DMGetCoordinatesLocalized(dm, &localized);
3457: DMPlexRefineUniform(dm, cr, dmRefined);
3458: DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);
3459: DMCopyDisc(dm, *dmRefined);
3460: DMGetCoordinateDM(dm, &cdm);
3461: DMGetCoordinateDM(*dmRefined, &rcdm);
3462: DMCopyDisc(cdm, rcdm);
3463: RefineDiscLabels_Internal(cr, *dmRefined);
3464: DMCopyBoundary(dm, *dmRefined);
3465: DMPlexCellRefinerDestroy(&cr);
3466: } else {
3467: DMPlexRefine_Internal(dm, NULL, dmRefined);
3468: }
3469: return(0);
3470: }
3472: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
3473: {
3474: DM cdm = dm;
3475: PetscInt r;
3476: PetscBool isUniform, localized;
3480: DMPlexGetRefinementUniform(dm, &isUniform);
3481: DMGetCoordinatesLocalized(dm, &localized);
3482: if (isUniform) {
3483: for (r = 0; r < nlevels; ++r) {
3484: DMPlexCellRefiner cr;
3485: DM codm, rcodm;
3487: DMPlexCellRefinerCreate(cdm, &cr);
3488: DMPlexCellRefinerSetUp(cr);
3489: DMPlexRefineUniform(cdm, cr, &dmRefined[r]);
3490: DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);
3491: DMSetRefineLevel(dmRefined[r], cdm->levelup+1);
3492: DMCopyDisc(cdm, dmRefined[r]);
3493: DMGetCoordinateDM(dm, &codm);
3494: DMGetCoordinateDM(dmRefined[r], &rcodm);
3495: DMCopyDisc(codm, rcodm);
3496: RefineDiscLabels_Internal(cr, dmRefined[r]);
3497: DMCopyBoundary(cdm, dmRefined[r]);
3498: DMSetCoarseDM(dmRefined[r], cdm);
3499: DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
3500: cdm = dmRefined[r];
3501: DMPlexCellRefinerDestroy(&cr);
3502: }
3503: } else {
3504: for (r = 0; r < nlevels; ++r) {
3505: DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);
3506: DMCopyDisc(cdm, dmRefined[r]);
3507: DMCopyBoundary(cdm, dmRefined[r]);
3508: if (localized) {DMLocalizeCoordinates(dmRefined[r]);}
3509: DMSetCoarseDM(dmRefined[r], cdm);
3510: cdm = dmRefined[r];
3511: }
3512: }
3513: return(0);
3514: }