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;

     21:   PetscFunctionBeginHot;
     22:   PetscCall(DMPlexTransformGetDM(tr, &dm));
     23:   PetscCall(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:     PetscCall(DMGetCoordinateDM(dm, &cdm));
     33:     PetscCall(DMPlexGetCone(dm, edge, &cone));
     34:     PetscCall(DMPlexGetConeSize(dm, edge, &coneSize));
     35:     PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Edge %" PetscInt_FMT " cone size must be 2, not %" PetscInt_FMT, edge, coneSize);
     36:     PetscCall(DMGetCoordinateDim(dm, &cdim));
     37:     PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordsLocal));
     38:     PetscCall(VecGetArrayRead(coordsLocal, &coords));
     39:     PetscCall(DMPlexPointLocalRead(cdm, cone[0], coords, &cA));
     40:     PetscCall(DMPlexPointLocalRead(cdm, cone[1], coords, &cB));
     41:     sbr->edgeLen[off] = DMPlex_DistD_Internal(cdim, cA, cB);
     42:     PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
     43:   }
     44:   *len = sbr->edgeLen[off];
     45:   PetscFunctionReturn(PETSC_SUCCESS);
     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:   PetscFunctionBegin;
     56:   PetscCall(DMPlexTransformGetDM(tr, &dm));
     57:   while (!DMPlexPointQueueEmpty(queue)) {
     58:     PetscInt        p = -1;
     59:     const PetscInt *support;
     60:     PetscInt        supportSize, s;

     62:     PetscCall(DMPlexPointQueueDequeue(queue, &p));
     63:     PetscCall(DMPlexGetSupport(dm, p, &support));
     64:     PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
     65:     for (s = 0; s < supportSize; ++s) {
     66:       const PetscInt  cell = support[s];
     67:       const PetscInt *cone;
     68:       PetscInt        coneSize, c;
     69:       PetscInt        cval, eval, maxedge;
     70:       PetscReal       len, maxlen;

     72:       PetscCall(DMLabelGetValue(sbr->splitPoints, cell, &cval));
     73:       if (cval == 2) continue;
     74:       PetscCall(DMPlexGetCone(dm, cell, &cone));
     75:       PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
     76:       PetscCall(SBRGetEdgeLen_Private(tr, cone[0], &maxlen));
     77:       maxedge = cone[0];
     78:       for (c = 1; c < coneSize; ++c) {
     79:         PetscCall(SBRGetEdgeLen_Private(tr, cone[c], &len));
     80:         if (len > maxlen) {
     81:           maxlen  = len;
     82:           maxedge = cone[c];
     83:         }
     84:       }
     85:       PetscCall(DMLabelGetValue(sbr->splitPoints, maxedge, &eval));
     86:       if (eval != 1) {
     87:         PetscCall(DMLabelSetValue(sbr->splitPoints, maxedge, 1));
     88:         PetscCall(DMPlexPointQueueEnqueue(queue, maxedge));
     89:       }
     90:       PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 2));
     91:     }
     92:   }
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: static PetscErrorCode splitPoint(PETSC_UNUSED DMLabel label, PetscInt p, PETSC_UNUSED PetscInt val, PetscCtx ctx)
     97: {
     98:   DMPlexPointQueue queue = (DMPlexPointQueue)ctx;

    100:   PetscFunctionBegin;
    101:   PetscCall(DMPlexPointQueueEnqueue(queue, p));
    102:   PetscFunctionReturn(PETSC_SUCCESS);
    103: }

    105: /*
    106:   The 'splitPoints' label marks mesh points to be divided. It marks edges with 1, triangles with 2, and tetrahedra with 3.
    107:   Then the refinement type is calculated as follows:

    109:     vertex:                   0
    110:     edge unsplit:             1
    111:     edge split:               2
    112:     triangle unsplit:         3
    113:     triangle split all edges: 4
    114:     triangle split edges 0 1: 5
    115:     triangle split edges 1 0: 6
    116:     triangle split edges 1 2: 7
    117:     triangle split edges 2 1: 8
    118:     triangle split edges 2 0: 9
    119:     triangle split edges 0 2: 10
    120:     triangle split edge 0:    11
    121:     triangle split edge 1:    12
    122:     triangle split edge 2:    13
    123: */
    124: typedef enum {
    125:   RT_VERTEX,
    126:   RT_EDGE,
    127:   RT_EDGE_SPLIT,
    128:   RT_TRIANGLE,
    129:   RT_TRIANGLE_SPLIT,
    130:   RT_TRIANGLE_SPLIT_01,
    131:   RT_TRIANGLE_SPLIT_10,
    132:   RT_TRIANGLE_SPLIT_12,
    133:   RT_TRIANGLE_SPLIT_21,
    134:   RT_TRIANGLE_SPLIT_20,
    135:   RT_TRIANGLE_SPLIT_02,
    136:   RT_TRIANGLE_SPLIT_0,
    137:   RT_TRIANGLE_SPLIT_1,
    138:   RT_TRIANGLE_SPLIT_2
    139: } RefinementType;

    141: static PetscErrorCode DMPlexTransformSetUp_SBR(DMPlexTransform tr)
    142: {
    143:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;
    144:   DM                dm;
    145:   DMLabel           active;
    146:   PetscSF           pointSF;
    147:   DMPlexPointQueue  queue = NULL;
    148:   IS                refineIS;
    149:   const PetscInt   *refineCells;
    150:   PetscInt          pStart, pEnd, p, eStart, eEnd, e, edgeLenSize, Nc, c;
    151:   PetscBool         empty;

    153:   PetscFunctionBegin;
    154:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    155:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Split Points", &sbr->splitPoints));
    156:   /* Create edge lengths */
    157:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
    158:   PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
    159:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &sbr->secEdgeLen));
    160:   PetscCall(PetscSectionSetChart(sbr->secEdgeLen, eStart, eEnd));
    161:   for (e = eStart; e < eEnd; ++e) PetscCall(PetscSectionSetDof(sbr->secEdgeLen, e, 1));
    162:   PetscCall(PetscSectionSetUp(sbr->secEdgeLen));
    163:   PetscCall(PetscSectionGetStorageSize(sbr->secEdgeLen, &edgeLenSize));
    164:   PetscCall(PetscCalloc1(edgeLenSize, &sbr->edgeLen));
    165:   /* Add edges of cells that are marked for refinement to edge queue */
    166:   PetscCall(DMPlexTransformGetActive(tr, &active));
    167:   PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONGSTATE, "DMPlexTransform must have an adaptation label in order to use SBR algorithm");
    168:   PetscCall(DMPlexPointQueueCreate(1024, &queue));
    169:   PetscCall(DMLabelGetStratumIS(active, DM_ADAPT_REFINE, &refineIS));
    170:   PetscCall(DMLabelGetStratumSize(active, DM_ADAPT_REFINE, &Nc));
    171:   if (refineIS) PetscCall(ISGetIndices(refineIS, &refineCells));
    172:   for (c = 0; c < Nc; ++c) {
    173:     const PetscInt cell = refineCells[c];
    174:     PetscInt       depth;

    176:     PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
    177:     if (depth == 1) {
    178:       PetscCall(DMLabelSetValue(sbr->splitPoints, cell, 1));
    179:       PetscCall(DMPlexPointQueueEnqueue(queue, cell));
    180:     } else {
    181:       PetscInt *closure = NULL;
    182:       PetscInt  Ncl, cl;

    184:       PetscCall(DMLabelSetValue(sbr->splitPoints, cell, depth));
    185:       PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
    186:       for (cl = 0; cl < Ncl; cl += 2) {
    187:         const PetscInt edge = closure[cl];

    189:         if (edge >= eStart && edge < eEnd) {
    190:           PetscCall(DMLabelSetValue(sbr->splitPoints, edge, 1));
    191:           PetscCall(DMPlexPointQueueEnqueue(queue, edge));
    192:         }
    193:       }
    194:       PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
    195:     }
    196:   }
    197:   if (refineIS) PetscCall(ISRestoreIndices(refineIS, &refineCells));
    198:   PetscCall(ISDestroy(&refineIS));
    199:   /* Setup communication */
    200:   PetscCall(DMGetPointSF(dm, &pointSF));
    201:   PetscCall(DMLabelPropagateBegin(sbr->splitPoints, pointSF));
    202:   /* While edge queue is not empty: */
    203:   PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
    204:   while (!empty) {
    205:     PetscCall(SBRSplitLocalEdges_Private(tr, queue));
    206:     /* Communicate marked edges
    207:          An easy implementation is to allocate an array the size of the number of points. We put the splitPoints marks into the
    208:          array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.

    210:          TODO: We could use in-place communication with a different SF
    211:            We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
    212:            already been marked. If not, it might have been handled on the process in this round, but we add it anyway.

    214:            In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
    215:            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
    216:            edge to the queue.
    217:     */
    218:     PetscCall(DMLabelPropagatePush(sbr->splitPoints, pointSF, splitPoint, queue));
    219:     PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
    220:   }
    221:   PetscCall(DMLabelPropagateEnd(sbr->splitPoints, pointSF));
    222:   /* Calculate refineType for each cell */
    223:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
    224:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    225:   for (p = pStart; p < pEnd; ++p) {
    226:     DMLabel        trType = tr->trType;
    227:     DMPolytopeType ct;
    228:     PetscInt       val;

    230:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    231:     switch (ct) {
    232:     case DM_POLYTOPE_POINT:
    233:       PetscCall(DMLabelSetValue(trType, p, RT_VERTEX));
    234:       break;
    235:     case DM_POLYTOPE_SEGMENT:
    236:       PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
    237:       if (val == 1) PetscCall(DMLabelSetValue(trType, p, RT_EDGE_SPLIT));
    238:       else PetscCall(DMLabelSetValue(trType, p, RT_EDGE));
    239:       break;
    240:     case DM_POLYTOPE_TRIANGLE:
    241:       PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
    242:       if (val == 2) {
    243:         const PetscInt *cone;
    244:         PetscReal       lens[3];
    245:         PetscInt        vals[3], i;

    247:         PetscCall(DMPlexGetCone(dm, p, &cone));
    248:         for (i = 0; i < 3; ++i) {
    249:           PetscCall(DMLabelGetValue(sbr->splitPoints, cone[i], &vals[i]));
    250:           vals[i] = vals[i] < 0 ? 0 : vals[i];
    251:           PetscCall(SBRGetEdgeLen_Private(tr, cone[i], &lens[i]));
    252:         }
    253:         if (vals[0] && vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT));
    254:         else if (vals[0] && vals[1]) PetscCall(DMLabelSetValue(trType, p, lens[0] > lens[1] ? RT_TRIANGLE_SPLIT_01 : RT_TRIANGLE_SPLIT_10));
    255:         else if (vals[1] && vals[2]) PetscCall(DMLabelSetValue(trType, p, lens[1] > lens[2] ? RT_TRIANGLE_SPLIT_12 : RT_TRIANGLE_SPLIT_21));
    256:         else if (vals[2] && vals[0]) PetscCall(DMLabelSetValue(trType, p, lens[2] > lens[0] ? RT_TRIANGLE_SPLIT_20 : RT_TRIANGLE_SPLIT_02));
    257:         else if (vals[0]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_0));
    258:         else if (vals[1]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_1));
    259:         else if (vals[2]) PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE_SPLIT_2));
    260:         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]);
    261:       } else PetscCall(DMLabelSetValue(trType, p, RT_TRIANGLE));
    262:       break;
    263:     default:
    264:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
    265:     }
    266:     PetscCall(DMLabelGetValue(sbr->splitPoints, p, &val));
    267:   }
    268:   /* Cleanup */
    269:   PetscCall(DMPlexPointQueueDestroy(&queue));
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: static PetscErrorCode DMPlexTransformGetSubcellOrientation_SBR(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
    274: {
    275:   PetscInt rt;

    277:   PetscFunctionBeginHot;
    278:   PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
    279:   *rnew = r;
    280:   *onew = o;
    281:   switch (rt) {
    282:   case RT_TRIANGLE_SPLIT_01:
    283:   case RT_TRIANGLE_SPLIT_10:
    284:   case RT_TRIANGLE_SPLIT_12:
    285:   case RT_TRIANGLE_SPLIT_21:
    286:   case RT_TRIANGLE_SPLIT_20:
    287:   case RT_TRIANGLE_SPLIT_02:
    288:     switch (tct) {
    289:     case DM_POLYTOPE_SEGMENT:
    290:       break;
    291:     case DM_POLYTOPE_TRIANGLE:
    292:       break;
    293:     default:
    294:       break;
    295:     }
    296:     break;
    297:   case RT_TRIANGLE_SPLIT_0:
    298:   case RT_TRIANGLE_SPLIT_1:
    299:   case RT_TRIANGLE_SPLIT_2:
    300:     switch (tct) {
    301:     case DM_POLYTOPE_SEGMENT:
    302:       break;
    303:     case DM_POLYTOPE_TRIANGLE:
    304:       *onew = so < 0 ? -(o + 1) : o;
    305:       *rnew = so < 0 ? (r + 1) % 2 : r;
    306:       break;
    307:     default:
    308:       break;
    309:     }
    310:     break;
    311:   case RT_EDGE_SPLIT:
    312:   case RT_TRIANGLE_SPLIT:
    313:     PetscCall(DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew));
    314:     break;
    315:   default:
    316:     PetscCall(DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew));
    317:   }
    318:   PetscFunctionReturn(PETSC_SUCCESS);
    319: }

    321: /* Add 1 edge inside this triangle, making 2 new triangles.
    322:  2
    323:  |\
    324:  | \
    325:  |  \
    326:  |   \
    327:  |    1
    328:  |     \
    329:  |  B   \
    330:  2       1
    331:  |      / \
    332:  | ____/   0
    333:  |/    A    \
    334:  0-----0-----1
    335: */
    336: static PetscErrorCode SBRGetTriangleSplitSingle(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    337: {
    338:   const PetscInt       *arr     = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
    339:   static DMPolytopeType triT1[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
    340:   static PetscInt       triS1[] = {1, 2};
    341:   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,
    342:                                    DM_POLYTOPE_SEGMENT, 0, 0};
    343:   static PetscInt       triO1[] = {0, 0, 0, 0, -1, 0, 0, 0};

    345:   PetscFunctionBeginHot;
    346:   /* To get the other divisions, we reorient the triangle */
    347:   triC1[2]  = arr[0 * 2];
    348:   triC1[7]  = arr[1 * 2];
    349:   triC1[11] = arr[0 * 2];
    350:   triC1[15] = arr[1 * 2];
    351:   triC1[22] = arr[1 * 2];
    352:   triC1[26] = arr[2 * 2];
    353:   *Nt       = 2;
    354:   *target   = triT1;
    355:   *size     = triS1;
    356:   *cone     = triC1;
    357:   *ornt     = triO1;
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: /* Add 2 edges inside this triangle, making 3 new triangles.
    362:  RT_TRIANGLE_SPLIT_12
    363:  2
    364:  |\
    365:  | \
    366:  |  \
    367:  0   \
    368:  |    1
    369:  |     \
    370:  |  B   \
    371:  2-------1
    372:  |   C  / \
    373:  1 ____/   0
    374:  |/    A    \
    375:  0-----0-----1
    376:  RT_TRIANGLE_SPLIT_10
    377:  2
    378:  |\
    379:  | \
    380:  |  \
    381:  0   \
    382:  |    1
    383:  |     \
    384:  |  A   \
    385:  2       1
    386:  |      /|\
    387:  1 ____/ / 0
    388:  |/ C   / B \
    389:  0-----0-----1
    390:  RT_TRIANGLE_SPLIT_20
    391:  2
    392:  |\
    393:  | \
    394:  |  \
    395:  0   \
    396:  |    \
    397:  |     \
    398:  |      \
    399:  2   A   1
    400:  |\       \
    401:  1 ---\    \
    402:  |B \_C----\\
    403:  0-----0-----1
    404:  RT_TRIANGLE_SPLIT_21
    405:  2
    406:  |\
    407:  | \
    408:  |  \
    409:  0   \
    410:  |    \
    411:  |  B  \
    412:  |      \
    413:  2-------1
    414:  |\     C \
    415:  1 ---\    \
    416:  |  A  ----\\
    417:  0-----0-----1
    418:  RT_TRIANGLE_SPLIT_01
    419:  2
    420:  |\
    421:  |\\
    422:  || \
    423:  | \ \
    424:  |  | \
    425:  |  |  \
    426:  |  |   \
    427:  2   \ C 1
    428:  |  A | / \
    429:  |    | |B \
    430:  |     \/   \
    431:  0-----0-----1
    432:  RT_TRIANGLE_SPLIT_02
    433:  2
    434:  |\
    435:  |\\
    436:  || \
    437:  | \ \
    438:  |  | \
    439:  |  |  \
    440:  |  |   \
    441:  2 C \   1
    442:  |\   |   \
    443:  | \__|  A \
    444:  | B  \\    \
    445:  0-----0-----1
    446: */
    447: static PetscErrorCode SBRGetTriangleSplitDouble(PetscInt o, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    448: {
    449:   PetscInt              e0, e1;
    450:   const PetscInt       *arr     = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, o);
    451:   static DMPolytopeType triT2[] = {DM_POLYTOPE_SEGMENT, DM_POLYTOPE_TRIANGLE};
    452:   static PetscInt       triS2[] = {2, 3};
    453:   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};
    454:   static PetscInt triO2[] = {0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0};

    456:   PetscFunctionBeginHot;
    457:   /* To get the other divisions, we reorient the triangle */
    458:   triC2[2]  = arr[0 * 2];
    459:   triC2[3]  = arr[0 * 2 + 1] ? 1 : 0;
    460:   triC2[7]  = arr[1 * 2];
    461:   triC2[11] = arr[1 * 2];
    462:   triC2[15] = arr[2 * 2];
    463:   /* Swap the first two edges if the triangle is reversed */
    464:   e0            = o < 0 ? 23 : 19;
    465:   e1            = o < 0 ? 19 : 23;
    466:   triC2[e0]     = arr[0 * 2];
    467:   triC2[e0 + 1] = 0;
    468:   triC2[e1]     = arr[1 * 2];
    469:   triC2[e1 + 1] = o < 0 ? 1 : 0;
    470:   triO2[6]      = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
    471:   /* Swap the first two edges if the triangle is reversed */
    472:   e0            = o < 0 ? 34 : 30;
    473:   e1            = o < 0 ? 30 : 34;
    474:   triC2[e0]     = arr[1 * 2];
    475:   triC2[e0 + 1] = o < 0 ? 0 : 1;
    476:   triC2[e1]     = arr[2 * 2];
    477:   triC2[e1 + 1] = o < 0 ? 1 : 0;
    478:   triO2[9]      = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, -1, arr[2 * 2 + 1]);
    479:   /* Swap the last two edges if the triangle is reversed */
    480:   triC2[41] = arr[2 * 2];
    481:   triC2[42] = o < 0 ? 0 : 1;
    482:   triC2[45] = o < 0 ? 1 : 0;
    483:   triC2[48] = o < 0 ? 0 : 1;
    484:   triO2[11] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[1 * 2 + 1]);
    485:   triO2[12] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, 0, arr[2 * 2 + 1]);
    486:   *Nt       = 2;
    487:   *target   = triT2;
    488:   *size     = triS2;
    489:   *cone     = triC2;
    490:   *ornt     = triO2;
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: static PetscErrorCode DMPlexTransformCellTransform_SBR(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    495: {
    496:   DMLabel  trType = tr->trType;
    497:   PetscInt val;

    499:   PetscFunctionBeginHot;
    500:   PetscCheck(p >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point argument is invalid");
    501:   PetscCall(DMLabelGetValue(trType, p, &val));
    502:   if (rt) *rt = val;
    503:   switch (source) {
    504:   case DM_POLYTOPE_POINT:
    505:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
    506:   case DM_POLYTOPE_QUADRILATERAL:
    507:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
    508:   case DM_POLYTOPE_TETRAHEDRON:
    509:   case DM_POLYTOPE_HEXAHEDRON:
    510:   case DM_POLYTOPE_TRI_PRISM:
    511:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
    512:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    513:   case DM_POLYTOPE_PYRAMID:
    514:     PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
    515:     break;
    516:   case DM_POLYTOPE_SEGMENT:
    517:     if (val == RT_EDGE) PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
    518:     else PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
    519:     break;
    520:   case DM_POLYTOPE_TRIANGLE:
    521:     switch (val) {
    522:     case RT_TRIANGLE_SPLIT_0:
    523:       PetscCall(SBRGetTriangleSplitSingle(2, Nt, target, size, cone, ornt));
    524:       break;
    525:     case RT_TRIANGLE_SPLIT_1:
    526:       PetscCall(SBRGetTriangleSplitSingle(0, Nt, target, size, cone, ornt));
    527:       break;
    528:     case RT_TRIANGLE_SPLIT_2:
    529:       PetscCall(SBRGetTriangleSplitSingle(1, Nt, target, size, cone, ornt));
    530:       break;
    531:     case RT_TRIANGLE_SPLIT_21:
    532:       PetscCall(SBRGetTriangleSplitDouble(-3, Nt, target, size, cone, ornt));
    533:       break;
    534:     case RT_TRIANGLE_SPLIT_10:
    535:       PetscCall(SBRGetTriangleSplitDouble(-2, Nt, target, size, cone, ornt));
    536:       break;
    537:     case RT_TRIANGLE_SPLIT_02:
    538:       PetscCall(SBRGetTriangleSplitDouble(-1, Nt, target, size, cone, ornt));
    539:       break;
    540:     case RT_TRIANGLE_SPLIT_12:
    541:       PetscCall(SBRGetTriangleSplitDouble(0, Nt, target, size, cone, ornt));
    542:       break;
    543:     case RT_TRIANGLE_SPLIT_20:
    544:       PetscCall(SBRGetTriangleSplitDouble(1, Nt, target, size, cone, ornt));
    545:       break;
    546:     case RT_TRIANGLE_SPLIT_01:
    547:       PetscCall(SBRGetTriangleSplitDouble(2, Nt, target, size, cone, ornt));
    548:       break;
    549:     case RT_TRIANGLE_SPLIT:
    550:       PetscCall(DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt));
    551:       break;
    552:     default:
    553:       PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
    554:     }
    555:     break;
    556:   default:
    557:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
    558:   }
    559:   PetscFunctionReturn(PETSC_SUCCESS);
    560: }

    562: static PetscErrorCode DMPlexTransformSetFromOptions_SBR(DMPlexTransform tr, PetscOptionItems PetscOptionsObject)
    563: {
    564:   PetscInt  cells[256], n = 256, i;
    565:   PetscBool flg;

    567:   PetscFunctionBegin;
    568:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
    569:   PetscCall(PetscOptionsIntArray("-dm_plex_transform_sbr_ref_cell", "Mark cells for refinement", "", cells, &n, &flg));
    570:   if (flg) {
    571:     DMLabel active;

    573:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active));
    574:     for (i = 0; i < n; ++i) PetscCall(DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE));
    575:     PetscCall(DMPlexTransformSetActive(tr, active));
    576:     PetscCall(DMLabelDestroy(&active));
    577:   }
    578:   PetscOptionsHeadEnd();
    579:   PetscFunctionReturn(PETSC_SUCCESS);
    580: }

    582: static PetscErrorCode DMPlexTransformView_SBR(DMPlexTransform tr, PetscViewer viewer)
    583: {
    584:   PetscBool isascii;

    586:   PetscFunctionBegin;
    589:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    590:   if (isascii) {
    591:     PetscViewerFormat format;
    592:     const char       *name;

    594:     PetscCall(PetscObjectGetName((PetscObject)tr, &name));
    595:     PetscCall(PetscViewerASCIIPrintf(viewer, "SBR refinement %s\n", name ? name : ""));
    596:     PetscCall(PetscViewerGetFormat(viewer, &format));
    597:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(DMLabelView(tr->trType, viewer));
    598:   } else {
    599:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
    600:   }
    601:   PetscFunctionReturn(PETSC_SUCCESS);
    602: }

    604: static PetscErrorCode DMPlexTransformDestroy_SBR(DMPlexTransform tr)
    605: {
    606:   DMPlexRefine_SBR *sbr = (DMPlexRefine_SBR *)tr->data;

    608:   PetscFunctionBegin;
    609:   PetscCall(PetscFree(sbr->edgeLen));
    610:   PetscCall(PetscSectionDestroy(&sbr->secEdgeLen));
    611:   PetscCall(DMLabelDestroy(&sbr->splitPoints));
    612:   PetscCall(PetscFree(tr->data));
    613:   PetscFunctionReturn(PETSC_SUCCESS);
    614: }

    616: static PetscErrorCode DMPlexTransformInitialize_SBR(DMPlexTransform tr)
    617: {
    618:   PetscFunctionBegin;
    619:   tr->ops->view                  = DMPlexTransformView_SBR;
    620:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_SBR;
    621:   tr->ops->setup                 = DMPlexTransformSetUp_SBR;
    622:   tr->ops->destroy               = DMPlexTransformDestroy_SBR;
    623:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Internal;
    624:   tr->ops->celltransform         = DMPlexTransformCellTransform_SBR;
    625:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_SBR;
    626:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
    627:   PetscFunctionReturn(PETSC_SUCCESS);
    628: }

    630: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform tr)
    631: {
    632:   DMPlexRefine_SBR *f;

    634:   PetscFunctionBegin;
    636:   PetscCall(PetscNew(&f));
    637:   tr->data = f;

    639:   PetscCall(DMPlexTransformInitialize_SBR(tr));
    640:   PetscCall(PetscCitationsRegister(SBRCitation, &SBRcite));
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }