1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/petscfeimpl.h> /* For PetscFEInterpolate_Static() */
3 #include <petscsf.h>
4
5 const char * const DMPlexCellRefinerTypes[] = {"Regular", "ToBox", "ToSimplex", "Alfeld2D", "Alfeld3D", "PowellSabin", "BoundaryLayer", "DMPlexCellRefinerTypes", "DM_REFINER_", NULL};
6
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 */
DMPlexCellRefinerGetAffineFaceTransforms_Regular(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nf,PetscReal * v0[],PetscReal * J[],PetscReal * invJ[],PetscReal * detJ[])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
47
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};
57
58 PetscFunctionBegin;
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 PetscFunctionReturn(0);
66 }
67
68 /* Gets the affine map from the original cell to each subcell */
DMPlexCellRefinerGetAffineTransforms_Regular(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nc,PetscReal * v0[],PetscReal * J[],PetscReal * invJ[])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,
89
90 0.5, 0.0,
91 0.0, 0.5,
92
93 0.5, 0.0,
94 0.0, 0.5,
95
96 0.0, -0.5,
97 0.5, 0.5};
98 static PetscReal tri_invJ[] = {2.0, 0.0,
99 0.0, 2.0,
100
101 2.0, 0.0,
102 0.0, 2.0,
103
104 2.0, 0.0,
105 0.0, 2.0,
106
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,
123
124 0.5, 0.0,
125 0.0, 0.5,
126
127 0.5, 0.0,
128 0.0, 0.5,
129
130 0.5, 0.0,
131 0.0, 0.5};
132 static PetscReal quad_invJ[] = {2.0, 0.0,
133 0.0, 2.0,
134
135 2.0, 0.0,
136 0.0, 2.0,
137
138 2.0, 0.0,
139 0.0, 2.0,
140
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,
160
161 0.5, 0.0, 0.0,
162 0.0, 0.5, 0.0,
163 0.0, 0.0, 0.5,
164
165 0.5, 0.0, 0.0,
166 0.0, 0.5, 0.0,
167 0.0, 0.0, 0.5,
168
169 0.5, 0.0, 0.0,
170 0.0, 0.5, 0.0,
171 0.0, 0.0, 0.5,
172
173 0.5, 0.0, 0.0,
174 0.0, 0.5, 0.0,
175 0.0, 0.0, 0.5,
176
177 0.5, 0.0, 0.0,
178 0.0, 0.5, 0.0,
179 0.0, 0.0, 0.5,
180
181 0.5, 0.0, 0.0,
182 0.0, 0.5, 0.0,
183 0.0, 0.0, 0.5,
184
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,
191
192 2.0, 0.0, 0.0,
193 0.0, 2.0, 0.0,
194 0.0, 0.0, 2.0,
195
196 2.0, 0.0, 0.0,
197 0.0, 2.0, 0.0,
198 0.0, 0.0, 2.0,
199
200 2.0, 0.0, 0.0,
201 0.0, 2.0, 0.0,
202 0.0, 0.0, 2.0,
203
204 2.0, 0.0, 0.0,
205 0.0, 2.0, 0.0,
206 0.0, 0.0, 2.0,
207
208 2.0, 0.0, 0.0,
209 0.0, 2.0, 0.0,
210 0.0, 0.0, 2.0,
211
212 2.0, 0.0, 0.0,
213 0.0, 2.0, 0.0,
214 0.0, 0.0, 2.0,
215
216 2.0, 0.0, 0.0,
217 0.0, 2.0, 0.0,
218 0.0, 0.0, 2.0};
219
220 PetscFunctionBegin;
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 PetscFunctionReturn(0);
229 }
230
231 /* Should this be here or in the DualSpace somehow? */
CellRefinerInCellTest_Internal(DMPolytopeType ct,const PetscReal point[],PetscBool * inside)232 PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType ct, const PetscReal point[], PetscBool *inside)
233 {
234 PetscReal sum = 0.0;
235 PetscInt d;
236
237 PetscFunctionBegin;
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 PetscFunctionReturn(0);
257 }
258
259 /* Regular Refinment of Hybrid Meshes
260
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.
264
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.
267
268 We need the group multiplication table for group actions from the dihedral group for each cell type.
269
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 */
274
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
280
281 Output Parameters:
282 . cpnew - The new cone point, taking into account the orientation co
283 */
DMPolytopeMapCell(DMPolytopeType ct,PetscInt co,PetscInt cp,PetscInt * cpnew)284 PETSC_STATIC_INLINE PetscErrorCode DMPolytopeMapCell(DMPolytopeType ct, PetscInt co, PetscInt cp, PetscInt *cpnew)
285 {
286 const PetscInt csize = DMPolytopeTypeGetConeSize(ct);
287
288 PetscFunctionBeginHot;
289 if (ct == DM_POLYTOPE_POINT) {*cpnew = cp;}
290 else {*cpnew = (co < 0 ? -(co+1)-cp+csize : co+cp)%csize;}
291 PetscFunctionReturn(0);
292 }
293
DMPlexCellRefinerGetCellVertices_Regular(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nv,PetscReal * subcellV[])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};
311
312 PetscFunctionBegin;
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 PetscFunctionReturn(0);
322 }
323
DMPlexCellRefinerGetCellVertices_ToBox(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nv,PetscReal * subcellV[])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;
333
334 PetscFunctionBegin;
335 switch (ct) {
336 case DM_POLYTOPE_SEGMENT:
337 case DM_POLYTOPE_QUADRILATERAL:
338 case DM_POLYTOPE_HEXAHEDRON:
339 ierr = DMPlexCellRefinerGetCellVertices_Regular(cr, ct, Nv, subcellV);CHKERRQ(ierr);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 PetscFunctionReturn(0);
345 }
346
347 /*
348 DMPlexCellRefinerGetCellVertices - Get the set of refined vertices lying in the closure of a reference cell of given type
349
350 Input Parameters:
351 + cr - The DMPlexCellRefiner object
352 - ct - The cell type
353
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
357
358 Level: developer
359
360 .seealso: DMPlexCellRefinerGetSubcellVertices()
361 */
DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nv,PetscReal * subcellV[])362 static PetscErrorCode DMPlexCellRefinerGetCellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nv, PetscReal *subcellV[])
363 {
364 PetscErrorCode ierr;
365
366 PetscFunctionBegin;
367 if (!cr->ops->getcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
368 ierr = (*cr->ops->getcellvertices)(cr, ct, Nv, subcellV);CHKERRQ(ierr);
369 PetscFunctionReturn(0);
370 }
371
DMPlexCellRefinerGetSubcellVertices_Regular(DMPlexCellRefiner cr,DMPolytopeType ct,DMPolytopeType rct,PetscInt r,PetscInt * Nv,PetscInt * subcellV[])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};
381
382 PetscFunctionBegin;
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 PetscFunctionReturn(0);
393 }
394
DMPlexCellRefinerGetSubcellVertices_ToBox(DMPlexCellRefiner cr,DMPolytopeType ct,DMPolytopeType rct,PetscInt r,PetscInt * Nv,PetscInt * subcellV[])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;
400
401 PetscFunctionBegin;
402 switch (ct) {
403 case DM_POLYTOPE_SEGMENT:
404 case DM_POLYTOPE_QUADRILATERAL:
405 case DM_POLYTOPE_HEXAHEDRON:
406 ierr = DMPlexCellRefinerGetSubcellVertices_Regular(cr, ct, rct, r, Nv, subcellV);CHKERRQ(ierr);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 PetscFunctionReturn(0);
416 }
417
418 /*
419 DMPlexCellRefinerGetSubcellVertices - Get the set of refined vertices defining a subcell in the reference cell of given type
420
421 Input Parameters:
422 + cr - The DMPlexCellRefiner object
423 . ct - The cell type
424 . rct - The type of subcell
425 - r - The subcell index
426
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()
430
431 Level: developer
432
433 .seealso: DMPlexCellRefinerGetCellVertices()
434 */
DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr,DMPolytopeType ct,DMPolytopeType rct,PetscInt r,PetscInt * Nv,PetscInt * subcellV[])435 static PetscErrorCode DMPlexCellRefinerGetSubcellVertices(DMPlexCellRefiner cr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *Nv, PetscInt *subcellV[])
436 {
437 PetscErrorCode ierr;
438
439 PetscFunctionBegin;
440 if (!cr->ops->getsubcellvertices) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
441 ierr = (*cr->ops->getsubcellvertices)(cr, ct, rct, r, Nv, subcellV);CHKERRQ(ierr);
442 PetscFunctionReturn(0);
443 }
444
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
454
455 Output Parameters:
456 . out - The coordinates of the new vertices
457 */
DMPlexCellRefinerMapCoordinates(DMPlexCellRefiner cr,DMPolytopeType pct,DMPolytopeType ct,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])458 static PetscErrorCode DMPlexCellRefinerMapCoordinates(DMPlexCellRefiner cr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
459 {
460 PetscErrorCode ierr;
461
462 PetscFunctionBeginHot;
463 if (!cr->ops->mapcoords) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
464 ierr = (*cr->ops->mapcoords)(cr, pct, ct, r, Nv, dE, in, out);CHKERRQ(ierr);
465 PetscFunctionReturn(0);
466 }
467
468 /* Computes new vertex as the barycenter */
DMPlexCellRefinerMapCoordinates_Barycenter(DMPlexCellRefiner cr,DMPolytopeType pct,DMPolytopeType ct,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])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;
472
473 PetscFunctionBeginHot;
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 PetscFunctionReturn(0);
479 }
480
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
489
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 */
DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)494 static PetscErrorCode DMPlexCellRefinerMapSubcells(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
495 {
496 PetscErrorCode ierr;
497
498 PetscFunctionBeginHot;
499 if (!cr->ops->mapsubcells) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
500 ierr = (*cr->ops->mapsubcells)(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr);
501 PetscFunctionReturn(0);
502 }
503
504 /*
505 This is the group multiplication table for the dihedral group of the cell.
506 */
ComposeOrientation_Private(PetscInt n,PetscInt o1,PetscInt o2,PetscInt * o)507 static PetscErrorCode ComposeOrientation_Private(PetscInt n, PetscInt o1, PetscInt o2, PetscInt *o)
508 {
509 PetscFunctionBeginHot;
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 PetscFunctionReturn(0);
516 }
517
DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)518 static PetscErrorCode DMPlexCellRefinerMapSubcells_None(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
519 {
520 PetscErrorCode ierr;
521
522 PetscFunctionBeginHot;
523 *rnew = r;
524 ierr = ComposeOrientation_Private(DMPolytopeTypeGetConeSize(ct), po, o, onew);CHKERRQ(ierr);
525 PetscFunctionReturn(0);
526 }
527
DMPlexCellRefinerMapSubcells_Regular(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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};
598
599 PetscFunctionBeginHot;
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 PetscFunctionReturn(0);
652 }
653
DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)654 static PetscErrorCode DMPlexCellRefinerMapSubcells_ToBox(DMPlexCellRefiner cr, DMPolytopeType pct, PetscInt po, DMPolytopeType ct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
655 {
656 PetscErrorCode ierr;
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};
698
699 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerMapSubcells_Regular(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr);
743 }
744 PetscFunctionReturn(0);
745 }
746
DMPlexCellRefinerMapSubcells_ToSimplex(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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 }
751
752 /*@
753 DMPlexCellRefinerRefine - Return a description of the refinement for a given cell type
754
755 Input Parameter:
756 . source - The cell type for a source point
757
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
764
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}
777
778 Level: developer
779
780 .seealso: DMPlexCellRefinerCreate(), DMPlexRefineUniform()
781 @*/
DMPlexCellRefinerRefine(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])782 PetscErrorCode DMPlexCellRefinerRefine(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
783 {
784 PetscErrorCode ierr;
785
786 PetscFunctionBeginHot;
787 if (!cr->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Not for refiner type %s",DMPlexCellRefinerTypes[cr->type]);
788 ierr = (*cr->ops->refine)(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
789 PetscFunctionReturn(0);
790 }
791
DMPlexCellRefinerRefine_None(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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};
842
843 PetscFunctionBegin;
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 PetscFunctionReturn(0);
859 }
860
DMPlexCellRefinerRefine_Regular(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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
1160
1161 2
1162
1163 0 1
1164
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};
1208
1209 PetscFunctionBegin;
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 PetscFunctionReturn(0);
1225 }
1226
DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])1227 static PetscErrorCode DMPlexCellRefinerRefine_ToBox(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1228 {
1229 PetscErrorCode ierr;
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
1387
1388 2
1389
1390 0 1
1391
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
1420
1421 2
1422
1423 0 1
1424
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;
1469
1470 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
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 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1507 }
1508
DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])1509 static PetscErrorCode DMPlexCellRefinerRefine_ToSimplex(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1510 {
1511 PetscErrorCode ierr;
1512
1513 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerRefine_Regular(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
1527 break;
1528 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1529 }
1530 PetscFunctionReturn(0);
1531 }
1532
DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])1533 static PetscErrorCode DMPlexCellRefinerRefine_Alfeld2D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1534 {
1535 PetscErrorCode ierr;
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};
1569
1570 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1588 }
1589
DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])1590 static PetscErrorCode DMPlexCellRefinerRefine_Alfeld3D(DMPlexCellRefiner cr, DMPolytopeType source, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
1591 {
1592 PetscErrorCode ierr;
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};
1641
1642 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerRefine_None(cr, source, Nt, target, size, cone, ornt);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1660 }
1661
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;
1672
DMPlexCellRefinerSetUp_BL(DMPlexCellRefiner cr)1673 static PetscErrorCode DMPlexCellRefinerSetUp_BL(DMPlexCellRefiner cr)
1674 {
1675 PlexRefiner_BL *crbl;
1676 PetscErrorCode ierr;
1677 PetscInt i,n;
1678 PetscReal r;
1679 PetscInt c1,c2,o1,o2;
1680
1681 PetscFunctionBegin;
1682 ierr = PetscNew(&crbl);CHKERRQ(ierr);
1683 cr->data = crbl;
1684 crbl->n = 1; /* 1 split -> 2 new cells */
1685 crbl->r = 1; /* linear progression */
1686
1687 /* TODO: add setfromoptions to the refiners? */
1688 ierr = PetscOptionsGetInt(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_splits", &crbl->n, NULL);CHKERRQ(ierr);
1689 if (crbl->n < 1) SETERRQ1(PetscObjectComm((PetscObject)cr),PETSC_ERR_SUP,"Number of splits %D must be positive",crbl->n);
1690 ierr = PetscOptionsGetReal(((PetscObject) cr->dm)->options,((PetscObject) cr->dm)->prefix, "-dm_plex_refine_boundarylayer_progression", &crbl->r, NULL);CHKERRQ(ierr);
1691 n = crbl->n;
1692 r = crbl->r;
1693
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 ierr = PetscMalloc5(4,&crbl->Nt,4,&crbl->target,4,&crbl->size,4,&crbl->cone,4,&crbl->ornt);CHKERRQ(ierr);
1696
1697 /* progression */
1698 ierr = PetscMalloc1(n,&crbl->h);CHKERRQ(ierr);
1699 if (r > 1) {
1700 PetscReal d = (r-1.)/(PetscPowRealInt(r,n+1)-1.);
1701
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 }
1710
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;
1715
1716 ierr = PetscMalloc4(crbl->Nt[0],&crbl->target[0],crbl->Nt[0],&crbl->size[0],c1,&crbl->cone[0],o1,&crbl->ornt[0]);CHKERRQ(ierr);
1717
1718 crbl->target[0][0] = DM_POLYTOPE_POINT;
1719 crbl->target[0][1] = DM_POLYTOPE_POINT_PRISM_TENSOR;
1720
1721 crbl->size[0][0] = n;
1722 crbl->size[0][1] = n+1;
1723
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;
1748
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;
1755
1756 ierr = PetscMalloc4(crbl->Nt[1],&crbl->target[1],crbl->Nt[1],&crbl->size[1],c1+c2,&crbl->cone[1],o1+o2,&crbl->ornt[1]);CHKERRQ(ierr);
1757
1758 crbl->target[1][0] = DM_POLYTOPE_SEGMENT;
1759 crbl->target[1][1] = DM_POLYTOPE_SEG_PRISM_TENSOR;
1760
1761 crbl->size[1][0] = n;
1762 crbl->size[1][1] = n+1;
1763
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 }
1775
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;
1824
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;
1831
1832 ierr = PetscMalloc4(crbl->Nt[2],&crbl->target[2],crbl->Nt[2],&crbl->size[2],c1+c2,&crbl->cone[2],o1+o2,&crbl->ornt[2]);CHKERRQ(ierr);
1833
1834 crbl->target[2][0] = DM_POLYTOPE_TRIANGLE;
1835 crbl->target[2][1] = DM_POLYTOPE_TRI_PRISM_TENSOR;
1836
1837 crbl->size[2][0] = n;
1838 crbl->size[2][1] = n+1;
1839
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 }
1855
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;
1916
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;
1923
1924 ierr = PetscMalloc4(crbl->Nt[3],&crbl->target[3],crbl->Nt[3],&crbl->size[3],c1+c2,&crbl->cone[3],o1+o2,&crbl->ornt[3]);CHKERRQ(ierr);
1925
1926 crbl->target[3][0] = DM_POLYTOPE_QUADRILATERAL;
1927 crbl->target[3][1] = DM_POLYTOPE_QUAD_PRISM_TENSOR;
1928
1929 crbl->size[3][0] = n;
1930 crbl->size[3][1] = n+1;
1931
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 }
1951
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 PetscFunctionReturn(0);
2025 }
2026
DMPlexCellRefinerDestroy_BL(DMPlexCellRefiner cr)2027 static PetscErrorCode DMPlexCellRefinerDestroy_BL(DMPlexCellRefiner cr)
2028 {
2029 PlexRefiner_BL *crbl = (PlexRefiner_BL *)cr->data;
2030 PetscErrorCode ierr;
2031
2032 PetscFunctionBegin;
2033 ierr = PetscFree4(crbl->target[0],crbl->size[0],crbl->cone[0],crbl->ornt[0]);CHKERRQ(ierr);
2034 ierr = PetscFree4(crbl->target[1],crbl->size[1],crbl->cone[1],crbl->ornt[1]);CHKERRQ(ierr);
2035 ierr = PetscFree4(crbl->target[2],crbl->size[2],crbl->cone[2],crbl->ornt[2]);CHKERRQ(ierr);
2036 ierr = PetscFree4(crbl->target[3],crbl->size[3],crbl->cone[3],crbl->ornt[3]);CHKERRQ(ierr);
2037 ierr = PetscFree5(crbl->Nt,crbl->target,crbl->size,crbl->cone,crbl->ornt);CHKERRQ(ierr);
2038 ierr = PetscFree(crbl->h);CHKERRQ(ierr);
2039 ierr = PetscFree(cr->data);CHKERRQ(ierr);
2040 PetscFunctionReturn(0);
2041 }
2042
DMPlexCellRefinerRefine_BL(DMPlexCellRefiner cr,DMPolytopeType source,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])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;
2047
2048 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerRefine_None(cr,source,Nt,target,size,cone,ornt);CHKERRQ(ierr);
2080 }
2081 PetscFunctionReturn(0);
2082 }
2083
DMPlexCellRefinerMapSubcells_BL(DMPlexCellRefiner cr,DMPolytopeType pct,PetscInt po,DMPolytopeType ct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)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;
2101 PetscErrorCode ierr;
2102
2103 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerMapSubcells_None(cr, pct, po, ct, r, o, rnew, onew);CHKERRQ(ierr);
2129 }
2130 PetscFunctionReturn(0);
2131 }
2132
DMPlexCellRefinerMapCoordinates_BL(DMPlexCellRefiner cr,DMPolytopeType pct,DMPolytopeType ct,PetscInt r,PetscInt Nv,PetscInt dE,const PetscScalar in[],PetscScalar out[])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;
2138
2139 PetscFunctionBeginHot;
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 ierr = DMPlexCellRefinerMapCoordinates_Barycenter(cr,pct,ct,r,Nv,dE,in,out);CHKERRQ(ierr);
2149 }
2150 PetscFunctionReturn(0);
2151 }
2152
CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr,PetscInt ctOrder[],PetscInt ctStart[],PetscInt ** offset)2153 static PetscErrorCode CellRefinerCreateOffset_Internal(DMPlexCellRefiner cr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
2154 {
2155 PetscInt c, cN, *off;
2156 PetscErrorCode ierr;
2157
2158 PetscFunctionBegin;
2159 ierr = PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);CHKERRQ(ierr);
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;
2167
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];
2173
2174 ierr = DMPlexCellRefinerRefine(cr, ict, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
2186 }
2187
DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr,const PetscInt ctStart[],const PetscInt ctStartNew[])2188 static PetscErrorCode DMPlexCellRefinerSetStarts(DMPlexCellRefiner cr, const PetscInt ctStart[], const PetscInt ctStartNew[])
2189 {
2190 const PetscInt ctSize = DM_NUM_POLYTOPES+1;
2191 PetscErrorCode ierr;
2192
2193 PetscFunctionBegin;
2194 if (cr->setupcalled) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_ARG_WRONGSTATE, "Must call this function before DMPlexCellRefinerSetUp()");
2195 ierr = PetscCalloc2(ctSize, &cr->ctStart, ctSize, &cr->ctStartNew);CHKERRQ(ierr);
2196 ierr = PetscArraycpy(cr->ctStart, ctStart, ctSize);CHKERRQ(ierr);
2197 ierr = PetscArraycpy(cr->ctStartNew, ctStartNew, ctSize);CHKERRQ(ierr);
2198 PetscFunctionReturn(0);
2199 }
2200
2201 /* Construct cell type order since we must loop over cell types in depth order */
DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell,PetscInt * ctOrder[],PetscInt * ctOrderInv[])2202 PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
2203 {
2204 PetscInt *ctO, *ctOInv;
2205 PetscInt c, d, off = 0;
2206 PetscErrorCode ierr;
2207
2208 PetscFunctionBegin;
2209 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
2239 }
2240
DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)2241 PetscErrorCode DMPlexCellRefinerSetUp(DMPlexCellRefiner cr)
2242 {
2243 DM dm = cr->dm;
2244 DMPolytopeType ctCell;
2245 PetscInt pStart, pEnd, p, c;
2246 PetscErrorCode ierr;
2247
2248 PetscFunctionBegin;
2249 PetscValidHeader(cr, 1);
2250 if (cr->setupcalled) PetscFunctionReturn(0);
2251 if (cr->ops->setup) {
2252 ierr = (*cr->ops->setup)(cr);CHKERRQ(ierr);
2253 }
2254 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
2255 if (pEnd > pStart) {
2256 ierr = DMPlexGetCellType(dm, 0, &ctCell);CHKERRQ(ierr);
2257 } else {
2258 PetscInt dim;
2259 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
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 ierr = DMPlexCreateCellTypeOrder_Internal(ctCell, &cr->ctOrder, &cr->ctOrderInv);CHKERRQ(ierr);
2269 /* Construct sizes and offsets for each cell type */
2270 if (!cr->ctStart) {
2271 PetscInt *ctS, *ctSN, *ctC, *ctCN;
2272
2273 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);CHKERRQ(ierr);
2274 ierr = PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);CHKERRQ(ierr);
2275 for (p = pStart; p < pEnd; ++p) {
2276 DMPolytopeType ct;
2277 DMPolytopeType *rct;
2278 PetscInt *rsize, *cone, *ornt;
2279 PetscInt Nct, n;
2280
2281 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
2282 if ((PetscInt) ct < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
2283 ++ctC[ct];
2284 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr);
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];
2290
2291 ctS[ctn] = ctS[ct] + ctC[ct];
2292 ctSN[ctn] = ctSN[ct] + ctCN[ct];
2293 }
2294 ierr = PetscFree2(ctC, ctCN);CHKERRQ(ierr);
2295 cr->ctStart = ctS;
2296 cr->ctStartNew = ctSN;
2297 }
2298 ierr = CellRefinerCreateOffset_Internal(cr, cr->ctOrder, cr->ctStart, &cr->offset);CHKERRQ(ierr);
2299 cr->setupcalled = PETSC_TRUE;
2300 PetscFunctionReturn(0);
2301 }
2302
DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr,PetscViewer v)2303 static PetscErrorCode DMPlexCellRefinerView_Ascii(DMPlexCellRefiner cr, PetscViewer v)
2304 {
2305 PetscErrorCode ierr;
2306
2307 PetscFunctionBegin;
2308 ierr = PetscViewerASCIIPrintf(v, "Cell Refiner: %s\n", DMPlexCellRefinerTypes[cr->type]);CHKERRQ(ierr);
2309 PetscFunctionReturn(0);
2310 }
2311
2312 /*
2313 DMPlexCellRefinerView - Views a DMPlexCellRefiner object
2314
2315 Collective on cr
2316
2317 Input Parameters:
2318 + cr - The DMPlexCellRefiner object
2319 - viewer - The PetscViewer object
2320
2321 Level: beginner
2322
2323 .seealso: DMPlexCellRefinerCreate()
2324 */
DMPlexCellRefinerView(DMPlexCellRefiner cr,PetscViewer viewer)2325 static PetscErrorCode DMPlexCellRefinerView(DMPlexCellRefiner cr, PetscViewer viewer)
2326 {
2327 PetscBool iascii;
2328 PetscErrorCode ierr;
2329
2330 PetscFunctionBegin;
2331 PetscValidHeader(cr, 1);
2332 if (viewer) PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
2333 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) cr), &viewer);CHKERRQ(ierr);}
2334 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
2335 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
2336 if (iascii) {ierr = DMPlexCellRefinerView_Ascii(cr, viewer);CHKERRQ(ierr);}
2337 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
2338 PetscFunctionReturn(0);
2339 }
2340
DMPlexCellRefinerDestroy(DMPlexCellRefiner * cr)2341 PetscErrorCode DMPlexCellRefinerDestroy(DMPlexCellRefiner *cr)
2342 {
2343 PetscInt c;
2344 PetscErrorCode ierr;
2345
2346 PetscFunctionBegin;
2347 if (!*cr) PetscFunctionReturn(0);
2348 PetscValidHeaderSpecific(*cr, DM_CLASSID, 1);
2349 if ((*cr)->ops->destroy) {
2350 ierr = ((*cr)->ops->destroy)(*cr);CHKERRQ(ierr);
2351 }
2352 ierr = PetscObjectDereference((PetscObject) (*cr)->dm);CHKERRQ(ierr);
2353 ierr = PetscFree2((*cr)->ctOrder, (*cr)->ctOrderInv);CHKERRQ(ierr);
2354 ierr = PetscFree2((*cr)->ctStart, (*cr)->ctStartNew);CHKERRQ(ierr);
2355 ierr = PetscFree((*cr)->offset);CHKERRQ(ierr);
2356 for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
2357 ierr = PetscFEDestroy(&(*cr)->coordFE[c]);CHKERRQ(ierr);
2358 ierr = PetscFEGeomDestroy(&(*cr)->refGeom[c]);CHKERRQ(ierr);
2359 }
2360 ierr = PetscFree2((*cr)->coordFE, (*cr)->refGeom);CHKERRQ(ierr);
2361 ierr = PetscHeaderDestroy(cr);CHKERRQ(ierr);
2362 PetscFunctionReturn(0);
2363 }
2364
DMPlexCellRefinerCreate(DM dm,DMPlexCellRefiner * cr)2365 PetscErrorCode DMPlexCellRefinerCreate(DM dm, DMPlexCellRefiner *cr)
2366 {
2367 DMPlexCellRefiner tmp;
2368 PetscErrorCode ierr;
2369
2370 PetscFunctionBegin;
2371 PetscValidPointer(cr, 2);
2372 *cr = NULL;
2373 ierr = PetscHeaderCreate(tmp, DM_CLASSID, "DMPlexCellRefiner", "Cell Refiner", "DMPlexCellRefiner", PETSC_COMM_SELF, DMPlexCellRefinerDestroy, DMPlexCellRefinerView);CHKERRQ(ierr);
2374 tmp->setupcalled = PETSC_FALSE;
2375
2376 tmp->dm = dm;
2377 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
2378 ierr = DMPlexGetCellRefinerType(dm, &tmp->type);CHKERRQ(ierr);
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 ierr = PetscCalloc2(DM_NUM_POLYTOPES, &tmp->coordFE, DM_NUM_POLYTOPES, &tmp->refGeom);CHKERRQ(ierr);
2421 *cr = tmp;
2422 PetscFunctionReturn(0);
2423 }
2424
2425 /*@
2426 DMPlexCellRefinerGetAffineTransforms - Gets the affine map from the reference cell to each subcell
2427
2428 Input Parameters:
2429 + cr - The DMPlexCellRefiner object
2430 - ct - The cell type
2431
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
2437
2438 Level: developer
2439
2440 .seealso: DMPlexCellRefinerGetAffineFaceTransforms(), Create()
2441 @*/
DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nc,PetscReal * v0[],PetscReal * J[],PetscReal * invJ[])2442 PetscErrorCode DMPlexCellRefinerGetAffineTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nc, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[])
2443 {
2444 PetscErrorCode ierr;
2445
2446 PetscFunctionBegin;
2447 if (!cr->ops->getaffinetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine transforms from this refiner");
2448 ierr = (*cr->ops->getaffinetransforms)(cr, ct, Nc, v0, J, invJ);CHKERRQ(ierr);
2449 PetscFunctionReturn(0);
2450 }
2451
2452 /*@
2453 DMPlexCellRefinerGetAffineFaceTransforms - Gets the affine map from the reference face cell to each face in the given cell
2454
2455 Input Parameters:
2456 + cr - The DMPlexCellRefiner object
2457 - ct - The cell type
2458
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
2465
2466 Note: The Jacobian and inverse Jacboian will be rectangular, and the inverse is really a generalized inverse.
2467
2468 Level: developer
2469
2470 .seealso: DMPlexCellRefinerGetAffineTransforms(), Create()
2471 @*/
DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr,DMPolytopeType ct,PetscInt * Nf,PetscReal * v0[],PetscReal * J[],PetscReal * invJ[],PetscReal * detJ[])2472 PetscErrorCode DMPlexCellRefinerGetAffineFaceTransforms(DMPlexCellRefiner cr, DMPolytopeType ct, PetscInt *Nf, PetscReal *v0[], PetscReal *J[], PetscReal *invJ[], PetscReal *detJ[])
2473 {
2474 PetscErrorCode ierr;
2475
2476 PetscFunctionBegin;
2477 if (!cr->ops->getaffinefacetransforms) SETERRQ(PetscObjectComm((PetscObject) cr), PETSC_ERR_SUP, "No support for affine face transforms from this refiner");
2478 ierr = (*cr->ops->getaffinefacetransforms)(cr, ct, Nf, v0, J, invJ, detJ);CHKERRQ(ierr);
2479 PetscFunctionReturn(0);
2480 }
2481
2482 /* Numbering regularly refined meshes
2483
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.
2487
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.
2492
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 */
DMPlexCellRefinerGetNewPoint(DMPlexCellRefiner cr,DMPolytopeType ct,DMPolytopeType ctNew,PetscInt p,PetscInt r,PetscInt * pNew)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;
2506 PetscErrorCode ierr;
2507
2508 PetscFunctionBeginHot;
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]);
2511
2512 newp += off;
2513 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &cone, &ornt);CHKERRQ(ierr);
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 }
2522
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 PetscFunctionReturn(0);
2526 }
2527
DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr,DM rdm)2528 static PetscErrorCode DMPlexCellRefinerSetConeSizes(DMPlexCellRefiner cr, DM rdm)
2529 {
2530 DM dm = cr->dm;
2531 PetscInt pStart, pEnd, p, pNew;
2532 PetscErrorCode ierr;
2533
2534 PetscFunctionBegin;
2535 /* Must create the celltype label here so that we do not automatically try to compute the types */
2536 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr);
2537 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
2538 for (p = pStart; p < pEnd; ++p) {
2539 DMPolytopeType ct;
2540 DMPolytopeType *rct;
2541 PetscInt *rsize, *rcone, *rornt;
2542 PetscInt Nct, n, r;
2543
2544 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
2545 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
2546 for (n = 0; n < Nct; ++n) {
2547 for (r = 0; r < rsize[n]; ++r) {
2548 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr);
2549 ierr = DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));CHKERRQ(ierr);
2550 ierr = DMPlexSetCellType(rdm, pNew, rct[n]);CHKERRQ(ierr);
2551 }
2552 }
2553 }
2554 {
2555 DMLabel ctLabel;
2556 DM_Plex *plex = (DM_Plex *) rdm->data;
2557
2558 ierr = DMPlexGetCellTypeLabel(rdm, &ctLabel);CHKERRQ(ierr);
2559 ierr = PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);CHKERRQ(ierr);
2560 }
2561 PetscFunctionReturn(0);
2562 }
2563
DMPlexCellRefinerSetCones(DMPlexCellRefiner cr,DM rdm)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;
2570 PetscErrorCode ierr;
2571
2572 PetscFunctionBegin;
2573 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
2574 ierr = PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);CHKERRQ(ierr);
2575 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
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 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
2583 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
2584 ierr = DMPlexGetConeOrientation(dm, p, &ornt);CHKERRQ(ierr);
2585 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
2586 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
2587 const DMPolytopeType ctNew = rct[n];
2588 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
2589
2590 for (r = 0; r < rsize[n]; ++r) {
2591 /* pNew is a subcell produced by subdividing p */
2592 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr);
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;
2604
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;
2609
2610 ierr = DMPolytopeMapCell(pct, po, rcone[coff++], &pcp);CHKERRQ(ierr);
2611 ppp = pp;
2612 pp = pcone[pcp];
2613 ierr = DMPlexGetCellType(dm, pp, &pct);CHKERRQ(ierr);
2614 ierr = DMPlexGetCone(dm, pp, &pcone);CHKERRQ(ierr);
2615 ierr = DMPlexGetConeOrientation(dm, ppp, &ppornt);CHKERRQ(ierr);
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 ierr = DMPlexCellRefinerMapSubcells(cr, pct, po, ft, pr, fo, &pr, &fo);CHKERRQ(ierr);
2629 ierr = DMPlexCellRefinerGetNewPoint(cr, pct, ft, pp, pr, &coneNew[c]);CHKERRQ(ierr);
2630 orntNew[c] = fo;
2631 }
2632 ierr = DMPlexSetCone(rdm, pNew, coneNew);CHKERRQ(ierr);
2633 ierr = DMPlexSetConeOrientation(rdm, pNew, orntNew);CHKERRQ(ierr);
2634 }
2635 }
2636 }
2637 ierr = PetscFree2(coneNew, orntNew);CHKERRQ(ierr);
2638 ierr = DMPlexSymmetrize(rdm);CHKERRQ(ierr);
2639 ierr = DMPlexStratify(rdm);CHKERRQ(ierr);
2640 PetscFunctionReturn(0);
2641 }
2642
DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr,DMPolytopeType ct,PetscFE * fe)2643 static PetscErrorCode DMPlexCellRefinerGetCoordinateFE(DMPlexCellRefiner cr, DMPolytopeType ct, PetscFE *fe)
2644 {
2645 PetscErrorCode ierr;
2646
2647 PetscFunctionBegin;
2648 if (!cr->coordFE[ct]) {
2649 PetscInt dim, cdim;
2650 PetscBool isSimplex;
2651
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 ierr = DMGetCoordinateDim(cr->dm, &cdim);CHKERRQ(ierr);
2661 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &cr->coordFE[ct]);CHKERRQ(ierr);
2662 {
2663 PetscDualSpace dsp;
2664 PetscQuadrature quad;
2665 DM K;
2666 PetscFEGeom *cg;
2667 PetscReal *Xq, *xq, *wq;
2668 PetscInt Nq, q;
2669
2670 ierr = DMPlexCellRefinerGetCellVertices(cr, ct, &Nq, &Xq);CHKERRQ(ierr);
2671 ierr = PetscMalloc1(Nq*cdim, &xq);CHKERRQ(ierr);
2672 for (q = 0; q < Nq*cdim; ++q) xq[q] = Xq[q];
2673 ierr = PetscMalloc1(Nq, &wq);CHKERRQ(ierr);
2674 for (q = 0; q < Nq; ++q) wq[q] = 1.0;
2675 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &quad);CHKERRQ(ierr);
2676 ierr = PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);CHKERRQ(ierr);
2677 ierr = PetscFESetQuadrature(cr->coordFE[ct], quad);CHKERRQ(ierr);
2678
2679 ierr = PetscFEGetDualSpace(cr->coordFE[ct], &dsp);CHKERRQ(ierr);
2680 ierr = PetscDualSpaceGetDM(dsp, &K);CHKERRQ(ierr);
2681 ierr = PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &cr->refGeom[ct]);CHKERRQ(ierr);
2682 cg = cr->refGeom[ct];
2683 ierr = DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);CHKERRQ(ierr);
2684 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr);
2685 }
2686 }
2687 *fe = cr->coordFE[ct];
2688 PetscFunctionReturn(0);
2689 }
2690
2691 /*
2692 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
2693
2694 Not collective
2695
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
2702
2703 Output Parameter:
2704 . xr - The localized coordinates for the produced cell
2705
2706 Level: developer
2707
2708 .seealso: DMPlexCellRefinerSetCoordinates()
2709 */
DMPlexCellRefinerMapLocalizedCoordinates(DMPlexCellRefiner cr,DMPolytopeType ct,DMPolytopeType rct,PetscInt r,const PetscScalar x[],PetscScalar xr[])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;
2714 PetscErrorCode ierr;
2715
2716 PetscFunctionBegin;
2717 ierr = DMPlexCellRefinerGetCoordinateFE(cr, ct, &fe);CHKERRQ(ierr);
2718 ierr = DMPlexCellRefinerGetSubcellVertices(cr, ct, rct, r, &Nv, &subcellV);CHKERRQ(ierr);
2719 ierr = PetscFEGetNumComponents(fe, &cdim);CHKERRQ(ierr);
2720 for (v = 0; v < Nv; ++v) {
2721 ierr = PetscFEInterpolate_Static(fe, x, cr->refGeom[ct], subcellV[v], &xr[v*cdim]);CHKERRQ(ierr);
2722 }
2723 PetscFunctionReturn(0);
2724 }
2725
DMPlexCellRefinerSetCoordinates(DMPlexCellRefiner cr,DM rdm)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;
2738
2739 PetscFunctionBegin;
2740 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
2741 ierr = DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);CHKERRQ(ierr);
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 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
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 }
2752
2753 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
2754 ierr = PetscSectionGetFieldComponents(coordSection, 0, &dE);CHKERRQ(ierr);
2755 if (maxCell) {
2756 PetscReal maxCellNew[3];
2757
2758 for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
2759 ierr = DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);CHKERRQ(ierr);
2760 } else {
2761 ierr = DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);CHKERRQ(ierr);
2762 }
2763 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);CHKERRQ(ierr);
2764 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr);
2765 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dE);CHKERRQ(ierr);
2766 ierr = DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);CHKERRQ(ierr);
2767 if (localizeCells) {ierr = PetscSectionSetChart(coordSectionNew, 0, vEndNew);CHKERRQ(ierr);}
2768 else {ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);CHKERRQ(ierr);}
2769
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 ierr = PetscSectionSetDof(coordSectionNew, v, dE);CHKERRQ(ierr);
2775 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);CHKERRQ(ierr);
2776 }
2777 if (localizeCells) {
2778 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
2779 for (c = cStart; c < cEnd; ++c) {
2780 PetscInt dof;
2781
2782 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr);
2783 if (dof) {
2784 DMPolytopeType ct;
2785 DMPolytopeType *rct;
2786 PetscInt *rsize, *rcone, *rornt;
2787 PetscInt dim, cNew, Nct, n, r;
2788
2789 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr);
2790 dim = DMPolytopeTypeGetDim(ct);
2791 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
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;
2798
2799 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], c, r, &cNew);CHKERRQ(ierr);
2800 ierr = DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr);
2801 for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
2802 ierr = DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr);
2803 ierr = PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);CHKERRQ(ierr);
2804 ierr = PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);CHKERRQ(ierr);
2805 }
2806 }
2807 }
2808 }
2809 }
2810 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr);
2811 ierr = DMViewFromOptions(dm, NULL, "-coarse_dm_view");CHKERRQ(ierr);
2812 ierr = DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);CHKERRQ(ierr);
2813 {
2814 VecType vtype;
2815 PetscInt coordSizeNew, bs;
2816 const char *name;
2817
2818 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
2819 ierr = VecCreate(PETSC_COMM_SELF, &coordsLocalNew);CHKERRQ(ierr);
2820 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr);
2821 ierr = VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr);
2822 ierr = PetscObjectGetName((PetscObject) coordsLocal, &name);CHKERRQ(ierr);
2823 ierr = PetscObjectSetName((PetscObject) coordsLocalNew, name);CHKERRQ(ierr);
2824 ierr = VecGetBlockSize(coordsLocal, &bs);CHKERRQ(ierr);
2825 ierr = VecSetBlockSize(coordsLocalNew, bs);CHKERRQ(ierr);
2826 ierr = VecGetType(coordsLocal, &vtype);CHKERRQ(ierr);
2827 ierr = VecSetType(coordsLocalNew, vtype);CHKERRQ(ierr);
2828 }
2829 ierr = VecGetArrayRead(coordsLocal, &coords);CHKERRQ(ierr);
2830 ierr = VecGetArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr);
2831 ierr = PetscSectionGetChart(coordSection, &ocStart, &ocEnd);CHKERRQ(ierr);
2832 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
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;
2840
2841 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
2842 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
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 ierr = PetscSectionGetDof(coordSection, p, &dof);CHKERRQ(ierr);
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;
2855
2856 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr);
2857
2858 icoords = pcoords;
2859 Nv = Nc/dE;
2860 if (ct != DM_POLYTOPE_POINT) {
2861 if (localizeVertices) {
2862 PetscScalar anchor[3];
2863
2864 for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
2865 if (!isLocalized) {
2866 for (v = 0; v < Nv; ++v) {ierr = DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);CHKERRQ(ierr);}
2867 } else {
2868 Nv = Nc/(2*dE);
2869 icoords = pcoords + Nv*dE;
2870 for (v = Nv; v < Nv*2; ++v) {ierr = DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);CHKERRQ(ierr);}
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;
2879
2880 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &vNew);CHKERRQ(ierr);
2881 ierr = PetscSectionGetOffset(coordSectionNew, vNew, &off);CHKERRQ(ierr);
2882 ierr = DMPlexCellRefinerMapCoordinates(cr, ct, rct[n], r, Nv, dE, icoords, vcoords);CHKERRQ(ierr);
2883 ierr = DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);CHKERRQ(ierr);
2884 }
2885 }
2886 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr);
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;
2896
2897 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
2898 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
2899 if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
2900 PetscInt dof;
2901 ierr = PetscSectionGetDof(coordSection, p, &dof);CHKERRQ(ierr);
2902 if (dof) isLocalized = PETSC_TRUE;
2903 }
2904 if (isLocalized) {
2905 const PetscScalar *pcoords;
2906
2907 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr);
2908 for (n = 0; n < Nct; ++n) {
2909 const PetscInt Nr = rsize[n];
2910
2911 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
2912 for (r = 0; r < Nr; ++r) {
2913 PetscInt pNew, offNew;
2914
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 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr);
2919 ierr = PetscSectionGetOffset(coordSectionNew, pNew, &offNew);CHKERRQ(ierr);
2920 ierr = DMPlexCellRefinerMapLocalizedCoordinates(cr, ct, rct[n], r, pcoords, &coordsNew[offNew]);CHKERRQ(ierr);
2921 }
2922 }
2923 }
2924 }
2925 ierr = VecRestoreArrayRead(coordsLocal, &coords);CHKERRQ(ierr);
2926 ierr = VecRestoreArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr);
2927 ierr = DMSetCoordinatesLocal(rdm, coordsLocalNew);CHKERRQ(ierr);
2928 /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
2929 ierr = VecDestroy(&coordsLocalNew);CHKERRQ(ierr);
2930 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr);
2931 if (!localizeCells) {ierr = DMLocalizeCoordinates(rdm);CHKERRQ(ierr);}
2932 PetscFunctionReturn(0);
2933 }
2934
2935 /*@
2936 DMPlexCreateProcessSF - Create an SF which just has process connectivity
2937
2938 Collective on dm
2939
2940 Input Parameters:
2941 + dm - The DM
2942 - sfPoint - The PetscSF which encodes point connectivity
2943
2944 Output Parameters:
2945 + processRanks - A list of process neighbors, or NULL
2946 - sfProcess - An SF encoding the process connectivity, or NULL
2947
2948 Level: developer
2949
2950 .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF()
2951 @*/
DMPlexCreateProcessSF(DM dm,PetscSF sfPoint,IS * processRanks,PetscSF * sfProcess)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;
2962
2963 PetscFunctionBegin;
2964 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2965 PetscValidHeaderSpecific(sfPoint, PETSCSF_CLASSID, 2);
2966 if (processRanks) {PetscValidPointer(processRanks, 3);}
2967 if (sfProcess) {PetscValidPointer(sfProcess, 4);}
2968 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
2969 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr);
2970 ierr = PetscMalloc1(numLeaves, &ranks);CHKERRQ(ierr);
2971 for (l = 0; l < numLeaves; ++l) {
2972 ranks[l] = remotePoints[l].rank;
2973 }
2974 ierr = PetscSortRemoveDupsInt(&numLeaves, ranks);CHKERRQ(ierr);
2975 ierr = PetscMalloc1(numLeaves, &ranksNew);CHKERRQ(ierr);
2976 ierr = PetscMalloc1(numLeaves, &localPointsNew);CHKERRQ(ierr);
2977 ierr = PetscMalloc1(numLeaves, &remotePointsNew);CHKERRQ(ierr);
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 ierr = PetscFree(ranks);CHKERRQ(ierr);
2985 if (processRanks) {ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks);CHKERRQ(ierr);}
2986 else {ierr = PetscFree(ranksNew);CHKERRQ(ierr);}
2987 if (sfProcess) {
2988 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess);CHKERRQ(ierr);
2989 ierr = PetscObjectSetName((PetscObject) *sfProcess, "Process SF");CHKERRQ(ierr);
2990 ierr = PetscSFSetFromOptions(*sfProcess);CHKERRQ(ierr);
2991 ierr = PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr);
2992 }
2993 PetscFunctionReturn(0);
2994 }
2995
DMPlexCellRefinerCreateSF(DMPlexCellRefiner cr,DM rdm)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;
3011
3012 PetscFunctionBegin;
3013 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr);
3014 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
3015 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr);
3016 /* Calculate size of new SF */
3017 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr);
3018 if (numRoots < 0) PetscFunctionReturn(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;
3025
3026 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
3027 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
3028 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
3029 }
3030 /* Communicate ctStart and cStartNew for each remote rank */
3031 ierr = DMPlexCreateProcessSF(dm, sf, &processRanks, &sfProcess);CHKERRQ(ierr);
3032 ierr = ISGetLocalSize(processRanks, &numNeighbors);CHKERRQ(ierr);
3033 ierr = PetscMalloc2(ctSize*numNeighbors, &ctStartRem, ctSize*numNeighbors, &ctStartNewRem);CHKERRQ(ierr);
3034 ierr = MPI_Type_contiguous(ctSize, MPIU_INT, &ctType);CHKERRQ(ierr);
3035 ierr = MPI_Type_commit(&ctType);CHKERRQ(ierr);
3036 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr);
3037 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStart, ctStartRem);CHKERRQ(ierr);
3038 ierr = PetscSFBcastBegin(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr);
3039 ierr = PetscSFBcastEnd(sfProcess, ctType, cr->ctStartNew, ctStartNewRem);CHKERRQ(ierr);
3040 ierr = MPI_Type_free(&ctType);CHKERRQ(ierr);
3041 ierr = PetscSFDestroy(&sfProcess);CHKERRQ(ierr);
3042 ierr = PetscMalloc1(numNeighbors, &crRem);CHKERRQ(ierr);
3043 for (n = 0; n < numNeighbors; ++n) {
3044 ierr = DMPlexCellRefinerCreate(dm, &crRem[n]);CHKERRQ(ierr);
3045 ierr = DMPlexCellRefinerSetStarts(crRem[n], &ctStartRem[n*ctSize], &ctStartNewRem[n*ctSize]);
3046 ierr = DMPlexCellRefinerSetUp(crRem[n]);CHKERRQ(ierr);
3047 }
3048 ierr = PetscFree2(ctStartRem, ctStartNewRem);CHKERRQ(ierr);
3049 /* Calculate new point SF */
3050 ierr = PetscMalloc1(numLeavesNew, &localPointsNew);CHKERRQ(ierr);
3051 ierr = PetscMalloc1(numLeavesNew, &remotePointsNew);CHKERRQ(ierr);
3052 ierr = ISGetIndices(processRanks, &neighbors);CHKERRQ(ierr);
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;
3061
3062 ierr = PetscFindInt(rankRem, numNeighbors, neighbors, &neighbor);CHKERRQ(ierr);
3063 if (neighbor < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not locate remote rank %D", rankRem);
3064 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr);
3065 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
3066 for (n = 0; n < Nct; ++n) {
3067 for (r = 0; r < rsize[n]; ++r) {
3068 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr);
3069 ierr = DMPlexCellRefinerGetNewPoint(crRem[neighbor], ct, rct[n], pRem, r, &pNewRem);CHKERRQ(ierr);
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) {ierr = DMPlexCellRefinerDestroy(&crRem[n]);CHKERRQ(ierr);}
3078 ierr = PetscFree(crRem);CHKERRQ(ierr);
3079 if (m != numLeavesNew) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of leaf point %D should be %D", m, numLeavesNew);
3080 ierr = ISRestoreIndices(processRanks, &neighbors);CHKERRQ(ierr);
3081 ierr = ISDestroy(&processRanks);CHKERRQ(ierr);
3082 {
3083 PetscSFNode *rp, *rtmp;
3084 PetscInt *lp, *idx, *ltmp, i;
3085
3086 /* SF needs sorted leaves to correct calculate Gather */
3087 ierr = PetscMalloc1(numLeavesNew, &idx);CHKERRQ(ierr);
3088 ierr = PetscMalloc1(numLeavesNew, &lp);CHKERRQ(ierr);
3089 ierr = PetscMalloc1(numLeavesNew, &rp);CHKERRQ(ierr);
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 ierr = PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);CHKERRQ(ierr);
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 ierr = PetscFree(idx);CHKERRQ(ierr);
3104 ierr = PetscFree(ltmp);CHKERRQ(ierr);
3105 ierr = PetscFree(rtmp);CHKERRQ(ierr);
3106 }
3107 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr);
3108 PetscFunctionReturn(0);
3109 }
3110
RefineLabel_Internal(DMPlexCellRefiner cr,DMLabel label,DMLabel labelNew)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;
3118
3119 PetscFunctionBegin;
3120 ierr = DMLabelGetDefaultValue(label, &defVal);CHKERRQ(ierr);
3121 ierr = DMLabelSetDefaultValue(labelNew, defVal);CHKERRQ(ierr);
3122 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr);
3123 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr);
3124 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr);
3125 for (val = 0; val < Nv; ++val) {
3126 IS pointIS;
3127 const PetscInt *points;
3128 PetscInt numPoints, p;
3129
3130 /* Ensure refined label is created with same number of strata as
3131 * original (even if no entries here). */
3132 ierr = DMLabelAddStratum(labelNew, values[val]);CHKERRQ(ierr);
3133 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr);
3134 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr);
3135 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
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;
3142
3143 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr);
3144 ierr = DMPlexCellRefinerRefine(cr, ct, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr);
3145 for (n = 0; n < Nct; ++n) {
3146 for (r = 0; r < rsize[n]; ++r) {
3147 ierr = DMPlexCellRefinerGetNewPoint(cr, ct, rct[n], point, r, &pNew);CHKERRQ(ierr);
3148 ierr = DMLabelSetValue(labelNew, pNew, values[val]);CHKERRQ(ierr);
3149 }
3150 }
3151 }
3152 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
3153 ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
3154 }
3155 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr);
3156 ierr = ISDestroy(&valueIS);CHKERRQ(ierr);
3157 PetscFunctionReturn(0);
3158 }
3159
DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr,DM rdm)3160 static PetscErrorCode DMPlexCellRefinerCreateLabels(DMPlexCellRefiner cr, DM rdm)
3161 {
3162 DM dm = cr->dm;
3163 PetscInt numLabels, l;
3164 PetscErrorCode ierr;
3165
3166 PetscFunctionBegin;
3167 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr);
3168 for (l = 0; l < numLabels; ++l) {
3169 DMLabel label, labelNew;
3170 const char *lname;
3171 PetscBool isDepth, isCellType;
3172
3173 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr);
3174 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr);
3175 if (isDepth) continue;
3176 ierr = PetscStrcmp(lname, "celltype", &isCellType);CHKERRQ(ierr);
3177 if (isCellType) continue;
3178 ierr = DMCreateLabel(rdm, lname);CHKERRQ(ierr);
3179 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr);
3180 ierr = DMGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr);
3181 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr);
3182 }
3183 PetscFunctionReturn(0);
3184 }
3185
3186 /* This will only work for interpolated meshes */
DMPlexRefineUniform(DM dm,DMPlexCellRefiner cr,DM * dmRefined)3187 PetscErrorCode DMPlexRefineUniform(DM dm, DMPlexCellRefiner cr, DM *dmRefined)
3188 {
3189 DM rdm;
3190 PetscInt dim, embedDim, depth;
3191 PetscErrorCode ierr;
3192
3193 PetscFunctionBegin;
3194 PetscValidHeader(cr, 1);
3195 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr);
3196 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr);
3197 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
3198 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr);
3199 ierr = DMGetCoordinateDim(dm, &embedDim);CHKERRQ(ierr);
3200 ierr = DMSetCoordinateDim(rdm, embedDim);CHKERRQ(ierr);
3201 /* Calculate number of new points of each depth */
3202 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
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 ierr = DMPlexSetChart(rdm, 0, cr->ctStartNew[cr->ctOrder[DM_NUM_POLYTOPES]]);CHKERRQ(ierr);
3206 /* Step 2: Set cone/support sizes (automatically stratifies) */
3207 ierr = DMPlexCellRefinerSetConeSizes(cr, rdm);CHKERRQ(ierr);
3208 /* Step 3: Setup refined DM */
3209 ierr = DMSetUp(rdm);CHKERRQ(ierr);
3210 /* Step 4: Set cones and supports (automatically symmetrizes) */
3211 ierr = DMPlexCellRefinerSetCones(cr, rdm);CHKERRQ(ierr);
3212 /* Step 5: Create pointSF */
3213 ierr = DMPlexCellRefinerCreateSF(cr, rdm);CHKERRQ(ierr);
3214 /* Step 6: Create labels */
3215 ierr = DMPlexCellRefinerCreateLabels(cr, rdm);CHKERRQ(ierr);
3216 /* Step 7: Set coordinates */
3217 ierr = DMPlexCellRefinerSetCoordinates(cr, rdm);CHKERRQ(ierr);
3218 *dmRefined = rdm;
3219 PetscFunctionReturn(0);
3220 }
3221
3222 /*@
3223 DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data
3224
3225 Input Parameter:
3226 . dm - The coarse DM
3227
3228 Output Parameter:
3229 . fpointIS - The IS of all the fine points which exist in the original coarse mesh
3230
3231 Level: developer
3232
3233 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS()
3234 @*/
DMPlexCreateCoarsePointIS(DM dm,IS * fpointIS)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;
3241
3242 PetscFunctionBegin;
3243 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
3244 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
3245 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr);
3246 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr);
3247 ierr = PetscMalloc1(pEnd-pStart, &fpoints);CHKERRQ(ierr);
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 */
3251
3252 ierr = DMPlexCellRefinerGetNewPoint(cr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew);CHKERRQ(ierr);
3253 fpoints[v-pStart] = vNew;
3254 }
3255 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr);
3256 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS);CHKERRQ(ierr);
3257 PetscFunctionReturn(0);
3258 }
3259
3260 /*@
3261 DMPlexSetRefinementUniform - Set the flag for uniform refinement
3262
3263 Input Parameters:
3264 + dm - The DM
3265 - refinementUniform - The flag for uniform refinement
3266
3267 Level: developer
3268
3269 .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3270 @*/
DMPlexSetRefinementUniform(DM dm,PetscBool refinementUniform)3271 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
3272 {
3273 DM_Plex *mesh = (DM_Plex*) dm->data;
3274
3275 PetscFunctionBegin;
3276 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3277 mesh->refinementUniform = refinementUniform;
3278 PetscFunctionReturn(0);
3279 }
3280
3281 /*@
3282 DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement
3283
3284 Input Parameter:
3285 . dm - The DM
3286
3287 Output Parameter:
3288 . refinementUniform - The flag for uniform refinement
3289
3290 Level: developer
3291
3292 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3293 @*/
DMPlexGetRefinementUniform(DM dm,PetscBool * refinementUniform)3294 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
3295 {
3296 DM_Plex *mesh = (DM_Plex*) dm->data;
3297
3298 PetscFunctionBegin;
3299 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3300 PetscValidPointer(refinementUniform, 2);
3301 *refinementUniform = mesh->refinementUniform;
3302 PetscFunctionReturn(0);
3303 }
3304
3305 /*@
3306 DMPlexSetRefinementLimit - Set the maximum cell volume for refinement
3307
3308 Input Parameters:
3309 + dm - The DM
3310 - refinementLimit - The maximum cell volume in the refined mesh
3311
3312 Level: developer
3313
3314 .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3315 @*/
DMPlexSetRefinementLimit(DM dm,PetscReal refinementLimit)3316 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
3317 {
3318 DM_Plex *mesh = (DM_Plex*) dm->data;
3319
3320 PetscFunctionBegin;
3321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3322 mesh->refinementLimit = refinementLimit;
3323 PetscFunctionReturn(0);
3324 }
3325
3326 /*@
3327 DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement
3328
3329 Input Parameter:
3330 . dm - The DM
3331
3332 Output Parameter:
3333 . refinementLimit - The maximum cell volume in the refined mesh
3334
3335 Level: developer
3336
3337 .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform()
3338 @*/
DMPlexGetRefinementLimit(DM dm,PetscReal * refinementLimit)3339 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
3340 {
3341 DM_Plex *mesh = (DM_Plex*) dm->data;
3342
3343 PetscFunctionBegin;
3344 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3345 PetscValidPointer(refinementLimit, 2);
3346 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
3347 *refinementLimit = mesh->refinementLimit;
3348 PetscFunctionReturn(0);
3349 }
3350
3351 /*@
3352 DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement
3353
3354 Input Parameters:
3355 + dm - The DM
3356 - refinementFunc - Function giving the maximum cell volume in the refined mesh
3357
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
3361
3362 Level: developer
3363
3364 .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3365 @*/
DMPlexSetRefinementFunction(DM dm,PetscErrorCode (* refinementFunc)(const PetscReal[],PetscReal *))3366 PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *))
3367 {
3368 DM_Plex *mesh = (DM_Plex*) dm->data;
3369
3370 PetscFunctionBegin;
3371 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3372 mesh->refinementFunc = refinementFunc;
3373 PetscFunctionReturn(0);
3374 }
3375
3376 /*@
3377 DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement
3378
3379 Input Parameter:
3380 . dm - The DM
3381
3382 Output Parameter:
3383 . refinementFunc - Function giving the maximum cell volume in the refined mesh
3384
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
3388
3389 Level: developer
3390
3391 .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit()
3392 @*/
DMPlexGetRefinementFunction(DM dm,PetscErrorCode (** refinementFunc)(const PetscReal[],PetscReal *))3393 PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *))
3394 {
3395 DM_Plex *mesh = (DM_Plex*) dm->data;
3396
3397 PetscFunctionBegin;
3398 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
3399 PetscValidPointer(refinementFunc, 2);
3400 *refinementFunc = mesh->refinementFunc;
3401 PetscFunctionReturn(0);
3402 }
3403
RefineDiscLabels_Internal(DMPlexCellRefiner cr,DM rdm)3404 static PetscErrorCode RefineDiscLabels_Internal(DMPlexCellRefiner cr, DM rdm)
3405 {
3406 DM dm = cr->dm;
3407 PetscInt Nf, f, Nds, s;
3408 PetscErrorCode ierr;
3409
3410 PetscFunctionBegin;
3411 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
3412 for (f = 0; f < Nf; ++f) {
3413 DMLabel label, labelNew;
3414 PetscObject obj;
3415 const char *lname;
3416
3417 ierr = DMGetField(rdm, f, &label, &obj);CHKERRQ(ierr);
3418 if (!label) continue;
3419 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr);
3420 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr);
3421 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr);
3422 ierr = DMSetField_Internal(rdm, f, labelNew, obj);CHKERRQ(ierr);
3423 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr);
3424 }
3425 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr);
3426 for (s = 0; s < Nds; ++s) {
3427 DMLabel label, labelNew;
3428 const char *lname;
3429
3430 ierr = DMGetRegionNumDS(rdm, s, &label, NULL, NULL);CHKERRQ(ierr);
3431 if (!label) continue;
3432 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr);
3433 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr);
3434 ierr = RefineLabel_Internal(cr, label, labelNew);CHKERRQ(ierr);
3435 ierr = DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);CHKERRQ(ierr);
3436 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr);
3437 }
3438 PetscFunctionReturn(0);
3439 }
3440
DMRefine_Plex(DM dm,MPI_Comm comm,DM * dmRefined)3441 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
3442 {
3443 PetscBool isUniform;
3444 DMPlexCellRefiner cr;
3445 PetscErrorCode ierr;
3446
3447 PetscFunctionBegin;
3448 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
3449 ierr = DMViewFromOptions(dm, NULL, "-initref_dm_view");CHKERRQ(ierr);
3450 if (isUniform) {
3451 DM cdm, rcdm;
3452 PetscBool localized;
3453
3454 ierr = DMPlexCellRefinerCreate(dm, &cr);CHKERRQ(ierr);
3455 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr);
3456 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
3457 ierr = DMPlexRefineUniform(dm, cr, dmRefined);CHKERRQ(ierr);
3458 ierr = DMPlexSetRegularRefinement(*dmRefined, PETSC_TRUE);CHKERRQ(ierr);
3459 ierr = DMCopyDisc(dm, *dmRefined);CHKERRQ(ierr);
3460 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
3461 ierr = DMGetCoordinateDM(*dmRefined, &rcdm);CHKERRQ(ierr);
3462 ierr = DMCopyDisc(cdm, rcdm);CHKERRQ(ierr);
3463 ierr = RefineDiscLabels_Internal(cr, *dmRefined);CHKERRQ(ierr);
3464 ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
3465 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr);
3466 } else {
3467 ierr = DMPlexRefine_Internal(dm, NULL, dmRefined);CHKERRQ(ierr);
3468 }
3469 PetscFunctionReturn(0);
3470 }
3471
DMRefineHierarchy_Plex(DM dm,PetscInt nlevels,DM dmRefined[])3472 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
3473 {
3474 DM cdm = dm;
3475 PetscInt r;
3476 PetscBool isUniform, localized;
3477 PetscErrorCode ierr;
3478
3479 PetscFunctionBegin;
3480 ierr = DMPlexGetRefinementUniform(dm, &isUniform);CHKERRQ(ierr);
3481 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
3482 if (isUniform) {
3483 for (r = 0; r < nlevels; ++r) {
3484 DMPlexCellRefiner cr;
3485 DM codm, rcodm;
3486
3487 ierr = DMPlexCellRefinerCreate(cdm, &cr);CHKERRQ(ierr);
3488 ierr = DMPlexCellRefinerSetUp(cr);CHKERRQ(ierr);
3489 ierr = DMPlexRefineUniform(cdm, cr, &dmRefined[r]);CHKERRQ(ierr);
3490 ierr = DMSetCoarsenLevel(dmRefined[r], cdm->leveldown);CHKERRQ(ierr);
3491 ierr = DMSetRefineLevel(dmRefined[r], cdm->levelup+1);CHKERRQ(ierr);
3492 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr);
3493 ierr = DMGetCoordinateDM(dm, &codm);CHKERRQ(ierr);
3494 ierr = DMGetCoordinateDM(dmRefined[r], &rcodm);CHKERRQ(ierr);
3495 ierr = DMCopyDisc(codm, rcodm);CHKERRQ(ierr);
3496 ierr = RefineDiscLabels_Internal(cr, dmRefined[r]);CHKERRQ(ierr);
3497 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr);
3498 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
3499 ierr = DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);CHKERRQ(ierr);
3500 cdm = dmRefined[r];
3501 ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr);
3502 }
3503 } else {
3504 for (r = 0; r < nlevels; ++r) {
3505 ierr = DMRefine(cdm, PetscObjectComm((PetscObject) dm), &dmRefined[r]);CHKERRQ(ierr);
3506 ierr = DMCopyDisc(cdm, dmRefined[r]);CHKERRQ(ierr);
3507 ierr = DMCopyBoundary(cdm, dmRefined[r]);CHKERRQ(ierr);
3508 if (localized) {ierr = DMLocalizeCoordinates(dmRefined[r]);CHKERRQ(ierr);}
3509 ierr = DMSetCoarseDM(dmRefined[r], cdm);CHKERRQ(ierr);
3510 cdm = dmRefined[r];
3511 }
3512 }
3513 PetscFunctionReturn(0);
3514 }
3515