Actual source code: plexrefsbr.c
1: #include <petsc/private/dmplextransformimpl.h>
2: #include <petscsf.h>
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";
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;
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;
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: }
48: /* Mark local edges that should be split */
49: /* TODO This will not work in 3D */
50: static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, DMPlexPointQueue queue)
51: {
52: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
53: DM dm;
55: PetscFunctionBegin;
56: PetscCall(DMPlexTransformGetDM(tr, &dm));
57: while (!DMPlexPointQueueEmpty(queue)) {
58: PetscInt p = -1;
59: const PetscInt *support;
60: PetscInt supportSize, s;
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;
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: }
96: static PetscErrorCode splitPoint(PETSC_UNUSED DMLabel label, PetscInt p, PETSC_UNUSED PetscInt val, void *ctx)
97: {
98: DMPlexPointQueue queue = (DMPlexPointQueue)ctx;
100: PetscFunctionBegin;
101: PetscCall(DMPlexPointQueueEnqueue(queue, p));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
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:
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;
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;
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;
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;
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];
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.
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.
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;
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;
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: }
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;
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: }
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: */
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};
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: }
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: */
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};
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: }
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;
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: }
562: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
563: {
564: PetscInt cells[256], n = 256, i;
565: PetscBool flg;
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;
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: }
582: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
583: {
584: PetscBool isascii;
586: PetscFunctionBegin;
589: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
590: if (isascii) {
591: PetscViewerFormat format;
592: const char *name;
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: }
604: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
605: {
606: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
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: }
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: }
630: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
631: {
632: DMPlexRefine_SBR *f;
634: PetscFunctionBegin;
636: PetscCall(PetscNew(&f));
637: tr->data = f;
639: PetscCall(DMPlexTransformInitialize_SBR(tr));
640: PetscCall(PetscCitationsRegister(SBRCitation, &SBRcite));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }