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