xref: /petsc/src/dm/impls/plex/tutorials/ex11.c (revision 5f80ce2ab25dff0f4601e710601cbbcecf323266)
1 static char help[] = "Mesh Orientation Tutorial\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscdmplextransform.h>
5 
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;
15 
16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18   PetscInt       n = 2;
19   PetscBool      flg;
20   PetscErrorCode ierr;
21 
22   PetscFunctionBeginUser;
23   options->genArr        = PETSC_FALSE;
24   options->refArr        = PETSC_FALSE;
25   options->printTable    = PETSC_FALSE;
26   options->orntBounds[0] = PETSC_MIN_INT;
27   options->orntBounds[1] = PETSC_MAX_INT;
28   options->numOrnt       = -1;
29   options->initOrnt      = 0;
30 
31   ierr = PetscOptionsBegin(comm, "", "Mesh Orientation Tutorials Options", "DMPLEX");CHKERRQ(ierr);
32   CHKERRQ(PetscOptionsBool("-gen_arrangements", "Flag for generating all arrangements of the cell", "ex11.c", options->genArr, &options->genArr, NULL));
33   CHKERRQ(PetscOptionsBool("-ref_arrangements", "Flag for refining all arrangements of the cell", "ex11.c", options->refArr, &options->refArr, NULL));
34   CHKERRQ(PetscOptionsBool("-print_table", "Print the Cayley table", "ex11.c", options->printTable, &options->printTable, NULL));
35   CHKERRQ(PetscOptionsIntArray("-ornt_bounds", "Bounds for orientation checks", "ex11.c", options->orntBounds, &n, NULL));
36   n    = 48;
37   CHKERRQ(PetscOptionsIntArray("-ornts", "Specific orientations for checks", "ex11.c", options->ornts, &n, &flg));
38   if (flg) {
39     options->numOrnt = n;
40     CHKERRQ(PetscSortInt(n, options->ornts));
41   }
42   CHKERRQ(PetscOptionsInt("-init_ornt", "Initial orientation for starting mesh", "ex11.c", options->initOrnt, &options->initOrnt, NULL));
43   ierr = PetscOptionsEnd();
44   PetscFunctionReturn(0);
45 }
46 
47 static PetscBool ignoreOrnt(AppCtx *user, PetscInt o)
48 {
49   PetscInt       loc;
50   PetscErrorCode ierr;
51 
52   if (user->numOrnt < 0) return PETSC_FALSE;
53   ierr = PetscFindInt(o, user->numOrnt, user->ornts, &loc);
54   if (loc < 0 || ierr)   return PETSC_TRUE;
55   return PETSC_FALSE;
56 }
57 
58 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
59 {
60   PetscFunctionBeginUser;
61   CHKERRQ(DMCreate(comm, dm));
62   CHKERRQ(DMSetType(*dm, DMPLEX));
63   CHKERRQ(DMSetFromOptions(*dm));
64   CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view"));
65   PetscFunctionReturn(0);
66 }
67 
68 static PetscErrorCode CheckCellVertices(DM dm, PetscInt cell, PetscInt o)
69 {
70   DMPolytopeType  ct;
71   const PetscInt *arrVerts;
72   PetscInt       *closure = NULL;
73   PetscInt        Ncl, cl, Nv, vStart, vEnd, v;
74   MPI_Comm        comm;
75 
76   PetscFunctionBeginUser;
77   CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm));
78   CHKERRQ(DMPlexGetCellType(dm, cell, &ct));
79   CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
80   CHKERRQ(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
81   for (cl = 0, Nv = 0; cl < Ncl*2; cl += 2) {
82     const PetscInt vertex = closure[cl];
83 
84     if (vertex < vStart || vertex >= vEnd) continue;
85     closure[Nv++] = vertex;
86   }
87   PetscCheckFalse(Nv != DMPolytopeTypeGetNumVertices(ct),comm, PETSC_ERR_ARG_WRONG, "Cell %D has %D vertices != %D vertices in a %s", cell, Nv, DMPolytopeTypeGetNumVertices(ct), DMPolytopeTypes[ct]);
88   arrVerts = DMPolytopeTypeGetVertexArrangment(ct, o);
89   for (v = 0; v < Nv; ++v) {
90     PetscCheckFalse(closure[v] != arrVerts[v]+vStart,comm, PETSC_ERR_ARG_WRONG, "Cell %D vertex[%D]: %D should be %D for arrangement %D", cell, v, closure[v], arrVerts[v]+vStart, o);
91   }
92   CHKERRQ(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
93   PetscFunctionReturn(0);
94 }
95 
96 /* Transform cell with group operation o */
97 static PetscErrorCode ReorientCell(DM dm, PetscInt cell, PetscInt o, PetscBool swapCoords)
98 {
99   DM              cdm;
100   Vec             coordinates;
101   PetscScalar    *coords, *ccoords = NULL;
102   PetscInt       *closure = NULL;
103   PetscInt        cdim, d, Nc, Ncl, cl, vStart, vEnd, Nv;
104 
105   PetscFunctionBegin;
106   /* Change vertex coordinates so that it plots as we expect */
107   CHKERRQ(DMGetCoordinateDM(dm, &cdm));
108   CHKERRQ(DMGetCoordinateDim(dm, &cdim));
109   CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates));
110   CHKERRQ(DMPlexVecGetClosure(cdm, NULL, coordinates, cell, &Nc, &ccoords));
111   /* Reorient cone */
112   CHKERRQ(DMPlexOrientPoint(dm, cell, o));
113   /* Finish resetting coordinates */
114   if (swapCoords) {
115     CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
116     CHKERRQ(VecGetArrayWrite(coordinates, &coords));
117     CHKERRQ(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
118     for (cl = 0, Nv = 0; cl < Ncl*2; cl += 2) {
119       const PetscInt vertex = closure[cl];
120       PetscScalar   *vcoords;
121 
122       if (vertex < vStart || vertex >= vEnd) continue;
123       CHKERRQ(DMPlexPointLocalRef(cdm, vertex, coords, &vcoords));
124       for (d = 0; d < cdim; ++d) vcoords[d] = ccoords[Nv*cdim + d];
125       ++Nv;
126     }
127     CHKERRQ(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure));
128     CHKERRQ(VecRestoreArrayWrite(coordinates, &coords));
129   }
130   CHKERRQ(DMPlexVecRestoreClosure(cdm, NULL, coordinates, cell, &Nc, &ccoords));
131   PetscFunctionReturn(0);
132 }
133 
134 static PetscErrorCode GenerateArrangments(DM dm, AppCtx *user)
135 {
136   DM             odm;
137   DMPolytopeType ct;
138   PetscInt       No, o;
139   const char    *name;
140 
141   PetscFunctionBeginUser;
142   if (!user->genArr) PetscFunctionReturn(0);
143   CHKERRQ(PetscObjectGetName((PetscObject) dm, &name));
144   CHKERRQ(DMPlexGetCellType(dm, 0, &ct));
145   No   = DMPolytopeTypeGetNumArrangments(ct)/2;
146   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
147     if (ignoreOrnt(user, o)) continue;
148     CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &odm));
149     CHKERRQ(ReorientCell(odm, 0, o, PETSC_TRUE));
150     CHKERRQ(PetscPrintf(PetscObjectComm((PetscObject) dm), "%s orientation %D\n", name, o));
151     CHKERRQ(DMViewFromOptions(odm, NULL, "-gen_dm_view"));
152     CHKERRQ(CheckCellVertices(odm, 0, o));
153     CHKERRQ(DMDestroy(&odm));
154   }
155   PetscFunctionReturn(0);
156 }
157 
158 static PetscErrorCode VerifyCayleyTable(DM dm, AppCtx *user)
159 {
160   DM              dm1, dm2;
161   DMPolytopeType  ct;
162   const PetscInt *refcone, *cone;
163   PetscInt        No, o1, o2, o3, o4;
164   PetscBool       equal;
165   const char     *name;
166 
167   PetscFunctionBeginUser;
168   if (!user->genArr) PetscFunctionReturn(0);
169   CHKERRQ(PetscObjectGetName((PetscObject) dm, &name));
170   CHKERRQ(DMPlexGetCellType(dm, 0, &ct));
171   CHKERRQ(DMPlexGetCone(dm, 0, &refcone));
172   No   = DMPolytopeTypeGetNumArrangments(ct)/2;
173   if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "Cayley Table for %s\n", DMPolytopeTypes[ct]));
174   for (o1 = PetscMax(-No, user->orntBounds[0]); o1 < PetscMin(No, user->orntBounds[1]); ++o1) {
175     for (o2 = PetscMax(-No, user->orntBounds[0]); o2 < PetscMin(No, user->orntBounds[1]); ++o2) {
176       CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &dm1));
177       CHKERRQ(DMPlexOrientPoint(dm1, 0, o2));
178       CHKERRQ(DMPlexCheckFaces(dm1, 0));
179       CHKERRQ(DMPlexOrientPoint(dm1, 0, o1));
180       CHKERRQ(DMPlexCheckFaces(dm1, 0));
181       o3   = DMPolytopeTypeComposeOrientation(ct, o1, o2);
182       /* First verification */
183       CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &dm2));
184       CHKERRQ(DMPlexOrientPoint(dm2, 0, o3));
185       CHKERRQ(DMPlexCheckFaces(dm2, 0));
186       CHKERRQ(DMPlexEqual(dm1, dm2, &equal));
187       if (!equal) {
188         CHKERRQ(DMViewFromOptions(dm1, NULL, "-error_dm_view"));
189         CHKERRQ(DMViewFromOptions(dm2, NULL, "-error_dm_view"));
190         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cayley table error for %s: %D * %D != %D", DMPolytopeTypes[ct], o1, o2, o3);
191       }
192       /* Second verification */
193       CHKERRQ(DMPlexGetCone(dm1, 0, &cone));
194       CHKERRQ(DMPolytopeGetOrientation(ct, refcone, cone, &o4));
195       if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "%D, ", o4));
196       PetscCheckFalse(o3 != o4,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cayley table error for %s: %D * %D = %D != %D", DMPolytopeTypes[ct], o1, o2, o3, o4);
197       CHKERRQ(DMDestroy(&dm1));
198       CHKERRQ(DMDestroy(&dm2));
199     }
200     if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "\n"));
201   }
202   PetscFunctionReturn(0);
203 }
204 
205 static PetscErrorCode VerifyInverse(DM dm, AppCtx *user)
206 {
207   DM              dm1, dm2;
208   DMPolytopeType  ct;
209   const PetscInt *refcone, *cone;
210   PetscInt        No, o, oi, o2;
211   PetscBool       equal;
212   const char     *name;
213 
214   PetscFunctionBeginUser;
215   if (!user->genArr) PetscFunctionReturn(0);
216   CHKERRQ(PetscObjectGetName((PetscObject) dm, &name));
217   CHKERRQ(DMPlexGetCellType(dm, 0, &ct));
218   CHKERRQ(DMPlexGetCone(dm, 0, &refcone));
219   No   = DMPolytopeTypeGetNumArrangments(ct)/2;
220   if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "Inverse table for %s\n", DMPolytopeTypes[ct]));
221   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
222     if (ignoreOrnt(user, o)) continue;
223     oi   = DMPolytopeTypeComposeOrientationInv(ct, 0, o);
224     CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &dm1));
225     CHKERRQ(DMPlexOrientPoint(dm1, 0, o));
226     CHKERRQ(DMPlexCheckFaces(dm1, 0));
227     CHKERRQ(DMPlexOrientPoint(dm1, 0, oi));
228     CHKERRQ(DMPlexCheckFaces(dm1, 0));
229     /* First verification */
230     CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &dm2));
231     CHKERRQ(DMPlexEqual(dm1, dm2, &equal));
232     if (!equal) {
233       CHKERRQ(DMViewFromOptions(dm1, NULL, "-error_dm_view"));
234       CHKERRQ(DMViewFromOptions(dm2, NULL, "-error_dm_view"));
235       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inverse error for %s: %D * %D != 0", DMPolytopeTypes[ct], o, oi);
236     }
237     /* Second verification */
238     CHKERRQ(DMPlexGetCone(dm1, 0, &cone));
239     CHKERRQ(DMPolytopeGetOrientation(ct, refcone, cone, &o2));
240     if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "%D, ", oi));
241     PetscCheckFalse(o2 != 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inverse error for %s: %D * %D = %D != 0", DMPolytopeTypes[ct], o, oi, o2);
242     CHKERRQ(DMDestroy(&dm1));
243     CHKERRQ(DMDestroy(&dm2));
244   }
245   if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "\n"));
246   PetscFunctionReturn(0);
247 }
248 
249 /* Suppose that point p has the same arrangement as o from canonical, compare the subcells to canonical subcells */
250 static PetscErrorCode CheckSubcells(DM dm, DM odm, PetscInt p, PetscInt o, AppCtx *user)
251 {
252   DMPlexTransform tr, otr;
253   DMPolytopeType  ct;
254   DMPolytopeType *rct;
255   const PetscInt *cone, *ornt, *ocone, *oornt;
256   PetscInt       *rsize, *rcone, *rornt;
257   PetscInt        Nct, n, oi, debug = 0;
258 
259   PetscFunctionBeginUser;
260   CHKERRQ(DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr));
261   CHKERRQ(DMPlexTransformSetDM(tr, dm));
262   CHKERRQ(DMPlexTransformSetFromOptions(tr));
263   CHKERRQ(DMPlexTransformSetUp(tr));
264 
265   CHKERRQ(DMPlexTransformCreate(PetscObjectComm((PetscObject) odm), &otr));
266   CHKERRQ(DMPlexTransformSetDM(otr, odm));
267   CHKERRQ(DMPlexTransformSetFromOptions(otr));
268   CHKERRQ(DMPlexTransformSetUp(otr));
269 
270   CHKERRQ(DMPlexGetCellType(dm, p, &ct));
271   CHKERRQ(DMPlexGetCone(dm, p, &cone));
272   CHKERRQ(DMPlexGetConeOrientation(dm, p, &ornt));
273   CHKERRQ(DMPlexGetCone(odm, p, &ocone));
274   CHKERRQ(DMPlexGetConeOrientation(odm, p, &oornt));
275   oi   = DMPolytopeTypeComposeOrientationInv(ct, 0, o);
276   if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "Orientation %D\n", oi));
277 
278   CHKERRQ(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt));
279   for (n = 0; n < Nct; ++n) {
280     DMPolytopeType ctNew = rct[n];
281     PetscInt       r, ro;
282 
283     if (debug) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "  Checking type %s\n", DMPolytopeTypes[ctNew]));
284     for (r = 0; r < rsize[n]; ++r) {
285       const PetscInt *qcone, *qornt, *oqcone, *oqornt;
286       PetscInt        pNew, opNew, oo, pr, fo;
287       PetscBool       restore = PETSC_TRUE;
288 
289       CHKERRQ(DMPlexTransformGetTargetPoint(tr, ct, ctNew, p, r, &pNew));
290       CHKERRQ(DMPlexTransformGetCone(tr, pNew, &qcone, &qornt));
291       if (debug) {
292         PetscInt c;
293 
294         CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "    Checking replica %D (%D)\n      Original Cone", r, pNew));
295         for (c = 0; c < DMPolytopeTypeGetConeSize(ctNew); ++c) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " %D (%D)", qcone[c], qornt[c]));
296         CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "\n"));
297       }
298       for (ro = 0; ro < rsize[n]; ++ro) {
299         PetscBool found;
300 
301         CHKERRQ(DMPlexTransformGetTargetPoint(otr, ct, ctNew, p, ro, &opNew));
302         CHKERRQ(DMPlexTransformGetConeOriented(otr, opNew, o, &oqcone, &oqornt));
303         CHKERRQ(DMPolytopeMatchOrientation(ctNew, oqcone, qcone, &oo, &found));
304         if (found) break;
305         CHKERRQ(DMPlexTransformRestoreCone(otr, pNew, &oqcone, &oqornt));
306       }
307       if (debug) {
308         PetscInt c;
309 
310         CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "    Checking transform replica %D (%D) (%D)\n      Transform Cone", ro, opNew, o));
311         for (c = 0; c < DMPolytopeTypeGetConeSize(ctNew); ++c) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " %D (%D)", oqcone[c], oqornt[c]));
312         CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "\n"));
313         CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "    Matched %D\n", oo));
314       }
315       if (ro == rsize[n]) {
316         /* The tetrahedron has 3 pairs of opposing edges, and any pair can be connected by the interior segment */
317         if (ct == DM_POLYTOPE_TETRAHEDRON) {
318           /* The segment in a tetrahedron does not map into itself under the group action */
319           if (ctNew == DM_POLYTOPE_SEGMENT) {restore = PETSC_FALSE; ro = r; oo = 0;}
320           /* The last four interior faces do not map into themselves under the group action */
321           if (r > 3 && ctNew == DM_POLYTOPE_TRIANGLE) {restore = PETSC_FALSE; ro = r; oo = 0;}
322           /* The last four interior faces do not map into themselves under the group action */
323           if (r > 3 && ctNew == DM_POLYTOPE_TETRAHEDRON) {restore = PETSC_FALSE; ro = r; oo = 0;}
324         }
325         PetscCheckFalse(restore,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to find matching %s %D orientation for cell orientation %D", DMPolytopeTypes[ctNew], r, o);
326       }
327       if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "%D, %D, ", ro, oo));
328       CHKERRQ(DMPlexTransformGetSubcellOrientation(tr, ct, p, oi, ctNew, r, 0, &pr, &fo));
329       if (!user->printTable) {
330         PetscCheckFalse(pr != ro,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Choose wrong replica %D != %D", pr, ro);
331         PetscCheckFalse(fo != oo,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Choose wrong orientation %D != %D", fo, oo);
332       }
333       CHKERRQ(DMPlexTransformRestoreCone(tr, pNew, &qcone, &qornt));
334       if (restore) CHKERRQ(DMPlexTransformRestoreCone(otr, pNew, &oqcone, &oqornt));
335     }
336     if (user->printTable) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "\n"));
337   }
338   CHKERRQ(DMPlexTransformDestroy(&tr));
339   CHKERRQ(DMPlexTransformDestroy(&otr));
340   PetscFunctionReturn(0);
341 }
342 
343 static PetscErrorCode RefineArrangments(DM dm, AppCtx *user)
344 {
345   DM             odm, rdm;
346   DMPolytopeType ct;
347   PetscInt       No, o;
348   const char    *name;
349 
350   PetscFunctionBeginUser;
351   if (!user->refArr) PetscFunctionReturn(0);
352   CHKERRQ(PetscObjectGetName((PetscObject) dm, &name));
353   CHKERRQ(DMPlexGetCellType(dm, 0, &ct));
354   No   = DMPolytopeTypeGetNumArrangments(ct)/2;
355   for (o = PetscMax(-No, user->orntBounds[0]); o < PetscMin(No, user->orntBounds[1]); ++o) {
356     if (ignoreOrnt(user, o)) continue;
357     CHKERRQ(CreateMesh(PetscObjectComm((PetscObject) dm), user, &odm));
358     if (user->initOrnt) CHKERRQ(ReorientCell(odm, 0, user->initOrnt, PETSC_FALSE));
359     CHKERRQ(ReorientCell(odm, 0, o, PETSC_TRUE));
360     CHKERRQ(DMViewFromOptions(odm, NULL, "-orig_dm_view"));
361     CHKERRQ(DMRefine(odm, MPI_COMM_NULL, &rdm));
362     CHKERRQ(DMSetFromOptions(rdm));
363     CHKERRQ(PetscPrintf(PetscObjectComm((PetscObject) dm), "%s orientation %D\n", name, o));
364     CHKERRQ(DMViewFromOptions(rdm, NULL, "-ref_dm_view"));
365     CHKERRQ(CheckSubcells(dm, odm, 0, o, user));
366     CHKERRQ(DMDestroy(&odm));
367     CHKERRQ(DMDestroy(&rdm));
368   }
369   PetscFunctionReturn(0);
370 }
371 
372 int main(int argc, char **argv)
373 {
374   DM             dm;
375   AppCtx         user;
376   PetscErrorCode ierr;
377 
378   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
379   CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user));
380   CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
381   if (user.initOrnt) {
382     CHKERRQ(ReorientCell(dm, 0, user.initOrnt, PETSC_FALSE));
383     CHKERRQ(DMViewFromOptions(dm, NULL, "-ornt_dm_view"));
384   }
385   CHKERRQ(GenerateArrangments(dm, &user));
386   CHKERRQ(VerifyCayleyTable(dm, &user));
387   CHKERRQ(VerifyInverse(dm, &user));
388   CHKERRQ(RefineArrangments(dm, &user));
389   CHKERRQ(DMDestroy(&dm));
390   ierr = PetscFinalize();
391   return ierr;
392 }
393 
394 /*TEST
395 
396   testset:
397     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -gen_arrangements \
398           -gen_dm_view ::ascii_latex -dm_plex_view_tikzscale 0.5
399 
400     test:
401       suffix: segment
402       args: -dm_plex_cell segment \
403             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0
404 
405     test:
406       suffix: triangle
407       args: -dm_plex_cell triangle \
408             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
409 
410     test:
411       suffix: quadrilateral
412       args: -dm_plex_cell quadrilateral \
413             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
414 
415     test:
416       suffix: tensor_segment
417       args: -dm_plex_cell tensor_quad \
418             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
419 
420     test:
421       suffix: tetrahedron
422       args: -dm_plex_cell tetrahedron \
423             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
424 
425     test:
426       suffix: hexahedron
427       args: -dm_plex_cell hexahedron \
428             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0 -dm_plex_view_tikzscale 0.3
429 
430     test:
431       suffix: triangular_prism
432       args: -dm_plex_cell triangular_prism \
433             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
434 
435     test:
436       suffix: tensor_triangular_prism
437       args: -dm_plex_cell tensor_triangular_prism \
438             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
439 
440     test:
441       suffix: tensor_quadrilateral_prism
442       args: -dm_plex_cell tensor_quadrilateral_prism \
443             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
444 
445     test:
446       suffix: pyramid
447       args: -dm_plex_cell pyramid \
448             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
449 
450   testset:
451     args: -dm_coord_space 0 -dm_plex_reference_cell_domain -ref_arrangements -dm_plex_check_all \
452           -ref_dm_view ::ascii_latex -dm_plex_view_tikzscale 1.0
453 
454     test:
455       suffix: ref_segment
456       args: -dm_plex_cell segment \
457             -dm_plex_view_numbers_depth 1,0 -dm_plex_view_colors_depth 1,0
458 
459     test:
460       suffix: ref_triangle
461       args: -dm_plex_cell triangle \
462             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
463 
464     test:
465       suffix: ref_quadrilateral
466       args: -dm_plex_cell quadrilateral \
467             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
468 
469     test:
470       suffix: ref_tensor_segment
471       args: -dm_plex_cell tensor_quad \
472             -dm_plex_view_numbers_depth 1,1,0 -dm_plex_view_colors_depth 1,1,0
473 
474     test:
475       suffix: ref_tetrahedron
476       args: -dm_plex_cell tetrahedron \
477             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
478 
479     test:
480       suffix: ref_hexahedron
481       args: -dm_plex_cell hexahedron \
482             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
483 
484     test:
485       suffix: ref_triangular_prism
486       args: -dm_plex_cell triangular_prism \
487             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
488 
489     test:
490       suffix: ref_tensor_triangular_prism
491       args: -dm_plex_cell tensor_triangular_prism \
492             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
493 
494     test:
495       suffix: ref_tensor_quadrilateral_prism
496       args: -dm_plex_cell tensor_quadrilateral_prism \
497             -dm_plex_view_numbers_depth 1,0,0,0 -dm_plex_view_colors_depth 1,0,0,0
498 
499   # ToBox should recreate the coordinate space since the cell shape changes
500   testset:
501     args: -dm_coord_space 0 -dm_plex_transform_type refine_tobox -ref_arrangements -dm_plex_check_all
502 
503     test:
504       suffix: tobox_triangle
505       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
506             -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
507 
508     test:
509       suffix: tobox_tensor_segment
510       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
511             -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
512 
513     test:
514       suffix: tobox_tetrahedron
515       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
516             -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
517 
518     test:
519       suffix: tobox_triangular_prism
520       args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
521             -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
522 
523     test:
524       suffix: tobox_tensor_triangular_prism
525       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
526             -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
527 
528     test:
529       suffix: tobox_tensor_quadrilateral_prism
530       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quadrilateral_prism \
531             -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
532 
533   testset:
534     args: -dm_coord_space 0 -dm_plex_transform_type refine_alfeld -ref_arrangements -dm_plex_check_all
535 
536     test:
537       suffix: alfeld_triangle
538       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle \
539             -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
540 
541     test:
542       suffix: alfeld_tetrahedron
543       args: -dm_plex_reference_cell_domain -dm_plex_cell tetrahedron \
544             -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
545 
546   testset:
547     args: -dm_plex_transform_type refine_sbr -ref_arrangements -dm_plex_check_all
548 
549     # This splits edge 1 of the triangle, and reflects about the added edge
550     test:
551       suffix: sbr_triangle_0
552       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -ornts -2,0 \
553             -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
554 
555     # This splits edge 0 of the triangle, and reflects about the added edge
556     test:
557       suffix: sbr_triangle_1
558       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 1 -ornts -3,0 \
559             -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
560 
561     # This splits edge 2 of the triangle, and reflects about the added edge
562     test:
563       suffix: sbr_triangle_2
564       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5 -init_ornt 2 -ornts -1,0 \
565             -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 
567     # This splits edges 1 and 2 of the triangle
568     test:
569       suffix: sbr_triangle_3
570       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -ornts 0 \
571             -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 
573     # This splits edges 1 and 0 of the triangle
574     test:
575       suffix: sbr_triangle_4
576       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -ornts 0 \
577             -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 
579     # This splits edges 0 and 1 of the triangle
580     test:
581       suffix: sbr_triangle_5
582       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 1 -ornts 0 \
583             -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 
585     # This splits edges 0 and 2 of the triangle
586     test:
587       suffix: sbr_triangle_6
588       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 1 -ornts 0 \
589             -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 
591     # This splits edges 2 and 0 of the triangle
592     test:
593       suffix: sbr_triangle_7
594       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 5,6 -init_ornt 2 -ornts 0 \
595             -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 
597     # This splits edges 2 and 1 of the triangle
598     test:
599       suffix: sbr_triangle_8
600       args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -dm_plex_transform_sbr_ref_cell 4,5 -init_ornt 2 -ornts 0 \
601             -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 
603   testset:
604     args: -dm_plex_transform_type refine_boundary_layer -dm_plex_transform_bl_splits 2 -ref_arrangements -dm_plex_check_all
605 
606     test:
607       suffix: bl_tensor_segment
608       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_quad \
609             -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 
611     # The subcell check is broken because at orientation 3, the internal triangles do not get properly permuted for the check
612     test:
613       suffix: bl_tensor_triangular_prism
614       requires: TODO
615       args: -dm_plex_reference_cell_domain -dm_plex_cell tensor_triangular_prism \
616             -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
617 
618 TEST*/
619