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       ro;

282:     if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Checking type %s\n", DMPolytopeTypes[ctNew]));
283:     for (PetscInt 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:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Checking replica %" PetscInt_FMT " (%" PetscInt_FMT ")\n      Original Cone", r, pNew));
292:         for (PetscInt c = 0; c < DMPolytopeTypeGetConeSize(ctNew); ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscInt_FMT " (%" PetscInt_FMT ")", qcone[c], qornt[c]));
293:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
294:       }
295:       for (ro = 0; ro < rsize[n]; ++ro) {
296:         PetscBool found;

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

350: static PetscErrorCode RefineArrangements(DM dm, AppCtx *user)
351: {
352:   DM             odm, rdm;
353:   DMPolytopeType ct;
354:   PetscInt       No, o;
355:   const char    *name;

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

379: int main(int argc, char **argv)
380: {
381:   DM     dm;
382:   AppCtx user;

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

401: /*TEST

403:   testset:
404:     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -gen_arrangements \
405:           -gen_dm_view ::ascii_latex -dm_plex_view_tikzscale 0.5

407:     test:
408:       suffix: segment
409:       args: -dm_plex_cell segment \
410:             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0

412:     test:
413:       suffix: triangle
414:       args: -dm_plex_cell triangle \
415:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

417:     test:
418:       suffix: quadrilateral
419:       args: -dm_plex_cell quadrilateral \
420:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

422:     test:
423:       suffix: tensor_segment
424:       args: -dm_plex_cell tensor_quad \
425:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

427:     test:
428:       suffix: tetrahedron
429:       args: -dm_plex_cell tetrahedron \
430:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

432:     test:
433:       suffix: hexahedron
434:       args: -dm_plex_cell hexahedron \
435:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 0.3

437:     test:
438:       suffix: triangular_prism
439:       args: -dm_plex_cell triangular_prism \
440:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

442:     test:
443:       suffix: tensor_triangular_prism
444:       args: -dm_plex_cell tensor_triangular_prism \
445:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

447:     test:
448:       suffix: tensor_quadrilateral_prism
449:       args: -dm_plex_cell tensor_quadrilateral_prism \
450:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

452:     test:
453:       suffix: pyramid
454:       args: -dm_plex_cell pyramid \
455:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

457:   testset:
458:     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -ref_arrangements -dm_plex_check_all \
459:           -ref_dm_view ::ascii_latex -dm_plex_view_tikzscale 1.0

461:     test:
462:       suffix: ref_segment
463:       args: -dm_plex_cell segment \
464:             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0

466:     test:
467:       suffix: ref_triangle
468:       args: -dm_plex_cell triangle \
469:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

471:     test:
472:       suffix: ref_quadrilateral
473:       args: -dm_plex_cell quadrilateral \
474:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

476:     test:
477:       suffix: ref_tensor_segment
478:       args: -dm_plex_cell tensor_quad \
479:             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0

481:     test:
482:       suffix: ref_tetrahedron
483:       args: -dm_plex_cell tetrahedron \
484:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

486:     test:
487:       suffix: ref_hexahedron
488:       args: -dm_plex_cell hexahedron \
489:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

491:     test:
492:       suffix: ref_triangular_prism
493:       args: -dm_plex_cell triangular_prism \
494:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

496:     test:
497:       suffix: ref_tensor_triangular_prism
498:       args: -dm_plex_cell tensor_triangular_prism \
499:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

501:     test:
502:       suffix: ref_tensor_quadrilateral_prism
503:       args: -dm_plex_cell tensor_quadrilateral_prism \
504:             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0

506:   # ToBox should recreate the coordinate space since the cell shape changes
507:   testset:
508:     args: -dm_coord_space 0 -dm_plex_transform_type refine_tobox -ref_arrangements -dm_plex_check_all

510:     test:
511:       suffix: tobox_triangle
512:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
513:             -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

515:     test:
516:       suffix: tobox_tensor_segment
517:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
518:             -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

520:     test:
521:       suffix: tobox_tetrahedron
522:       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
523:             -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

525:     test:
526:       suffix: tobox_triangular_prism
527:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
528:             -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

530:     test:
531:       suffix: tobox_tensor_triangular_prism
532:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
533:             -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

535:     test:
536:       suffix: tobox_tensor_quadrilateral_prism
537:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quadrilateral_prism \
538:             -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

540:   testset:
541:     args: -dm_coord_space 0 -dm_plex_transform_type refine_alfeld -ref_arrangements -dm_plex_check_all

543:     test:
544:       suffix: alfeld_triangle
545:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
546:             -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

548:     test:
549:       suffix: alfeld_tetrahedron
550:       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
551:             -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

553:   testset:
554:     args: -dm_plex_transform_type refine_sbr -ref_arrangements -dm_plex_check_all

556:     # This splits edge 1 of the triangle, and reflects about the added edge
557:     test:
558:       suffix: sbr_triangle_0
559:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -ornts -2,0 \
560:             -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

562:     # This splits edge 0 of the triangle, and reflects about the added edge
563:     test:
564:       suffix: sbr_triangle_1
565:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 1 -ornts -3,0 \
566:             -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

568:     # This splits edge 2 of the triangle, and reflects about the added edge
569:     test:
570:       suffix: sbr_triangle_2
571:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 2 -ornts -1,0 \
572:             -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

574:     # This splits edges 1 and 2 of the triangle
575:     test:
576:       suffix: sbr_triangle_3
577:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -ornts 0 \
578:             -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

580:     # This splits edges 1 and 0 of the triangle
581:     test:
582:       suffix: sbr_triangle_4
583:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -ornts 0 \
584:             -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

586:     # This splits edges 0 and 1 of the triangle
587:     test:
588:       suffix: sbr_triangle_5
589:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 1 -ornts 0 \
590:             -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

592:     # This splits edges 0 and 2 of the triangle
593:     test:
594:       suffix: sbr_triangle_6
595:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 1 -ornts 0 \
596:             -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

598:     # This splits edges 2 and 0 of the triangle
599:     test:
600:       suffix: sbr_triangle_7
601:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 2 -ornts 0 \
602:             -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

604:     # This splits edges 2 and 1 of the triangle
605:     test:
606:       suffix: sbr_triangle_8
607:       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 2 -ornts 0 \
608:             -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

610:   testset:
611:     args: -dm_plex_transform_type refine_boundary_layer -dm_plex_transform_bl_splits 2 -ref_arrangements -dm_plex_check_all

613:     test:
614:       suffix: bl_tensor_segment
615:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
616:             -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

618:     # The subcell check is broken because at orientation 3, the internal triangles do not get properly permuted for the check
619:     test:
620:       suffix: bl_tensor_triangular_prism
621:       requires: TODO
622:       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
623:             -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

625: TEST*/