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