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