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_MIN_INT;
 26:   options->orntBounds[1] = PETSC_MAX_INT;
 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*/