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
ProcessOptions(MPI_Comm comm,AppCtx * options)16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18 PetscInt n = 2;
19 PetscBool flg;
20
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;
29
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 }
45
ignoreOrnt(AppCtx * user,PetscInt o)46 static PetscBool ignoreOrnt(AppCtx *user, PetscInt o)
47 {
48 PetscInt loc;
49 PetscErrorCode ierr;
50
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 }
56
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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 }
66
CheckCellVertices(DM dm,PetscInt cell,PetscInt o)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;
74
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];
82
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 }
94
95 /* Transform cell with group operation o */
ReorientCell(DM dm,PetscInt cell,PetscInt o,PetscBool swapCoords)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;
103
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;
120
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 }
132
GenerateArrangements(DM dm,AppCtx * user)133 static PetscErrorCode GenerateArrangements(DM dm, AppCtx *user)
134 {
135 DM odm;
136 DMPolytopeType ct;
137 PetscInt No, o;
138 const char *name;
139
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 }
156
VerifyCayleyTable(DM dm,AppCtx * user)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;
165
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 }
203
VerifyInverse(DM dm,AppCtx * user)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;
212
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 }
247
248 /* Suppose that point p has the same arrangement as o from canonical, compare the subcells to canonical subcells */
CheckSubcells(DM dm,DM odm,PetscInt p,PetscInt o,AppCtx * user)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;
257
258 PetscFunctionBeginUser;
259 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
260 PetscCall(DMPlexTransformSetDM(tr, dm));
261 PetscCall(DMPlexTransformSetFromOptions(tr));
262 PetscCall(DMPlexTransformSetUp(tr));
263
264 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)odm), &otr));
265 PetscCall(DMPlexTransformSetDM(otr, odm));
266 PetscCall(DMPlexTransformSetFromOptions(otr));
267 PetscCall(DMPlexTransformSetUp(otr));
268
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));
276
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;
281
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;
287
288 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, ctNew, p, r, &pNew));
289 PetscCall(DMPlexTransformGetCone(tr, pNew, &qcone, &qornt));
290 if (debug) {
291 PetscInt c;
292
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;
299
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;
308
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 }
353
RefineArrangements(DM dm,AppCtx * user)354 static PetscErrorCode RefineArrangements(DM dm, AppCtx *user)
355 {
356 DM odm, rdm;
357 DMPolytopeType ct;
358 PetscInt No, o;
359 const char *name;
360
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 }
382
main(int argc,char ** argv)383 int main(int argc, char **argv)
384 {
385 DM dm;
386 AppCtx user;
387
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 }
404
405 /*TEST
406
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
410
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
415
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
420
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
425
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
430
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
435
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
440
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
445
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
450
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
455
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
460
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
464
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
469
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
474
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
479
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
484
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
489
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
494
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
499
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
504
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
509
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
513
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
518
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
523
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
528
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
533
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
538
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
543
544 testset:
545 args: -dm_coord_space 0 -dm_plex_transform_type refine_alfeld -ref_arrangements -dm_plex_check_all
546
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
551
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
556
557 testset:
558 args: -dm_plex_transform_type refine_sbr -ref_arrangements -dm_plex_check_all
559
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
565
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
571
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
577
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
583
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
589
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
595
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
601
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
607
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
613
614 testset:
615 args: -dm_plex_transform_type refine_boundary_layer -dm_plex_transform_bl_splits 2 -ref_arrangements -dm_plex_check_all
616
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
621
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
628
629 TEST*/
630