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 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; 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 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 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 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 */ 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 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 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 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 */ 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 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 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