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: }