Actual source code: plextransform.c

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

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

      5: PetscClassId DMPLEXTRANSFORM_CLASSID;

      7: PetscFunctionList DMPlexTransformList              = NULL;
      8: PetscBool         DMPlexTransformRegisterAllCalled = PETSC_FALSE;

     10: PetscLogEvent DMPLEXTRANSFORM_SetUp, DMPLEXTRANSFORM_Apply, DMPLEXTRANSFORM_SetConeSizes, DMPLEXTRANSFORM_SetCones, DMPLEXTRANSFORM_CreateSF, DMPLEXTRANSFORM_CreateLabels, DMPLEXTRANSFORM_SetCoordinates;

     12: /* Construct cell type order since we must loop over cell types in the same dimensional order they are stored in the plex if dm != NULL
     13:         OR in standard plex ordering if dm == NULL */
     14: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DM dm, PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
     15: {
     16:   PetscInt *ctO, *ctOInv;
     17:   PetscInt  d, c, off = 0;
     18:   PetscInt  dimOrder[5] = {3, 2, 1, 0, -1};

     20:   PetscFunctionBegin;
     21:   PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv));
     22:   if (dm) { // Order the dimensions by their starting location
     23:     PetscInt hStart[4] = {-1, -1, -1, -1};
     24:     for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, dim - d, &hStart[d], NULL));
     25:     PetscCall(PetscSortIntWithArray(dim + 1, hStart, &dimOrder[3 - dim]));
     26:   } else if (dim > 1) { // Standard plex ordering. dimOrder is in correct order if dim > 1
     27:     off             = 4 - dim;
     28:     dimOrder[off++] = 0;
     29:     for (d = dim - 1; d > 0; --d) dimOrder[off++] = d;
     30:   }

     32:   off = 0;
     33:   for (d = 0; d < 5; ++d) {
     34:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
     35:       if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue;
     36:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) == dimOrder[d]) ctO[off++] = c;
     37:     }
     38:   }
     39:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
     40:     if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) ctO[off++] = c;
     41:   }
     42:   ctO[off++] = DM_NUM_POLYTOPES;
     43:   PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off);

     45:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;

     47:   *ctOrder    = ctO;
     48:   *ctOrderInv = ctOInv;
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: /*@C
     53:   DMPlexTransformRegister - Adds a new transform component implementation

     55:   Not Collective

     57:   Input Parameters:
     58: + name        - The name of a new user-defined creation routine
     59: - create_func - The creation routine

     61:   Example Usage:
     62: .vb
     63:   DMPlexTransformRegister("my_transform", MyTransformCreate);
     64: .ve

     66:   Then, your transform type can be chosen with the procedural interface via
     67: .vb
     68:   DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
     69:   DMPlexTransformSetType(DMPlexTransform, "my_transform");
     70: .ve
     71:   or at runtime via the option
     72: .vb
     73:   -dm_plex_transform_type my_transform
     74: .ve

     76:   Level: advanced

     78:   Note:
     79:   `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms

     81: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
     82: @*/
     83: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
     84: {
     85:   PetscFunctionBegin;
     86:   PetscCall(DMInitializePackage());
     87:   PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func));
     88:   PetscFunctionReturn(PETSC_SUCCESS);
     89: }

     91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
     92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
     93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
     94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToSimplex(DMPlexTransform);
     95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
     96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
     97: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
     98: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
     99: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
    100: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform);

    102: /*@C
    103:   DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package.

    105:   Not Collective

    107:   Level: advanced

    109: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
    110: @*/
    111: PetscErrorCode DMPlexTransformRegisterAll(void)
    112: {
    113:   PetscFunctionBegin;
    114:   if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    115:   DMPlexTransformRegisterAllCalled = PETSC_TRUE;

    117:   PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter));
    118:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular));
    119:   PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox));
    120:   PetscCall(DMPlexTransformRegister(DMPLEXREFINETOSIMPLEX, DMPlexTransformCreate_ToSimplex));
    121:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld));
    122:   PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL));
    123:   PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR));
    124:   PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D));
    125:   PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude));
    126:   PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive));
    127:   PetscFunctionReturn(PETSC_SUCCESS);
    128: }

    130: /*@C
    131:   DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`.

    133:   Not collective

    135:   Level: developer

    137: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()`
    138: @*/
    139: PetscErrorCode DMPlexTransformRegisterDestroy(void)
    140: {
    141:   PetscFunctionBegin;
    142:   PetscCall(PetscFunctionListDestroy(&DMPlexTransformList));
    143:   DMPlexTransformRegisterAllCalled = PETSC_FALSE;
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }

    147: /*@
    148:   DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`.

    150:   Collective

    152:   Input Parameter:
    153: . comm - The communicator for the transform object

    155:   Output Parameter:
    156: . tr - The transform object

    158:   Level: beginner

    160: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
    161: @*/
    162: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
    163: {
    164:   DMPlexTransform t;

    166:   PetscFunctionBegin;
    167:   PetscAssertPointer(tr, 2);
    168:   *tr = NULL;
    169:   PetscCall(DMInitializePackage());

    171:   PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView));
    172:   t->setupcalled = PETSC_FALSE;
    173:   t->redFactor   = 2.0;
    174:   PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom));
    175:   *tr = t;
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: /*@
    180:   DMPlexTransformSetType - Sets the particular implementation for a transform.

    182:   Collective

    184:   Input Parameters:
    185: + tr     - The transform
    186: - method - The name of the transform type

    188:   Options Database Key:
    189: . -dm_plex_transform_type type - Sets the transform type; see `DMPlexTransformType`

    191:   Level: intermediate

    193: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
    194: @*/
    195: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
    196: {
    197:   PetscErrorCode (*r)(DMPlexTransform);
    198:   PetscBool match;

    200:   PetscFunctionBegin;
    202:   PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match));
    203:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    205:   PetscCall(DMPlexTransformRegisterAll());
    206:   PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r));
    207:   PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);

    209:   PetscTryTypeMethod(tr, destroy);
    210:   PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops)));
    211:   PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method));
    212:   PetscCall((*r)(tr));
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: /*@
    217:   DMPlexTransformGetType - Gets the type name (as a string) from the transform.

    219:   Not Collective

    221:   Input Parameter:
    222: . tr - The `DMPlexTransform`

    224:   Output Parameter:
    225: . type - The `DMPlexTransformType` name

    227:   Level: intermediate

    229: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
    230: @*/
    231: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
    232: {
    233:   PetscFunctionBegin;
    235:   PetscAssertPointer(type, 2);
    236:   PetscCall(DMPlexTransformRegisterAll());
    237:   *type = ((PetscObject)tr)->type_name;
    238:   PetscFunctionReturn(PETSC_SUCCESS);
    239: }

    241: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
    242: {
    243:   PetscViewerFormat format;

    245:   PetscFunctionBegin;
    246:   PetscCall(PetscViewerGetFormat(v, &format));
    247:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    248:     const PetscInt *trTypes = NULL;
    249:     IS              trIS;
    250:     PetscInt        cols = 8;
    251:     PetscInt        Nrt  = 8, f, g;

    253:     if (tr->trType) PetscCall(DMLabelView(tr->trType, v));
    254:     PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n"));
    255:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
    256:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
    257:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]));
    258:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
    259:     PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n"));
    260:     for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
    261:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
    262:     for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]));
    263:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));

    265:     if (tr->trType) {
    266:       PetscCall(DMLabelGetNumValues(tr->trType, &Nrt));
    267:       PetscCall(DMLabelGetValueIS(tr->trType, &trIS));
    268:       PetscCall(ISGetIndices(trIS, &trTypes));
    269:     }
    270:     PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n"));
    271:     PetscCall(PetscViewerASCIIPrintf(v, "     "));
    272:     for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]));
    273:     PetscCall(PetscViewerASCIIPrintf(v, "\n"));
    274:     for (f = 0; f < Nrt; ++f) {
    275:       PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT "  |", trTypes ? trTypes[f] : f));
    276:       for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]));
    277:       PetscCall(PetscViewerASCIIPrintf(v, " |\n"));
    278:     }
    279:     if (tr->trType) {
    280:       PetscCall(ISRestoreIndices(trIS, &trTypes));
    281:       PetscCall(ISDestroy(&trIS));
    282:     }
    283:   }
    284:   PetscFunctionReturn(PETSC_SUCCESS);
    285: }

    287: /*@
    288:   DMPlexTransformView - Views a `DMPlexTransform`

    290:   Collective

    292:   Input Parameters:
    293: + tr - the `DMPlexTransform` object to view
    294: - v  - the viewer

    296:   Level: beginner

    298: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
    299: @*/
    300: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
    301: {
    302:   PetscBool isascii;

    304:   PetscFunctionBegin;
    306:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v));
    308:   PetscCheckSameComm(tr, 1, v, 2);
    309:   PetscCall(PetscViewerCheckWritable(v));
    310:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v));
    311:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
    312:   if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v));
    313:   PetscTryTypeMethod(tr, view, v);
    314:   PetscFunctionReturn(PETSC_SUCCESS);
    315: }

    317: /*@
    318:   DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database

    320:   Collective

    322:   Input Parameter:
    323: . tr - the `DMPlexTransform` object to set options for

    325:   Options Database Keys:
    326: + -dm_plex_transform_type type               - Set the transform type, e.g. refine_regular
    327: . -dm_plex_transform_label_match_strata      - Only label points of the same stratum as the producing point
    328: . -dm_plex_transform_label_replica_inc inc   - Increment for the label value to be multiplied by the replica number, so that the new label value is oldValue + r * inc
    329: . -dm_plex_transform_active name             - Name for active mesh label
    330: - -dm_plex_transform_active_values v0,v1,... - Values in the active label

    332:   Level: intermediate

    334: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
    335: @*/
    336: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
    337: {
    338:   char        typeName[1024], active[PETSC_MAX_PATH_LEN];
    339:   const char *defName = DMPLEXREFINEREGULAR;
    340:   PetscBool   flg, match;

    342:   PetscFunctionBegin;
    344:   PetscObjectOptionsBegin((PetscObject)tr);
    345:   PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg));
    346:   if (flg) PetscCall(DMPlexTransformSetType(tr, typeName));
    347:   else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName));
    348:   PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg));
    349:   if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match));
    350:   PetscCall(PetscOptionsInt("-dm_plex_transform_label_replica_inc", "Increment for the label value to be multiplied by the replica number", "", tr->labelReplicaInc, &tr->labelReplicaInc, NULL));
    351:   PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg));
    352:   if (flg) {
    353:     DM       dm;
    354:     DMLabel  label;
    355:     PetscInt values[16];
    356:     PetscInt n = 16;

    358:     PetscCall(DMPlexTransformGetDM(tr, &dm));
    359:     PetscCall(DMGetLabel(dm, active, &label));
    360:     PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg));
    361:     if (flg && n) {
    362:       DMLabel newlabel;

    364:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel));
    365:       for (PetscInt i = 0; i < n; ++i) {
    366:         IS is;

    368:         PetscCall(DMLabelGetStratumIS(label, values[i], &is));
    369:         PetscCall(DMLabelInsertIS(newlabel, is, values[i]));
    370:         PetscCall(ISDestroy(&is));
    371:       }
    372:       PetscCall(DMPlexTransformSetActive(tr, newlabel));
    373:       PetscCall(DMLabelDestroy(&newlabel));
    374:     } else {
    375:       PetscCall(DMPlexTransformSetActive(tr, label));
    376:     }
    377:   }
    378:   PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
    379:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
    380:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject));
    381:   PetscOptionsEnd();
    382:   PetscFunctionReturn(PETSC_SUCCESS);
    383: }

    385: /*@
    386:   DMPlexTransformDestroy - Destroys a `DMPlexTransform`

    388:   Collective

    390:   Input Parameter:
    391: . tr - the transform object to destroy

    393:   Level: beginner

    395: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()`
    396: @*/
    397: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
    398: {
    399:   PetscInt c;

    401:   PetscFunctionBegin;
    402:   if (!*tr) PetscFunctionReturn(PETSC_SUCCESS);
    404:   if (--((PetscObject)*tr)->refct > 0) {
    405:     *tr = NULL;
    406:     PetscFunctionReturn(PETSC_SUCCESS);
    407:   }

    409:   PetscTryTypeMethod(*tr, destroy);
    410:   PetscCall(DMDestroy(&(*tr)->dm));
    411:   PetscCall(DMLabelDestroy(&(*tr)->active));
    412:   PetscCall(DMLabelDestroy(&(*tr)->trType));
    413:   PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld));
    414:   PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew));
    415:   PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew));
    416:   PetscCall(PetscFree((*tr)->offset));
    417:   PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd));
    418:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    419:     PetscCall(PetscFEDestroy(&(*tr)->coordFE[c]));
    420:     PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c]));
    421:   }
    422:   if ((*tr)->trVerts) {
    423:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    424:       DMPolytopeType *rct;
    425:       PetscInt       *rsize, *rcone, *rornt, Nct, n, r;

    427:       if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) {
    428:         PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    429:         for (n = 0; n < Nct; ++n) {
    430:           if (rct[n] == DM_POLYTOPE_POINT) continue;
    431:           for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r]));
    432:           PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]]));
    433:         }
    434:       }
    435:       PetscCall(PetscFree((*tr)->trSubVerts[c]));
    436:       PetscCall(PetscFree((*tr)->trVerts[c]));
    437:     }
    438:   }
    439:   PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts));
    440:   PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom));
    441:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
    442:   PetscCall(PetscHeaderDestroy(tr));
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
    447: {
    448:   DMLabel  trType = tr->trType;
    449:   PetscInt c, cN, *off;

    451:   PetscFunctionBegin;
    452:   if (trType) {
    453:     DM              dm;
    454:     IS              rtIS;
    455:     const PetscInt *reftypes;
    456:     PetscInt        Nrt, r;

    458:     PetscCall(DMPlexTransformGetDM(tr, &dm));
    459:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
    460:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
    461:     PetscCall(ISGetIndices(rtIS, &reftypes));
    462:     PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off));
    463:     for (r = 0; r < Nrt; ++r) {
    464:       const PetscInt  rt = reftypes[r];
    465:       IS              rtIS;
    466:       const PetscInt *points;
    467:       DMPolytopeType  ct;
    468:       PetscInt        np, p;

    470:       PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS));
    471:       PetscCall(ISGetLocalSize(rtIS, &np));
    472:       PetscCall(ISGetIndices(rtIS, &points));
    473:       if (!np) continue;
    474:       p = points[0];
    475:       PetscCall(ISRestoreIndices(rtIS, &points));
    476:       PetscCall(ISDestroy(&rtIS));
    477:       PetscCall(DMPlexGetCellType(dm, p, &ct));
    478:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
    479:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
    480:         DMPolytopeType      *rct;
    481:         PetscInt            *rsize, *cone, *ornt;
    482:         PetscInt             Nct, n, s;

    484:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
    485:           off[r * DM_NUM_POLYTOPES + ctNew] = -1;
    486:           break;
    487:         }
    488:         off[r * DM_NUM_POLYTOPES + ctNew] = 0;
    489:         for (s = 0; s <= r; ++s) {
    490:           const PetscInt st = reftypes[s];
    491:           DMPolytopeType sct;
    492:           PetscInt       q, qrt;

    494:           PetscCall(DMLabelGetStratumIS(trType, st, &rtIS));
    495:           PetscCall(ISGetLocalSize(rtIS, &np));
    496:           PetscCall(ISGetIndices(rtIS, &points));
    497:           if (!np) continue;
    498:           q = points[0];
    499:           PetscCall(ISRestoreIndices(rtIS, &points));
    500:           PetscCall(ISDestroy(&rtIS));
    501:           PetscCall(DMPlexGetCellType(dm, q, &sct));
    502:           PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt));
    503:           PetscCheck(st == qrt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %" PetscInt_FMT " of point %" PetscInt_FMT " does not match predicted type %" PetscInt_FMT, qrt, q, st);
    504:           if (st == rt) {
    505:             for (n = 0; n < Nct; ++n)
    506:               if (rct[n] == ctNew) break;
    507:             if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
    508:             break;
    509:           }
    510:           for (n = 0; n < Nct; ++n) {
    511:             if (rct[n] == ctNew) {
    512:               PetscInt sn;

    514:               PetscCall(DMLabelGetStratumSize(trType, st, &sn));
    515:               off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
    516:             }
    517:           }
    518:         }
    519:       }
    520:     }
    521:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
    522:     PetscCall(ISDestroy(&rtIS));
    523:   } else {
    524:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off));
    525:     for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
    526:       const DMPolytopeType ct = (DMPolytopeType)c;
    527:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
    528:         const DMPolytopeType ctNew = (DMPolytopeType)cN;
    529:         DMPolytopeType      *rct;
    530:         PetscInt            *rsize, *cone, *ornt;
    531:         PetscInt             Nct, n, i;

    533:         if (DMPolytopeTypeGetDim(ct) < 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE || DMPolytopeTypeGetDim(ctNew) < 0 || ctNew == DM_POLYTOPE_UNKNOWN_CELL || ctNew == DM_POLYTOPE_UNKNOWN_FACE) {
    534:           off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
    535:           continue;
    536:         }
    537:         off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
    538:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
    539:           const DMPolytopeType ict  = (DMPolytopeType)ctOrderOld[i];
    540:           const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];

    542:           PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt));
    543:           if (ict == ct) {
    544:             for (n = 0; n < Nct; ++n)
    545:               if (rct[n] == ctNew) break;
    546:             if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
    547:             break;
    548:           }
    549:           for (n = 0; n < Nct; ++n)
    550:             if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
    551:         }
    552:       }
    553:     }
    554:   }
    555:   *offset = off;
    556:   PetscFunctionReturn(PETSC_SUCCESS);
    557: }

    559: /*@
    560:   DMPlexTransformSetUp - Create the tables that drive the transform

    562:   Input Parameter:
    563: . tr - The `DMPlexTransform` object

    565:   Level: intermediate

    567: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    568: @*/
    569: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
    570: {
    571:   DMPolytopeType ctCell;
    572:   DM             dm;
    573:   PetscInt       pStart, pEnd, p, c, celldim = 0;

    575:   PetscFunctionBegin;
    577:   if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
    578:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    579:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
    580:   PetscTryTypeMethod(tr, setup);
    581:   PetscCall(DMSetSnapToGeomModel(dm, NULL));
    582:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));

    584:   if (pEnd > pStart) {
    585:     // Ignore cells hanging off of embedded surfaces
    586:     PetscInt c = pStart;

    588:     ctCell = DM_POLYTOPE_FV_GHOST;
    589:     while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell));
    590:   } else {
    591:     PetscInt dim;

    593:     PetscCall(DMGetDimension(dm, &dim));
    594:     switch (dim) {
    595:     case 0:
    596:       ctCell = DM_POLYTOPE_POINT;
    597:       break;
    598:     case 1:
    599:       ctCell = DM_POLYTOPE_SEGMENT;
    600:       break;
    601:     case 2:
    602:       ctCell = DM_POLYTOPE_TRIANGLE;
    603:       break;
    604:     case 3:
    605:       ctCell = DM_POLYTOPE_TETRAHEDRON;
    606:       break;
    607:     default:
    608:       ctCell = DM_POLYTOPE_UNKNOWN;
    609:     }
    610:   }
    611:   PetscCall(DMPlexCreateCellTypeOrder_Internal(dm, DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld));
    612:   for (p = pStart; p < pEnd; ++p) {
    613:     DMPolytopeType  ct;
    614:     DMPolytopeType *rct;
    615:     PetscInt       *rsize, *cone, *ornt;
    616:     PetscInt        Nct, n;

    618:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    619:     PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
    620:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
    621:     for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
    622:   }
    623:   PetscCall(DMPlexCreateCellTypeOrder_Internal(NULL, celldim, &tr->ctOrderNew, &tr->ctOrderInvNew));
    624:   /* Construct sizes and offsets for each cell type */
    625:   if (!tr->ctStart) {
    626:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

    628:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN));
    629:     PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN));
    630:     for (p = pStart; p < pEnd; ++p) {
    631:       DMPolytopeType  ct;
    632:       DMPolytopeType *rct;
    633:       PetscInt       *rsize, *cone, *ornt;
    634:       PetscInt        Nct, n;

    636:       PetscCall(DMPlexGetCellType(dm, p, &ct));
    637:       PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p);
    638:       ++ctC[ct];
    639:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt));
    640:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
    641:     }
    642:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    643:       const PetscInt cto  = tr->ctOrderOld[c];
    644:       const PetscInt cton = tr->ctOrderOld[c + 1];
    645:       const PetscInt ctn  = tr->ctOrderNew[c];
    646:       const PetscInt ctnn = tr->ctOrderNew[c + 1];

    648:       ctS[cton]  = ctS[cto] + ctC[cto];
    649:       ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
    650:     }
    651:     PetscCall(PetscFree2(ctC, ctCN));
    652:     tr->ctStart    = ctS;
    653:     tr->ctStartNew = ctSN;
    654:   }
    655:   PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset));
    656:   // Compute depth information
    657:   tr->depth = -1;
    658:   for (c = 0; c < DM_NUM_POLYTOPES; ++c)
    659:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]));
    660:   PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd));
    661:   for (PetscInt d = 0; d <= tr->depth; ++d) {
    662:     tr->depthStart[d] = PETSC_INT_MAX;
    663:     tr->depthEnd[d]   = -1;
    664:   }
    665:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    666:     const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]);

    668:     if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue;
    669:     tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]);
    670:     tr->depthEnd[dep]   = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]);
    671:   }
    672:   tr->setupcalled = PETSC_TRUE;
    673:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0));
    674:   PetscFunctionReturn(PETSC_SUCCESS);
    675: }

    677: /*@
    678:   DMPlexTransformGetDM - Get the base `DM` for the transform

    680:   Input Parameter:
    681: . tr - The `DMPlexTransform` object

    683:   Output Parameter:
    684: . dm - The original `DM` which will be transformed

    686:   Level: intermediate

    688: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    689: @*/
    690: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
    691: {
    692:   PetscFunctionBegin;
    694:   PetscAssertPointer(dm, 2);
    695:   *dm = tr->dm;
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: /*@
    700:   DMPlexTransformSetDM - Set the base `DM` for the transform

    702:   Input Parameters:
    703: + tr - The `DMPlexTransform` object
    704: - dm - The original `DM` which will be transformed

    706:   Level: intermediate

    708:   Note:
    709:   The user does not typically call this, as it is called by `DMPlexTransformApply()`.

    711: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    712: @*/
    713: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
    714: {
    715:   PetscFunctionBegin;
    718:   PetscCall(PetscObjectReference((PetscObject)dm));
    719:   PetscCall(DMDestroy(&tr->dm));
    720:   tr->dm = dm;
    721:   PetscFunctionReturn(PETSC_SUCCESS);
    722: }

    724: /*@
    725:   DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform

    727:   Input Parameter:
    728: . tr - The `DMPlexTransform` object

    730:   Output Parameter:
    731: . active - The `DMLabel` indicating which points will be transformed

    733:   Level: intermediate

    735: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    736: @*/
    737: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
    738: {
    739:   PetscFunctionBegin;
    741:   PetscAssertPointer(active, 2);
    742:   *active = tr->active;
    743:   PetscFunctionReturn(PETSC_SUCCESS);
    744: }

    746: /*@
    747:   DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform

    749:   Input Parameters:
    750: + tr     - The `DMPlexTransform` object
    751: - active - The `DMLabel` indicating which points will be transformed

    753:   Level: intermediate

    755:   Note:
    756:   This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh.

    758: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    759: @*/
    760: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
    761: {
    762:   PetscFunctionBegin;
    765:   PetscCall(PetscObjectReference((PetscObject)active));
    766:   PetscCall(DMLabelDestroy(&tr->active));
    767:   tr->active = active;
    768:   PetscFunctionReturn(PETSC_SUCCESS);
    769: }

    771: /*@
    772:   DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform

    774:   Input Parameter:
    775: . tr - The `DMPlexTransform` object

    777:   Output Parameter:
    778: . trType - The `DMLabel` indicating the transform type for each point

    780:   Level: intermediate

    782: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    783: @*/
    784: PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType)
    785: {
    786:   PetscFunctionBegin;
    788:   PetscAssertPointer(trType, 2);
    789:   *trType = tr->trType;
    790:   PetscFunctionReturn(PETSC_SUCCESS);
    791: }

    793: /*@
    794:   DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform

    796:   Input Parameters:
    797: + tr     - The `DMPlexTransform` object
    798: - trType - The original `DM` which will be transformed

    800:   Level: intermediate

    802: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive())`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    803: @*/
    804: PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType)
    805: {
    806:   PetscFunctionBegin;
    809:   PetscCall(PetscObjectReference((PetscObject)trType));
    810:   PetscCall(DMLabelDestroy(&tr->trType));
    811:   tr->trType = trType;
    812:   PetscFunctionReturn(PETSC_SUCCESS);
    813: }

    815: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
    816: {
    817:   PetscFunctionBegin;
    818:   if (!tr->coordFE[ct]) {
    819:     PetscInt dim, cdim;

    821:     dim = DMPolytopeTypeGetDim(ct);
    822:     PetscCall(DMGetCoordinateDim(tr->dm, &cdim));
    823:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]));
    824:     {
    825:       PetscDualSpace  dsp;
    826:       PetscQuadrature quad;
    827:       DM              K;
    828:       PetscFEGeom    *cg;
    829:       PetscScalar    *Xq;
    830:       PetscReal      *xq, *wq;
    831:       PetscInt        Nq, q;

    833:       PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq));
    834:       PetscCall(PetscMalloc1(Nq * cdim, &xq));
    835:       for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
    836:       PetscCall(PetscMalloc1(Nq, &wq));
    837:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
    838:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad));
    839:       PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq));
    840:       PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad));

    842:       PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp));
    843:       PetscCall(PetscDualSpaceGetDM(dsp, &K));
    844:       PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct]));
    845:       cg = tr->refGeom[ct];
    846:       PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ));
    847:       PetscCall(PetscQuadratureDestroy(&quad));
    848:     }
    849:   }
    850:   *fe = tr->coordFE[ct];
    851:   PetscFunctionReturn(PETSC_SUCCESS);
    852: }

    854: PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm)
    855: {
    856:   PetscInt dim, cdim;

    858:   PetscFunctionBegin;
    859:   PetscCall(DMGetDimension(dm, &dim));
    860:   PetscCall(DMSetDimension(tdm, dim));
    861:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    862:   PetscCall(DMSetCoordinateDim(tdm, cdim));
    863:   PetscFunctionReturn(PETSC_SUCCESS);
    864: }

    866: /*@
    867:   DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM`

    869:   Input Parameters:
    870: + tr - The `DMPlexTransform` object
    871: - dm - The original `DM`

    873:   Output Parameter:
    874: . trdm - The transformed `DM`

    876:   Level: advanced

    878: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    879: @*/
    880: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM trdm)
    881: {
    882:   PetscFunctionBegin;
    883:   PetscUseTypeMethod(tr, setdimensions, dm, trdm);
    884:   PetscFunctionReturn(PETSC_SUCCESS);
    885: }

    887: PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd)
    888: {
    889:   PetscFunctionBegin;
    890:   if (pStart) *pStart = 0;
    891:   if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]];
    892:   PetscFunctionReturn(PETSC_SUCCESS);
    893: }

    895: PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype)
    896: {
    897:   PetscInt ctNew;

    899:   PetscFunctionBegin;
    901:   PetscAssertPointer(celltype, 3);
    902:   /* TODO Can do bisection since everything is sorted */
    903:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
    904:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];

    906:     if (cell >= ctSN && cell < ctEN) break;
    907:   }
    908:   PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell);
    909:   *celltype = (DMPolytopeType)ctNew;
    910:   PetscFunctionReturn(PETSC_SUCCESS);
    911: }

    913: PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end)
    914: {
    915:   PetscFunctionBegin;
    917:   if (start) *start = tr->ctStartNew[celltype];
    918:   if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]];
    919:   PetscFunctionReturn(PETSC_SUCCESS);
    920: }

    922: PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth)
    923: {
    924:   PetscFunctionBegin;
    926:   *depth = tr->depth;
    927:   PetscFunctionReturn(PETSC_SUCCESS);
    928: }

    930: PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end)
    931: {
    932:   PetscFunctionBegin;
    934:   if (start) *start = tr->depthStart[depth];
    935:   if (end) *end = tr->depthEnd[depth];
    936:   PetscFunctionReturn(PETSC_SUCCESS);
    937: }

    939: /*@
    940:   DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels

    942:   Not Collective

    944:   Input Parameter:
    945: . tr - The `DMPlexTransform`

    947:   Output Parameter:
    948: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels

    950:   Level: intermediate

    952: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()`
    953: @*/
    954: PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match)
    955: {
    956:   PetscFunctionBegin;
    958:   PetscAssertPointer(match, 2);
    959:   *match = tr->labelMatchStrata;
    960:   PetscFunctionReturn(PETSC_SUCCESS);
    961: }

    963: /*@
    964:   DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels

    966:   Not Collective

    968:   Input Parameters:
    969: + tr    - The `DMPlexTransform`
    970: - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels

    972:   Level: intermediate

    974: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()`
    975: @*/
    976: PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match)
    977: {
    978:   PetscFunctionBegin;
    980:   tr->labelMatchStrata = match;
    981:   PetscFunctionReturn(PETSC_SUCCESS);
    982: }

    984: /*@
    985:   DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.

    987:   Not Collective

    989:   Input Parameters:
    990: + tr    - The `DMPlexTransform`
    991: . ct    - The type of the original point which produces the new point
    992: . ctNew - The type of the new point
    993: . p     - The original point which produces the new point
    994: - r     - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p`

    996:   Output Parameter:
    997: . pNew - The new point number

    999:   Level: developer

    1001: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
    1002: @*/
    1003: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
    1004: {
    1005:   DMPolytopeType *rct;
    1006:   PetscInt       *rsize, *cone, *ornt;
    1007:   PetscInt        rt, Nct, n, off, rp;
    1008:   DMLabel         trType = tr->trType;
    1009:   PetscInt        ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
    1010:   PetscInt        ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
    1011:   PetscInt        newp = ctSN, cind;

    1013:   PetscFunctionBeginHot;
    1014:   PetscCheck(p >= ctS && p < ctE, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], ctS, ctE);
    1015:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt));
    1016:   if (trType) {
    1017:     PetscCall(DMLabelGetValueIndex(trType, rt, &cind));
    1018:     PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp));
    1019:     PetscCheck(rp >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %" PetscInt_FMT " does not have refine type %" PetscInt_FMT, DMPolytopeTypes[ct], p, rt);
    1020:   } else {
    1021:     cind = ct;
    1022:     rp   = p - ctS;
    1023:   }
    1024:   off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
    1025:   PetscCheck(off >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%" PetscInt_FMT ") of point %" PetscInt_FMT " does not produce type %s for transform %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew], tr->hdr.type_name);
    1026:   newp += off;
    1027:   for (n = 0; n < Nct; ++n) {
    1028:     if (rct[n] == ctNew) {
    1029:       if (rsize[n] && r >= rsize[n])
    1030:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " for point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, p, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
    1031:       newp += rp * rsize[n] + r;
    1032:       break;
    1033:     }
    1034:   }

    1036:   PetscCheck(!(newp < ctSN) && !(newp >= ctEN), PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
    1037:   *pNew = newp;
    1038:   PetscFunctionReturn(PETSC_SUCCESS);
    1039: }

    1041: /*@
    1042:   DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.

    1044:   Not Collective

    1046:   Input Parameters:
    1047: + tr   - The `DMPlexTransform`
    1048: - pNew - The new point number

    1050:   Output Parameters:
    1051: + ct    - The type of the original point which produces the new point
    1052: . ctNew - The type of the new point
    1053: . p     - The original point which produces the new point
    1054: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

    1056:   Level: developer

    1058: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
    1059: @*/
    1060: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
    1061: {
    1062:   DMLabel         trType = tr->trType;
    1063:   DMPolytopeType *rct, ctN;
    1064:   PetscInt       *rsize, *cone, *ornt;
    1065:   PetscInt        rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO;
    1066:   PetscInt        offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS;

    1068:   PetscFunctionBegin;
    1069:   PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN));
    1070:   if (trType) {
    1071:     DM              dm;
    1072:     IS              rtIS;
    1073:     const PetscInt *reftypes;
    1074:     PetscInt        Nrt, r, rtStart;

    1076:     PetscCall(DMPlexTransformGetDM(tr, &dm));
    1077:     PetscCall(DMLabelGetNumValues(trType, &Nrt));
    1078:     PetscCall(DMLabelGetValueIS(trType, &rtIS));
    1079:     PetscCall(ISGetIndices(rtIS, &reftypes));
    1080:     for (r = 0; r < Nrt; ++r) {
    1081:       const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];

    1083:       if (tr->ctStartNew[ctN] + off > pNew) continue;
    1084:       /* Check that any of this refinement type exist */
    1085:       /* TODO Actually keep track of the number produced here instead */
    1086:       if (off > offset) {
    1087:         rt     = reftypes[r];
    1088:         offset = off;
    1089:       }
    1090:     }
    1091:     PetscCall(ISRestoreIndices(rtIS, &reftypes));
    1092:     PetscCall(ISDestroy(&rtIS));
    1093:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
    1094:     /* TODO Map refinement types to cell types */
    1095:     PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL));
    1096:     PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt);
    1097:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
    1098:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];

    1100:       if ((rtStart >= ctS) && (rtStart < ctE)) break;
    1101:     }
    1102:     PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt);
    1103:   } else {
    1104:     for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
    1105:       const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];

    1107:       if (tr->ctStartNew[ctN] + off > pNew) continue;
    1108:       if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
    1109:       /* TODO Actually keep track of the number produced here instead */
    1110:       if (off > offset) {
    1111:         ctO    = ctTmp;
    1112:         offset = off;
    1113:       }
    1114:     }
    1115:     rt = -1;
    1116:     PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew);
    1117:   }
    1118:   ctS = tr->ctStart[ctO];
    1119:   ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
    1120:   if (trType) {
    1121:     for (rtS = ctS; rtS < ctE; ++rtS) {
    1122:       PetscInt val;
    1123:       PetscCall(DMLabelGetValue(trType, rtS, &val));
    1124:       if (val == rt) break;
    1125:     }
    1126:     PetscCheck(rtS < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find point of type %s with refine type %" PetscInt_FMT, DMPolytopeTypes[ctO], rt);
    1127:   } else rtS = ctS;
    1128:   PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt));
    1129:   PetscCheck(!trType || rt == rtTmp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " has refine type %" PetscInt_FMT " != %" PetscInt_FMT " refine type which produced point %" PetscInt_FMT, rtS, rtTmp, rt, pNew);
    1130:   for (n = 0; n < Nct; ++n) {
    1131:     if (rct[n] == ctN) {
    1132:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c;

    1134:       if (trType) {
    1135:         for (c = ctS; c < ctE; ++c) {
    1136:           PetscCall(DMLabelGetValue(trType, c, &val));
    1137:           if (val == rt) {
    1138:             if (tmp < rsize[n]) break;
    1139:             tmp -= rsize[n];
    1140:           }
    1141:         }
    1142:         PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew);
    1143:         rp = c - ctS;
    1144:         rO = tmp;
    1145:       } else {
    1146:         // This assumes that all points of type ctO transform the same way
    1147:         rp = tmp / rsize[n];
    1148:         rO = tmp % rsize[n];
    1149:       }
    1150:       break;
    1151:     }
    1152:   }
    1153:   PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew);
    1154:   pO = rp + ctS;
    1155:   PetscCheck(!(pO < ctS) && !(pO >= ctE), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Source point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", pO, DMPolytopeTypes[ctO], ctS, ctE);
    1156:   if (ct) *ct = (DMPolytopeType)ctO;
    1157:   if (ctNew) *ctNew = ctN;
    1158:   if (p) *p = pO;
    1159:   if (r) *r = rO;
    1160:   PetscFunctionReturn(PETSC_SUCCESS);
    1161: }

    1163: /*@
    1164:   DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.

    1166:   Input Parameters:
    1167: + tr     - The `DMPlexTransform` object
    1168: . source - The source cell type
    1169: - p      - The source point, which can also determine the refine type

    1171:   Output Parameters:
    1172: + rt     - The refine type for this point
    1173: . Nt     - The number of types produced by this point
    1174: . target - An array of length `Nt` giving the types produced
    1175: . size   - An array of length `Nt` giving the number of cells of each type produced
    1176: . cone   - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
    1177: - ornt   - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

    1179:   Level: advanced

    1181:   Notes:
    1182:   The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
    1183:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
    1184:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
    1185: .vb
    1186:    the number of cones to be taken, or 0 for the current cell
    1187:    the cell cone point number at each level from which it is subdivided
    1188:    the replica number r of the subdivision.
    1189: .ve
    1190:   The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
    1191: .vb
    1192:    Nt     = 2
    1193:    target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
    1194:    size   = {1, 2}
    1195:    cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
    1196:    ornt   = {                         0,                       0,                        0,                          0}
    1197: .ve

    1199: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
    1200: @*/
    1201: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    1202: {
    1203:   PetscFunctionBegin;
    1204:   PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
    1205:   PetscFunctionReturn(PETSC_SUCCESS);
    1206: }

    1208: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
    1209: {
    1210:   PetscFunctionBegin;
    1211:   *rnew = r;
    1212:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
    1213:   PetscFunctionReturn(PETSC_SUCCESS);
    1214: }

    1216: /* Returns the same thing */
    1217: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
    1218: {
    1219:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
    1220:   static PetscInt       vertexS[] = {1};
    1221:   static PetscInt       vertexC[] = {0};
    1222:   static PetscInt       vertexO[] = {0};
    1223:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
    1224:   static PetscInt       edgeS[]   = {1};
    1225:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
    1226:   static PetscInt       edgeO[]   = {0, 0};
    1227:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
    1228:   static PetscInt       tedgeS[]  = {1};
    1229:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
    1230:   static PetscInt       tedgeO[]  = {0, 0};
    1231:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
    1232:   static PetscInt       triS[]    = {1};
    1233:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
    1234:   static PetscInt       triO[]    = {0, 0, 0};
    1235:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
    1236:   static PetscInt       quadS[]   = {1};
    1237:   static PetscInt       quadC[]   = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
    1238:   static PetscInt       quadO[]   = {0, 0, 0, 0};
    1239:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
    1240:   static PetscInt       tquadS[]  = {1};
    1241:   static PetscInt       tquadC[]  = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
    1242:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
    1243:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
    1244:   static PetscInt       tetS[]    = {1};
    1245:   static PetscInt       tetC[]    = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0};
    1246:   static PetscInt       tetO[]    = {0, 0, 0, 0};
    1247:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
    1248:   static PetscInt       hexS[]    = {1};
    1249:   static PetscInt       hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
    1250:   static PetscInt       hexO[] = {0, 0, 0, 0, 0, 0};
    1251:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
    1252:   static PetscInt       tripS[]   = {1};
    1253:   static PetscInt       tripC[]   = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
    1254:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
    1255:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
    1256:   static PetscInt       ttripS[]  = {1};
    1257:   static PetscInt       ttripC[]  = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
    1258:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
    1259:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
    1260:   static PetscInt       tquadpS[] = {1};
    1261:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL,    1, 0, 0, DM_POLYTOPE_QUADRILATERAL,    1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
    1262:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
    1263:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
    1264:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
    1265:   static PetscInt       pyrS[]    = {1};
    1266:   static PetscInt       pyrC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
    1267:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

    1269:   PetscFunctionBegin;
    1270:   if (rt) *rt = 0;
    1271:   switch (source) {
    1272:   case DM_POLYTOPE_POINT:
    1273:     *Nt     = 1;
    1274:     *target = vertexT;
    1275:     *size   = vertexS;
    1276:     *cone   = vertexC;
    1277:     *ornt   = vertexO;
    1278:     break;
    1279:   case DM_POLYTOPE_SEGMENT:
    1280:     *Nt     = 1;
    1281:     *target = edgeT;
    1282:     *size   = edgeS;
    1283:     *cone   = edgeC;
    1284:     *ornt   = edgeO;
    1285:     break;
    1286:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
    1287:     *Nt     = 1;
    1288:     *target = tedgeT;
    1289:     *size   = tedgeS;
    1290:     *cone   = tedgeC;
    1291:     *ornt   = tedgeO;
    1292:     break;
    1293:   case DM_POLYTOPE_TRIANGLE:
    1294:     *Nt     = 1;
    1295:     *target = triT;
    1296:     *size   = triS;
    1297:     *cone   = triC;
    1298:     *ornt   = triO;
    1299:     break;
    1300:   case DM_POLYTOPE_QUADRILATERAL:
    1301:     *Nt     = 1;
    1302:     *target = quadT;
    1303:     *size   = quadS;
    1304:     *cone   = quadC;
    1305:     *ornt   = quadO;
    1306:     break;
    1307:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
    1308:     *Nt     = 1;
    1309:     *target = tquadT;
    1310:     *size   = tquadS;
    1311:     *cone   = tquadC;
    1312:     *ornt   = tquadO;
    1313:     break;
    1314:   case DM_POLYTOPE_TETRAHEDRON:
    1315:     *Nt     = 1;
    1316:     *target = tetT;
    1317:     *size   = tetS;
    1318:     *cone   = tetC;
    1319:     *ornt   = tetO;
    1320:     break;
    1321:   case DM_POLYTOPE_HEXAHEDRON:
    1322:     *Nt     = 1;
    1323:     *target = hexT;
    1324:     *size   = hexS;
    1325:     *cone   = hexC;
    1326:     *ornt   = hexO;
    1327:     break;
    1328:   case DM_POLYTOPE_TRI_PRISM:
    1329:     *Nt     = 1;
    1330:     *target = tripT;
    1331:     *size   = tripS;
    1332:     *cone   = tripC;
    1333:     *ornt   = tripO;
    1334:     break;
    1335:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
    1336:     *Nt     = 1;
    1337:     *target = ttripT;
    1338:     *size   = ttripS;
    1339:     *cone   = ttripC;
    1340:     *ornt   = ttripO;
    1341:     break;
    1342:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    1343:     *Nt     = 1;
    1344:     *target = tquadpT;
    1345:     *size   = tquadpS;
    1346:     *cone   = tquadpC;
    1347:     *ornt   = tquadpO;
    1348:     break;
    1349:   case DM_POLYTOPE_PYRAMID:
    1350:     *Nt     = 1;
    1351:     *target = pyrT;
    1352:     *size   = pyrS;
    1353:     *cone   = pyrC;
    1354:     *ornt   = pyrO;
    1355:     break;
    1356:   default:
    1357:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
    1358:   }
    1359:   PetscFunctionReturn(PETSC_SUCCESS);
    1360: }

    1362: /*@
    1363:   DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point

    1365:   Not Collective

    1367:   Input Parameters:
    1368: + tr  - The `DMPlexTransform`
    1369: . sct - The source point cell type, from whom the new cell is being produced
    1370: . sp  - The source point
    1371: . so  - The orientation of the source point in its enclosing parent
    1372: . tct - The target point cell type
    1373: . r   - The replica number requested for the produced cell type
    1374: - o   - The orientation of the replica

    1376:   Output Parameters:
    1377: + rnew - The replica number, given the orientation of the parent
    1378: - onew - The replica orientation, given the orientation of the parent

    1380:   Level: advanced

    1382: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
    1383: @*/
    1384: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
    1385: {
    1386:   PetscFunctionBeginHot;
    1387:   PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
    1388:   PetscFunctionReturn(PETSC_SUCCESS);
    1389: }

    1391: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
    1392: {
    1393:   DM       dm;
    1394:   PetscInt pStart, pEnd, pNew;

    1396:   PetscFunctionBegin;
    1397:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1398:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
    1399:   /* Must create the celltype label here so that we do not automatically try to compute the types */
    1400:   PetscCall(DMCreateLabel(rdm, "celltype"));
    1401:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1402:   for (PetscInt p = pStart; p < pEnd; ++p) {
    1403:     DMPolytopeType  ct;
    1404:     DMPolytopeType *rct;
    1405:     PetscInt       *rsize, *rcone, *rornt;
    1406:     PetscInt        Nct, n, r;

    1408:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    1409:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1410:     for (n = 0; n < Nct; ++n) {
    1411:       for (r = 0; r < rsize[n]; ++r) {
    1412:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
    1413:         PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n])));
    1414:         PetscCall(DMPlexSetCellType(rdm, pNew, rct[n]));
    1415:       }
    1416:     }
    1417:   }
    1418:   /* Let the DM know we have set all the cell types */
    1419:   {
    1420:     DMLabel  ctLabel;
    1421:     DM_Plex *plex = (DM_Plex *)rdm->data;

    1423:     PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel));
    1424:     PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState));
    1425:   }
    1426:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0));
    1427:   PetscFunctionReturn(PETSC_SUCCESS);
    1428: }

    1430: PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
    1431: {
    1432:   DMPolytopeType ctNew;

    1434:   PetscFunctionBegin;
    1436:   PetscAssertPointer(coneSize, 3);
    1437:   PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew));
    1438:   *coneSize = DMPolytopeTypeGetConeSize(ctNew);
    1439:   PetscFunctionReturn(PETSC_SUCCESS);
    1440: }

    1442: /* The orientation o is for the interior of the cell p */
    1443: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew, const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff, PetscInt coneNew[], PetscInt orntNew[])
    1444: {
    1445:   DM              dm;
    1446:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
    1447:   const PetscInt *cone;
    1448:   DMPolytopeType *newft = NULL;
    1449:   PetscInt        c, coff = *coneoff, ooff = *orntoff;
    1450:   PetscInt        dim, cr = 0, co = 0, nr, no;

    1452:   PetscFunctionBegin;
    1453:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1454:   PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL));
    1455:   // Check if we have to permute this cell
    1456:   PetscCall(DMGetDimension(dm, &dim));
    1457:   if (DMPolytopeTypeGetDim(ctNew) == dim && DMPolytopeTypeGetDim(ct) == dim - 1) {
    1458:     PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, o, ctNew, cr, co, &nr, &no));
    1459:     if (cr != nr || co != no) PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newft));
    1460:   }
    1461:   for (c = 0; c < csizeNew; ++c) {
    1462:     PetscInt             ppp   = -1;                            /* Parent Parent point: Parent of point pp */
    1463:     PetscInt             pp    = p;                             /* Parent point: Point in the original mesh producing new cone point */
    1464:     PetscInt             po    = 0;                             /* Orientation of parent point pp in parent parent point ppp */
    1465:     DMPolytopeType       pct   = ct;                            /* Parent type: Cell type for parent of new cone point */
    1466:     const PetscInt      *pcone = cone;                          /* Parent cone: Cone of parent point pp */
    1467:     PetscInt             pr    = -1;                            /* Replica number of pp that produces new cone point  */
    1468:     const DMPolytopeType ft    = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
    1469:     const PetscInt       fn    = rcone[coff++];                 /* Number of cones of p that need to be taken when producing new cone point */
    1470:     PetscInt             fo    = rornt[ooff++];                 /* Orientation of new cone point in pNew */
    1471:     PetscInt             lc;

    1473:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
    1474:     for (lc = 0; lc < fn; ++lc) {
    1475:       const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po);
    1476:       const PetscInt  acp  = rcone[coff++];
    1477:       const PetscInt  pcp  = parr[acp * 2];
    1478:       const PetscInt  pco  = parr[acp * 2 + 1];
    1479:       const PetscInt *ppornt;

    1481:       ppp = pp;
    1482:       pp  = pcone[pcp];
    1483:       PetscCall(DMPlexGetCellType(dm, pp, &pct));
    1484:       // Restore the parent cone from the last iterate
    1485:       if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL));
    1486:       PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL));
    1487:       PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt));
    1488:       po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
    1489:       PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt));
    1490:     }
    1491:     if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL));
    1492:     pr = rcone[coff++];
    1493:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
    1494:     PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo));
    1495:     PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]));
    1496:     orntNew[c] = fo;
    1497:     if (newft) newft[c] = ft;
    1498:   }
    1499:   PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL));
    1500:   if (newft) {
    1501:     const PetscInt *arr;
    1502:     PetscInt       *newcone, *newornt;

    1504:     arr = DMPolytopeTypeGetArrangement(ctNew, no);
    1505:     PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newcone));
    1506:     PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newornt));
    1507:     for (PetscInt c = 0; c < csizeNew; ++c) {
    1508:       DMPolytopeType ft = newft[c];
    1509:       PetscInt       nO;

    1511:       nO         = DMPolytopeTypeGetNumArrangements(ft) / 2;
    1512:       newcone[c] = coneNew[arr[c * 2 + 0]];
    1513:       newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], orntNew[arr[c * 2 + 0]]);
    1514:       PetscCheck(!newornt[c] || !(newornt[c] >= nO || newornt[c] < -nO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid orientation %" PetscInt_FMT " not in [%" PetscInt_FMT ",%" PetscInt_FMT ") for %s %" PetscInt_FMT, newornt[c], -nO, nO, DMPolytopeTypes[ft], coneNew[c]);
    1515:     }
    1516:     for (PetscInt c = 0; c < csizeNew; ++c) {
    1517:       coneNew[c] = newcone[c];
    1518:       orntNew[c] = newornt[c];
    1519:     }
    1520:     PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newcone));
    1521:     PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newornt));
    1522:     PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newft));
    1523:   }
    1524:   *coneoff = coff;
    1525:   *orntoff = ooff;
    1526:   PetscFunctionReturn(PETSC_SUCCESS);
    1527: }

    1529: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
    1530: {
    1531:   DM             dm;
    1532:   DMPolytopeType ct;
    1533:   PetscInt      *coneNew, *orntNew;
    1534:   PetscInt       maxConeSize = 0, pStart, pEnd, p, pNew;

    1536:   PetscFunctionBegin;
    1537:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1538:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
    1539:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
    1540:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
    1541:   PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
    1542:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1543:   for (p = pStart; p < pEnd; ++p) {
    1544:     PetscInt        coff, ooff;
    1545:     DMPolytopeType *rct;
    1546:     PetscInt       *rsize, *rcone, *rornt;
    1547:     PetscInt        Nct, n, r;

    1549:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    1550:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1551:     for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
    1552:       const DMPolytopeType ctNew = rct[n];

    1554:       for (r = 0; r < rsize[n]; ++r) {
    1555:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
    1556:         PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew));
    1557:         PetscCall(DMPlexSetCone(rdm, pNew, coneNew));
    1558:         PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew));
    1559:       }
    1560:     }
    1561:   }
    1562:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew));
    1563:   PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew));
    1564:   PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view"));
    1565:   PetscCall(DMPlexSymmetrize(rdm));
    1566:   PetscCall(DMPlexStratify(rdm));
    1567:   PetscTryTypeMethod(tr, ordersupports, dm, rdm);
    1568:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0));
    1569:   PetscFunctionReturn(PETSC_SUCCESS);
    1570: }

    1572: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
    1573: {
    1574:   DM              dm;
    1575:   DMPolytopeType  ct, qct;
    1576:   DMPolytopeType *rct;
    1577:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
    1578:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

    1580:   PetscFunctionBegin;
    1582:   PetscAssertPointer(cone, 4);
    1583:   PetscAssertPointer(ornt, 5);
    1584:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
    1585:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1586:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
    1587:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
    1588:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
    1589:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1590:   for (n = 0; n < Nct; ++n) {
    1591:     const DMPolytopeType ctNew    = rct[n];
    1592:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
    1593:     PetscInt             Nr       = rsize[n], fn, c;

    1595:     if (ctNew == qct) Nr = r;
    1596:     for (nr = 0; nr < Nr; ++nr) {
    1597:       for (c = 0; c < csizeNew; ++c) {
    1598:         ++coff;             /* Cell type of new cone point */
    1599:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
    1600:         coff += fn;
    1601:         ++coff; /* Replica number of new cone point */
    1602:         ++ooff; /* Orientation of new cone point */
    1603:       }
    1604:     }
    1605:     if (ctNew == qct) break;
    1606:   }
    1607:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
    1608:   *cone = qcone;
    1609:   *ornt = qornt;
    1610:   PetscFunctionReturn(PETSC_SUCCESS);
    1611: }

    1613: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
    1614: {
    1615:   DM              dm;
    1616:   DMPolytopeType  ct, qct;
    1617:   DMPolytopeType *rct;
    1618:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
    1619:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;

    1621:   PetscFunctionBegin;
    1623:   if (cone) PetscAssertPointer(cone, 3);
    1624:   if (ornt) PetscAssertPointer(ornt, 4);
    1625:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
    1626:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1627:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
    1628:   PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
    1629:   PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r));
    1630:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1631:   for (n = 0; n < Nct; ++n) {
    1632:     const DMPolytopeType ctNew    = rct[n];
    1633:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
    1634:     PetscInt             Nr       = rsize[n], fn, c;

    1636:     if (ctNew == qct) Nr = r;
    1637:     for (nr = 0; nr < Nr; ++nr) {
    1638:       for (c = 0; c < csizeNew; ++c) {
    1639:         ++coff;             /* Cell type of new cone point */
    1640:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
    1641:         coff += fn;
    1642:         ++coff; /* Replica number of new cone point */
    1643:         ++ooff; /* Orientation of new cone point */
    1644:       }
    1645:     }
    1646:     if (ctNew == qct) break;
    1647:   }
    1648:   PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt));
    1649:   if (cone) *cone = qcone;
    1650:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone));
    1651:   if (ornt) *ornt = qornt;
    1652:   else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt));
    1653:   PetscFunctionReturn(PETSC_SUCCESS);
    1654: }

    1656: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
    1657: {
    1658:   DM dm;

    1660:   PetscFunctionBegin;
    1662:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1663:   if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone));
    1664:   if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt));
    1665:   PetscFunctionReturn(PETSC_SUCCESS);
    1666: }

    1668: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
    1669: {
    1670:   PetscInt ict;

    1672:   PetscFunctionBegin;
    1673:   PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts));
    1674:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
    1675:     const DMPolytopeType ct = (DMPolytopeType)ict;
    1676:     DMPlexTransform      reftr;
    1677:     DM                   refdm, trdm;
    1678:     Vec                  coordinates;
    1679:     const PetscScalar   *coords;
    1680:     DMPolytopeType      *rct;
    1681:     PetscInt            *rsize, *rcone, *rornt;
    1682:     PetscInt             Nct, n, r, pNew = 0;
    1683:     PetscInt             trdim, vStart, vEnd, Nc;
    1684:     const PetscInt       debug = 0;
    1685:     const char          *typeName;

    1687:     /* Since points are 0-dimensional, coordinates make no sense */
    1688:     if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue;
    1689:     PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm));
    1690:     PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr));
    1691:     PetscCall(DMPlexTransformSetDM(reftr, refdm));
    1692:     PetscCall(DMPlexTransformGetType(tr, &typeName));
    1693:     PetscCall(DMPlexTransformSetType(reftr, typeName));
    1694:     PetscCall(DMPlexTransformSetUp(reftr));
    1695:     PetscCall(DMPlexTransformApply(reftr, refdm, &trdm));

    1697:     PetscCall(DMGetDimension(trdm, &trdim));
    1698:     PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd));
    1699:     tr->trNv[ct] = vEnd - vStart;
    1700:     PetscCall(DMGetCoordinatesLocal(trdm, &coordinates));
    1701:     PetscCall(VecGetLocalSize(coordinates, &Nc));
    1702:     PetscCheck(tr->trNv[ct] * trdim == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell type %s, transformed coordinate size %" PetscInt_FMT " != %" PetscInt_FMT " size of coordinate storage", DMPolytopeTypes[ct], tr->trNv[ct] * trdim, Nc);
    1703:     PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct]));
    1704:     PetscCall(VecGetArrayRead(coordinates, &coords));
    1705:     PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc));
    1706:     PetscCall(VecRestoreArrayRead(coordinates, &coords));

    1708:     PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]));
    1709:     PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1710:     for (n = 0; n < Nct; ++n) {
    1711:       /* Since points are 0-dimensional, coordinates make no sense */
    1712:       if (rct[n] == DM_POLYTOPE_POINT) continue;
    1713:       PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]));
    1714:       for (r = 0; r < rsize[n]; ++r) {
    1715:         PetscInt *closure = NULL;
    1716:         PetscInt  clSize, cl, Nv = 0;

    1718:         PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]));
    1719:         PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew));
    1720:         PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
    1721:         for (cl = 0; cl < clSize * 2; cl += 2) {
    1722:           const PetscInt sv = closure[cl];

    1724:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
    1725:         }
    1726:         PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure));
    1727:         PetscCheck(Nv == DMPolytopeTypeGetNumVertices(rct[n]), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of vertices %" PetscInt_FMT " != %" PetscInt_FMT " for %s subcell %" PetscInt_FMT " from cell %s", Nv, DMPolytopeTypeGetNumVertices(rct[n]), DMPolytopeTypes[rct[n]], r, DMPolytopeTypes[ct]);
    1728:       }
    1729:     }
    1730:     if (debug) {
    1731:       DMPolytopeType *rct;
    1732:       PetscInt       *rsize, *rcone, *rornt;
    1733:       PetscInt        v, dE = trdim, d, off = 0;

    1735:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]));
    1736:       for (v = 0; v < tr->trNv[ct]; ++v) {
    1737:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
    1738:         for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++])));
    1739:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1740:       }

    1742:       PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1743:       for (n = 0; n < Nct; ++n) {
    1744:         if (rct[n] == DM_POLYTOPE_POINT) continue;
    1745:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]));
    1746:         for (r = 0; r < rsize[n]; ++r) {
    1747:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  "));
    1748:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]));
    1749:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1750:         }
    1751:       }
    1752:     }
    1753:     PetscCall(DMDestroy(&refdm));
    1754:     PetscCall(DMDestroy(&trdm));
    1755:     PetscCall(DMPlexTransformDestroy(&reftr));
    1756:   }
    1757:   PetscFunctionReturn(PETSC_SUCCESS);
    1758: }

    1760: /*@C
    1761:   DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type

    1763:   Input Parameters:
    1764: + tr - The `DMPlexTransform` object
    1765: - ct - The cell type

    1767:   Output Parameters:
    1768: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
    1769: - trVerts - The coordinates of these vertices in the reference cell

    1771:   Level: developer

    1773: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()`
    1774: @*/
    1775: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
    1776: {
    1777:   PetscFunctionBegin;
    1778:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
    1779:   if (Nv) *Nv = tr->trNv[ct];
    1780:   if (trVerts) *trVerts = tr->trVerts[ct];
    1781:   PetscFunctionReturn(PETSC_SUCCESS);
    1782: }

    1784: /*@C
    1785:   DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type

    1787:   Input Parameters:
    1788: + tr  - The `DMPlexTransform` object
    1789: . ct  - The cell type
    1790: . rct - The subcell type
    1791: - r   - The subcell index

    1793:   Output Parameter:
    1794: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()`

    1796:   Level: developer

    1798: .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()`
    1799: @*/
    1800: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
    1801: {
    1802:   PetscFunctionBegin;
    1803:   if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr));
    1804:   PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
    1805:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
    1806:   PetscFunctionReturn(PETSC_SUCCESS);
    1807: }

    1809: /* Computes new vertex as the barycenter, or centroid */
    1810: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
    1811: {
    1812:   PetscInt v, d;

    1814:   PetscFunctionBeginHot;
    1815:   PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]);
    1816:   for (d = 0; d < dE; ++d) out[d] = 0.0;
    1817:   for (v = 0; v < Nv; ++v)
    1818:     for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
    1819:   for (d = 0; d < dE; ++d) out[d] /= Nv;
    1820:   PetscFunctionReturn(PETSC_SUCCESS);
    1821: }

    1823: /*@
    1824:   DMPlexTransformMapCoordinates - Calculate new coordinates for produced points

    1826:   Not collective

    1828:   Input Parameters:
    1829: + tr  - The `DMPlexTransform`
    1830: . pct - The cell type of the parent, from whom the new cell is being produced
    1831: . ct  - The type being produced
    1832: . p   - The original point
    1833: . r   - The replica number requested for the produced cell type
    1834: . Nv  - Number of vertices in the closure of the parent cell
    1835: . dE  - Spatial dimension
    1836: - in  - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

    1838:   Output Parameter:
    1839: . out - The coordinates of the new vertices

    1841:   Level: intermediate

    1843: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`
    1844: @*/
    1845: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
    1846: {
    1847:   PetscFunctionBeginHot;
    1848:   if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
    1849:   PetscFunctionReturn(PETSC_SUCCESS);
    1850: }

    1852: /*
    1853:   DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label

    1855:   Not Collective

    1857:   Input Parameters:
    1858: + tr    - The `DMPlexTransform`
    1859: . label - The label in the transformed mesh
    1860: . pp    - The parent point in the original mesh
    1861: . pct   - The cell type of the parent point
    1862: . p     - The point in the transformed mesh
    1863: . ct    - The cell type of the point
    1864: . r     - The replica number of the point
    1865: - val   - The label value of the parent point

    1867:   Level: developer

    1869: .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()`
    1870: */
    1871: static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val)
    1872: {
    1873:   PetscFunctionBeginHot;
    1874:   if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS);
    1875:   PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r));
    1876:   PetscFunctionReturn(PETSC_SUCCESS);
    1877: }

    1879: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
    1880: {
    1881:   DM              dm;
    1882:   IS              valueIS;
    1883:   const PetscInt *values;
    1884:   PetscInt        defVal, Nv, val;

    1886:   PetscFunctionBegin;
    1887:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1888:   PetscCall(DMLabelGetDefaultValue(label, &defVal));
    1889:   PetscCall(DMLabelSetDefaultValue(labelNew, defVal));
    1890:   PetscCall(DMLabelGetValueIS(label, &valueIS));
    1891:   PetscCall(ISGetLocalSize(valueIS, &Nv));
    1892:   PetscCall(ISGetIndices(valueIS, &values));
    1893:   for (val = 0; val < Nv; ++val) {
    1894:     IS              pointIS;
    1895:     const PetscInt *points;
    1896:     PetscInt        numPoints, p;

    1898:     /* Ensure refined label is created with same number of strata as
    1899:      * original (even if no entries here). */
    1900:     PetscCall(DMLabelAddStratum(labelNew, values[val]));
    1901:     PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
    1902:     PetscCall(ISGetLocalSize(pointIS, &numPoints));
    1903:     PetscCall(ISGetIndices(pointIS, &points));
    1904:     for (p = 0; p < numPoints; ++p) {
    1905:       const PetscInt  point = points[p];
    1906:       DMPolytopeType  ct;
    1907:       DMPolytopeType *rct;
    1908:       PetscInt       *rsize, *rcone, *rornt;
    1909:       PetscInt        Nct, n, r, pNew = 0;

    1911:       PetscCall(DMPlexGetCellType(dm, point, &ct));
    1912:       PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    1913:       for (n = 0; n < Nct; ++n) {
    1914:         for (r = 0; r < rsize[n]; ++r) {
    1915:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew));
    1916:           PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val]));
    1917:         }
    1918:       }
    1919:     }
    1920:     PetscCall(ISRestoreIndices(pointIS, &points));
    1921:     PetscCall(ISDestroy(&pointIS));
    1922:   }
    1923:   PetscCall(ISRestoreIndices(valueIS, &values));
    1924:   PetscCall(ISDestroy(&valueIS));
    1925:   PetscFunctionReturn(PETSC_SUCCESS);
    1926: }

    1928: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
    1929: {
    1930:   DM       dm;
    1931:   PetscInt numLabels, l;

    1933:   PetscFunctionBegin;
    1934:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1935:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
    1936:   PetscCall(DMGetNumLabels(dm, &numLabels));
    1937:   for (l = 0; l < numLabels; ++l) {
    1938:     DMLabel     label, labelNew;
    1939:     const char *lname;
    1940:     PetscBool   isDepth, isCellType;

    1942:     PetscCall(DMGetLabelName(dm, l, &lname));
    1943:     PetscCall(PetscStrcmp(lname, "depth", &isDepth));
    1944:     if (isDepth) continue;
    1945:     PetscCall(PetscStrcmp(lname, "celltype", &isCellType));
    1946:     if (isCellType) continue;
    1947:     PetscCall(DMCreateLabel(rdm, lname));
    1948:     PetscCall(DMGetLabel(dm, lname, &label));
    1949:     PetscCall(DMGetLabel(rdm, lname, &labelNew));
    1950:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
    1951:   }
    1952:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0));
    1953:   PetscFunctionReturn(PETSC_SUCCESS);
    1954: }

    1956: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
    1957: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
    1958: {
    1959:   DM       dm;
    1960:   PetscInt Nf, f, Nds, s;

    1962:   PetscFunctionBegin;
    1963:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    1964:   PetscCall(DMGetNumFields(dm, &Nf));
    1965:   for (f = 0; f < Nf; ++f) {
    1966:     DMLabel     label, labelNew;
    1967:     PetscObject obj;
    1968:     const char *lname;

    1970:     PetscCall(DMGetField(rdm, f, &label, &obj));
    1971:     if (!label) continue;
    1972:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
    1973:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
    1974:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
    1975:     PetscCall(DMSetField_Internal(rdm, f, labelNew, obj));
    1976:     PetscCall(DMLabelDestroy(&labelNew));
    1977:   }
    1978:   PetscCall(DMGetNumDS(dm, &Nds));
    1979:   for (s = 0; s < Nds; ++s) {
    1980:     DMLabel     label, labelNew;
    1981:     const char *lname;

    1983:     PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL));
    1984:     if (!label) continue;
    1985:     PetscCall(PetscObjectGetName((PetscObject)label, &lname));
    1986:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew));
    1987:     PetscCall(RefineLabel_Internal(tr, label, labelNew));
    1988:     PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL));
    1989:     PetscCall(DMLabelDestroy(&labelNew));
    1990:   }
    1991:   PetscFunctionReturn(PETSC_SUCCESS);
    1992: }

    1994: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
    1995: {
    1996:   DM                 dm;
    1997:   PetscSF            sf, sfNew;
    1998:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
    1999:   const PetscInt    *localPoints;
    2000:   const PetscSFNode *remotePoints;
    2001:   PetscInt          *localPointsNew;
    2002:   PetscSFNode       *remotePointsNew;
    2003:   PetscInt           pStartNew, pEndNew, pNew;
    2004:   /* Brute force algorithm */
    2005:   PetscSF         rsf;
    2006:   PetscSection    s;
    2007:   const PetscInt *rootdegree;
    2008:   PetscInt       *rootPointsNew, *remoteOffsets;
    2009:   PetscInt        numPointsNew, pStart, pEnd, p;

    2011:   PetscFunctionBegin;
    2012:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    2013:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
    2014:   PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew));
    2015:   PetscCall(DMGetPointSF(dm, &sf));
    2016:   PetscCall(DMGetPointSF(rdm, &sfNew));
    2017:   /* Calculate size of new SF */
    2018:   PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
    2019:   if (numRoots < 0) {
    2020:     PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
    2021:     PetscFunctionReturn(PETSC_SUCCESS);
    2022:   }
    2023:   for (l = 0; l < numLeaves; ++l) {
    2024:     const PetscInt  p = localPoints[l];
    2025:     DMPolytopeType  ct;
    2026:     DMPolytopeType *rct;
    2027:     PetscInt       *rsize, *rcone, *rornt;
    2028:     PetscInt        Nct, n;

    2030:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    2031:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2032:     for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
    2033:   }
    2034:   /* Send new root point numbers
    2035:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
    2036:   */
    2037:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    2038:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s));
    2039:   PetscCall(PetscSectionSetChart(s, pStart, pEnd));
    2040:   for (p = pStart; p < pEnd; ++p) {
    2041:     DMPolytopeType  ct;
    2042:     DMPolytopeType *rct;
    2043:     PetscInt       *rsize, *rcone, *rornt;
    2044:     PetscInt        Nct, n;

    2046:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    2047:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2048:     for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n]));
    2049:   }
    2050:   PetscCall(PetscSectionSetUp(s));
    2051:   PetscCall(PetscSectionGetStorageSize(s, &numPointsNew));
    2052:   PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets));
    2053:   PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf));
    2054:   PetscCall(PetscFree(remoteOffsets));
    2055:   PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
    2056:   PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
    2057:   PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew));
    2058:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
    2059:   for (p = pStart; p < pEnd; ++p) {
    2060:     DMPolytopeType  ct;
    2061:     DMPolytopeType *rct;
    2062:     PetscInt       *rsize, *rcone, *rornt;
    2063:     PetscInt        Nct, n, r, off;

    2065:     if (!rootdegree[p - pStart]) continue;
    2066:     PetscCall(PetscSectionGetOffset(s, p, &off));
    2067:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    2068:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2069:     for (n = 0, m = 0; n < Nct; ++n) {
    2070:       for (r = 0; r < rsize[n]; ++r, ++m) {
    2071:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
    2072:         rootPointsNew[off + m] = pNew;
    2073:       }
    2074:     }
    2075:   }
    2076:   PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
    2077:   PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE));
    2078:   PetscCall(PetscSFDestroy(&rsf));
    2079:   PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew));
    2080:   PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew));
    2081:   for (l = 0, m = 0; l < numLeaves; ++l) {
    2082:     const PetscInt  p = localPoints[l];
    2083:     DMPolytopeType  ct;
    2084:     DMPolytopeType *rct;
    2085:     PetscInt       *rsize, *rcone, *rornt;
    2086:     PetscInt        Nct, n, r, q, off;

    2088:     PetscCall(PetscSectionGetOffset(s, p, &off));
    2089:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    2090:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2091:     for (n = 0, q = 0; n < Nct; ++n) {
    2092:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
    2093:         PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
    2094:         localPointsNew[m]        = pNew;
    2095:         remotePointsNew[m].index = rootPointsNew[off + q];
    2096:         remotePointsNew[m].rank  = remotePoints[l].rank;
    2097:       }
    2098:     }
    2099:   }
    2100:   PetscCall(PetscSectionDestroy(&s));
    2101:   PetscCall(PetscFree(rootPointsNew));
    2102:   /* SF needs sorted leaves to correctly calculate Gather */
    2103:   {
    2104:     PetscSFNode *rp, *rtmp;
    2105:     PetscInt    *lp, *idx, *ltmp, i;

    2107:     PetscCall(PetscMalloc1(numLeavesNew, &idx));
    2108:     PetscCall(PetscMalloc1(numLeavesNew, &lp));
    2109:     PetscCall(PetscMalloc1(numLeavesNew, &rp));
    2110:     for (i = 0; i < numLeavesNew; ++i) {
    2111:       PetscCheck(!(localPointsNew[i] < pStartNew) && !(localPointsNew[i] >= pEndNew), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %" PetscInt_FMT " (%" PetscInt_FMT ") not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[i], i, pStartNew, pEndNew);
    2112:       idx[i] = i;
    2113:     }
    2114:     PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx));
    2115:     for (i = 0; i < numLeavesNew; ++i) {
    2116:       lp[i] = localPointsNew[idx[i]];
    2117:       rp[i] = remotePointsNew[idx[i]];
    2118:     }
    2119:     ltmp            = localPointsNew;
    2120:     localPointsNew  = lp;
    2121:     rtmp            = remotePointsNew;
    2122:     remotePointsNew = rp;
    2123:     PetscCall(PetscFree(idx));
    2124:     PetscCall(PetscFree(ltmp));
    2125:     PetscCall(PetscFree(rtmp));
    2126:   }
    2127:   PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
    2128:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0));
    2129:   PetscFunctionReturn(PETSC_SUCCESS);
    2130: }

    2132: /*
    2133:   DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`.

    2135:   Not Collective

    2137:   Input Parameters:
    2138: + tr  - The `DMPlexTransform`
    2139: . ct  - The type of the parent cell
    2140: . rct - The type of the produced cell
    2141: . r   - The index of the produced cell
    2142: - x   - The localized coordinates for the parent cell

    2144:   Output Parameter:
    2145: . xr  - The localized coordinates for the produced cell

    2147:   Level: developer

    2149: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
    2150: */
    2151: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
    2152: {
    2153:   PetscFE  fe = NULL;
    2154:   PetscInt cdim, v, *subcellV;

    2156:   PetscFunctionBegin;
    2157:   PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe));
    2158:   PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV));
    2159:   PetscCall(PetscFEGetNumComponents(fe, &cdim));
    2160:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]));
    2161:   PetscFunctionReturn(PETSC_SUCCESS);
    2162: }

    2164: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
    2165: {
    2166:   DM                 dm, cdm, cdmCell, cdmNew, cdmCellNew;
    2167:   PetscSection       coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
    2168:   Vec                coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
    2169:   const PetscScalar *coords;
    2170:   PetscScalar       *coordsNew;
    2171:   const PetscReal   *maxCell, *Lstart, *L;
    2172:   PetscBool          localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE, sparseLocalize;
    2173:   PetscInt           dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;

    2175:   PetscFunctionBegin;
    2176:   // Need to clear the DMField for coordinates
    2177:   PetscCall(DMSetCoordinateField(rdm, NULL));
    2178:   PetscCall(DMPlexTransformGetDM(tr, &dm));
    2179:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
    2180:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    2181:   PetscCall(DMGetCellCoordinateDM(dm, &cdmCell));
    2182:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
    2183:   PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
    2184:   PetscCall(DMSetSparseLocalize(rdm, sparseLocalize));
    2185:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
    2186:   if (localized) {
    2187:     /* Localize coordinates of new vertices */
    2188:     localizeVertices = PETSC_TRUE;
    2189:     /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
    2190:     if (!maxCell) localizeCells = PETSC_TRUE;
    2191:   }
    2192:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    2193:   PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo));
    2194:   PetscCall(DMGetCoordinateDim(rdm, &dE));
    2195:   if (maxCell) {
    2196:     PetscReal *LstartNew, *LNew, *maxCellNew;

    2198:     PetscCall(PetscMalloc3(dE, &LstartNew, dE, &LNew, dE, &maxCellNew));
    2199:     for (d = 0; d < dEo; ++d) {
    2200:       LstartNew[d]  = Lstart[d];
    2201:       LNew[d]       = L[d];
    2202:       maxCellNew[d] = maxCell[d] / tr->redFactor;
    2203:     }
    2204:     for (d = dEo; d < dE; ++d) {
    2205:       LstartNew[d]  = 0.;
    2206:       LNew[d]       = -1.;
    2207:       maxCellNew[d] = -1.;
    2208:     }
    2209:     PetscCall(DMSetPeriodicity(rdm, maxCellNew, LstartNew, LNew));
    2210:     PetscCall(PetscFree3(LstartNew, LNew, maxCellNew));
    2211:   }
    2212:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew));
    2213:   PetscCall(PetscSectionSetNumFields(coordSectionNew, 1));
    2214:   PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE));
    2215:   PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew));
    2216:   PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew));
    2217:   /* Localization should be inherited */
    2218:   /*   Stefano calculates parent cells for each new cell for localization */
    2219:   /*   Localized cells need coordinates of closure */
    2220:   for (v = vStartNew; v < vEndNew; ++v) {
    2221:     PetscCall(PetscSectionSetDof(coordSectionNew, v, dE));
    2222:     PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE));
    2223:   }
    2224:   PetscCall(PetscSectionSetUp(coordSectionNew));
    2225:   PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew));

    2227:   if (localizeCells) {
    2228:     PetscCall(DMGetCoordinateDM(rdm, &cdmNew));
    2229:     PetscCall(DMClone(cdmNew, &cdmCellNew));
    2230:     PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew));
    2231:     PetscCall(DMDestroy(&cdmCellNew));

    2233:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew));
    2234:     PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1));
    2235:     PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE));
    2236:     PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew));
    2237:     PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew));

    2239:     PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell));
    2240:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    2241:     for (c = cStart; c < cEnd; ++c) {
    2242:       PetscInt dof;

    2244:       PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof));
    2245:       if (dof) {
    2246:         DMPolytopeType  ct;
    2247:         DMPolytopeType *rct;
    2248:         PetscInt       *rsize, *rcone, *rornt;
    2249:         PetscInt        dim, cNew, Nct, n, r;

    2251:         PetscCall(DMPlexGetCellType(dm, c, &ct));
    2252:         dim = DMPolytopeTypeGetDim(ct);
    2253:         PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2254:         /* This allows for different cell types */
    2255:         for (n = 0; n < Nct; ++n) {
    2256:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
    2257:           for (r = 0; r < rsize[n]; ++r) {
    2258:             PetscInt *closure = NULL;
    2259:             PetscInt  clSize, cl, Nv = 0;

    2261:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew));
    2262:             PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
    2263:             for (cl = 0; cl < clSize * 2; cl += 2) {
    2264:               if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
    2265:             }
    2266:             PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure));
    2267:             PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE));
    2268:             PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE));
    2269:           }
    2270:         }
    2271:       }
    2272:     }
    2273:     PetscCall(PetscSectionSetUp(coordSectionCellNew));
    2274:     PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew));
    2275:   }
    2276:   PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view"));
    2277:   {
    2278:     VecType     vtype;
    2279:     PetscInt    coordSizeNew, bs;
    2280:     const char *name;

    2282:     PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
    2283:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew));
    2284:     PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew));
    2285:     PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE));
    2286:     PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name));
    2287:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name));
    2288:     PetscCall(VecGetBlockSize(coordsLocal, &bs));
    2289:     PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE));
    2290:     PetscCall(VecGetType(coordsLocal, &vtype));
    2291:     PetscCall(VecSetType(coordsLocalNew, vtype));
    2292:   }
    2293:   PetscCall(VecGetArrayRead(coordsLocal, &coords));
    2294:   PetscCall(VecGetArray(coordsLocalNew, &coordsNew));
    2295:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    2296:   /* First set coordinates for vertices */
    2297:   for (p = pStart; p < pEnd; ++p) {
    2298:     DMPolytopeType  ct;
    2299:     DMPolytopeType *rct;
    2300:     PetscInt       *rsize, *rcone, *rornt;
    2301:     PetscInt        Nct, n, r;
    2302:     PetscBool       hasVertex = PETSC_FALSE;

    2304:     PetscCall(DMPlexGetCellType(dm, p, &ct));
    2305:     PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2306:     for (n = 0; n < Nct; ++n) {
    2307:       if (rct[n] == DM_POLYTOPE_POINT) {
    2308:         hasVertex = PETSC_TRUE;
    2309:         break;
    2310:       }
    2311:     }
    2312:     if (hasVertex) {
    2313:       const PetscScalar *icoords = NULL;
    2314:       const PetscScalar *array   = NULL;
    2315:       PetscScalar       *pcoords = NULL;
    2316:       PetscBool          isDG;
    2317:       PetscInt           Nc, Nv, v, d;

    2319:       PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));

    2321:       icoords = pcoords;
    2322:       Nv      = Nc / dEo;
    2323:       if (ct != DM_POLYTOPE_POINT) {
    2324:         if (localizeVertices && maxCell) {
    2325:           PetscScalar anchor[3];

    2327:           for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
    2328:           for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]));
    2329:         }
    2330:       }
    2331:       for (n = 0; n < Nct; ++n) {
    2332:         if (rct[n] != DM_POLYTOPE_POINT) continue;
    2333:         for (r = 0; r < rsize[n]; ++r) {
    2334:           PetscScalar vcoords[3];
    2335:           PetscInt    vNew, off;

    2337:           PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew));
    2338:           PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off));
    2339:           PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords));
    2340:           PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]));
    2341:         }
    2342:       }
    2343:       PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords));
    2344:     }
    2345:   }
    2346:   PetscCall(VecRestoreArrayRead(coordsLocal, &coords));
    2347:   PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew));
    2348:   PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew));
    2349:   PetscCall(VecDestroy(&coordsLocalNew));
    2350:   PetscCall(PetscSectionDestroy(&coordSectionNew));
    2351:   /* Then set coordinates for cells by localizing */
    2352:   if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm));
    2353:   else {
    2354:     VecType     vtype;
    2355:     PetscInt    coordSizeNew, bs;
    2356:     const char *name;

    2358:     PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell));
    2359:     PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew));
    2360:     PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew));
    2361:     PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE));
    2362:     PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name));
    2363:     PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name));
    2364:     PetscCall(VecGetBlockSize(coordsLocalCell, &bs));
    2365:     PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE));
    2366:     PetscCall(VecGetType(coordsLocalCell, &vtype));
    2367:     PetscCall(VecSetType(coordsLocalCellNew, vtype));
    2368:     PetscCall(VecGetArrayRead(coordsLocalCell, &coords));
    2369:     PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew));

    2371:     for (p = pStart; p < pEnd; ++p) {
    2372:       DMPolytopeType  ct;
    2373:       DMPolytopeType *rct;
    2374:       PetscInt       *rsize, *rcone, *rornt;
    2375:       PetscInt        dof = 0, Nct, n, r;

    2377:       PetscCall(DMPlexGetCellType(dm, p, &ct));
    2378:       PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    2379:       if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof));
    2380:       if (dof) {
    2381:         const PetscScalar *pcoords;

    2383:         PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords));
    2384:         for (n = 0; n < Nct; ++n) {
    2385:           const PetscInt Nr = rsize[n];

    2387:           if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
    2388:           for (r = 0; r < Nr; ++r) {
    2389:             PetscInt pNew, offNew;

    2391:             /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
    2392:                DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
    2393:                cell to the ones it produces. */
    2394:             PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew));
    2395:             PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew));
    2396:             PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]));
    2397:           }
    2398:         }
    2399:       }
    2400:     }
    2401:     PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords));
    2402:     PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew));
    2403:     PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew));
    2404:     PetscCall(VecDestroy(&coordsLocalCellNew));
    2405:     PetscCall(PetscSectionDestroy(&coordSectionCellNew));
    2406:   }
    2407:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0));
    2408:   PetscFunctionReturn(PETSC_SUCCESS);
    2409: }

    2411: /*@
    2412:   DMPlexTransformApply - Execute the transformation, producing another `DM`

    2414:   Collective

    2416:   Input Parameters:
    2417: + tr - The `DMPlexTransform` object
    2418: - dm - The original `DM`

    2420:   Output Parameter:
    2421: . trdm - The transformed `DM`

    2423:   Level: intermediate

    2425:   Options Database Keys:
    2426: + -dm_plex_transform_label_match_strata    - Only label points of the same stratum as the producing point
    2427: . -dm_plex_transform_label_replica_inc num - Increment for the label value to be multiplied by the replica number
    2428: - -dm_plex_transform_active name           - Name for active mesh label

    2430: .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()`
    2431: @*/
    2432: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *trdm)
    2433: {
    2434:   DM                     rdm;
    2435:   DMPlexInterpolatedFlag interp;
    2436:   PetscInt               pStart, pEnd;

    2438:   PetscFunctionBegin;
    2441:   PetscAssertPointer(trdm, 3);
    2442:   PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
    2443:   PetscCall(DMPlexTransformSetDM(tr, dm));

    2445:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm));
    2446:   PetscCall(DMSetType(rdm, DMPLEX));
    2447:   PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm));
    2448:   /* Calculate number of new points of each depth */
    2449:   PetscCall(DMPlexIsInterpolatedCollective(dm, &interp));
    2450:   PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
    2451:   /* Step 1: Set chart */
    2452:   PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd));
    2453:   PetscCall(DMPlexSetChart(rdm, pStart, pEnd));
    2454:   /* Step 2: Set cone/support sizes (automatically stratifies) */
    2455:   PetscCall(DMPlexTransformSetConeSizes(tr, rdm));
    2456:   /* Step 3: Setup refined DM */
    2457:   PetscCall(DMSetUp(rdm));
    2458:   /* Step 4: Set cones and supports (automatically symmetrizes) */
    2459:   PetscCall(DMPlexTransformSetCones(tr, rdm));
    2460:   /* Step 5: Create pointSF */
    2461:   PetscCall(DMPlexTransformCreateSF(tr, rdm));
    2462:   /* Step 6: Create labels */
    2463:   PetscCall(DMPlexTransformCreateLabels(tr, rdm));
    2464:   /* Step 7: Set coordinates */
    2465:   PetscCall(DMPlexTransformSetCoordinates(tr, rdm));
    2466:   //   Do not copy periodicity, which was handled in DMPlexTransformSetCoordinates()
    2467:   PetscCall(DMPlexCopy_Internal(dm, PETSC_FALSE, PETSC_TRUE, rdm));
    2468:   // If the original DM was configured from options, the transformed DM should be as well
    2469:   rdm->setfromoptionscalled = dm->setfromoptionscalled;
    2470:   PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0));
    2471:   *trdm = rdm;
    2472:   PetscFunctionReturn(PETSC_SUCCESS);
    2473: }

    2475: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
    2476: {
    2477:   DMPlexTransform tr;
    2478:   DM              cdm, rcdm;
    2479:   const char     *prefix;
    2480:   PetscBool       save;

    2482:   PetscFunctionBegin;
    2483:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
    2484:   PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform"));
    2485:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    2486:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
    2487:   PetscCall(DMPlexTransformSetDM(tr, dm));
    2488:   PetscCall(DMPlexTransformSetFromOptions(tr));
    2489:   if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
    2490:   PetscCall(DMPlexTransformSetUp(tr));
    2491:   PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
    2492:   PetscCall(DMPlexTransformApply(tr, dm, rdm));
    2493:   PetscCall(DMCopyDisc(dm, *rdm));
    2494:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    2495:   PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
    2496:   PetscCall(DMCopyDisc(cdm, rcdm));
    2497:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
    2498:   PetscCall(DMCopyDisc(dm, *rdm));
    2499:   PetscCall(DMPlexGetSaveTransform(dm, &save));
    2500:   if (save) PetscCall(DMPlexSetTransform(*rdm, tr));
    2501:   PetscCall(DMPlexTransformDestroy(&tr));
    2502:   ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
    2503:   PetscFunctionReturn(PETSC_SUCCESS);
    2504: }