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