Actual source code: plextrcohesive.c

      1: #include <petsc/private/dmplextransformimpl.h>

      3: #include <petsc/private/dmlabelimpl.h>

      5: /*
      6:   The cohesive transformation extrudes cells into a mesh from faces along an internal boundary.

      8:   Orientation:

     10:   We will say that a face has a positive and negative side. The positive side is defined by the cell which attaches the face with a positive orientation, and the negative side cell attaches it with a negative orientation (a reflection). However, this means that the positive side is in the opposite direction of the face normal, and the negative side is in the direction of the face normal, since all cells have outward facing normals. For clarity, in 2D the cross product of the normal and the edge is in the positive z direction.

     12:   Labeling:

     14:   We require an active label on input, which marks all points on the internal surface. Each point is
     15:   labeled with its depth. This label is passed to DMPlexLabelCohesiveComplete(), which adds all points
     16:   which ``impinge'' on the surface, meaning a point has a face on the surface. These points are labeled
     17:   with celltype + 100 on the positive side, and -(celltype + 100) on the negative side.

     19:   Point Creation:

     21:   We split points on the fault surface, creating a new partner point for each one. The negative side
     22:   receives the old point, while the positive side receives the new partner. In addition, points are
     23:   created with the two split points as boundaries. For example, split vertices have a segment between
     24:   them, split edges a quadrilaterial, split triangles a prism, and split quads a hexahedron. By
     25:   default, these spanning points have tensor ordering, but the user can choose to have them use the
     26:   outward normal convention instead.

     28: */

     30: static PetscErrorCode DMPlexTransformView_Cohesive(DMPlexTransform tr, PetscViewer viewer)
     31: {
     32:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
     33:   PetscBool                 isascii;

     35:   PetscFunctionBegin;
     38:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     39:   if (isascii) {
     40:     const char *name;

     42:     PetscCall(PetscObjectGetName((PetscObject)tr, &name));
     43:     PetscCall(PetscViewerASCIIPrintf(viewer, "Cohesive extrusion transformation %s\n", name ? name : ""));
     44:     PetscCall(PetscViewerASCIIPrintf(viewer, "  create tensor cells: %s\n", ex->useTensor ? "YES" : "NO"));
     45:   } else {
     46:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
     47:   }
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: static PetscErrorCode DMPlexTransformSetFromOptions_Cohesive(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
     52: {
     53:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
     54:   PetscReal                 width;
     55:   PetscBool                 tensor, flg;

     57:   PetscFunctionBegin;
     58:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Cohesive Extrusion Options");
     59:   PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg));
     60:   if (flg) PetscCall(DMPlexTransformCohesiveExtrudeSetTensor(tr, tensor));
     61:   PetscCall(PetscOptionsReal("-dm_plex_transform_cohesive_width", "Width of a cohesive cell", "", ex->width, &width, &flg));
     62:   if (flg) PetscCall(DMPlexTransformCohesiveExtrudeSetWidth(tr, width));
     63:   PetscCall(PetscOptionsInt("-dm_plex_transform_cohesive_debug", "Det debugging level", "", ex->debug, &ex->debug, NULL));
     64:   PetscOptionsHeadEnd();
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: /*
     69:   ComputeSplitFaceNumber - Compute an encoding describing which faces of p are split by the surface

     71:   Not collective

     73:   Input Parameters:
     74:   + dm    - The `DM`
     75:   . label - `DMLabel` marking the surface and adjacent points
     76:   - p     - Impinging point, adjacent to the surface

     78:   Output Parameter:
     79:   . fsplit - A number encoding the faces which are split by the surface

     81:   Level: developer

     83:   Note: We will use a bit encoding, where bit k is 1 if face k is split.

     85: .seealso: ComputeUnsplitFaceNumber()
     86: */
     87: static PetscErrorCode ComputeSplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *fsplit)
     88: {
     89:   const PetscInt *cone;
     90:   PetscInt        coneSize, val;

     92:   PetscFunctionBegin;
     93:   *fsplit = 0;
     94:   PetscCall(DMPlexGetCone(dm, p, &cone));
     95:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
     96:   PetscCheck(coneSize < (PetscInt)sizeof(*fsplit) * 8, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cone of size %" PetscInt_FMT " is too large to be contained in an integer", coneSize);
     97:   for (PetscInt c = 0; c < coneSize; ++c) {
     98:     PetscCall(DMLabelGetValue(label, cone[c], &val));
     99:     if (val >= 0 && val < 100) *fsplit |= 1 << c;
    100:   }
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: /*
    105:   ComputeUnsplitFaceNumber - Compute an encoding describing which faces of p are unsplit by the surface

    107:   Not collective

    109:   Input Parameters:
    110:   + dm    - The `DM`
    111:   . label - `DMLabel` marking the surface and adjacent points
    112:   - p     - Split point, on the surface

    114:   Output Parameter:
    115:   . funsplit - A number encoding the faces which are split by the surface

    117:   Level: developer

    119:   Note: We will use a bit encoding, where bit k is 1 if face k is unsplit.

    121: .seealso: ComputeSplitFaceNumber()
    122: */
    123: static PetscErrorCode ComputeUnsplitFaceNumber(DM dm, DMLabel label, PetscInt p, PetscInt *funsplit)
    124: {
    125:   const PetscInt *cone;
    126:   PetscInt        coneSize, val;

    128:   PetscFunctionBegin;
    129:   *funsplit = 0;
    130:   PetscCall(DMPlexGetCone(dm, p, &cone));
    131:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    132:   PetscCheck(coneSize < (PetscInt)sizeof(*funsplit) * 8, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cone of size %" PetscInt_FMT " is too large to be contained in an integer", coneSize);
    133:   for (PetscInt c = 0; c < coneSize; ++c) {
    134:     PetscCall(DMLabelGetValue(label, cone[c], &val));
    135:     if (val >= 200) *funsplit |= 1 << c;
    136:   }
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: /* DM_POLYTOPE_POINT produces
    141:      2 points when split, or 1 point when unsplit, and
    142:      1 segment, or tensor segment
    143: */
    144: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Point(DMPlexTransform_Cohesive *ex)
    145: {
    146:   PetscInt rt, Nc, No;

    148:   PetscFunctionBegin;
    149:   // Unsplit vertex
    150:   rt         = DM_POLYTOPE_POINT * 2 + 1;
    151:   ex->Nt[rt] = 2;
    152:   Nc         = 6;
    153:   No         = 2;
    154:   PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    155:   ex->target[rt][0] = DM_POLYTOPE_POINT;
    156:   ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
    157:   ex->size[rt][0]   = 1;
    158:   ex->size[rt][1]   = 1;
    159:   //   cone for segment/tensor segment
    160:   ex->cone[rt][0] = DM_POLYTOPE_POINT;
    161:   ex->cone[rt][1] = 0;
    162:   ex->cone[rt][2] = 0;
    163:   ex->cone[rt][3] = DM_POLYTOPE_POINT;
    164:   ex->cone[rt][4] = 0;
    165:   ex->cone[rt][5] = 0;
    166:   for (PetscInt i = 0; i < No; ++i) ex->ornt[rt][i] = 0;
    167:   // Split vertex
    168:   rt         = (DM_POLYTOPE_POINT * 2 + 1) * 100 + 0;
    169:   ex->Nt[rt] = 2;
    170:   Nc         = 6;
    171:   No         = 2;
    172:   PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    173:   ex->target[rt][0] = DM_POLYTOPE_POINT;
    174:   ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT;
    175:   ex->size[rt][0]   = 2;
    176:   ex->size[rt][1]   = 1;
    177:   //   cone for segment/tensor segment
    178:   ex->cone[rt][0] = DM_POLYTOPE_POINT;
    179:   ex->cone[rt][1] = 0;
    180:   ex->cone[rt][2] = 0;
    181:   ex->cone[rt][3] = DM_POLYTOPE_POINT;
    182:   ex->cone[rt][4] = 0;
    183:   ex->cone[rt][5] = 1;
    184:   for (PetscInt i = 0; i < No; ++i) ex->ornt[rt][i] = 0;
    185:   PetscFunctionReturn(PETSC_SUCCESS);
    186: }

    188: /* DM_POLYTOPE_SEGMENT produces
    189:      2 segments when split, or 1 segment when unsplit, and
    190:      1 quad, or tensor quad
    191: */
    192: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Segment(DMPlexTransform_Cohesive *ex)
    193: {
    194:   PetscInt rt, Nc, No, coff, ooff;

    196:   PetscFunctionBegin;
    197:   // Unsplit segment
    198:   rt         = DM_POLYTOPE_SEGMENT * 2 + 1;
    199:   ex->Nt[rt] = 2;
    200:   Nc         = 8 + 14;
    201:   No         = 2 + 4;
    202:   PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    203:   ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
    204:   ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
    205:   ex->size[rt][0]   = 1;
    206:   ex->size[rt][1]   = 1;
    207:   //   cones for segment
    208:   ex->cone[rt][0] = DM_POLYTOPE_POINT;
    209:   ex->cone[rt][1] = 1;
    210:   ex->cone[rt][2] = 0;
    211:   ex->cone[rt][3] = 0;
    212:   ex->cone[rt][4] = DM_POLYTOPE_POINT;
    213:   ex->cone[rt][5] = 1;
    214:   ex->cone[rt][6] = 1;
    215:   ex->cone[rt][7] = 0;
    216:   for (PetscInt i = 0; i < 2; ++i) ex->ornt[rt][i] = 0;
    217:   //   cone for quad/tensor quad
    218:   coff = 8;
    219:   ooff = 2;
    220:   if (ex->useTensor) {
    221:     ex->cone[rt][coff + 0]  = DM_POLYTOPE_SEGMENT;
    222:     ex->cone[rt][coff + 1]  = 0;
    223:     ex->cone[rt][coff + 2]  = 0;
    224:     ex->cone[rt][coff + 3]  = DM_POLYTOPE_SEGMENT;
    225:     ex->cone[rt][coff + 4]  = 0;
    226:     ex->cone[rt][coff + 5]  = 0;
    227:     ex->cone[rt][coff + 6]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
    228:     ex->cone[rt][coff + 7]  = 1;
    229:     ex->cone[rt][coff + 8]  = 0;
    230:     ex->cone[rt][coff + 9]  = 0;
    231:     ex->cone[rt][coff + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
    232:     ex->cone[rt][coff + 11] = 1;
    233:     ex->cone[rt][coff + 12] = 1;
    234:     ex->cone[rt][coff + 13] = 0;
    235:     ex->ornt[rt][ooff + 0]  = 0;
    236:     ex->ornt[rt][ooff + 1]  = 0;
    237:     ex->ornt[rt][ooff + 2]  = 0;
    238:     ex->ornt[rt][ooff + 3]  = 0;
    239:   } else {
    240:     ex->cone[rt][coff + 0]  = DM_POLYTOPE_SEGMENT;
    241:     ex->cone[rt][coff + 1]  = 0;
    242:     ex->cone[rt][coff + 2]  = 0;
    243:     ex->cone[rt][coff + 3]  = DM_POLYTOPE_SEGMENT;
    244:     ex->cone[rt][coff + 4]  = 1;
    245:     ex->cone[rt][coff + 5]  = 1;
    246:     ex->cone[rt][coff + 6]  = 0;
    247:     ex->cone[rt][coff + 7]  = DM_POLYTOPE_SEGMENT;
    248:     ex->cone[rt][coff + 8]  = 0;
    249:     ex->cone[rt][coff + 9]  = 0;
    250:     ex->cone[rt][coff + 10] = DM_POLYTOPE_SEGMENT;
    251:     ex->cone[rt][coff + 11] = 1;
    252:     ex->cone[rt][coff + 12] = 0;
    253:     ex->cone[rt][coff + 13] = 0;
    254:     ex->ornt[rt][ooff + 0]  = 0;
    255:     ex->ornt[rt][ooff + 1]  = 0;
    256:     ex->ornt[rt][ooff + 2]  = -1;
    257:     ex->ornt[rt][ooff + 3]  = -1;
    258:   }
    259:   // Split segment
    260:   //   0: no unsplit vertex
    261:   //   1: unsplit vertex 0
    262:   //   2: unsplit vertex 1
    263:   //   3: both vertices unsplit (impossible)
    264:   for (PetscInt s = 0; s < 3; ++s) {
    265:     rt         = (DM_POLYTOPE_SEGMENT * 2 + 1) * 100 + s;
    266:     ex->Nt[rt] = 2;
    267:     Nc         = 8 * 2 + 14;
    268:     No         = 2 * 2 + 4;
    269:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    270:     ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
    271:     ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL;
    272:     ex->size[rt][0]   = 2;
    273:     ex->size[rt][1]   = 1;
    274:     //   cones for segments
    275:     for (PetscInt i = 0; i < 2; ++i) {
    276:       ex->cone[rt][8 * i + 0] = DM_POLYTOPE_POINT;
    277:       ex->cone[rt][8 * i + 1] = 1;
    278:       ex->cone[rt][8 * i + 2] = 0;
    279:       ex->cone[rt][8 * i + 3] = s == 1 ? 0 : i;
    280:       ex->cone[rt][8 * i + 4] = DM_POLYTOPE_POINT;
    281:       ex->cone[rt][8 * i + 5] = 1;
    282:       ex->cone[rt][8 * i + 6] = 1;
    283:       ex->cone[rt][8 * i + 7] = s == 2 ? 0 : i;
    284:     }
    285:     for (PetscInt i = 0; i < 2 * 2; ++i) ex->ornt[rt][i] = 0;
    286:     //   cone for quad/tensor quad
    287:     coff = 8 * 2;
    288:     ooff = 2 * 2;
    289:     if (ex->useTensor) {
    290:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_SEGMENT;
    291:       ex->cone[rt][coff + 1]  = 0;
    292:       ex->cone[rt][coff + 2]  = 0;
    293:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_SEGMENT;
    294:       ex->cone[rt][coff + 4]  = 0;
    295:       ex->cone[rt][coff + 5]  = 1;
    296:       ex->cone[rt][coff + 6]  = DM_POLYTOPE_POINT_PRISM_TENSOR;
    297:       ex->cone[rt][coff + 7]  = 1;
    298:       ex->cone[rt][coff + 8]  = 0;
    299:       ex->cone[rt][coff + 9]  = 0;
    300:       ex->cone[rt][coff + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR;
    301:       ex->cone[rt][coff + 11] = 1;
    302:       ex->cone[rt][coff + 12] = 1;
    303:       ex->cone[rt][coff + 13] = 0;
    304:       ex->ornt[rt][ooff + 0]  = 0;
    305:       ex->ornt[rt][ooff + 1]  = 0;
    306:       ex->ornt[rt][ooff + 2]  = 0;
    307:       ex->ornt[rt][ooff + 3]  = 0;
    308:     } else {
    309:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_SEGMENT;
    310:       ex->cone[rt][coff + 1]  = 0;
    311:       ex->cone[rt][coff + 2]  = 0;
    312:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_SEGMENT;
    313:       ex->cone[rt][coff + 4]  = 1;
    314:       ex->cone[rt][coff + 5]  = 1;
    315:       ex->cone[rt][coff + 6]  = 0;
    316:       ex->cone[rt][coff + 7]  = DM_POLYTOPE_SEGMENT;
    317:       ex->cone[rt][coff + 8]  = 0;
    318:       ex->cone[rt][coff + 9]  = 1;
    319:       ex->cone[rt][coff + 10] = DM_POLYTOPE_SEGMENT;
    320:       ex->cone[rt][coff + 11] = 1;
    321:       ex->cone[rt][coff + 12] = 0;
    322:       ex->cone[rt][coff + 13] = 0;
    323:       ex->ornt[rt][ooff + 0]  = 0;
    324:       ex->ornt[rt][ooff + 1]  = 0;
    325:       ex->ornt[rt][ooff + 2]  = -1;
    326:       ex->ornt[rt][ooff + 3]  = -1;
    327:     }
    328:   }
    329:   // Impinging segment
    330:   //   0: no splits (impossible)
    331:   //   1: split vertex 0
    332:   //   2: split vertex 1
    333:   //   3: split both vertices (impossible)
    334:   for (PetscInt s = 1; s < 3; ++s) {
    335:     rt         = (DM_POLYTOPE_SEGMENT * 2 + 0) * 100 + s;
    336:     ex->Nt[rt] = 1;
    337:     Nc         = 8;
    338:     No         = 2;
    339:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    340:     ex->target[rt][0] = DM_POLYTOPE_SEGMENT;
    341:     ex->size[rt][0]   = 1;
    342:     //   cone for segment
    343:     ex->cone[rt][0] = DM_POLYTOPE_POINT;
    344:     ex->cone[rt][1] = 1;
    345:     ex->cone[rt][2] = 0;
    346:     ex->cone[rt][3] = s == 1 ? 1 : 0;
    347:     ex->cone[rt][4] = DM_POLYTOPE_POINT;
    348:     ex->cone[rt][5] = 1;
    349:     ex->cone[rt][6] = 1;
    350:     ex->cone[rt][7] = s == 2 ? 1 : 0;
    351:     for (PetscInt i = 0; i < 2; ++i) ex->ornt[rt][i] = 0;
    352:   }
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    356: /* DM_POLYTOPE_TRIANGLE produces
    357:      2 triangles, and
    358:      1 triangular prism/tensor triangular prism
    359: */
    360: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Triangle(DMPlexTransform_Cohesive *ex)
    361: {
    362:   PetscInt rt, Nc, No, coff, ooff;

    364:   PetscFunctionBegin;
    365:   // No unsplit triangles
    366:   // Split triangles
    367:   //   0: no unsplit edge
    368:   //   1: unsplit edge 0
    369:   //   2: unsplit edge 1
    370:   //   3: unsplit edge 0 1
    371:   //   4: unsplit edge 2
    372:   //   5: unsplit edge 0 2
    373:   //   6: unsplit edge 1 2
    374:   //   7: all edges unsplit (impossible)
    375:   for (PetscInt s = 0; s < 7; ++s) {
    376:     rt         = (DM_POLYTOPE_TRIANGLE * 2 + 1) * 100 + s;
    377:     ex->Nt[rt] = 2;
    378:     Nc         = 12 * 2 + 18;
    379:     No         = 3 * 2 + 5;
    380:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    381:     ex->target[rt][0] = DM_POLYTOPE_TRIANGLE;
    382:     ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM;
    383:     ex->size[rt][0]   = 2;
    384:     ex->size[rt][1]   = 1;
    385:     //   cones for triangles
    386:     for (PetscInt i = 0; i < 2; ++i) {
    387:       ex->cone[rt][12 * i + 0]  = DM_POLYTOPE_SEGMENT;
    388:       ex->cone[rt][12 * i + 1]  = 1;
    389:       ex->cone[rt][12 * i + 2]  = 0;
    390:       ex->cone[rt][12 * i + 3]  = s & 1 ? 0 : i;
    391:       ex->cone[rt][12 * i + 4]  = DM_POLYTOPE_SEGMENT;
    392:       ex->cone[rt][12 * i + 5]  = 1;
    393:       ex->cone[rt][12 * i + 6]  = 1;
    394:       ex->cone[rt][12 * i + 7]  = s & 2 ? 0 : i;
    395:       ex->cone[rt][12 * i + 8]  = DM_POLYTOPE_SEGMENT;
    396:       ex->cone[rt][12 * i + 9]  = 1;
    397:       ex->cone[rt][12 * i + 10] = 2;
    398:       ex->cone[rt][12 * i + 11] = s & 4 ? 0 : i;
    399:     }
    400:     for (PetscInt i = 0; i < 3 * 2; ++i) ex->ornt[rt][i] = 0;
    401:     //   cone for triangular prism/tensor triangular prism
    402:     coff = 12 * 2;
    403:     ooff = 3 * 2;
    404:     if (ex->useTensor) {
    405:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_TRIANGLE;
    406:       ex->cone[rt][coff + 1]  = 0;
    407:       ex->cone[rt][coff + 2]  = 0;
    408:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_TRIANGLE;
    409:       ex->cone[rt][coff + 4]  = 0;
    410:       ex->cone[rt][coff + 5]  = 1;
    411:       ex->cone[rt][coff + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
    412:       ex->cone[rt][coff + 7]  = 1;
    413:       ex->cone[rt][coff + 8]  = 0;
    414:       ex->cone[rt][coff + 9]  = 0;
    415:       ex->cone[rt][coff + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    416:       ex->cone[rt][coff + 11] = 1;
    417:       ex->cone[rt][coff + 12] = 1;
    418:       ex->cone[rt][coff + 13] = 0;
    419:       ex->cone[rt][coff + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    420:       ex->cone[rt][coff + 15] = 1;
    421:       ex->cone[rt][coff + 16] = 2;
    422:       ex->cone[rt][coff + 17] = 0;
    423:       ex->ornt[rt][ooff + 0]  = 0;
    424:       ex->ornt[rt][ooff + 1]  = 0;
    425:       ex->ornt[rt][ooff + 2]  = 0;
    426:       ex->ornt[rt][ooff + 3]  = 0;
    427:       ex->ornt[rt][ooff + 4]  = 0;
    428:     } else {
    429:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_TRIANGLE;
    430:       ex->cone[rt][coff + 1]  = 0;
    431:       ex->cone[rt][coff + 2]  = 0;
    432:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_TRIANGLE;
    433:       ex->cone[rt][coff + 4]  = 0;
    434:       ex->cone[rt][coff + 5]  = 1;
    435:       ex->cone[rt][coff + 6]  = DM_POLYTOPE_QUADRILATERAL;
    436:       ex->cone[rt][coff + 7]  = 1;
    437:       ex->cone[rt][coff + 8]  = 0;
    438:       ex->cone[rt][coff + 9]  = 0;
    439:       ex->cone[rt][coff + 10] = DM_POLYTOPE_QUADRILATERAL;
    440:       ex->cone[rt][coff + 11] = 1;
    441:       ex->cone[rt][coff + 12] = 1;
    442:       ex->cone[rt][coff + 13] = 0;
    443:       ex->cone[rt][coff + 14] = DM_POLYTOPE_QUADRILATERAL;
    444:       ex->cone[rt][coff + 15] = 1;
    445:       ex->cone[rt][coff + 16] = 2;
    446:       ex->cone[rt][coff + 17] = 0;
    447:       ex->ornt[rt][ooff + 0]  = -2;
    448:       ex->ornt[rt][ooff + 1]  = 0;
    449:       ex->ornt[rt][ooff + 2]  = 0;
    450:       ex->ornt[rt][ooff + 3]  = 0;
    451:       ex->ornt[rt][ooff + 4]  = 0;
    452:     }
    453:   }
    454:   // Impinging triangles
    455:   //   0: no splits (impossible)
    456:   //   1: split edge 0
    457:   //   2: split edge 1
    458:   //   3: split edges 0 and 1
    459:   //   4: split edge 2
    460:   //   5: split edges 0 and 2
    461:   //   6: split edges 1 and 2
    462:   //   7: split all edges (impossible)
    463:   for (PetscInt s = 1; s < 7; ++s) {
    464:     rt         = (DM_POLYTOPE_TRIANGLE * 2 + 0) * 100 + s;
    465:     ex->Nt[rt] = 1;
    466:     Nc         = 12;
    467:     No         = 3;
    468:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    469:     ex->target[rt][0] = DM_POLYTOPE_TRIANGLE;
    470:     ex->size[rt][0]   = 1;
    471:     //   cone for triangle
    472:     ex->cone[rt][0]  = DM_POLYTOPE_SEGMENT;
    473:     ex->cone[rt][1]  = 1;
    474:     ex->cone[rt][2]  = 0;
    475:     ex->cone[rt][3]  = s & 1 ? 1 : 0;
    476:     ex->cone[rt][4]  = DM_POLYTOPE_SEGMENT;
    477:     ex->cone[rt][5]  = 1;
    478:     ex->cone[rt][6]  = 1;
    479:     ex->cone[rt][7]  = s & 2 ? 1 : 0;
    480:     ex->cone[rt][8]  = DM_POLYTOPE_SEGMENT;
    481:     ex->cone[rt][9]  = 1;
    482:     ex->cone[rt][10] = 2;
    483:     ex->cone[rt][11] = s & 4 ? 1 : 0;
    484:     for (PetscInt i = 0; i < 3; ++i) ex->ornt[rt][i] = 0;
    485:   }
    486:   PetscFunctionReturn(PETSC_SUCCESS);
    487: }

    489: /* DM_POLYTOPE_QUADRILATERAL produces
    490:      2 quads, and
    491:      1 hex/tensor hex
    492: */
    493: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(DMPlexTransform_Cohesive *ex)
    494: {
    495:   PetscInt rt, Nc, No, coff, ooff;

    497:   PetscFunctionBegin;
    498:   // No unsplit quadrilaterals
    499:   // Split quadrilateral
    500:   //   0: no unsplit edge
    501:   //   1: unsplit edge 0
    502:   //   2: unsplit edge 1
    503:   //   3: unsplit edge 0 1
    504:   //   4: unsplit edge 2
    505:   //   5: unsplit edge 0 2
    506:   //   6: unsplit edge 1 2
    507:   //   7: unsplit edge 0 1 2
    508:   //   8: unsplit edge 3
    509:   //   9: unsplit edge 0 3
    510:   //  10: unsplit edge 1 3
    511:   //  11: unsplit edge 0 1 3
    512:   //  12: unsplit edge 2 3
    513:   //  13: unsplit edge 0 2 3
    514:   //  14: unsplit edge 1 2 3
    515:   //  15: all edges unsplit (impossible)
    516:   for (PetscInt s = 0; s < 15; ++s) {
    517:     rt         = (DM_POLYTOPE_QUADRILATERAL * 2 + 1) * 100 + s;
    518:     ex->Nt[rt] = 2;
    519:     Nc         = 16 * 2 + 22;
    520:     No         = 4 * 2 + 6;
    521:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    522:     ex->target[rt][0] = DM_POLYTOPE_QUADRILATERAL;
    523:     ex->target[rt][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON;
    524:     ex->size[rt][0]   = 2;
    525:     ex->size[rt][1]   = 1;
    526:     //   cones for quads
    527:     for (PetscInt i = 0; i < 2; ++i) {
    528:       ex->cone[rt][16 * i + 0]  = DM_POLYTOPE_SEGMENT;
    529:       ex->cone[rt][16 * i + 1]  = 1;
    530:       ex->cone[rt][16 * i + 2]  = 0;
    531:       ex->cone[rt][16 * i + 3]  = s & 1 ? 0 : i;
    532:       ex->cone[rt][16 * i + 4]  = DM_POLYTOPE_SEGMENT;
    533:       ex->cone[rt][16 * i + 5]  = 1;
    534:       ex->cone[rt][16 * i + 6]  = 1;
    535:       ex->cone[rt][16 * i + 7]  = s & 2 ? 0 : i;
    536:       ex->cone[rt][16 * i + 8]  = DM_POLYTOPE_SEGMENT;
    537:       ex->cone[rt][16 * i + 9]  = 1;
    538:       ex->cone[rt][16 * i + 10] = 2;
    539:       ex->cone[rt][16 * i + 11] = s & 4 ? 0 : i;
    540:       ex->cone[rt][16 * i + 12] = DM_POLYTOPE_SEGMENT;
    541:       ex->cone[rt][16 * i + 13] = 1;
    542:       ex->cone[rt][16 * i + 14] = 3;
    543:       ex->cone[rt][16 * i + 15] = s & 8 ? 0 : i;
    544:     }
    545:     for (PetscInt i = 0; i < 4 * 2; ++i) ex->ornt[rt][i] = 0;
    546:     //   cones for hexes/tensor hexes
    547:     coff = 16 * 2;
    548:     ooff = 4 * 2;
    549:     if (ex->useTensor) {
    550:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_QUADRILATERAL;
    551:       ex->cone[rt][coff + 1]  = 0;
    552:       ex->cone[rt][coff + 2]  = 0;
    553:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_QUADRILATERAL;
    554:       ex->cone[rt][coff + 4]  = 0;
    555:       ex->cone[rt][coff + 5]  = 1;
    556:       ex->cone[rt][coff + 6]  = DM_POLYTOPE_SEG_PRISM_TENSOR;
    557:       ex->cone[rt][coff + 7]  = 1;
    558:       ex->cone[rt][coff + 8]  = 0;
    559:       ex->cone[rt][coff + 9]  = 0;
    560:       ex->cone[rt][coff + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    561:       ex->cone[rt][coff + 11] = 1;
    562:       ex->cone[rt][coff + 12] = 1;
    563:       ex->cone[rt][coff + 13] = 0;
    564:       ex->cone[rt][coff + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    565:       ex->cone[rt][coff + 15] = 1;
    566:       ex->cone[rt][coff + 16] = 2;
    567:       ex->cone[rt][coff + 17] = 0;
    568:       ex->cone[rt][coff + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR;
    569:       ex->cone[rt][coff + 19] = 1;
    570:       ex->cone[rt][coff + 20] = 3;
    571:       ex->cone[rt][coff + 21] = 0;
    572:       ex->ornt[rt][ooff + 0]  = 0;
    573:       ex->ornt[rt][ooff + 1]  = 0;
    574:       ex->ornt[rt][ooff + 2]  = 0;
    575:       ex->ornt[rt][ooff + 3]  = 0;
    576:       ex->ornt[rt][ooff + 4]  = 0;
    577:       ex->ornt[rt][ooff + 5]  = 0;
    578:     } else {
    579:       ex->cone[rt][coff + 0]  = DM_POLYTOPE_QUADRILATERAL;
    580:       ex->cone[rt][coff + 1]  = 0;
    581:       ex->cone[rt][coff + 2]  = 0;
    582:       ex->cone[rt][coff + 3]  = DM_POLYTOPE_QUADRILATERAL;
    583:       ex->cone[rt][coff + 4]  = 0;
    584:       ex->cone[rt][coff + 5]  = 1;
    585:       ex->cone[rt][coff + 6]  = DM_POLYTOPE_QUADRILATERAL;
    586:       ex->cone[rt][coff + 7]  = 1;
    587:       ex->cone[rt][coff + 8]  = 0;
    588:       ex->cone[rt][coff + 9]  = 0;
    589:       ex->cone[rt][coff + 10] = DM_POLYTOPE_QUADRILATERAL;
    590:       ex->cone[rt][coff + 11] = 1;
    591:       ex->cone[rt][coff + 12] = 2;
    592:       ex->cone[rt][coff + 13] = 0;
    593:       ex->cone[rt][coff + 14] = DM_POLYTOPE_QUADRILATERAL;
    594:       ex->cone[rt][coff + 15] = 1;
    595:       ex->cone[rt][coff + 16] = 1;
    596:       ex->cone[rt][coff + 17] = 0;
    597:       ex->cone[rt][coff + 18] = DM_POLYTOPE_QUADRILATERAL;
    598:       ex->cone[rt][coff + 19] = 1;
    599:       ex->cone[rt][coff + 20] = 3;
    600:       ex->cone[rt][coff + 21] = 0;
    601:       ex->ornt[rt][ooff + 0]  = -2;
    602:       ex->ornt[rt][ooff + 1]  = 0;
    603:       ex->ornt[rt][ooff + 2]  = 0;
    604:       ex->ornt[rt][ooff + 3]  = 0;
    605:       ex->ornt[rt][ooff + 4]  = 0;
    606:       ex->ornt[rt][ooff + 5]  = 1;
    607:     }
    608:   }
    609:   // Impinging quadrilaterals
    610:   //   0:  no splits (impossible)
    611:   //   1:  split edge 0
    612:   //   2:  split edge 1
    613:   //   3:  split edges 0 and 1
    614:   //   4:  split edge 2
    615:   //   5:  split edges 0 and 2
    616:   //   6:  split edges 1 and 2
    617:   //   7:  split edges 0, 1, and 2
    618:   //   8:  split edge 3
    619:   //   9:  split edges 0 and 3
    620:   //   10: split edges 1 and 3
    621:   //   11: split edges 0, 1, and 3
    622:   //   12: split edges 2 and 3
    623:   //   13: split edges 0, 2, and 3
    624:   //   14: split edges 1, 2, and 3
    625:   //   15: split all edges (impossible)
    626:   for (PetscInt s = 1; s < 15; ++s) {
    627:     rt         = (DM_POLYTOPE_QUADRILATERAL * 2 + 0) * 100 + s;
    628:     ex->Nt[rt] = 1;
    629:     Nc         = 16;
    630:     No         = 4;
    631:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    632:     ex->target[rt][0] = DM_POLYTOPE_QUADRILATERAL;
    633:     ex->size[rt][0]   = 1;
    634:     //   cone for quadrilateral
    635:     ex->cone[rt][0]  = DM_POLYTOPE_SEGMENT;
    636:     ex->cone[rt][1]  = 1;
    637:     ex->cone[rt][2]  = 0;
    638:     ex->cone[rt][3]  = s & 1 ? 1 : 0;
    639:     ex->cone[rt][4]  = DM_POLYTOPE_SEGMENT;
    640:     ex->cone[rt][5]  = 1;
    641:     ex->cone[rt][6]  = 1;
    642:     ex->cone[rt][7]  = s & 2 ? 1 : 0;
    643:     ex->cone[rt][8]  = DM_POLYTOPE_SEGMENT;
    644:     ex->cone[rt][9]  = 1;
    645:     ex->cone[rt][10] = 2;
    646:     ex->cone[rt][11] = s & 4 ? 1 : 0;
    647:     ex->cone[rt][12] = DM_POLYTOPE_SEGMENT;
    648:     ex->cone[rt][13] = 1;
    649:     ex->cone[rt][14] = 3;
    650:     ex->cone[rt][15] = s & 8 ? 1 : 0;
    651:     for (PetscInt i = 0; i < 4; ++i) ex->ornt[rt][i] = 0;
    652:   }
    653:   PetscFunctionReturn(PETSC_SUCCESS);
    654: }

    656: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(DMPlexTransform_Cohesive *ex)
    657: {
    658:   PetscInt rt, Nc, No;

    660:   PetscFunctionBegin;
    661:   // Impinging tetrahedra
    662:   //   0:  no splits (impossible)
    663:   //   1:  split face 0
    664:   //   2:  split face 1
    665:   //   3:  split faces 0 and 1
    666:   //   4:  split face 2
    667:   //   5:  split faces 0 and 2
    668:   //   6:  split faces 1 and 2
    669:   //   7:  split faces 0, 1, and 2
    670:   //   8:  split face 3
    671:   //   9:  split faces 0 and 3
    672:   //   10: split faces 1 and 3
    673:   //   11: split faces 0, 1, and 3
    674:   //   12: split faces 2 and 3
    675:   //   13: split faces 0, 2, and 3
    676:   //   14: split faces 1, 2, and 3
    677:   //   15: split all faces (impossible)
    678:   for (PetscInt s = 1; s < 15; ++s) {
    679:     rt         = (DM_POLYTOPE_TETRAHEDRON * 2 + 0) * 100 + s;
    680:     ex->Nt[rt] = 1;
    681:     Nc         = 16;
    682:     No         = 4;
    683:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    684:     ex->target[rt][0] = DM_POLYTOPE_TETRAHEDRON;
    685:     ex->size[rt][0]   = 1;
    686:     //   cone for triangle
    687:     ex->cone[rt][0]  = DM_POLYTOPE_TRIANGLE;
    688:     ex->cone[rt][1]  = 1;
    689:     ex->cone[rt][2]  = 0;
    690:     ex->cone[rt][3]  = s & 1 ? 1 : 0;
    691:     ex->cone[rt][4]  = DM_POLYTOPE_TRIANGLE;
    692:     ex->cone[rt][5]  = 1;
    693:     ex->cone[rt][6]  = 1;
    694:     ex->cone[rt][7]  = s & 2 ? 1 : 0;
    695:     ex->cone[rt][8]  = DM_POLYTOPE_TRIANGLE;
    696:     ex->cone[rt][9]  = 1;
    697:     ex->cone[rt][10] = 2;
    698:     ex->cone[rt][11] = s & 4 ? 1 : 0;
    699:     ex->cone[rt][12] = DM_POLYTOPE_TRIANGLE;
    700:     ex->cone[rt][13] = 1;
    701:     ex->cone[rt][14] = 3;
    702:     ex->cone[rt][15] = s & 8 ? 1 : 0;
    703:     for (PetscInt i = 0; i < 4; ++i) ex->ornt[rt][i] = 0;
    704:   }
    705:   PetscFunctionReturn(PETSC_SUCCESS);
    706: }

    708: static PetscErrorCode DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(DMPlexTransform_Cohesive *ex)
    709: {
    710:   PetscInt rt, Nc, No;

    712:   PetscFunctionBegin;
    713:   // Impinging hexahedra
    714:   //   0:  no splits (impossible)
    715:   //   bit is set if the face is split
    716:   //   63: split all faces (impossible)
    717:   for (PetscInt s = 1; s < 63; ++s) {
    718:     rt         = (DM_POLYTOPE_HEXAHEDRON * 2 + 0) * 100 + s;
    719:     ex->Nt[rt] = 1;
    720:     Nc         = 24;
    721:     No         = 6;
    722:     PetscCall(PetscMalloc4(ex->Nt[rt], &ex->target[rt], ex->Nt[rt], &ex->size[rt], Nc, &ex->cone[rt], No, &ex->ornt[rt]));
    723:     ex->target[rt][0] = DM_POLYTOPE_HEXAHEDRON;
    724:     ex->size[rt][0]   = 1;
    725:     //   cone for hexahedron
    726:     ex->cone[rt][0]  = DM_POLYTOPE_QUADRILATERAL;
    727:     ex->cone[rt][1]  = 1;
    728:     ex->cone[rt][2]  = 0;
    729:     ex->cone[rt][3]  = s & 1 ? 1 : 0;
    730:     ex->cone[rt][4]  = DM_POLYTOPE_QUADRILATERAL;
    731:     ex->cone[rt][5]  = 1;
    732:     ex->cone[rt][6]  = 1;
    733:     ex->cone[rt][7]  = s & 2 ? 1 : 0;
    734:     ex->cone[rt][8]  = DM_POLYTOPE_QUADRILATERAL;
    735:     ex->cone[rt][9]  = 1;
    736:     ex->cone[rt][10] = 2;
    737:     ex->cone[rt][11] = s & 4 ? 1 : 0;
    738:     ex->cone[rt][12] = DM_POLYTOPE_QUADRILATERAL;
    739:     ex->cone[rt][13] = 1;
    740:     ex->cone[rt][14] = 3;
    741:     ex->cone[rt][15] = s & 8 ? 1 : 0;
    742:     ex->cone[rt][16] = DM_POLYTOPE_QUADRILATERAL;
    743:     ex->cone[rt][17] = 1;
    744:     ex->cone[rt][18] = 4;
    745:     ex->cone[rt][19] = s & 16 ? 1 : 0;
    746:     ex->cone[rt][20] = DM_POLYTOPE_QUADRILATERAL;
    747:     ex->cone[rt][21] = 1;
    748:     ex->cone[rt][22] = 5;
    749:     ex->cone[rt][23] = s & 32 ? 1 : 0;
    750:     for (PetscInt i = 0; i < 6; ++i) ex->ornt[rt][i] = 0;
    751:   }
    752:   PetscFunctionReturn(PETSC_SUCCESS);
    753: }

    755: /*
    756:   The refine types for cohesive extrusion are:

    758:   ct * 2 + 0:                  For any point which should just return itself
    759:   ct * 2 + 1:                  For unsplit points
    760:   (ct * 2 + 0) * 100 + fsplit: For impinging points, one type for each combination of split faces
    761:   (ct * 2 + 1) * 100 + fsplit: For split points, one type for each combination of unsplit faces
    762: */
    763: static PetscErrorCode DMPlexTransformSetUp_Cohesive(DMPlexTransform tr)
    764: {
    765:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
    766:   DM                        dm;
    767:   DMLabel                   active, celltype;
    768:   PetscInt                  numRt, pStart, pEnd, ict;

    770:   PetscFunctionBegin;
    771:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    772:   PetscCall(DMPlexTransformGetActive(tr, &active));
    773:   PetscCheck(active, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cohesive extrusion requires an active label");
    774:   PetscCall(DMPlexGetCellTypeLabel(dm, &celltype));
    775:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType));
    776:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    777:   PetscCall(DMLabelMakeAllInvalid_Internal(active));
    778:   for (PetscInt p = pStart; p < pEnd; ++p) {
    779:     PetscInt ct, val;

    781:     PetscCall(DMLabelGetValue(celltype, p, &ct));
    782:     PetscCall(DMLabelGetValue(active, p, &val));
    783:     if (val < 0) {
    784:       // Also negative size impinging points
    785:       // ct * 2 + 0 is the identity transform
    786:       PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
    787:     } else {
    788:       PetscInt fsplit = -1, funsplit = -1;

    790:       // Unsplit points ct * 2 + 1
    791:       if (val >= 200) {
    792:         // Cohesive cells cannot be unsplit
    793:         //   This is faulty inheritance through the label
    794:         if (ct == DM_POLYTOPE_POINT_PRISM_TENSOR || ct == DM_POLYTOPE_SEG_PRISM_TENSOR || ct == DM_POLYTOPE_TRI_PRISM_TENSOR || ct == DM_POLYTOPE_QUAD_PRISM_TENSOR) PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
    795:         else PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 1));
    796:       } else if (val >= 100) {
    797:         // Impinging points: (ct * 2 + 0) * 100 + fsplit
    798:         PetscCall(ComputeSplitFaceNumber(dm, active, p, &fsplit));
    799:         if (!fsplit) PetscCall(DMLabelSetValue(tr->trType, p, ct * 2 + 0));
    800:         else PetscCall(DMLabelSetValue(tr->trType, p, (ct * 2 + 0) * 100 + fsplit));
    801:       } else {
    802:         // Split points: (ct * 2 + 1) * 100 + funsplit
    803:         PetscCall(ComputeUnsplitFaceNumber(dm, active, p, &funsplit));
    804:         PetscCall(DMLabelSetValue(tr->trType, p, (ct * 2 + 1) * 100 + funsplit));
    805:       }
    806:     }
    807:   }
    808:   if (ex->debug) {
    809:     PetscCall(DMLabelView(active, NULL));
    810:     PetscCall(DMLabelView(tr->trType, NULL));
    811:   }
    812:   numRt = DM_NUM_POLYTOPES * 2 * 100;
    813:   PetscCall(PetscMalloc5(numRt, &ex->Nt, numRt, &ex->target, numRt, &ex->size, numRt, &ex->cone, numRt, &ex->ornt));
    814:   for (ict = 0; ict < numRt; ++ict) {
    815:     ex->Nt[ict]     = -1;
    816:     ex->target[ict] = NULL;
    817:     ex->size[ict]   = NULL;
    818:     ex->cone[ict]   = NULL;
    819:     ex->ornt[ict]   = NULL;
    820:   }
    821:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Point(ex));
    822:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Segment(ex));
    823:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Triangle(ex));
    824:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Quadrilateral(ex));
    825:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Tetrahedron(ex));
    826:   PetscCall(DMPlexTransformCohesiveExtrudeSetUp_Hexahedron(ex));
    827:   PetscFunctionReturn(PETSC_SUCCESS);
    828: }

    830: static PetscErrorCode DMPlexTransformDestroy_Cohesive(DMPlexTransform tr)
    831: {
    832:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;
    833:   PetscInt                  ct;

    835:   PetscFunctionBegin;
    836:   if (ex->target) {
    837:     for (ct = 0; ct < DM_NUM_POLYTOPES * 2 * 100; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]));
    838:   }
    839:   PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt));
    840:   PetscCall(PetscFree(ex));
    841:   PetscFunctionReturn(PETSC_SUCCESS);
    842: }

    844: static PetscErrorCode DMPlexTransformGetSubcellOrientation_Cohesive(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
    845: {
    846:   DMPlexTransform_Cohesive *ex     = (DMPlexTransform_Cohesive *)tr->data;
    847:   DMLabel                   trType = tr->trType;
    848:   PetscInt                  rt;

    850:   PetscFunctionBeginHot;
    851:   *rnew = r;
    852:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
    853:   if (!so) PetscFunctionReturn(PETSC_SUCCESS);
    854:   if (trType) {
    855:     PetscCall(DMLabelGetValue(tr->trType, sp, &rt));
    856:     if (rt < 100 && !(rt % 2)) PetscFunctionReturn(PETSC_SUCCESS);
    857:   }
    858:   if (ex->useTensor) {
    859:     switch (sct) {
    860:     case DM_POLYTOPE_POINT:
    861:       break;
    862:     case DM_POLYTOPE_SEGMENT:
    863:       switch (tct) {
    864:       case DM_POLYTOPE_SEGMENT:
    865:         break;
    866:       case DM_POLYTOPE_SEG_PRISM_TENSOR:
    867:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0);
    868:         break;
    869:       default:
    870:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
    871:       }
    872:       break;
    873:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
    874:     case DM_POLYTOPE_TRIANGLE:
    875:       break;
    876:     case DM_POLYTOPE_QUADRILATERAL:
    877:       break;
    878:     default:
    879:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
    880:     }
    881:   } else {
    882:     switch (sct) {
    883:     case DM_POLYTOPE_POINT:
    884:       break;
    885:     case DM_POLYTOPE_SEGMENT:
    886:       switch (tct) {
    887:       case DM_POLYTOPE_SEGMENT:
    888:         break;
    889:       case DM_POLYTOPE_QUADRILATERAL:
    890:         *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0);
    891:         break;
    892:       default:
    893:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]);
    894:       }
    895:       break;
    896:     // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded
    897:     case DM_POLYTOPE_TRIANGLE:
    898:       break;
    899:     case DM_POLYTOPE_QUADRILATERAL:
    900:       break;
    901:     default:
    902:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]);
    903:     }
    904:   }
    905:   PetscFunctionReturn(PETSC_SUCCESS);
    906: }

    908: static PetscErrorCode DMPlexTransformCellTransform_Cohesive(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    909: {
    910:   DMPlexTransform_Cohesive *ex       = (DMPlexTransform_Cohesive *)tr->data;
    911:   DMLabel                   trType   = tr->trType;
    912:   PetscBool                 identity = PETSC_FALSE;
    913:   PetscInt                  val      = 0;

    915:   PetscFunctionBegin;
    916:   PetscCheck(trType, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing transform type label");
    917:   PetscCall(DMLabelGetValue(trType, p, &val));
    918:   identity = val < 100 && !(val % 2) ? PETSC_TRUE : PETSC_FALSE;
    919:   if (rt) *rt = val;
    920:   if (identity) {
    921:     PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt));
    922:   } else {
    923:     *Nt     = ex->Nt[val];
    924:     *target = ex->target[val];
    925:     *size   = ex->size[val];
    926:     *cone   = ex->cone[val];
    927:     *ornt   = ex->ornt[val];
    928:   }
    929:   PetscFunctionReturn(PETSC_SUCCESS);
    930: }

    932: /* New vertices have the same coordinates */
    933: static PetscErrorCode DMPlexTransformMapCoordinates_Cohesive(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
    934: {
    935:   PetscReal width;
    936:   PetscInt  pval;

    938:   PetscFunctionBeginHot;
    939:   PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]);
    940:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
    941:   PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv);
    942:   PetscCheck(r < 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should only have two replicas, not %" PetscInt_FMT, r);

    944:   PetscCall(DMPlexTransformCohesiveExtrudeGetWidth(tr, &width));
    945:   PetscCall(DMLabelGetValue(tr->trType, p, &pval));
    946:   if (width == 0. || pval < 100) {
    947:     for (PetscInt d = 0; d < dE; ++d) out[d] = in[d];
    948:   } else {
    949:     DM        dm;
    950:     PetscReal avgNormal[3] = {0., 0., 0.}, norm = 0.;
    951:     PetscInt *star = NULL;
    952:     PetscInt  Nst, fStart, fEnd, Nf = 0;

    954:     PetscCall(DMPlexTransformGetDM(tr, &dm));
    955:     PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    956:     PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
    957:     // Get support faces that are split, refine type (ct * 2 + 1) * 100 + fsplit
    958:     for (PetscInt st = 0; st < Nst * 2; st += 2) {
    959:       DMPolytopeType ct;
    960:       PetscInt       val;

    962:       if (star[st] < fStart || star[st] >= fEnd) continue;
    963:       PetscCall(DMPlexGetCellType(dm, star[st], &ct));
    964:       PetscCall(DMLabelGetValue(tr->trType, star[st], &val));
    965:       if (val < (PetscInt)(ct * 2 + 1) * 100) continue;
    966:       star[Nf++] = star[st];
    967:     }
    968:     PetscCheck(Nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Split vertex %" PetscInt_FMT " must be connected to at least one split face", p);
    969:     // Average normals
    970:     for (PetscInt f = 0; f < Nf; ++f) {
    971:       PetscReal normal[3], vol;

    973:       PetscCall(DMPlexComputeCellGeometryFVM(dm, star[f], &vol, NULL, normal));
    974:       for (PetscInt d = 0; d < dE; ++d) avgNormal[d] += normal[d];
    975:     }
    976:     PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &Nst, &star));
    977:     // Normalize normal
    978:     for (PetscInt d = 0; d < dE; ++d) norm += PetscSqr(avgNormal[d]);
    979:     norm = PetscSqrtReal(norm);
    980:     for (PetscInt d = 0; d < dE; ++d) avgNormal[d] /= norm;
    981:     // Symmetrically push vertices along normal
    982:     for (PetscInt d = 0; d < dE; ++d) out[d] = in[d] + width * avgNormal[d] * (r ? -0.5 : 0.5);
    983:   }
    984:   PetscFunctionReturn(PETSC_SUCCESS);
    985: }

    987: static PetscErrorCode DMPlexTransformInitialize_Cohesive(DMPlexTransform tr)
    988: {
    989:   PetscFunctionBegin;
    990:   tr->ops->view                  = DMPlexTransformView_Cohesive;
    991:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_Cohesive;
    992:   tr->ops->setup                 = DMPlexTransformSetUp_Cohesive;
    993:   tr->ops->destroy               = DMPlexTransformDestroy_Cohesive;
    994:   tr->ops->setdimensions         = DMPlexTransformSetDimensions_Internal;
    995:   tr->ops->celltransform         = DMPlexTransformCellTransform_Cohesive;
    996:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Cohesive;
    997:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinates_Cohesive;
    998:   PetscFunctionReturn(PETSC_SUCCESS);
    999: }

    1001: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform tr)
    1002: {
    1003:   DMPlexTransform_Cohesive *ex;

    1005:   PetscFunctionBegin;
    1007:   PetscCall(PetscNew(&ex));
    1008:   tr->data      = ex;
    1009:   ex->useTensor = PETSC_TRUE;
    1010:   PetscCall(DMPlexTransformInitialize_Cohesive(tr));
    1011:   PetscFunctionReturn(PETSC_SUCCESS);
    1012: }

    1014: /*@
    1015:   DMPlexTransformCohesiveExtrudeGetTensor - Get the flag to use tensor cells

    1017:   Not Collective

    1019:   Input Parameter:
    1020: . tr - The `DMPlexTransform`

    1022:   Output Parameter:
    1023: . useTensor - The flag to use tensor cells

    1025:   Note:
    1026:   This flag determines the orientation behavior of the created points.

    1028:   For example, if tensor is `PETSC_TRUE`, then
    1029: .vb
    1030:   DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
    1031:   DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
    1032:   DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
    1033:   DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
    1034: .ve

    1036:   Level: intermediate

    1038: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetTensor()`, `DMPlexTransformExtrudeGetTensor()`
    1039: @*/
    1040: PetscErrorCode DMPlexTransformCohesiveExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor)
    1041: {
    1042:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;

    1044:   PetscFunctionBegin;
    1046:   PetscAssertPointer(useTensor, 2);
    1047:   *useTensor = ex->useTensor;
    1048:   PetscFunctionReturn(PETSC_SUCCESS);
    1049: }

    1051: /*@
    1052:   DMPlexTransformCohesiveExtrudeSetTensor - Set the flag to use tensor cells

    1054:   Not Collective

    1056:   Input Parameters:
    1057: + tr        - The `DMPlexTransform`
    1058: - useTensor - The flag for tensor cells

    1060:   Note:
    1061:   This flag determines the orientation behavior of the created points
    1062:   For example, if tensor is `PETSC_TRUE`, then
    1063: .vb
    1064:   DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT,
    1065:   DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL,
    1066:   DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and
    1067:   DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON.
    1068: .ve

    1070:   Level: intermediate

    1072: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetTensor()`, `DMPlexTransformExtrudeSetTensor()`
    1073: @*/
    1074: PetscErrorCode DMPlexTransformCohesiveExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor)
    1075: {
    1076:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;

    1078:   PetscFunctionBegin;
    1080:   ex->useTensor = useTensor;
    1081:   PetscFunctionReturn(PETSC_SUCCESS);
    1082: }

    1084: /*@
    1085:   DMPlexTransformCohesiveExtrudeGetWidth - Get the width of extruded cells

    1087:   Not Collective

    1089:   Input Parameter:
    1090: . tr - The `DMPlexTransform`

    1092:   Output Parameter:
    1093: . width - The width of extruded cells, or 0.

    1095:   Level: intermediate

    1097: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeSetWidth()`
    1098: @*/
    1099: PetscErrorCode DMPlexTransformCohesiveExtrudeGetWidth(DMPlexTransform tr, PetscReal *width)
    1100: {
    1101:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;

    1103:   PetscFunctionBegin;
    1105:   PetscAssertPointer(width, 2);
    1106:   *width = ex->width;
    1107:   PetscFunctionReturn(PETSC_SUCCESS);
    1108: }

    1110: /*@
    1111:   DMPlexTransformCohesiveExtrudeSetWidth - Set the width of extruded cells

    1113:   Not Collective

    1115:   Input Parameters:
    1116: + tr    - The `DMPlexTransform`
    1117: - width - The width of the extruded cells, or 0.

    1119:   Level: intermediate

    1121: .seealso: `DMPlexTransform`, `DMPlexTransformCohesiveExtrudeGetWidth()`
    1122: @*/
    1123: PetscErrorCode DMPlexTransformCohesiveExtrudeSetWidth(DMPlexTransform tr, PetscReal width)
    1124: {
    1125:   DMPlexTransform_Cohesive *ex = (DMPlexTransform_Cohesive *)tr->data;

    1127:   PetscFunctionBegin;
    1129:   ex->width = width;
    1130:   PetscFunctionReturn(PETSC_SUCCESS);
    1131: }