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: typedef struct _p_PointQueue *PointQueue;
16: struct _p_PointQueue {
17: PetscInt size; /* Size of the storage array */
18: PetscInt *points; /* Array of mesh points */
19: PetscInt front; /* Index of the front of the queue */
20: PetscInt back; /* Index of the back of the queue */
21: PetscInt num; /* Number of enqueued points */
22: };
24: static PetscErrorCode PointQueueCreate(PetscInt size, PointQueue *queue)
25: {
26: PointQueue q;
30: if (size < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Queue size %D must be non-negative", size);
31: PetscCalloc1(1, &q);
32: q->size = size;
33: PetscMalloc1(q->size, &q->points);
34: q->num = 0;
35: q->front = 0;
36: q->back = q->size-1;
37: *queue = q;
38: return(0);
39: }
41: static PetscErrorCode PointQueueDestroy(PointQueue *queue)
42: {
43: PointQueue q = *queue;
47: PetscFree(q->points);
48: PetscFree(q);
49: *queue = NULL;
50: return(0);
51: }
53: static PetscErrorCode PointQueueEnsureSize(PointQueue queue)
54: {
58: if (queue->num < queue->size) return(0);
59: queue->size *= 2;
60: PetscRealloc(queue->size * sizeof(PetscInt), &queue->points);
61: return(0);
62: }
64: static PetscErrorCode PointQueueEnqueue(PointQueue queue, PetscInt p)
65: {
69: PointQueueEnsureSize(queue);
70: queue->back = (queue->back + 1) % queue->size;
71: queue->points[queue->back] = p;
72: ++queue->num;
73: return(0);
74: }
76: static PetscErrorCode PointQueueDequeue(PointQueue queue, PetscInt *p)
77: {
79: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot dequeue from an empty queue");
80: *p = queue->points[queue->front];
81: queue->front = (queue->front + 1) % queue->size;
82: --queue->num;
83: return(0);
84: }
86: #if 0
87: static PetscErrorCode PointQueueFront(PointQueue queue, PetscInt *p)
88: {
90: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the front of an empty queue");
91: *p = queue->points[queue->front];
92: return(0);
93: }
95: static PetscErrorCode PointQueueBack(PointQueue queue, PetscInt *p)
96: {
98: if (!queue->num) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot get the back of an empty queue");
99: *p = queue->points[queue->back];
100: return(0);
101: }
102: #endif
104: PETSC_STATIC_INLINE PetscBool PointQueueEmpty(PointQueue queue)
105: {
106: if (!queue->num) return PETSC_TRUE;
107: return PETSC_FALSE;
108: }
110: static PetscErrorCode SBRGetEdgeLen_Private(DMPlexTransform tr, PetscInt edge, PetscReal *len)
111: {
112: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
113: DM dm;
114: PetscInt off;
115: PetscErrorCode ierr;
118: DMPlexTransformGetDM(tr, &dm);
119: PetscSectionGetOffset(sbr->secEdgeLen, edge, &off);
120: if (sbr->edgeLen[off] <= 0.0) {
121: DM cdm;
122: Vec coordsLocal;
123: const PetscScalar *coords;
124: const PetscInt *cone;
125: PetscScalar *cA, *cB;
126: PetscInt coneSize, cdim;
128: DMGetCoordinateDM(dm, &cdm);
129: DMPlexGetCone(dm, edge, &cone);
130: DMPlexGetConeSize(dm, edge, &coneSize);
131: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %D cone size must be 2, not %D", edge, coneSize);
132: DMGetCoordinateDim(dm, &cdim);
133: DMGetCoordinatesLocal(dm, &coordsLocal);
134: VecGetArrayRead(coordsLocal, &coords);
135: DMPlexPointLocalRead(cdm, cone[0], coords, &cA);
136: DMPlexPointLocalRead(cdm, cone[1], coords, &cB);
137: sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
138: VecRestoreArrayRead(coordsLocal, &coords);
139: }
140: *len = sbr->edgeLen[off];
141: return(0);
142: }
144: /* Mark local edges that should be split */
145: /* TODO This will not work in 3D */
146: static PetscErrorCode SBRSplitLocalEdges_Private(DMPlexTransform tr, PointQueue queue)
147: {
148: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
149: DM dm;
150: PetscErrorCode ierr;
153: DMPlexTransformGetDM(tr, &dm);
154: while (!PointQueueEmpty(queue)) {
155: PetscInt p = -1;
156: const PetscInt *support;
157: PetscInt supportSize, s;
159: PointQueueDequeue(queue, &p);
160: DMPlexGetSupport(dm, p, &support);
161: DMPlexGetSupportSize(dm, p, &supportSize);
162: for (s = 0; s < supportSize; ++s) {
163: const PetscInt cell = support[s];
164: const PetscInt *cone;
165: PetscInt coneSize, c;
166: PetscInt cval, eval, maxedge;
167: PetscReal len, maxlen;
169: DMLabelGetValue(sbr->splitPoints, cell, &cval);
170: if (cval == 2) continue;
171: DMPlexGetCone(dm, cell, &cone);
172: DMPlexGetConeSize(dm, cell, &coneSize);
173: SBRGetEdgeLen_Private(tr, cone[0], &maxlen);
174: maxedge = cone[0];
175: for (c = 1; c < coneSize; ++c) {
176: SBRGetEdgeLen_Private(tr, cone[c], &len);
177: if (len > maxlen) {maxlen = len; maxedge = cone[c];}
178: }
179: DMLabelGetValue(sbr->splitPoints, maxedge, &eval);
180: if (eval != 1) {
181: DMLabelSetValue(sbr->splitPoints, maxedge, 1);
182: PointQueueEnqueue(queue, maxedge);
183: }
184: DMLabelSetValue(sbr->splitPoints, cell, 2);
185: }
186: }
187: return(0);
188: }
190: static PetscErrorCode SBRInitializeComm(DMPlexTransform tr, PetscSF pointSF)
191: {
192: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
193: DM dm;
194: DMLabel splitPoints = sbr->splitPoints;
195: PetscInt *splitArray = sbr->splitArray;
196: const PetscInt *degree;
197: const PetscInt *points;
198: PetscInt Nl, l, pStart, pEnd, p, val;
199: PetscErrorCode ierr;
202: DMPlexTransformGetDM(tr, &dm);
203: /* Add in leaves */
204: PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
205: for (l = 0; l < Nl; ++l) {
206: DMLabelGetValue(splitPoints, points[l], &val);
207: if (val > 0) splitArray[points[l]] = val;
208: }
209: /* Add in shared roots */
210: PetscSFComputeDegreeBegin(pointSF, °ree);
211: PetscSFComputeDegreeEnd(pointSF, °ree);
212: DMPlexGetChart(dm, &pStart, &pEnd);
213: for (p = pStart; p < pEnd; ++p) {
214: if (degree[p]) {
215: DMLabelGetValue(splitPoints, p, &val);
216: if (val > 0) splitArray[p] = val;
217: }
218: }
219: return(0);
220: }
222: static PetscErrorCode SBRFinalizeComm(DMPlexTransform tr, PetscSF pointSF, PointQueue queue)
223: {
224: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
225: DM dm;
226: DMLabel splitPoints = sbr->splitPoints;
227: PetscInt *splitArray = sbr->splitArray;
228: const PetscInt *degree;
229: const PetscInt *points;
230: PetscInt Nl, l, pStart, pEnd, p, val;
231: PetscErrorCode ierr;
234: DMPlexTransformGetDM(tr, &dm);
235: /* Read out leaves */
236: PetscSFGetGraph(pointSF, NULL, &Nl, &points, NULL);
237: for (l = 0; l < Nl; ++l) {
238: const PetscInt p = points[l];
239: const PetscInt cval = splitArray[p];
241: if (cval) {
242: DMLabelGetValue(splitPoints, p, &val);
243: if (val <= 0) {
244: DMLabelSetValue(splitPoints, p, cval);
245: PointQueueEnqueue(queue, p);
246: }
247: }
248: }
249: /* Read out shared roots */
250: PetscSFComputeDegreeBegin(pointSF, °ree);
251: PetscSFComputeDegreeEnd(pointSF, °ree);
252: DMPlexGetChart(dm, &pStart, &pEnd);
253: for (p = pStart; p < pEnd; ++p) {
254: if (degree[p]) {
255: const PetscInt cval = splitArray[p];
257: if (cval) {
258: DMLabelGetValue(splitPoints, p, &val);
259: if (val <= 0) {
260: DMLabelSetValue(splitPoints, p, cval);
261: PointQueueEnqueue(queue, p);
262: }
263: }
264: }
265: }
266: return(0);
267: }
269: /*
270: The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3.
271: Then the refinement type is calculated as follows:
273: vertex: 0
274: edge unsplit: 1
275: edge split: 2
276: triangle unsplit: 3
277: triangle split all edges: 4
278: triangle split edges 0 1: 5
279: triangle split edges 1 0: 6
280: triangle split edges 1 2: 7
281: triangle split edges 2 1: 8
282: triangle split edges 2 0: 9
283: triangle split edges 0 2: 10
284: triangle split edge 0: 11
285: triangle split edge 1: 12
286: triangle split edge 2: 13
287: */
288: typedef enum {RT_VERTEX, RT_EDGE, RT_EDGE_SPLIT, RT_TRIANGLE, RT_TRIANGLE_SPLIT, RT_TRIANGLE_SPLIT_01, RT_TRIANGLE_SPLIT_10, RT_TRIANGLE_SPLIT_12, RT_TRIANGLE_SPLIT_21, RT_TRIANGLE_SPLIT_20, RT_TRIANGLE_SPLIT_02, RT_TRIANGLE_SPLIT_0, RT_TRIANGLE_SPLIT_1, RT_TRIANGLE_SPLIT_2} RefinementType;
290: static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr)
291: {
292: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
293: DM dm;
294: DMLabel active;
295: PetscSF pointSF;
296: PointQueue queue = NULL;
297: IS refineIS;
298: const PetscInt *refineCells;
299: PetscMPIInt size;
300: PetscInt pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
301: PetscBool empty;
302: PetscErrorCode ierr;
305: DMPlexTransformGetDM(tr, &dm);
306: DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints);
307: /* Create edge lengths */
308: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
309: PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen);
310: PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd);
311: for (e = eStart; e < eEnd; ++e) {
312: PetscSectionSetDof(sbr->secEdgeLen, e, 1);
313: }
314: PetscSectionSetUp(sbr->secEdgeLen);
315: PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize);
316: PetscCalloc1(edgeLenSize, &sbr->edgeLen);
317: /* Add edges of cells that are marked for refinement to edge queue */
318: DMPlexTransformGetActive(tr, &active);
319: if (!active) SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm");
320: PointQueueCreate(1024, &queue);
321: DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS);
322: DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc);
323: if (refineIS) {ISGetIndices(refineIS, &refineCells);}
324: for (c = 0; c < Nc; ++c) {
325: const PetscInt cell = refineCells[c];
326: PetscInt depth;
328: DMPlexGetPointDepth(dm, cell, &depth);
329: if (depth == 1) {
330: DMLabelSetValue(sbr->splitPoints, cell, 1);
331: PointQueueEnqueue(queue, cell);
332: } else {
333: PetscInt *closure = NULL;
334: PetscInt Ncl, cl;
336: DMLabelSetValue(sbr->splitPoints, cell, depth);
337: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
338: for (cl = 0; cl < Ncl; cl += 2) {
339: const PetscInt edge = closure[cl];
341: if (edge >= eStart && edge < eEnd) {
342: DMLabelSetValue(sbr->splitPoints, edge, 1);
343: PointQueueEnqueue(queue, edge);
344: }
345: }
346: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
347: }
348: }
349: if (refineIS) {ISRestoreIndices(refineIS, &refineCells);}
350: ISDestroy(&refineIS);
351: /* Setup communication */
352: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
353: DMGetPointSF(dm, &pointSF);
354: if (size > 1) {
355: PetscInt pStart, pEnd;
357: DMPlexGetChart(dm, &pStart, &pEnd);
358: PetscCalloc1(pEnd-pStart, &sbr->splitArray);
359: }
360: /* While edge queue is not empty: */
361: empty = PointQueueEmpty(queue);
362: MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
363: while (!empty) {
364: SBRSplitLocalEdges_Private(tr, queue);
365: /* Communicate marked edges
366: An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
367: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
369: TODO: We could use in-place communication with a different SF
370: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
371: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
373: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
374: 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
375: edge to the queue.
376: */
377: if (size > 1) {
378: SBRInitializeComm(tr, pointSF);
379: PetscSFReduceBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
380: PetscSFReduceEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray, MPI_MAX);
381: PetscSFBcastBegin(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
382: PetscSFBcastEnd(pointSF, MPIU_INT, sbr->splitArray, sbr->splitArray,MPI_REPLACE);
383: SBRFinalizeComm(tr, pointSF, queue);
384: }
385: empty = PointQueueEmpty(queue);
386: MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject) dm));
387: }
388: PetscFree(sbr->splitArray);
389: /* Calculate refineType for each cell */
390: DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
391: DMPlexGetChart(dm, &pStart, &pEnd);
392: for (p = pStart; p < pEnd; ++p) {
393: DMLabel trType = tr->trType;
394: DMPolytopeType ct;
395: PetscInt val;
397: DMPlexGetCellType(dm, p, &ct);
398: switch (ct) {
399: case DM_POLYTOPE_POINT:
400: DMLabelSetValue(trType, p, RT_VERTEX);break;
401: case DM_POLYTOPE_SEGMENT:
402: DMLabelGetValue(sbr->splitPoints, p, &val);
403: if (val == 1) {DMLabelSetValue(trType, p, RT_EDGE_SPLIT);}
404: else {DMLabelSetValue(trType, p, RT_EDGE);}
405: break;
406: case DM_POLYTOPE_TRIANGLE:
407: DMLabelGetValue(sbr->splitPoints, p, &val);
408: if (val == 2) {
409: const PetscInt *cone;
410: PetscReal lens[3];
411: PetscInt vals[3], i;
413: DMPlexGetCone(dm, p, &cone);
414: for (i = 0; i < 3; ++i) {
415: DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]);
416: vals[i] = vals[i] < 0 ? 0 : vals[i];
417: SBRGetEdgeLen_Private(tr, cone[i], &lens[i]);
418: }
419: if (vals[0] && vals[1] && vals[2]) {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT);}
420: else if (vals[0] && vals[1]) {DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10);}
421: else if (vals[1] && vals[2]) {DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21);}
422: else if (vals[2] && vals[0]) {DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02);}
423: else if (vals[0]) {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0);}
424: else if (vals[1]) {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1);}
425: else if (vals[2]) {DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2);}
426: else SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D does not fit any refinement type (%D, %D, %D)", p, vals[0], vals[1], vals[2]);
427: } else {DMLabelSetValue(trType, p, RT_TRIANGLE);}
428: break;
429: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
430: }
431: DMLabelGetValue(sbr->splitPoints, p, &val);
432: }
433: /* Cleanup */
434: PointQueueDestroy(&queue);
435: return(0);
436: }
438: static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
439: {
440: PetscInt rt;
441: PetscErrorCode ierr;
444: DMLabelGetValue(tr->trType, sp, &rt);
445: *rnew = r;
446: *onew = o;
447: switch (rt) {
448: case RT_TRIANGLE_SPLIT_01:
449: case RT_TRIANGLE_SPLIT_10:
450: case RT_TRIANGLE_SPLIT_12:
451: case RT_TRIANGLE_SPLIT_21:
452: case RT_TRIANGLE_SPLIT_20:
453: case RT_TRIANGLE_SPLIT_02:
454: switch (tct) {
455: case DM_POLYTOPE_SEGMENT: break;
456: case DM_POLYTOPE_TRIANGLE: break;
457: default: break;
458: }
459: break;
460: case RT_TRIANGLE_SPLIT_0:
461: case RT_TRIANGLE_SPLIT_1:
462: case RT_TRIANGLE_SPLIT_2:
463: switch (tct) {
464: case DM_POLYTOPE_SEGMENT:
465: break;
466: case DM_POLYTOPE_TRIANGLE:
467: *onew = so < 0 ? -(o+1) : o;
468: *rnew = so < 0 ? (r+1)%2 : r;
469: break;
470: default: break;
471: }
472: break;
473: case RT_EDGE_SPLIT:
474: case RT_TRIANGLE_SPLIT:
475: DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew);
476: break;
477: default: DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
478: }
479: return(0);
480: }
482: /* Add 1 edge inside this triangle, making 2 new triangles.
483: 2
484: |\
485: | \
486: | \
487: | \
488: | 1
489: | \
490: | B \
491: 2 1
492: | / \
493: | ____/ 0
494: |/ A \
495: 0-----0-----1
496: */
497: static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
498: {
499: const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
500: static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
501: static PetscInt triS1[] = {1, 2};
502: static PetscInt triC1[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
503: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
504: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 0};
505: static PetscInt triO1[] = {0, 0,
506: 0, 0, -1,
507: 0, 0, 0};
510: /* To get the other divisions, we reorient the triangle */
511: triC1[2] = arr[0*2];
512: triC1[7] = arr[1*2];
513: triC1[11] = arr[0*2];
514: triC1[15] = arr[1*2];
515: triC1[22] = arr[1*2];
516: triC1[26] = arr[2*2];
517: *Nt = 2; *target = triT1; *size = triS1; *cone = triC1; *ornt = triO1;
518: return(0);
519: }
521: /* Add 2 edges inside this triangle, making 3 new triangles.
522: RT_TRIANGLE_SPLIT_12
523: 2
524: |\
525: | \
526: | \
527: 0 \
528: | 1
529: | \
530: | B \
531: 2-------1
532: | C / \
533: 1 ____/ 0
534: |/ A \
535: 0-----0-----1
536: RT_TRIANGLE_SPLIT_10
537: 2
538: |\
539: | \
540: | \
541: 0 \
542: | 1
543: | \
544: | A \
545: 2 1
546: | /|\
547: 1 ____/ / 0
548: |/ C / B \
549: 0-----0-----1
550: RT_TRIANGLE_SPLIT_20
551: 2
552: |\
553: | \
554: | \
555: 0 \
556: | \
557: | \
558: | \
559: 2 A 1
560: |\ \
561: 1 ---\ \
562: |B \_C----\\
563: 0-----0-----1
564: RT_TRIANGLE_SPLIT_21
565: 2
566: |\
567: | \
568: | \
569: 0 \
570: | \
571: | B \
572: | \
573: 2-------1
574: |\ C \
575: 1 ---\ \
576: | A ----\\
577: 0-----0-----1
578: RT_TRIANGLE_SPLIT_01
579: 2
580: |\
581: |\\
582: || \
583: | \ \
584: | | \
585: | | \
586: | | \
587: 2 \ C 1
588: | A | / \
589: | | |B \
590: | \/ \
591: 0-----0-----1
592: RT_TRIANGLE_SPLIT_02
593: 2
594: |\
595: |\\
596: || \
597: | \ \
598: | | \
599: | | \
600: | | \
601: 2 C \ 1
602: |\ | \
603: | \__| A \
604: | B \\ \
605: 0-----0-----1
606: */
607: static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
608: {
609: PetscInt e0, e1;
610: const PetscInt *arr = DMPolytopeTypeGetArrangment(DM_POLYTOPE_TRIANGLE, o);
611: static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
612: static PetscInt triS2[] = {2, 3};
613: static PetscInt triC2[] = {DM_POLYTOPE_POINT, 2, 0, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0,
614: DM_POLYTOPE_POINT, 1, 1, 0, DM_POLYTOPE_POINT, 1, 2, 0,
615: DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 0, 0,
616: DM_POLYTOPE_SEGMENT, 1, 1, 1, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 0, 1,
617: DM_POLYTOPE_SEGMENT, 1, 2, 1, DM_POLYTOPE_SEGMENT, 0, 0, DM_POLYTOPE_SEGMENT, 0, 1};
618: static PetscInt triO2[] = {0, 0,
619: 0, 0,
620: 0, 0, -1,
621: 0, 0, -1,
622: 0, 0, 0};
625: /* To get the other divisions, we reorient the triangle */
626: triC2[2] = arr[0*2]; triC2[3] = arr[0*2+1] ? 1 : 0;
627: triC2[7] = arr[1*2];
628: triC2[11] = arr[1*2];
629: triC2[15] = arr[2*2];
630: /* Swap the first two edges if the triangle is reversed */
631: e0 = o < 0 ? 23: 19;
632: e1 = o < 0 ? 19: 23;
633: triC2[e0] = arr[0*2]; triC2[e0+1] = 0;
634: triC2[e1] = arr[1*2]; triC2[e1+1] = o < 0 ? 1 : 0;
635: triO2[6] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2*2+1]);
636: /* Swap the first two edges if the triangle is reversed */
637: e0 = o < 0 ? 34: 30;
638: e1 = o < 0 ? 30: 34;
639: triC2[e0] = arr[1*2]; triC2[e0+1] = o < 0 ? 0 : 1;
640: triC2[e1] = arr[2*2]; triC2[e1+1] = o < 0 ? 1 : 0;
641: triO2[9] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2*2+1]);
642: /* Swap the last two edges if the triangle is reversed */
643: triC2[41] = arr[2*2]; triC2[42] = o < 0 ? 0 : 1;
644: triC2[45] = o < 0 ? 1 : 0;
645: triC2[48] = o < 0 ? 0 : 1;
646: triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1*2+1]);
647: triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2*2+1]);
648: *Nt = 2; *target = triT2; *size = triS2; *cone = triC2; *ornt = triO2;
649: return(0);
650: }
652: static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
653: {
654: DMLabel trType = tr->trType;
655: PetscInt val;
659: if (p < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
660: DMLabelGetValue(trType, p, &val);
661: if (rt) *rt = val;
662: switch (source) {
663: case DM_POLYTOPE_POINT:
664: case DM_POLYTOPE_POINT_PRISM_TENSOR:
665: case DM_POLYTOPE_QUADRILATERAL:
666: case DM_POLYTOPE_SEG_PRISM_TENSOR:
667: case DM_POLYTOPE_TETRAHEDRON:
668: case DM_POLYTOPE_HEXAHEDRON:
669: case DM_POLYTOPE_TRI_PRISM:
670: case DM_POLYTOPE_TRI_PRISM_TENSOR:
671: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
672: case DM_POLYTOPE_PYRAMID:
673: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
674: break;
675: case DM_POLYTOPE_SEGMENT:
676: if (val == RT_EDGE) {DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);}
677: else {DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);}
678: break;
679: case DM_POLYTOPE_TRIANGLE:
680: switch (val) {
681: case RT_TRIANGLE_SPLIT_0: SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt);break;
682: case RT_TRIANGLE_SPLIT_1: SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt);break;
683: case RT_TRIANGLE_SPLIT_2: SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt);break;
684: case RT_TRIANGLE_SPLIT_21: SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt);break;
685: case RT_TRIANGLE_SPLIT_10: SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt);break;
686: case RT_TRIANGLE_SPLIT_02: SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt);break;
687: case RT_TRIANGLE_SPLIT_12: SBRGetTriangleSplitDouble( 0, Nt, target, size, cone, ornt);break;
688: case RT_TRIANGLE_SPLIT_20: SBRGetTriangleSplitDouble( 1, Nt, target, size, cone, ornt);break;
689: case RT_TRIANGLE_SPLIT_01: SBRGetTriangleSplitDouble( 2, Nt, target, size, cone, ornt);break;
690: case RT_TRIANGLE_SPLIT: DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt); break;
691: default: DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
692: }
693: break;
694: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
695: }
696: return(0);
697: }
699: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
700: {
701: PetscInt cells[256], n = 256, i;
702: PetscBool flg;
707: PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
708: PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
709: if (flg) {
710: DMLabel active;
712: DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
713: for (i = 0; i < n; ++i) {DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);}
714: DMPlexTransformSetActive(tr, active);
715: DMLabelDestroy(&active);
716: }
717: PetscOptionsTail();
718: return(0);
719: }
721: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
722: {
723: PetscBool isascii;
729: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
730: if (isascii) {
731: PetscViewerFormat format;
732: const char *name;
734: PetscObjectGetName((PetscObject) tr, &name);
735: PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : "");
736: PetscViewerGetFormat(viewer, &format);
737: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
738: DMLabelView(tr->trType, viewer);
739: }
740: } else {
741: SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name);
742: }
743: return(0);
744: }
746: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
747: {
748: DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *) tr->data;
749: PetscErrorCode ierr;
752: PetscFree(sbr->edgeLen);
753: PetscSectionDestroy(&sbr->secEdgeLen);
754: DMLabelDestroy(&sbr->splitPoints);
755: PetscFree(tr->data);
756: return(0);
757: }
759: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
760: {
762: tr->ops->view = DMPlexTransformView_SBR;
763: tr->ops->setfromoptions = DMPlexTransformSetFromOptions_SBR;
764: tr->ops->setup = DMPlexTransformSetUp_SBR;
765: tr->ops->destroy = DMPlexTransformDestroy_SBR;
766: tr->ops->celltransform = DMPlexTransformCellTransform_SBR;
767: tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
768: tr->ops->mapcoordinates = DMPlexTransformMapCoordinatesBarycenter_Internal;
769: return(0);
770: }
772: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
773: {
774: DMPlexRefine_SBR *f;
775: PetscErrorCode ierr;
779: PetscNewLog(tr, &f);
780: tr->data = f;
782: DMPlexTransformInitialize_SBR(tr);
783: PetscCitationsRegister(SBRCitation, &SBRcite);
784: return(0);
785: }