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: DMGetCoordinatesLocalNoncollective(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: DMGetCoordinatesLocalSetUp(dm);
155: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
156: PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen);
157: PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd);
158: for (e = eStart; e < eEnd; ++e) PetscSectionSetDof(sbr->secEdgeLen, e, 1);
159: PetscSectionSetUp(sbr->secEdgeLen);
160: PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize);
161: PetscCalloc1(edgeLenSize, &sbr->edgeLen);
162: /* Add edges of cells that are marked for refinement to edge queue */
163: DMPlexTransformGetActive(tr, &active);
165: DMPlexPointQueueCreate(1024, &queue);
166: DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
167: DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
168: if (refineIS) ISGetIndices(refineIS, &refineCells);
169: for (c = 0; c < Nc; ++c) {
170: const PetscInt cell = refineCells[c];
171: PetscInt depth;
173: DMPlexGetPointDepth(dm, cell, &depth);
174: if (depth == 1) {
175: DMLabelSetValue(sbr->splitPoints, cell, 1);
176: DMPlexPointQueueEnqueue(queue, cell);
177: } else {
178: PetscInt *closure = NULL;
179: PetscInt Ncl, cl;
181: DMLabelSetValue(sbr->splitPoints, cell, depth);
182: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
183: for (cl = 0; cl < Ncl; cl += 2) {
184: const PetscInt edge = closure[cl];
186: if (edge >= eStart && edge < eEnd) {
187: DMLabelSetValue(sbr->splitPoints, edge, 1);
188: DMPlexPointQueueEnqueue(queue, edge);
189: }
190: }
191: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
192: }
193: }
194: if (refineIS) ISRestoreIndices(refineIS, &refineCells);
195: ISDestroy(&refineIS);
196: /* Setup communication */
197: DMGetPointSF(dm, &pointSF);
198: DMLabelPropagateBegin(sbr->splitPoints, pointSF);
199: /* While edge queue is not empty: */
200: DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
201: while (!empty) {
202: SBRSplitLocalEdges_Private(tr, queue);
203: /* Communicate marked edges
204: An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
205: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
207: TODO: We could use in-place communication with a different SF
208: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
209: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
211: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
212: 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
213: edge to the queue.
214: */
215: DMLabelPropagatePush(sbr->splitPoints, pointSF, splitPoint, queue);
216: DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
217: }
218: DMLabelPropagateEnd(sbr->splitPoints, pointSF);
219: /* Calculate refineType for each cell */
220: DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
221: DMPlexGetChart(dm, &pStart, &pEnd);
222: for (p = pStart; p < pEnd; ++p) {
223: DMLabel trType = tr->trType;
224: DMPolytopeType ct;
225: PetscInt val;
227: DMPlexGetCellType(dm, p, &ct);
228: switch (ct) {
229: case DM_POLYTOPE_POINT:
230: DMLabelSetValue(trType, p, RT_VERTEX);
231: break;
232: case DM_POLYTOPE_SEGMENT:
233: DMLabelGetValue(sbr->splitPoints, p, &val);
234: if (val == 1) DMLabelSetValue(trType, p, RT_EDGE_SPLIT);
235: else DMLabelSetValue(trType, p, RT_EDGE);
236: break;
237: case DM_POLYTOPE_TRIANGLE:
238: DMLabelGetValue(sbr->splitPoints, p, &val);
239: if (val == 2) {
240: const PetscInt *cone;
241: PetscReal lens[3];
242: PetscInt vals[3], i;
244: DMPlexGetCone(dm, p, &cone);
245: for (i = 0; i < 3; ++i) {
246: DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]);
247: vals[i] = vals[i] < 0 ? 0 : vals[i];
248: SBRGetEdgeLen_Private(tr, cone[i], &lens[i]);
249: }
250: if (vals[0] && vals[1] && vals[2]) DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT);
251: else if (vals[0] && vals[1]) DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10);
252: else if (vals[1] && vals[2]) DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21);
253: else if (vals[2] && vals[0]) DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02);
254: else if (vals[0]) DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0);
255: else if (vals[1]) DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1);
256: else if (vals[2]) DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2);
257: 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]);
258: } else DMLabelSetValue(trType, p, RT_TRIANGLE);
259: break;
260: default:
261: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
262: }
263: DMLabelGetValue(sbr->splitPoints, p, &val);
264: }
265: /* Cleanup */
266: DMPlexPointQueueDestroy(&queue);
267: return 0;
268: }
270: static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
271: {
272: PetscInt rt;
275: DMLabelGetValue(tr->trType, sp, &rt);
276: *rnew = r;
277: *onew = o;
278: switch (rt) {
279: case RT_TRIANGLE_SPLIT_01:
280: case RT_TRIANGLE_SPLIT_10:
281: case RT_TRIANGLE_SPLIT_12:
282: case RT_TRIANGLE_SPLIT_21:
283: case RT_TRIANGLE_SPLIT_20:
284: case RT_TRIANGLE_SPLIT_02:
285: switch (tct) {
286: case DM_POLYTOPE_SEGMENT:
287: break;
288: case DM_POLYTOPE_TRIANGLE:
289: break;
290: default:
291: break;
292: }
293: break;
294: case RT_TRIANGLE_SPLIT_0:
295: case RT_TRIANGLE_SPLIT_1:
296: case RT_TRIANGLE_SPLIT_2:
297: switch (tct) {
298: case DM_POLYTOPE_SEGMENT:
299: break;
300: case DM_POLYTOPE_TRIANGLE:
301: *onew = so < 0 ? -(o + 1) : o;
302: *rnew = so < 0 ? (r + 1) % 2 : r;
303: break;
304: default:
305: break;
306: }
307: break;
308: case RT_EDGE_SPLIT:
309: case RT_TRIANGLE_SPLIT:
310: DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew);
311: break;
312: default:
313: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
314: }
315: return 0;
316: }
318: /* Add 1 edge inside this triangle, making 2 new triangles.
319: 2
320: |\
321: | \
322: | \
323: | \
324: | 1
325: | \
326: | B \
327: 2 1
328: | / \
329: | ____/ 0
330: |/ A \
331: 0-----0-----1
332: */
333: static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
334: {
335: const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
336: static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
337: static PetscInt triS1[] = {1, 2};
338: 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,
339: DM_POLYTOPE_SEGMENT, 0, 0};
340: static PetscInt triO1[] = {0, 0, 0, 0, -1, 0, 0, 0};
343: /* To get the other divisions, we reorient the triangle */
344: triC1[2] = arr[0 * 2];
345: triC1[7] = arr[1 * 2];
346: triC1[11] = arr[0 * 2];
347: triC1[15] = arr[1 * 2];
348: triC1[22] = arr[1 * 2];
349: triC1[26] = arr[2 * 2];
350: *Nt = 2;
351: *target = triT1;
352: *size = triS1;
353: *cone = triC1;
354: *ornt = triO1;
355: return 0;
356: }
358: /* Add 2 edges inside this triangle, making 3 new triangles.
359: RT_TRIANGLE_SPLIT_12
360: 2
361: |\
362: | \
363: | \
364: 0 \
365: | 1
366: | \
367: | B \
368: 2-------1
369: | C / \
370: 1 ____/ 0
371: |/ A \
372: 0-----0-----1
373: RT_TRIANGLE_SPLIT_10
374: 2
375: |\
376: | \
377: | \
378: 0 \
379: | 1
380: | \
381: | A \
382: 2 1
383: | /|\
384: 1 ____/ / 0
385: |/ C / B \
386: 0-----0-----1
387: RT_TRIANGLE_SPLIT_20
388: 2
389: |\
390: | \
391: | \
392: 0 \
393: | \
394: | \
395: | \
396: 2 A 1
397: |\ \
398: 1 ---\ \
399: |B \_C----\\
400: 0-----0-----1
401: RT_TRIANGLE_SPLIT_21
402: 2
403: |\
404: | \
405: | \
406: 0 \
407: | \
408: | B \
409: | \
410: 2-------1
411: |\ C \
412: 1 ---\ \
413: | A ----\\
414: 0-----0-----1
415: RT_TRIANGLE_SPLIT_01
416: 2
417: |\
418: |\\
419: || \
420: | \ \
421: | | \
422: | | \
423: | | \
424: 2 \ C 1
425: | A | / \
426: | | |B \
427: | \/ \
428: 0-----0-----1
429: RT_TRIANGLE_SPLIT_02
430: 2
431: |\
432: |\\
433: || \
434: | \ \
435: | | \
436: | | \
437: | | \
438: 2 C \ 1
439: |\ | \
440: | \__| A \
441: | B \\ \
442: 0-----0-----1
443: */
444: static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
445: {
446: PetscInt e0, e1;
447: const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
448: static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
449: static PetscInt triS2[] = {2, 3};
450: 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};
451: static PetscInt triO2[] = {0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0};
454: /* To get the other divisions, we reorient the triangle */
455: triC2[2] = arr[0 * 2];
456: triC2[3] = arr[0 * 2 + 1] ? 1 : 0;
457: triC2[7] = arr[1 * 2];
458: triC2[11] = arr[1 * 2];
459: triC2[15] = arr[2 * 2];
460: /* Swap the first two edges if the triangle is reversed */
461: e0 = o < 0 ? 23 : 19;
462: e1 = o < 0 ? 19 : 23;
463: triC2[e0] = arr[0 * 2];
464: triC2[e0 + 1] = 0;
465: triC2[e1] = arr[1 * 2];
466: triC2[e1 + 1] = o < 0 ? 1 : 0;
467: triO2[6] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
468: /* Swap the first two edges if the triangle is reversed */
469: e0 = o < 0 ? 34 : 30;
470: e1 = o < 0 ? 30 : 34;
471: triC2[e0] = arr[1 * 2];
472: triC2[e0 + 1] = o < 0 ? 0 : 1;
473: triC2[e1] = arr[2 * 2];
474: triC2[e1 + 1] = o < 0 ? 1 : 0;
475: triO2[9] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
476: /* Swap the last two edges if the triangle is reversed */
477: triC2[41] = arr[2 * 2];
478: triC2[42] = o < 0 ? 0 : 1;
479: triC2[45] = o < 0 ? 1 : 0;
480: triC2[48] = o < 0 ? 0 : 1;
481: triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1 * 2 + 1]);
482: triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2 * 2 + 1]);
483: *Nt = 2;
484: *target = triT2;
485: *size = triS2;
486: *cone = triC2;
487: *ornt = triO2;
488: return 0;
489: }
491: static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
492: {
493: DMLabel trType = tr->trType;
494: PetscInt val;
498: DMLabelGetValue(trType, p, &val);
499: if (rt) *rt = val;
500: switch (source) {
501: case DM_POLYTOPE_POINT:
502: case DM_POLYTOPE_POINT_PRISM_TENSOR:
503: case DM_POLYTOPE_QUADRILATERAL:
504: case DM_POLYTOPE_SEG_PRISM_TENSOR:
505: case DM_POLYTOPE_TETRAHEDRON:
506: case DM_POLYTOPE_HEXAHEDRON:
507: case DM_POLYTOPE_TRI_PRISM:
508: case DM_POLYTOPE_TRI_PRISM_TENSOR:
509: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
510: case DM_POLYTOPE_PYRAMID:
511: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
512: break;
513: case DM_POLYTOPE_SEGMENT:
514: if (val == RT_EDGE) DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
515: else DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);
516: break;
517: case DM_POLYTOPE_TRIANGLE:
518: switch (val) {
519: case RT_TRIANGLE_SPLIT_0:
520: SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt);
521: break;
522: case RT_TRIANGLE_SPLIT_1:
523: SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt);
524: break;
525: case RT_TRIANGLE_SPLIT_2:
526: SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt);
527: break;
528: case RT_TRIANGLE_SPLIT_21:
529: SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt);
530: break;
531: case RT_TRIANGLE_SPLIT_10:
532: SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt);
533: break;
534: case RT_TRIANGLE_SPLIT_02:
535: SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt);
536: break;
537: case RT_TRIANGLE_SPLIT_12:
538: SBRGetTriangleSplitDouble(0, Nt, target, size, cone, ornt);
539: break;
540: case RT_TRIANGLE_SPLIT_20:
541: SBRGetTriangleSplitDouble(1, Nt, target, size, cone, ornt);
542: break;
543: case RT_TRIANGLE_SPLIT_01:
544: SBRGetTriangleSplitDouble(2, Nt, target, size, cone, ornt);
545: break;
546: case RT_TRIANGLE_SPLIT:
547: DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);
548: break;
549: default:
550: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
551: }
552: break;
553: default:
554: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
555: }
556: return 0;
557: }
559: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
560: {
561: PetscInt cells[256], n = 256, i;
562: PetscBool flg;
564: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
565: PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
566: if (flg) {
567: DMLabel active;
569: DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
570: for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
571: DMPlexTransformSetActive(tr, active);
572: DMLabelDestroy(&active);
573: }
574: PetscOptionsHeadEnd();
575: return 0;
576: }
578: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
579: {
580: PetscBool isascii;
584: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
585: if (isascii) {
586: PetscViewerFormat format;
587: const char *name;
589: PetscObjectGetName((PetscObject)tr, &name);
590: PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "");
591: PetscViewerGetFormat(viewer, &format);
592: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) DMLabelView(tr->trType, viewer);
593: } else {
594: SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
595: }
596: return 0;
597: }
599: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
600: {
601: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
603: PetscFree(sbr->edgeLen);
604: PetscSectionDestroy(&sbr->secEdgeLen);
605: DMLabelDestroy(&sbr->splitPoints);
606: PetscFree(tr->data);
607: return 0;
608: }
610: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
611: {
612: tr->ops->view = DMPlexTransformView_SBR;
613: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR;
614: tr->ops->setup = DMPlexTransformSetUp_SBR;
615: tr->ops->destroy = DMPlexTransformDestroy_SBR;
616: tr->ops->celltransform = DMPlexTransformCellTransform_SBR;
617: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
618: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
619: return 0;
620: }
622: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
623: {
624: DMPlexRefine_SBR *f;
627: PetscNew(&f);
628: tr->data = f;
630: DMPlexTransformInitialize_SBR(tr);
631: PetscCitationsRegister(SBRCitation, &SBRcite);
632: return 0;
633: }