Actual source code: ex11.c

      1: static char help[] = "Mesh Orientation Tutorial\n\n";

      3: #include <petscdmplex.h>
      4: #include <petscdmplextransform.h>

      6: typedef struct {
      7:   PetscBool genArr;        /* Generate all possible cell arrangements */
      8:   PetscBool refArr;        /* Refine all possible cell arrangements */
      9:   PetscBool printTable;    /* Print the CAyley table */
     10:   PetscInt  orntBounds[2]; /* Bounds for the orientation check */
     11:   PetscInt  numOrnt;       /* Number of specific orientations specified, or -1 for all orientations */
     12:   PetscInt  ornts[48];     /* Specific orientations if specified */
     13:   PetscInt  initOrnt;      /* Initial orientation for starting mesh */
     14: } AppCtx;

     16: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
     17: {
     18:   PetscInt  n = 2;
     19:   PetscBool flg;

     21:   PetscFunctionBeginUser;
     22:   options->genArr        = PETSC_FALSE;
     23:   options->refArr        = PETSC_FALSE;
     24:   options->printTable    = PETSC_FALSE;
     25:   options->orntBounds[0] = PETSC_INT_MIN;
     26:   options->orntBounds[1] = PETSC_INT_MAX;
     27:   options->numOrnt       = -1;
     28:   options->initOrnt      = 0;

     30:   PetscOptionsBegin(comm, "", "Mesh Orientation Tutorials Options", "DMPLEX");
     31:   PetscCall(PetscOptionsBool("-gen_arrangements", "Flag for generating all arrangements of the cell", "ex11.c", options->genArr, &options->genArr, NULL));
     32:   PetscCall(PetscOptionsBool("-ref_arrangements", "Flag for refining all arrangements of the cell", "ex11.c", options->refArr, &options->refArr, NULL));
     33:   PetscCall(PetscOptionsBool("-print_table", "Print the Cayley table", "ex11.c", options->printTable, &options->printTable, NULL));
     34:   PetscCall(PetscOptionsIntArray("-ornt_bounds", "Bounds for orientation checks", "ex11.c", options->orntBounds, &n, NULL));
     35:   n = 48;
     36:   PetscCall(PetscOptionsIntArray("-ornts", "Specific orientations for checks", "ex11.c", options->ornts, &n, &flg));
     37:   if (flg) {
     38:     options->numOrnt = n;
     39:     PetscCall(PetscSortInt(n, options->ornts));
     40:   }
     41:   PetscCall(PetscOptionsInt("-init_ornt", "Initial orientation for starting mesh", "ex11.c", options->initOrnt, &options->initOrnt, NULL));
     42:   PetscOptionsEnd();
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: static PetscBool ignoreOrnt(AppCtx *user, PetscInt o)
     47: {
     48:   PetscInt       loc;
     49:   PetscErrorCode ierr;

     51:   if (user->numOrnt < 0) return PETSC_FALSE;
     52:   ierr = PetscFindInt(o, user->numOrnt, user->ornts, &loc);
     53:   if (loc < 0 || ierr) return PETSC_TRUE;
     54:   return PETSC_FALSE;
     55: }

     57: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
     58: {
     59:   PetscFunctionBeginUser;
     60:   PetscCall(DMCreate(comm, dm));
     61:   PetscCall(DMSetType(*dm, DMPLEX));
     62:   PetscCall(DMSetFromOptions(*dm));
     63:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
     64:   PetscFunctionReturn(PETSC_SUCCESS);
     65: }

     67: static PetscErrorCode CheckCellVertices(DM dm, PetscInt cell, PetscInt o)
     68: {
     69:   DMPolytopeType  ct;
     70:   const PetscInt *arrVerts;
     71:   PetscInt       *closure = NULL;
     72:   PetscInt        Ncl, cl, Nv, vStart, vEnd, v;
     73:   MPI_Comm        comm;

     75:   PetscFunctionBeginUser;
     76:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
     77:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
     78:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
     79:   PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
     80:   for (cl = 0, Nv = 0; cl < Ncl * 2; cl += 2) {
     81:     const PetscInt vertex = closure[cl];

     83:     if (vertex < vStart || vertex >= vEnd) continue;
     84:     closure[Nv++] = vertex;
     85:   }
     86:   PetscCheck(Nv == DMPolytopeTypeGetNumVertices(ct), comm, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " has %" PetscInt_FMT " vertices != %" PetscInt_FMT " vertices in a %s", cell, Nv, DMPolytopeTypeGetNumVertices(ct), DMPolytopeTypes[ct]);
     87:   arrVerts = DMPolytopeTypeGetVertexArrangement(ct, o);
     88:   for (v = 0; v < Nv; ++v) {
     89:     PetscCheck(closure[v] == arrVerts[v] + vStart, comm, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " vertex[%" PetscInt_FMT "]: %" PetscInt_FMT " should be %" PetscInt_FMT " for arrangement %" PetscInt_FMT, cell, v, closure[v], arrVerts[v] + vStart, o);
     90:   }
     91:   PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }

     95: /* Transform cell with group operation o */
     96: static PetscErrorCode ReorientCell(DM dm, PetscInt cell, PetscInt o, PetscBool swapCoords)
     97: {
     98:   DM           cdm;
     99:   Vec          coordinates;
    100:   PetscScalar *coords, *ccoords = NULL;
    101:   PetscInt    *closure = NULL;
    102:   PetscInt     cdim, d, Nc, Ncl, cl, vStart, vEnd, Nv;

    104:   PetscFunctionBegin;
    105:   /* Change vertex coordinates so that it plots as we expect */
    106:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    107:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    108:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    109:   PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinates, cell, &Nc, &ccoords));
    110:   /* Reorient cone */
    111:   PetscCall(DMPlexOrientPoint(dm, cell, o));
    112:   /* Finish resetting coordinates */
    113:   if (swapCoords) {
    114:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    115:     PetscCall(VecGetArrayWrite(coordinates, &coords));
    116:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
    117:     for (cl = 0, Nv = 0; cl < Ncl * 2; cl += 2) {
    118:       const PetscInt vertex = closure[cl];
    119:       PetscScalar   *vcoords;

    121:       if (vertex < vStart || vertex >= vEnd) continue;
    122:       PetscCall(DMPlexPointLocalRef(cdm, vertex, coords, &vcoords));
    123:       for (d = 0; d < cdim; ++d) vcoords[d] = ccoords[Nv * cdim + d];
    124:       ++Nv;
    125:     }
    126:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
    127:     PetscCall(VecRestoreArrayWrite(coordinates, &coords));
    128:   }
    129:   PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinates, cell, &Nc, &ccoords));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: static PetscErrorCode GenerateArrangements(DM dm, AppCtx *user)
    134: {
    135:   DM             odm;
    136:   DMPolytopeType ct;
    137:   PetscInt       No, o;
    138:   const char    *name;

    140:   PetscFunctionBeginUser;
    141:   if (!user->genArr) PetscFunctionReturn(PETSC_SUCCESS);
    142:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    143:   PetscCall(DMPlexGetCellType(dm, 0, &ct));
    144:   No = DMPolytopeTypeGetNumArrangements(ct) / 2;
    145:   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
    146:     if (ignoreOrnt(user, o)) continue;
    147:     PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &odm));
    148:     PetscCall(ReorientCell(odm, 0, o, PETSC_TRUE));
    149:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "%s orientation %" PetscInt_FMT "\n", name, o));
    150:     PetscCall(DMViewFromOptions(odm, NULL, "-gen_dm_view"));
    151:     PetscCall(CheckCellVertices(odm, 0, o));
    152:     PetscCall(DMDestroy(&odm));
    153:   }
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: static PetscErrorCode VerifyCayleyTable(DM dm, AppCtx *user)
    158: {
    159:   DM              dm1, dm2;
    160:   DMPolytopeType  ct;
    161:   const PetscInt *refcone, *cone;
    162:   PetscInt        No, o1, o2, o3, o4;
    163:   PetscBool       equal;
    164:   const char     *name;

    166:   PetscFunctionBeginUser;
    167:   if (!user->genArr) PetscFunctionReturn(PETSC_SUCCESS);
    168:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    169:   PetscCall(DMPlexGetCellType(dm, 0, &ct));
    170:   PetscCall(DMPlexGetCone(dm, 0, &refcone));
    171:   No = DMPolytopeTypeGetNumArrangements(ct) / 2;
    172:   if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cayley Table for %s\n", DMPolytopeTypes[ct]));
    173:   for (o1 = PetscMax(-No, user->orntBounds[0]); o1 < PetscMin(No, user->orntBounds[1]); ++o1) {
    174:     for (o2 = PetscMax(-No, user->orntBounds[0]); o2 < PetscMin(No, user->orntBounds[1]); ++o2) {
    175:       PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &dm1));
    176:       PetscCall(DMPlexOrientPoint(dm1, 0, o2));
    177:       PetscCall(DMPlexCheckFaces(dm1, 0));
    178:       PetscCall(DMPlexOrientPoint(dm1, 0, o1));
    179:       PetscCall(DMPlexCheckFaces(dm1, 0));
    180:       o3 = DMPolytopeTypeComposeOrientation(ct, o1, o2);
    181:       /* First verification */
    182:       PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &dm2));
    183:       PetscCall(DMPlexOrientPoint(dm2, 0, o3));
    184:       PetscCall(DMPlexCheckFaces(dm2, 0));
    185:       PetscCall(DMPlexEqual(dm1, dm2, &equal));
    186:       if (!equal) {
    187:         PetscCall(DMViewFromOptions(dm1, NULL, "-error_dm_view"));
    188:         PetscCall(DMViewFromOptions(dm2, NULL, "-error_dm_view"));
    189:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cayley table error for %s: %" PetscInt_FMT " * %" PetscInt_FMT " != %" PetscInt_FMT, DMPolytopeTypes[ct], o1, o2, o3);
    190:       }
    191:       /* Second verification */
    192:       PetscCall(DMPlexGetCone(dm1, 0, &cone));
    193:       PetscCall(DMPolytopeGetOrientation(ct, refcone, cone, &o4));
    194:       if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT ", ", o4));
    195:       PetscCheck(o3 == o4, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cayley table error for %s: %" PetscInt_FMT " * %" PetscInt_FMT " = %" PetscInt_FMT " != %" PetscInt_FMT, DMPolytopeTypes[ct], o1, o2, o3, o4);
    196:       PetscCall(DMDestroy(&dm1));
    197:       PetscCall(DMDestroy(&dm2));
    198:     }
    199:     if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    200:   }
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: static PetscErrorCode VerifyInverse(DM dm, AppCtx *user)
    205: {
    206:   DM              dm1, dm2;
    207:   DMPolytopeType  ct;
    208:   const PetscInt *refcone, *cone;
    209:   PetscInt        No, o, oi, o2;
    210:   PetscBool       equal;
    211:   const char     *name;

    213:   PetscFunctionBeginUser;
    214:   if (!user->genArr) PetscFunctionReturn(PETSC_SUCCESS);
    215:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    216:   PetscCall(DMPlexGetCellType(dm, 0, &ct));
    217:   PetscCall(DMPlexGetCone(dm, 0, &refcone));
    218:   No = DMPolytopeTypeGetNumArrangements(ct) / 2;
    219:   if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Inverse table for %s\n", DMPolytopeTypes[ct]));
    220:   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
    221:     if (ignoreOrnt(user, o)) continue;
    222:     oi = DMPolytopeTypeComposeOrientationInv(ct, 0, o);
    223:     PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &dm1));
    224:     PetscCall(DMPlexOrientPoint(dm1, 0, o));
    225:     PetscCall(DMPlexCheckFaces(dm1, 0));
    226:     PetscCall(DMPlexOrientPoint(dm1, 0, oi));
    227:     PetscCall(DMPlexCheckFaces(dm1, 0));
    228:     /* First verification */
    229:     PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &dm2));
    230:     PetscCall(DMPlexEqual(dm1, dm2, &equal));
    231:     if (!equal) {
    232:       PetscCall(DMViewFromOptions(dm1, NULL, "-error_dm_view"));
    233:       PetscCall(DMViewFromOptions(dm2, NULL, "-error_dm_view"));
    234:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inverse error for %s: %" PetscInt_FMT " * %" PetscInt_FMT " != 0", DMPolytopeTypes[ct], o, oi);
    235:     }
    236:     /* Second verification */
    237:     PetscCall(DMPlexGetCone(dm1, 0, &cone));
    238:     PetscCall(DMPolytopeGetOrientation(ct, refcone, cone, &o2));
    239:     if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT ", ", oi));
    240:     PetscCheck(o2 == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inverse error for %s: %" PetscInt_FMT " * %" PetscInt_FMT " = %" PetscInt_FMT " != 0", DMPolytopeTypes[ct], o, oi, o2);
    241:     PetscCall(DMDestroy(&dm1));
    242:     PetscCall(DMDestroy(&dm2));
    243:   }
    244:   if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: /* Suppose that point p has the same arrangement as o from canonical, compare the subcells to canonical subcells */
    249: static PetscErrorCode CheckSubcells(DM dm, DM odm, PetscInt p, PetscInt o, AppCtx *user)
    250: {
    251:   DMPlexTransform tr, otr;
    252:   DMPolytopeType  ct;
    253:   DMPolytopeType *rct;
    254:   const PetscInt *cone, *ornt, *ocone, *oornt;
    255:   PetscInt       *rsize, *rcone, *rornt;
    256:   PetscInt        Nct, n, oi, debug = 0;

    258:   PetscFunctionBeginUser;
    259:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
    260:   PetscCall(DMPlexTransformSetDM(tr, dm));
    261:   PetscCall(DMPlexTransformSetFromOptions(tr));
    262:   PetscCall(DMPlexTransformSetUp(tr));

    264:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)odm), &otr));
    265:   PetscCall(DMPlexTransformSetDM(otr, odm));
    266:   PetscCall(DMPlexTransformSetFromOptions(otr));
    267:   PetscCall(DMPlexTransformSetUp(otr));

    269:   PetscCall(DMPlexGetCellType(dm, p, &ct));
    270:   PetscCall(DMPlexGetCone(dm, p, &cone));
    271:   PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
    272:   PetscCall(DMPlexGetCone(odm, p, &ocone));
    273:   PetscCall(DMPlexGetConeOrientation(odm, p, &oornt));
    274:   oi = DMPolytopeTypeComposeOrientationInv(ct, 0, o);
    275:   if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Orientation %" PetscInt_FMT "\n", oi));

    277:   PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
    278:   for (n = 0; n < Nct; ++n) {
    279:     DMPolytopeType ctNew = rct[n];
    280:     PetscInt       r, ro;

    282:     if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Checking type %s\n", DMPolytopeTypes[ctNew]));
    283:     for (r = 0; r < rsize[n]; ++r) {
    284:       const PetscInt *qcone, *qornt, *oqcone, *oqornt;
    285:       PetscInt        pNew, opNew, oo, pr, fo;
    286:       PetscBool       restore = PETSC_TRUE;

    288:       PetscCall(DMPlexTransformGetTargetPoint(tr, ct, ctNew, p, r, &pNew));
    289:       PetscCall(DMPlexTransformGetCone(tr, pNew, &qcone, &qornt));
    290:       if (debug) {
    291:         PetscInt c;

    293:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Checking replica %" PetscInt_FMT " (%" PetscInt_FMT ")\n      Original Cone", r, pNew));
    294:         for (c = 0; c < DMPolytopeTypeGetConeSize(ctNew); ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscInt_FMT " (%" PetscInt_FMT ")", qcone[c], qornt[c]));
    295:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    296:       }
    297:       for (ro = 0; ro < rsize[n]; ++ro) {
    298:         PetscBool found;

    300:         PetscCall(DMPlexTransformGetTargetPoint(otr, ct, ctNew, p, ro, &opNew));
    301:         PetscCall(DMPlexTransformGetConeOriented(otr, opNew, o, &oqcone, &oqornt));
    302:         PetscCall(DMPolytopeMatchOrientation(ctNew, oqcone, qcone, &oo, &found));
    303:         if (found) break;
    304:         PetscCall(DMPlexTransformRestoreCone(otr, pNew, &oqcone, &oqornt));
    305:       }
    306:       if (debug) {
    307:         PetscInt c;

    309:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Checking transform replica %" PetscInt_FMT " (%" PetscInt_FMT ") (%" PetscInt_FMT ")\n      Transform Cone", ro, opNew, o));
    310:         for (c = 0; c < DMPolytopeTypeGetConeSize(ctNew); ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscInt_FMT " (%" PetscInt_FMT ")", oqcone[c], oqornt[c]));
    311:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    312:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Matched %" PetscInt_FMT "\n", oo));
    313:       }
    314:       if (ro == rsize[n]) {
    315:         /* The tetrahedron has 3 pairs of opposing edges, and any pair can be connected by the interior segment */
    316:         if (ct == DM_POLYTOPE_TETRAHEDRON) {
    317:           /* The segment in a tetrahedron does not map into itself under the group action */
    318:           if (ctNew == DM_POLYTOPE_SEGMENT) {
    319:             restore = PETSC_FALSE;
    320:             ro      = r;
    321:             oo      = 0;
    322:           }
    323:           /* The last four interior faces do not map into themselves under the group action */
    324:           if (r > 3 && ctNew == DM_POLYTOPE_TRIANGLE) {
    325:             restore = PETSC_FALSE;
    326:             ro      = r;
    327:             oo      = 0;
    328:           }
    329:           /* The last four interior faces do not map into themselves under the group action */
    330:           if (r > 3 && ctNew == DM_POLYTOPE_TETRAHEDRON) {
    331:             restore = PETSC_FALSE;
    332:             ro      = r;
    333:             oo      = 0;
    334:           }
    335:         }
    336:         PetscCheck(!restore, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to find matching %s %" PetscInt_FMT " orientation for cell orientation %" PetscInt_FMT, DMPolytopeTypes[ctNew], r, o);
    337:       }
    338:       if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT ", %" PetscInt_FMT ", ", ro, oo));
    339:       PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, oi, ctNew, r, 0, &pr, &fo));
    340:       if (!user->printTable) {
    341:         PetscCheck(pr == ro, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Choose wrong replica %" PetscInt_FMT " != %" PetscInt_FMT, pr, ro);
    342:         PetscCheck(fo == oo, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Choose wrong orientation %" PetscInt_FMT " != %" PetscInt_FMT, fo, oo);
    343:       }
    344:       PetscCall(DMPlexTransformRestoreCone(tr, pNew, &qcone, &qornt));
    345:       if (restore) PetscCall(DMPlexTransformRestoreCone(otr, pNew, &oqcone, &oqornt));
    346:     }
    347:     if (user->printTable) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    348:   }
    349:   PetscCall(DMPlexTransformDestroy(&tr));
    350:   PetscCall(DMPlexTransformDestroy(&otr));
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: static PetscErrorCode RefineArrangements(DM dm, AppCtx *user)
    355: {
    356:   DM             odm, rdm;
    357:   DMPolytopeType ct;
    358:   PetscInt       No, o;
    359:   const char    *name;

    361:   PetscFunctionBeginUser;
    362:   if (!user->refArr) PetscFunctionReturn(PETSC_SUCCESS);
    363:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    364:   PetscCall(DMPlexGetCellType(dm, 0, &ct));
    365:   No = DMPolytopeTypeGetNumArrangements(ct) / 2;
    366:   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
    367:     if (ignoreOrnt(user, o)) continue;
    368:     PetscCall(CreateMesh(PetscObjectComm((PetscObject)dm), user, &odm));
    369:     if (user->initOrnt) PetscCall(ReorientCell(odm, 0, user->initOrnt, PETSC_FALSE));
    370:     PetscCall(ReorientCell(odm, 0, o, PETSC_TRUE));
    371:     PetscCall(DMViewFromOptions(odm, NULL, "-orig_dm_view"));
    372:     PetscCall(DMRefine(odm, MPI_COMM_NULL, &rdm));
    373:     PetscCall(DMSetFromOptions(rdm));
    374:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "%s orientation %" PetscInt_FMT "\n", name, o));
    375:     PetscCall(DMViewFromOptions(rdm, NULL, "-ref_dm_view"));
    376:     PetscCall(CheckSubcells(dm, odm, 0, o, user));
    377:     PetscCall(DMDestroy(&odm));
    378:     PetscCall(DMDestroy(&rdm));
    379:   }
    380:   PetscFunctionReturn(PETSC_SUCCESS);
    381: }

    383: int main(int argc, char **argv)
    384: {
    385:   DM     dm;
    386:   AppCtx user;

    388:   PetscFunctionBeginUser;
    389:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    390:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    391:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    392:   if (user.initOrnt) {
    393:     PetscCall(ReorientCell(dm, 0, user.initOrnt, PETSC_FALSE));
    394:     PetscCall(DMViewFromOptions(dm, NULL, "-ornt_dm_view"));
    395:   }
    396:   PetscCall(GenerateArrangements(dm, &user));
    397:   PetscCall(VerifyCayleyTable(dm, &user));
    398:   PetscCall(VerifyInverse(dm, &user));
    399:   PetscCall(RefineArrangements(dm, &user));
    400:   PetscCall(DMDestroy(&dm));
    401:   PetscCall(PetscFinalize());
    402:   return 0;
    403: }

    405: /*TEST

    407:   testset:
    408:     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -gen_arrangements \
    409:           -gen_dm_view ::ascii_latex -dm_plex_view_tikzscale 0.5

    411:     test:
    412:       suffix: segment
    413:       args: -dm_plex_cell segment \
    414:             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0

    416:     test:
    417:       suffix: triangle
    418:       args: -dm_plex_cell triangle \
    419:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    421:     test:
    422:       suffix: quadrilateral
    423:       args: -dm_plex_cell quadrilateral \
    424:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    426:     test:
    427:       suffix: tensor_segment
    428:       args: -dm_plex_cell tensor_quad \
    429:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    431:     test:
    432:       suffix: tetrahedron
    433:       args: -dm_plex_cell tetrahedron \
    434:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    436:     test:
    437:       suffix: hexahedron
    438:       args: -dm_plex_cell hexahedron \
    439:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 0.3

    441:     test:
    442:       suffix: triangular_prism
    443:       args: -dm_plex_cell triangular_prism \
    444:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    446:     test:
    447:       suffix: tensor_triangular_prism
    448:       args: -dm_plex_cell tensor_triangular_prism \
    449:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    451:     test:
    452:       suffix: tensor_quadrilateral_prism
    453:       args: -dm_plex_cell tensor_quadrilateral_prism \
    454:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    456:     test:
    457:       suffix: pyramid
    458:       args: -dm_plex_cell pyramid \
    459:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    461:   testset:
    462:     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -ref_arrangements -dm_plex_check_all \
    463:           -ref_dm_view ::ascii_latex -dm_plex_view_tikzscale 1.0

    465:     test:
    466:       suffix: ref_segment
    467:       args: -dm_plex_cell segment \
    468:             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0

    470:     test:
    471:       suffix: ref_triangle
    472:       args: -dm_plex_cell triangle \
    473:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    475:     test:
    476:       suffix: ref_quadrilateral
    477:       args: -dm_plex_cell quadrilateral \
    478:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    480:     test:
    481:       suffix: ref_tensor_segment
    482:       args: -dm_plex_cell tensor_quad \
    483:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

    485:     test:
    486:       suffix: ref_tetrahedron
    487:       args: -dm_plex_cell tetrahedron \
    488:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    490:     test:
    491:       suffix: ref_hexahedron
    492:       args: -dm_plex_cell hexahedron \
    493:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    495:     test:
    496:       suffix: ref_triangular_prism
    497:       args: -dm_plex_cell triangular_prism \
    498:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    500:     test:
    501:       suffix: ref_tensor_triangular_prism
    502:       args: -dm_plex_cell tensor_triangular_prism \
    503:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    505:     test:
    506:       suffix: ref_tensor_quadrilateral_prism
    507:       args: -dm_plex_cell tensor_quadrilateral_prism \
    508:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

    510:   # ToBox should recreate the coordinate space since the cell shape changes
    511:   testset:
    512:     args: -dm_coord_space 0 -dm_plex_transform_type refine_tobox -ref_arrangements -dm_plex_check_all

    514:     test:
    515:       suffix: tobox_triangle
    516:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
    517:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    519:     test:
    520:       suffix: tobox_tensor_segment
    521:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
    522:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    524:     test:
    525:       suffix: tobox_tetrahedron
    526:       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
    527:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    529:     test:
    530:       suffix: tobox_triangular_prism
    531:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
    532:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    534:     test:
    535:       suffix: tobox_tensor_triangular_prism
    536:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
    537:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    539:     test:
    540:       suffix: tobox_tensor_quadrilateral_prism
    541:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quadrilateral_prism \
    542:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    544:   testset:
    545:     args: -dm_coord_space 0 -dm_plex_transform_type refine_alfeld -ref_arrangements -dm_plex_check_all

    547:     test:
    548:       suffix: alfeld_triangle
    549:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
    550:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    552:     test:
    553:       suffix: alfeld_tetrahedron
    554:       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
    555:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    557:   testset:
    558:     args: -dm_plex_transform_type refine_sbr -ref_arrangements -dm_plex_check_all

    560:     # This splits edge 1 of the triangle, and reflects about the added edge
    561:     test:
    562:       suffix: sbr_triangle_0
    563:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -ornts -2,0 \
    564:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    566:     # This splits edge 0 of the triangle, and reflects about the added edge
    567:     test:
    568:       suffix: sbr_triangle_1
    569:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 1 -ornts -3,0 \
    570:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    572:     # This splits edge 2 of the triangle, and reflects about the added edge
    573:     test:
    574:       suffix: sbr_triangle_2
    575:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 2 -ornts -1,0 \
    576:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    578:     # This splits edges 1 and 2 of the triangle
    579:     test:
    580:       suffix: sbr_triangle_3
    581:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -ornts 0 \
    582:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    584:     # This splits edges 1 and 0 of the triangle
    585:     test:
    586:       suffix: sbr_triangle_4
    587:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -ornts 0 \
    588:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    590:     # This splits edges 0 and 1 of the triangle
    591:     test:
    592:       suffix: sbr_triangle_5
    593:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 1 -ornts 0 \
    594:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    596:     # This splits edges 0 and 2 of the triangle
    597:     test:
    598:       suffix: sbr_triangle_6
    599:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 1 -ornts 0 \
    600:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    602:     # This splits edges 2 and 0 of the triangle
    603:     test:
    604:       suffix: sbr_triangle_7
    605:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 2 -ornts 0 \
    606:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    608:     # This splits edges 2 and 1 of the triangle
    609:     test:
    610:       suffix: sbr_triangle_8
    611:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 2 -ornts 0 \
    612:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    614:   testset:
    615:     args: -dm_plex_transform_type refine_boundary_layer -dm_plex_transform_bl_splits 2 -ref_arrangements -dm_plex_check_all

    617:     test:
    618:       suffix: bl_tensor_segment
    619:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
    620:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0 -dm_plex_view_tikzscale 1.0

    622:     # The subcell check is broken because at orientation 3, the internal triangles do not get properly permuted for the check
    623:     test:
    624:       suffix: bl_tensor_triangular_prism
    625:       requires: TODO
    626:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
    627:             -ref_dm_view ::ascii_latex -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 1.0

    629: TEST*/