1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/
2 #include <petscsf.h>
3
4 PetscBool SBRcite = PETSC_FALSE;
5 const char SBRCitation[] = "@article{PlazaCarey2000,\n"
6 " title = {Local refinement of simplicial grids based on the skeleton},\n"
7 " journal = {Applied Numerical Mathematics},\n"
8 " author = {A. Plaza and Graham F. Carey},\n"
9 " volume = {32},\n"
10 " number = {3},\n"
11 " pages = {195--218},\n"
12 " doi = {10.1016/S0168-9274(99)00022-7},\n"
13 " year = {2000}\n}\n";
14
SBRGetEdgeLen_Private(DMPlexTransform tr,PetscInt edge,PetscReal * len)15 static PetscErrorCode SBRGetEdgeLen_Private(DMPlexTransform tr, PetscInt edge, PetscReal *len)
16 {
17 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
18 DM dm;
19 PetscInt off;
20
21 PetscFunctionBeginHot;
22 PetscCall(DMPlexTransformGetDM(tr, &dm));
23 PetscCall(PetscSectionGetOffset(sbr->secEdgeLen, edge, &off));
24 if (sbr->edgeLen[off] <= 0.0) {
25 DM cdm;
26 Vec coordsLocal;
27 const PetscScalar *coords;
28 const PetscInt *cone;
29 PetscScalar *cA, *cB;
30 PetscInt coneSize, cdim;
31
32 PetscCall(DMGetCoordinateDM(dm, &cdm));
33 PetscCall(DMPlexGetCone(dm, edge, &cone));
34 PetscCall(DMPlexGetConeSize(dm, edge, &coneSize));
35 PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %" PetscInt_FMT " cone size must be 2, not %" PetscInt_FMT, edge, coneSize);
36 PetscCall(DMGetCoordinateDim(dm, &cdim));
37 PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordsLocal));
38 PetscCall(VecGetArrayRead(coordsLocal, &coords));
39 PetscCall(DMPlexPointLocalRead(cdm, cone[0], coords, &cA));
40 PetscCall(DMPlexPointLocalRead(cdm, cone[1], coords, &cB));
41 sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
42 PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
43 }
44 *len = sbr->edgeLen[off];
45 PetscFunctionReturn(PETSC_SUCCESS);
46 }
47
48 /* Mark local edges that should be split */
49 /* TODO This will not work in 3D */
SBRSplitLocalEdges_Private(DMPlexTransform tr,DMPlexPointQueue queue)50 static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, DMPlexPointQueue queue)
51 {
52 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
53 DM dm;
54
55 PetscFunctionBegin;
56 PetscCall(DMPlexTransformGetDM(tr, &dm));
57 while (!DMPlexPointQueueEmpty(queue)) {
58 PetscInt p = -1;
59 const PetscInt *support;
60 PetscInt supportSize, s;
61
62 PetscCall(DMPlexPointQueueDequeue(queue, &p));
63 PetscCall(DMPlexGetSupport(dm, p, &support));
64 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
65 for (s = 0; s < supportSize; ++s) {
66 const PetscInt cell = support[s];
67 const PetscInt *cone;
68 PetscInt coneSize, c;
69 PetscInt cval, eval, maxedge;
70 PetscReal len, maxlen;
71
72 PetscCall(DMLabelGetValue(sbr->splitPoints, cell, &cval));
73 if (cval == 2) continue;
74 PetscCall(DMPlexGetCone(dm, cell, &cone));
75 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
76 PetscCall(SBRGetEdgeLen_Private(tr, cone[0], &maxlen));
77 maxedge = cone[0];
78 for (c = 1; c < coneSize; ++c) {
79 PetscCall(SBRGetEdgeLen_Private(tr, cone[c], &len));
80 if (len > maxlen) {
81 maxlen = len;
82 maxedge = cone[c];
83 }
84 }
85 PetscCall(DMLabelGetValue(sbr->splitPoints, maxedge, &eval));
86 if (eval != 1) {
87 PetscCall(DMLabelSetValue(sbr->splitPoints, maxedge, 1));
88 PetscCall(DMPlexPointQueueEnqueue(queue, maxedge));
89 }
90 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 2));
91 }
92 }
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
splitPoint(PETSC_UNUSED DMLabel label,PetscInt p,PETSC_UNUSED PetscInt val,PetscCtx ctx)96 static PetscErrorCode splitPoint(PETSC_UNUSED DMLabel label, PetscInt p, PETSC_UNUSED PetscInt val, PetscCtx ctx)
97 {
98 DMPlexPointQueue queue = (DMPlexPointQueue)ctx;
99
100 PetscFunctionBegin;
101 PetscCall(DMPlexPointQueueEnqueue(queue, p));
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
105 /*
106 The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3.
107 Then the refinement type is calculated as follows:
108
109 vertex: 0
110 edge unsplit: 1
111 edge split: 2
112 triangle unsplit: 3
113 triangle split all edges: 4
114 triangle split edges 0 1: 5
115 triangle split edges 1 0: 6
116 triangle split edges 1 2: 7
117 triangle split edges 2 1: 8
118 triangle split edges 2 0: 9
119 triangle split edges 0 2: 10
120 triangle split edge 0: 11
121 triangle split edge 1: 12
122 triangle split edge 2: 13
123 */
124 typedef enum {
125 RT_VERTEX,
126 RT_EDGE,
127 RT_EDGE_SPLIT,
128 RT_TRIANGLE,
129 RT_TRIANGLE_SPLIT,
130 RT_TRIANGLE_SPLIT_01,
131 RT_TRIANGLE_SPLIT_10,
132 RT_TRIANGLE_SPLIT_12,
133 RT_TRIANGLE_SPLIT_21,
134 RT_TRIANGLE_SPLIT_20,
135 RT_TRIANGLE_SPLIT_02,
136 RT_TRIANGLE_SPLIT_0,
137 RT_TRIANGLE_SPLIT_1,
138 RT_TRIANGLE_SPLIT_2
139 } RefinementType;
140
DMPlexTransformSetUp_SBR(DMPlexTransform tr)141 static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr)
142 {
143 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
144 DM dm;
145 DMLabel active;
146 PetscSF pointSF;
147 DMPlexPointQueue queue = NULL;
148 IS refineIS;
149 const PetscInt *refineCells;
150 PetscInt pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
151 PetscBool empty;
152
153 PetscFunctionBegin;
154 PetscCall(DMPlexTransformGetDM(tr, &dm));
155 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints));
156 /* Create edge lengths */
157 PetscCall(DMGetCoordinatesLocalSetUp(dm));
158 PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
159 PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen));
160 PetscCall(PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd));
161 for (e = eStart; e < eEnd; ++e) PetscCall(PetscSectionSetDof(sbr->secEdgeLen, e, 1));
162 PetscCall(PetscSectionSetUp(sbr->secEdgeLen));
163 PetscCall(PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize));
164 PetscCall(PetscCalloc1(edgeLenSize, &sbr->edgeLen));
165 /* Add edges of cells that are marked for refinement to edge queue */
166 PetscCall(DMPlexTransformGetActive(tr, &active));
167 PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm");
168 PetscCall(DMPlexPointQueueCreate(1024, &queue));
169 PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS));
170 PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc));
171 if (refineIS) PetscCall(ISGetIndices(refineIS, &refineCells));
172 for (c = 0; c < Nc; ++c) {
173 const PetscInt cell = refineCells[c];
174 PetscInt depth;
175
176 PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
177 if (depth == 1) {
178 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 1));
179 PetscCall(DMPlexPointQueueEnqueue(queue, cell));
180 } else {
181 PetscInt *closure = NULL;
182 PetscInt Ncl, cl;
183
184 PetscCall(DMLabelSetValue(sbr->splitPoints, cell, depth));
185 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
186 for (cl = 0; cl < Ncl; cl += 2) {
187 const PetscInt edge = closure[cl];
188
189 if (edge >= eStart && edge < eEnd) {
190 PetscCall(DMLabelSetValue(sbr->splitPoints, edge, 1));
191 PetscCall(DMPlexPointQueueEnqueue(queue, edge));
192 }
193 }
194 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
195 }
196 }
197 if (refineIS) PetscCall(ISRestoreIndices(refineIS, &refineCells));
198 PetscCall(ISDestroy(&refineIS));
199 /* Setup communication */
200 PetscCall(DMGetPointSF(dm, &pointSF));
201 PetscCall(DMLabelPropagateBegin(sbr->splitPoints, pointSF));
202 /* While edge queue is not empty: */
203 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
204 while (!empty) {
205 PetscCall(SBRSplitLocalEdges_Private(tr, queue));
206 /* Communicate marked edges
207 An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
208 array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
209
210 TODO: We could use in-place communication with a different SF
211 We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
212 already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
213
214 In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
215 values will have 1+0=1 and old values will have 1+1=2. Loop over these, resetting the values to 1, and adding any new
216 edge to the queue.
217 */
218 PetscCall(DMLabelPropagatePush(sbr->splitPoints, pointSF, splitPoint, queue));
219 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
220 }
221 PetscCall(DMLabelPropagateEnd(sbr->splitPoints, pointSF));
222 /* Calculate refineType for each cell */
223 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
224 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
225 for (p = pStart; p < pEnd; ++p) {
226 DMLabel trType = tr->trType;
227 DMPolytopeType ct;
228 PetscInt val;
229
230 PetscCall(DMPlexGetCellType(dm, p, &ct));
231 switch (ct) {
232 case DM_POLYTOPE_POINT:
233 PetscCall(DMLabelSetValue(trType, p, RT_VERTEX));
234 break;
235 case DM_POLYTOPE_SEGMENT:
236 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
237 if (val == 1) PetscCall(DMLabelSetValue(trType, p, RT_EDGE_SPLIT));
238 else PetscCall(DMLabelSetValue(trType, p, RT_EDGE));
239 break;
240 case DM_POLYTOPE_TRIANGLE:
241 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
242 if (val == 2) {
243 const PetscInt *cone;
244 PetscReal lens[3];
245 PetscInt vals[3], i;
246
247 PetscCall(DMPlexGetCone(dm, p, &cone));
248 for (i = 0; i < 3; ++i) {
249 PetscCall(DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]));
250 vals[i] = vals[i] < 0 ? 0 : vals[i];
251 PetscCall(SBRGetEdgeLen_Private(tr, cone[i], &lens[i]));
252 }
253 if (vals[0] && vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT));
254 else if (vals[0] && vals[1]) PetscCall(DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10));
255 else if (vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21));
256 else if (vals[2] && vals[0]) PetscCall(DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02));
257 else if (vals[0]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0));
258 else if (vals[1]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1));
259 else if (vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2));
260 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " does not fit any refinement type (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")", p, vals[0], vals[1], vals[2]);
261 } else PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE));
262 break;
263 default:
264 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
265 }
266 PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
267 }
268 /* Cleanup */
269 PetscCall(DMPlexPointQueueDestroy(&queue));
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr,DMPolytopeType sct,PetscInt sp,PetscInt so,DMPolytopeType tct,PetscInt r,PetscInt o,PetscInt * rnew,PetscInt * onew)273 static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
274 {
275 PetscInt rt;
276
277 PetscFunctionBeginHot;
278 PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
279 *rnew = r;
280 *onew = o;
281 switch (rt) {
282 case RT_TRIANGLE_SPLIT_01:
283 case RT_TRIANGLE_SPLIT_10:
284 case RT_TRIANGLE_SPLIT_12:
285 case RT_TRIANGLE_SPLIT_21:
286 case RT_TRIANGLE_SPLIT_20:
287 case RT_TRIANGLE_SPLIT_02:
288 switch (tct) {
289 case DM_POLYTOPE_SEGMENT:
290 break;
291 case DM_POLYTOPE_TRIANGLE:
292 break;
293 default:
294 break;
295 }
296 break;
297 case RT_TRIANGLE_SPLIT_0:
298 case RT_TRIANGLE_SPLIT_1:
299 case RT_TRIANGLE_SPLIT_2:
300 switch (tct) {
301 case DM_POLYTOPE_SEGMENT:
302 break;
303 case DM_POLYTOPE_TRIANGLE:
304 *onew = so < 0 ? -(o + 1) : o;
305 *rnew = so < 0 ? (r + 1) % 2 : r;
306 break;
307 default:
308 break;
309 }
310 break;
311 case RT_EDGE_SPLIT:
312 case RT_TRIANGLE_SPLIT:
313 PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
314 break;
315 default:
316 PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
317 }
318 PetscFunctionReturn(PETSC_SUCCESS);
319 }
320
321 /* Add 1 edge inside this triangle, making 2 new triangles.
322 2
323 |\
324 | \
325 | \
326 | \
327 | 1
328 | \
329 | B \
330 2 1
331 | / \
332 | ____/ 0
333 |/ A \
334 0-----0-----1
335 */
SBRGetTriangleSplitSingle(PetscInt o,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])336 static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
337 {
338 const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
339 static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
340 static PetscInt triS1[] = {1, 2};
341 static PetscInt triC1[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0,
342 DM_POLYTOPE_SEGMENT, 0, 0};
343 static PetscInt triO1[] = {0, 0, 0, 0, -1, 0, 0, 0};
344
345 PetscFunctionBeginHot;
346 /* To get the other divisions, we reorient the triangle */
347 triC1[2] = arr[0 * 2];
348 triC1[7] = arr[1 * 2];
349 triC1[11] = arr[0 * 2];
350 triC1[15] = arr[1 * 2];
351 triC1[22] = arr[1 * 2];
352 triC1[26] = arr[2 * 2];
353 *Nt = 2;
354 *target = triT1;
355 *size = triS1;
356 *cone = triC1;
357 *ornt = triO1;
358 PetscFunctionReturn(PETSC_SUCCESS);
359 }
360
361 /* Add 2 edges inside this triangle, making 3 new triangles.
362 RT_TRIANGLE_SPLIT_12
363 2
364 |\
365 | \
366 | \
367 0 \
368 | 1
369 | \
370 | B \
371 2-------1
372 | C / \
373 1 ____/ 0
374 |/ A \
375 0-----0-----1
376 RT_TRIANGLE_SPLIT_10
377 2
378 |\
379 | \
380 | \
381 0 \
382 | 1
383 | \
384 | A \
385 2 1
386 | /|\
387 1 ____/ / 0
388 |/ C / B \
389 0-----0-----1
390 RT_TRIANGLE_SPLIT_20
391 2
392 |\
393 | \
394 | \
395 0 \
396 | \
397 | \
398 | \
399 2 A 1
400 |\ \
401 1 ---\ \
402 |B \_C----\\
403 0-----0-----1
404 RT_TRIANGLE_SPLIT_21
405 2
406 |\
407 | \
408 | \
409 0 \
410 | \
411 | B \
412 | \
413 2-------1
414 |\ C \
415 1 ---\ \
416 | A ----\\
417 0-----0-----1
418 RT_TRIANGLE_SPLIT_01
419 2
420 |\
421 |\\
422 || \
423 | \ \
424 | | \
425 | | \
426 | | \
427 2 \ C 1
428 | A | / \
429 | | |B \
430 | \/ \
431 0-----0-----1
432 RT_TRIANGLE_SPLIT_02
433 2
434 |\
435 |\\
436 || \
437 | \ \
438 | | \
439 | | \
440 | | \
441 2 C \ 1
442 |\ | \
443 | \__| A \
444 | B \\ \
445 0-----0-----1
446 */
SBRGetTriangleSplitDouble(PetscInt o,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])447 static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
448 {
449 PetscInt e0, e1;
450 const PetscInt *arr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
451 static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
452 static PetscInt triS2[] = {2, 3};
453 static PetscInt triC2[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1, DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
454 static PetscInt triO2[] = {0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0};
455
456 PetscFunctionBeginHot;
457 /* To get the other divisions, we reorient the triangle */
458 triC2[2] = arr[0 * 2];
459 triC2[3] = arr[0 * 2 + 1] ? 1 : 0;
460 triC2[7] = arr[1 * 2];
461 triC2[11] = arr[1 * 2];
462 triC2[15] = arr[2 * 2];
463 /* Swap the first two edges if the triangle is reversed */
464 e0 = o < 0 ? 23 : 19;
465 e1 = o < 0 ? 19 : 23;
466 triC2[e0] = arr[0 * 2];
467 triC2[e0 + 1] = 0;
468 triC2[e1] = arr[1 * 2];
469 triC2[e1 + 1] = o < 0 ? 1 : 0;
470 triO2[6] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
471 /* Swap the first two edges if the triangle is reversed */
472 e0 = o < 0 ? 34 : 30;
473 e1 = o < 0 ? 30 : 34;
474 triC2[e0] = arr[1 * 2];
475 triC2[e0 + 1] = o < 0 ? 0 : 1;
476 triC2[e1] = arr[2 * 2];
477 triC2[e1 + 1] = o < 0 ? 1 : 0;
478 triO2[9] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
479 /* Swap the last two edges if the triangle is reversed */
480 triC2[41] = arr[2 * 2];
481 triC2[42] = o < 0 ? 0 : 1;
482 triC2[45] = o < 0 ? 1 : 0;
483 triC2[48] = o < 0 ? 0 : 1;
484 triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1 * 2 + 1]);
485 triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2 * 2 + 1]);
486 *Nt = 2;
487 *target = triT2;
488 *size = triS2;
489 *cone = triC2;
490 *ornt = triO2;
491 PetscFunctionReturn(PETSC_SUCCESS);
492 }
493
DMPlexTransformCellTransform_SBR(DMPlexTransform tr,DMPolytopeType source,PetscInt p,PetscInt * rt,PetscInt * Nt,DMPolytopeType * target[],PetscInt * size[],PetscInt * cone[],PetscInt * ornt[])494 static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
495 {
496 DMLabel trType = tr->trType;
497 PetscInt val;
498
499 PetscFunctionBeginHot;
500 PetscCheck(p >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
501 PetscCall(DMLabelGetValue(trType, p, &val));
502 if (rt) *rt = val;
503 switch (source) {
504 case DM_POLYTOPE_POINT:
505 case DM_POLYTOPE_POINT_PRISM_TENSOR:
506 case DM_POLYTOPE_QUADRILATERAL:
507 case DM_POLYTOPE_SEG_PRISM_TENSOR:
508 case DM_POLYTOPE_TETRAHEDRON:
509 case DM_POLYTOPE_HEXAHEDRON:
510 case DM_POLYTOPE_TRI_PRISM:
511 case DM_POLYTOPE_TRI_PRISM_TENSOR:
512 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
513 case DM_POLYTOPE_PYRAMID:
514 PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
515 break;
516 case DM_POLYTOPE_SEGMENT:
517 if (val == RT_EDGE) PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
518 else PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
519 break;
520 case DM_POLYTOPE_TRIANGLE:
521 switch (val) {
522 case RT_TRIANGLE_SPLIT_0:
523 PetscCall(SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt));
524 break;
525 case RT_TRIANGLE_SPLIT_1:
526 PetscCall(SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt));
527 break;
528 case RT_TRIANGLE_SPLIT_2:
529 PetscCall(SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt));
530 break;
531 case RT_TRIANGLE_SPLIT_21:
532 PetscCall(SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt));
533 break;
534 case RT_TRIANGLE_SPLIT_10:
535 PetscCall(SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt));
536 break;
537 case RT_TRIANGLE_SPLIT_02:
538 PetscCall(SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt));
539 break;
540 case RT_TRIANGLE_SPLIT_12:
541 PetscCall(SBRGetTriangleSplitDouble(0, Nt, target, size, cone, ornt));
542 break;
543 case RT_TRIANGLE_SPLIT_20:
544 PetscCall(SBRGetTriangleSplitDouble(1, Nt, target, size, cone, ornt));
545 break;
546 case RT_TRIANGLE_SPLIT_01:
547 PetscCall(SBRGetTriangleSplitDouble(2, Nt, target, size, cone, ornt));
548 break;
549 case RT_TRIANGLE_SPLIT:
550 PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
551 break;
552 default:
553 PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
554 }
555 break;
556 default:
557 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
558 }
559 PetscFunctionReturn(PETSC_SUCCESS);
560 }
561
DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr,PetscOptionItems PetscOptionsObject)562 static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
563 {
564 PetscInt cells[256], n = 256, i;
565 PetscBool flg;
566
567 PetscFunctionBegin;
568 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
569 PetscCall(PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg));
570 if (flg) {
571 DMLabel active;
572
573 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active));
574 for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE));
575 PetscCall(DMPlexTransformSetActive(tr, active));
576 PetscCall(DMLabelDestroy(&active));
577 }
578 PetscOptionsHeadEnd();
579 PetscFunctionReturn(PETSC_SUCCESS);
580 }
581
DMPlexTransformView_SBR(DMPlexTransform tr,PetscViewer viewer)582 static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
583 {
584 PetscBool isascii;
585
586 PetscFunctionBegin;
587 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
588 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
589 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
590 if (isascii) {
591 PetscViewerFormat format;
592 const char *name;
593
594 PetscCall(PetscObjectGetName((PetscObject)tr, &name));
595 PetscCall(PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : ""));
596 PetscCall(PetscViewerGetFormat(viewer, &format));
597 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer));
598 } else {
599 SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
600 }
601 PetscFunctionReturn(PETSC_SUCCESS);
602 }
603
DMPlexTransformDestroy_SBR(DMPlexTransform tr)604 static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
605 {
606 DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
607
608 PetscFunctionBegin;
609 PetscCall(PetscFree(sbr->edgeLen));
610 PetscCall(PetscSectionDestroy(&sbr->secEdgeLen));
611 PetscCall(DMLabelDestroy(&sbr->splitPoints));
612 PetscCall(PetscFree(tr->data));
613 PetscFunctionReturn(PETSC_SUCCESS);
614 }
615
DMPlexTransformInitialize_SBR(DMPlexTransform tr)616 static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
617 {
618 PetscFunctionBegin;
619 tr->ops->view = DMPlexTransformView_SBR;
620 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR;
621 tr->ops->setup = DMPlexTransformSetUp_SBR;
622 tr->ops->destroy = DMPlexTransformDestroy_SBR;
623 tr->ops->setdimensions = DMPlexTransformSetDimensions_Internal;
624 tr->ops->celltransform = DMPlexTransformCellTransform_SBR;
625 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
626 tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
627 PetscFunctionReturn(PETSC_SUCCESS);
628 }
629
DMPlexTransformCreate_SBR(DMPlexTransform tr)630 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
631 {
632 DMPlexRefine_SBR *f;
633
634 PetscFunctionBegin;
635 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1);
636 PetscCall(PetscNew(&f));
637 tr->data = f;
638
639 PetscCall(DMPlexTransformInitialize_SBR(tr));
640 PetscCall(PetscCitationsRegister(SBRCitation, &SBRcite));
641 PetscFunctionReturn(PETSC_SUCCESS);
642 }
643