1 static char help[] = "Tests for parallel mesh loading and parallel topological interpolation\n\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 /* List of test meshes 5 6 Network 7 ------- 8 Test 0 (2 ranks): 9 10 network=0: 11 --------- 12 cell 0 cell 1 cell 2 nCells-1 (edge) 13 0 ------ 1 ------ 2 ------ 3 -- -- v -- -- nCells (vertex) 14 15 vertex distribution: 16 rank 0: 0 1 17 rank 1: 2 3 ... nCells 18 cell(edge) distribution: 19 rank 0: 0 1 20 rank 1: 2 ... nCells-1 21 22 network=1: 23 --------- 24 v2 25 ^ 26 | 27 cell 2 28 | 29 v0 --cell 0--> v3--cell 1--> v1 30 31 vertex distribution: 32 rank 0: 0 1 3 33 rank 1: 2 34 cell(edge) distribution: 35 rank 0: 0 1 36 rank 1: 2 37 38 example: 39 mpiexec -n 2 ./ex18 -distribute 1 -dim 1 -orig_dm_view -dist_dm_view -dist_dm_view -petscpartitioner_type parmetis -ncells 50 40 41 Triangle 42 -------- 43 Test 0 (2 ranks): 44 Two triangles sharing a face 45 46 2 47 / | \ 48 / | \ 49 / | \ 50 0 0 | 1 3 51 \ | / 52 \ | / 53 \ | / 54 1 55 56 vertex distribution: 57 rank 0: 0 1 58 rank 1: 2 3 59 cell distribution: 60 rank 0: 0 61 rank 1: 1 62 63 Test 1 (3 ranks): 64 Four triangles partitioned across 3 ranks 65 66 0 _______ 3 67 | \ / | 68 | \ 1 / | 69 | \ / | 70 | 0 2 2 | 71 | / \ | 72 | / 3 \ | 73 | / \ | 74 1 ------- 4 75 76 vertex distribution: 77 rank 0: 0 1 78 rank 1: 2 3 79 rank 2: 4 80 cell distribution: 81 rank 0: 0 82 rank 1: 1 83 rank 2: 2 3 84 85 Test 2 (3 ranks): 86 Four triangles partitioned across 3 ranks 87 88 1 _______ 3 89 | \ / | 90 | \ 1 / | 91 | \ / | 92 | 0 0 2 | 93 | / \ | 94 | / 3 \ | 95 | / \ | 96 2 ------- 4 97 98 vertex distribution: 99 rank 0: 0 1 100 rank 1: 2 3 101 rank 2: 4 102 cell distribution: 103 rank 0: 0 104 rank 1: 1 105 rank 2: 2 3 106 107 Tetrahedron 108 ----------- 109 Test 0: 110 Two tets sharing a face 111 112 cell 3 _______ cell 113 0 / | \ \ 1 114 / | \ \ 115 / | \ \ 116 0----|----4-----2 117 \ | / / 118 \ | / / 119 \ | / / 120 1------- 121 y 122 | x 123 |/ 124 *----z 125 126 vertex distribution: 127 rank 0: 0 1 128 rank 1: 2 3 4 129 cell distribution: 130 rank 0: 0 131 rank 1: 1 132 133 Quadrilateral 134 ------------- 135 Test 0 (2 ranks): 136 Two quads sharing a face 137 138 3-------2-------5 139 | | | 140 | 0 | 1 | 141 | | | 142 0-------1-------4 143 144 vertex distribution: 145 rank 0: 0 1 2 146 rank 1: 3 4 5 147 cell distribution: 148 rank 0: 0 149 rank 1: 1 150 151 TODO Test 1: 152 A quad and a triangle sharing a face 153 154 5-------4 155 | | \ 156 | 0 | \ 157 | | 1 \ 158 2-------3----6 159 160 Hexahedron 161 ---------- 162 Test 0 (2 ranks): 163 Two hexes sharing a face 164 165 cell 7-------------6-------------11 cell 166 0 /| /| /| 1 167 / | F1 / | F7 / | 168 / | / | / | 169 4-------------5-------------10 | 170 | | F4 | | F10 | | 171 | | | | | | 172 |F5 | |F3 | |F9 | 173 | | F2 | | F8 | | 174 | 3---------|---2---------|---9 175 | / | / | / 176 | / F0 | / F6 | / 177 |/ |/ |/ 178 0-------------1-------------8 179 180 vertex distribution: 181 rank 0: 0 1 2 3 4 5 182 rank 1: 6 7 8 9 10 11 183 cell distribution: 184 rank 0: 0 185 rank 1: 1 186 187 */ 188 189 typedef enum {NONE, CREATE, AFTER_CREATE, AFTER_DISTRIBUTE} InterpType; 190 191 typedef struct { 192 PetscInt debug; /* The debugging level */ 193 PetscInt testNum; /* Indicates the mesh to create */ 194 PetscInt dim; /* The topological mesh dimension */ 195 PetscBool cellSimplex; /* Use simplices or hexes */ 196 PetscBool distribute; /* Distribute the mesh */ 197 InterpType interpolate; /* Interpolate the mesh before or after DMPlexDistribute() */ 198 PetscBool useGenerator; /* Construct mesh with a mesh generator */ 199 PetscBool testOrientIF; /* Test for different original interface orientations */ 200 PetscBool testHeavy; /* Run the heavy PointSF test */ 201 PetscBool customView; /* Show results of DMPlexIsInterpolated() etc. */ 202 PetscInt ornt[2]; /* Orientation of interface on rank 0 and rank 1 */ 203 PetscInt faces[3]; /* Number of faces per dimension for generator */ 204 PetscScalar coords[128]; 205 PetscReal coordsTol; 206 PetscInt ncoords; 207 PetscInt pointsToExpand[128]; 208 PetscInt nPointsToExpand; 209 PetscBool testExpandPointsEmpty; 210 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 211 } AppCtx; 212 213 struct _n_PortableBoundary { 214 Vec coordinates; 215 PetscInt depth; 216 PetscSection *sections; 217 }; 218 typedef struct _n_PortableBoundary * PortableBoundary; 219 220 #if defined(PETSC_USE_LOG) 221 static PetscLogStage stage[3]; 222 #endif 223 224 static PetscErrorCode DMPlexCheckPointSFHeavy(DM, PortableBoundary); 225 static PetscErrorCode DMPlexSetOrientInterface_Private(DM,PetscBool); 226 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM, PortableBoundary *); 227 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM, IS, PetscSection, IS *); 228 229 static PetscErrorCode PortableBoundaryDestroy(PortableBoundary *bnd) 230 { 231 PetscInt d; 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 if (!*bnd) PetscFunctionReturn(0); 236 ierr = VecDestroy(&(*bnd)->coordinates);CHKERRQ(ierr); 237 for (d=0; d < (*bnd)->depth; d++) { 238 ierr = PetscSectionDestroy(&(*bnd)->sections[d]);CHKERRQ(ierr); 239 } 240 ierr = PetscFree((*bnd)->sections);CHKERRQ(ierr); 241 ierr = PetscFree(*bnd);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 246 { 247 const char *interpTypes[4] = {"none", "create", "after_create", "after_distribute"}; 248 PetscInt interp=NONE, dim; 249 PetscBool flg1, flg2; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 options->debug = 0; 254 options->testNum = 0; 255 options->dim = 2; 256 options->cellSimplex = PETSC_TRUE; 257 options->distribute = PETSC_FALSE; 258 options->interpolate = NONE; 259 options->useGenerator = PETSC_FALSE; 260 options->testOrientIF = PETSC_FALSE; 261 options->testHeavy = PETSC_TRUE; 262 options->customView = PETSC_FALSE; 263 options->testExpandPointsEmpty = PETSC_FALSE; 264 options->ornt[0] = 0; 265 options->ornt[1] = 0; 266 options->faces[0] = 2; 267 options->faces[1] = 2; 268 options->faces[2] = 2; 269 options->filename[0] = '\0'; 270 options->coordsTol = PETSC_DEFAULT; 271 272 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 273 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 274 ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 275 ierr = PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 276 ierr = PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL);CHKERRQ(ierr); 277 ierr = PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL);CHKERRQ(ierr); 278 options->interpolate = (InterpType) interp; 279 PetscCheckFalse(!options->distribute && options->interpolate == AFTER_DISTRIBUTE,comm, PETSC_ERR_SUP, "-interpolate after_distribute needs -distribute 1"); 280 ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr); 281 options->ncoords = 128; 282 ierr = PetscOptionsScalarArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL);CHKERRQ(ierr); 283 ierr = PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL);CHKERRQ(ierr); 284 options->nPointsToExpand = 128; 285 ierr = PetscOptionsIntArray("-test_expand_points", "Expand given array of DAG point using DMPlexGetConeRecursive() and print results", "ex18.c", options->pointsToExpand, &options->nPointsToExpand, NULL);CHKERRQ(ierr); 286 if (options->nPointsToExpand) { 287 ierr = PetscOptionsBool("-test_expand_points_empty", "For -test_expand_points, rank 0 will have empty input array", "ex18.c", options->testExpandPointsEmpty, &options->testExpandPointsEmpty, NULL);CHKERRQ(ierr); 288 } 289 ierr = PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL);CHKERRQ(ierr); 290 ierr = PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL);CHKERRQ(ierr); 291 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3);CHKERRQ(ierr); 292 dim = 3; 293 ierr = PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2);CHKERRQ(ierr); 294 if (flg2) { 295 PetscCheckFalse(flg1 && dim != options->dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "specified -dim %D is not equal to length %D of -faces (note that -dim can be omitted)", options->dim, dim); 296 options->dim = dim; 297 } 298 ierr = PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 299 ierr = PetscOptionsBoundedInt("-rotate_interface_0", "Rotation (relative orientation) of interface on rank 0; implies -interpolate create -distribute 0", "ex18.c", options->ornt[0], &options->ornt[0], &options->testOrientIF,0);CHKERRQ(ierr); 300 ierr = PetscOptionsBoundedInt("-rotate_interface_1", "Rotation (relative orientation) of interface on rank 1; implies -interpolate create -distribute 0", "ex18.c", options->ornt[1], &options->ornt[1], &flg2,0);CHKERRQ(ierr); 301 PetscCheckFalse(flg2 != options->testOrientIF,comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set"); 302 if (options->testOrientIF) { 303 PetscInt i; 304 for (i=0; i<2; i++) { 305 if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10); /* 11 12 13 become -1 -2 -3 */ 306 } 307 options->filename[0] = 0; 308 options->useGenerator = PETSC_FALSE; 309 options->dim = 3; 310 options->cellSimplex = PETSC_TRUE; 311 options->interpolate = CREATE; 312 options->distribute = PETSC_FALSE; 313 } 314 ierr = PetscOptionsEnd();CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 319 { 320 PetscInt testNum = user->testNum; 321 PetscMPIInt rank,size; 322 PetscErrorCode ierr; 323 PetscInt numCorners=2,i; 324 PetscInt numCells,numVertices,network; 325 PetscInt *cells; 326 PetscReal *coords; 327 328 PetscFunctionBegin; 329 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 330 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 331 PetscCheckFalse(size > 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum); 332 333 numCells = 3; 334 ierr = PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL);CHKERRQ(ierr); 335 PetscCheckFalse(numCells < 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells); 336 337 if (size == 1) { 338 numVertices = numCells + 1; 339 ierr = PetscMalloc2(2*numCells,&cells,2*numVertices,&coords);CHKERRQ(ierr); 340 for (i=0; i<numCells; i++) { 341 cells[2*i] = i; cells[2*i+1] = i + 1; 342 coords[2*i] = i; coords[2*i+1] = i + 1; 343 } 344 345 ierr = DMPlexCreateFromCellListPetsc(comm, user->dim, numCells, numVertices, numCorners, PETSC_FALSE, cells, user->dim, coords, dm);CHKERRQ(ierr); 346 ierr = PetscFree2(cells,coords);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 350 network = 0; 351 ierr = PetscOptionsGetInt(NULL, NULL, "-network_case", &network, NULL);CHKERRQ(ierr); 352 if (network == 0) { 353 switch (rank) { 354 case 0: 355 { 356 numCells = 2; 357 numVertices = numCells; 358 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 359 cells[0] = 0; cells[1] = 1; 360 cells[2] = 1; cells[3] = 2; 361 coords[0] = 0.; coords[1] = 1.; 362 coords[2] = 1.; coords[3] = 2.; 363 } 364 break; 365 case 1: 366 { 367 numCells -= 2; 368 numVertices = numCells + 1; 369 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 370 for (i=0; i<numCells; i++) { 371 cells[2*i] = 2+i; cells[2*i+1] = 2 + i + 1; 372 coords[2*i] = 2+i; coords[2*i+1] = 2 + i + 1; 373 } 374 } 375 break; 376 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 377 } 378 } else { /* network_case = 1 */ 379 /* ----------------------- */ 380 switch (rank) { 381 case 0: 382 { 383 numCells = 2; 384 numVertices = 3; 385 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 386 cells[0] = 0; cells[1] = 3; 387 cells[2] = 3; cells[3] = 1; 388 } 389 break; 390 case 1: 391 { 392 numCells = 1; 393 numVertices = 1; 394 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 395 cells[0] = 3; cells[1] = 2; 396 } 397 break; 398 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 399 } 400 } 401 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, PETSC_FALSE, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 402 ierr = PetscFree2(cells,coords);CHKERRQ(ierr); 403 PetscFunctionReturn(0); 404 } 405 406 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 407 { 408 PetscInt testNum = user->testNum, p; 409 PetscMPIInt rank, size; 410 PetscErrorCode ierr; 411 412 PetscFunctionBegin; 413 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 414 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 415 switch (testNum) { 416 case 0: 417 PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 418 switch (rank) { 419 case 0: 420 { 421 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 422 const PetscInt cells[3] = {0, 1, 2}; 423 PetscReal coords[4] = {-0.5, 0.5, 0.0, 0.0}; 424 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 425 426 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 427 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 428 } 429 break; 430 case 1: 431 { 432 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 433 const PetscInt cells[3] = {1, 3, 2}; 434 PetscReal coords[4] = {0.0, 1.0, 0.5, 0.5}; 435 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 436 437 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 438 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 439 } 440 break; 441 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 442 } 443 break; 444 case 1: 445 PetscCheckFalse(size != 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); 446 switch (rank) { 447 case 0: 448 { 449 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 450 const PetscInt cells[3] = {0, 1, 2}; 451 PetscReal coords[4] = {0.0, 1.0, 0.0, 0.0}; 452 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 453 454 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 455 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 456 } 457 break; 458 case 1: 459 { 460 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 461 const PetscInt cells[3] = {0, 2, 3}; 462 PetscReal coords[4] = {0.5, 0.5, 1.0, 1.0}; 463 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 464 465 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 466 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 467 } 468 break; 469 case 2: 470 { 471 const PetscInt numCells = 2, numVertices = 1, numCorners = 3; 472 const PetscInt cells[6] = {2, 4, 3, 2, 1, 4}; 473 PetscReal coords[2] = {1.0, 0.0}; 474 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 475 476 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 477 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 478 } 479 break; 480 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 481 } 482 break; 483 case 2: 484 PetscCheckFalse(size != 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); 485 switch (rank) { 486 case 0: 487 { 488 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 489 const PetscInt cells[3] = {1, 2, 0}; 490 PetscReal coords[4] = {0.5, 0.5, 0.0, 1.0}; 491 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 492 493 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 494 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 495 } 496 break; 497 case 1: 498 { 499 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 500 const PetscInt cells[3] = {1, 0, 3}; 501 PetscReal coords[4] = {0.0, 0.0, 1.0, 1.0}; 502 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 503 504 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 505 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 506 } 507 break; 508 case 2: 509 { 510 const PetscInt numCells = 2, numVertices = 1, numCorners = 3; 511 const PetscInt cells[6] = {0, 4, 3, 0, 2, 4}; 512 PetscReal coords[2] = {1.0, 0.0}; 513 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 514 515 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 516 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 517 } 518 break; 519 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 520 } 521 break; 522 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 523 } 524 PetscFunctionReturn(0); 525 } 526 527 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 528 { 529 PetscInt testNum = user->testNum, p; 530 PetscMPIInt rank, size; 531 PetscErrorCode ierr; 532 533 PetscFunctionBegin; 534 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 535 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 536 switch (testNum) { 537 case 0: 538 PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 539 switch (rank) { 540 case 0: 541 { 542 const PetscInt numCells = 1, numVertices = 2, numCorners = 4; 543 const PetscInt cells[4] = {0, 2, 1, 3}; 544 PetscReal coords[6] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0}; 545 PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; 546 547 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 548 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 549 } 550 break; 551 case 1: 552 { 553 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 554 const PetscInt cells[4] = {1, 2, 4, 3}; 555 PetscReal coords[9] = {1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; 556 PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; 557 558 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 559 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 560 } 561 break; 562 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 563 } 564 break; 565 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 566 } 567 if (user->testOrientIF) { 568 PetscInt ifp[] = {8, 6}; 569 570 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh before orientation");CHKERRQ(ierr); 571 ierr = DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view");CHKERRQ(ierr); 572 /* rotate interface face ifp[rank] by given orientation ornt[rank] */ 573 ierr = DMPlexOrientPoint(*dm, ifp[rank], user->ornt[rank]);CHKERRQ(ierr); 574 ierr = DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view");CHKERRQ(ierr); 575 ierr = DMPlexCheckFaces(*dm, 0);CHKERRQ(ierr); 576 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 577 ierr = PetscPrintf(comm, "Orientation test PASSED\n");CHKERRQ(ierr); 578 } 579 PetscFunctionReturn(0); 580 } 581 582 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 583 { 584 PetscInt testNum = user->testNum, p; 585 PetscMPIInt rank, size; 586 PetscErrorCode ierr; 587 588 PetscFunctionBegin; 589 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 590 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 591 switch (testNum) { 592 case 0: 593 PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 594 switch (rank) { 595 case 0: 596 { 597 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 598 const PetscInt cells[4] = {0, 1, 2, 3}; 599 PetscReal coords[6] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0}; 600 PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; 601 602 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 603 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 604 } 605 break; 606 case 1: 607 { 608 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 609 const PetscInt cells[4] = {1, 4, 5, 2}; 610 PetscReal coords[6] = {-0.5, 1.0, 0.5, 0.0, 0.5, 1.0}; 611 PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; 612 613 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 614 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 615 } 616 break; 617 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 618 } 619 break; 620 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 621 } 622 PetscFunctionReturn(0); 623 } 624 625 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 626 { 627 PetscInt testNum = user->testNum, p; 628 PetscMPIInt rank, size; 629 PetscErrorCode ierr; 630 631 PetscFunctionBegin; 632 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 633 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 634 switch (testNum) { 635 case 0: 636 PetscCheckFalse(size != 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 637 switch (rank) { 638 case 0: 639 { 640 const PetscInt numCells = 1, numVertices = 6, numCorners = 8; 641 const PetscInt cells[8] = {0, 3, 2, 1, 4, 5, 6, 7}; 642 PetscReal coords[6*3] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, -0.5,0.0,1.0, 0.0,0.0,1.0}; 643 PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 644 645 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 646 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 647 } 648 break; 649 case 1: 650 { 651 const PetscInt numCells = 1, numVertices = 6, numCorners = 8; 652 const PetscInt cells[8] = {1, 2, 9, 8, 5, 10, 11, 6}; 653 PetscReal coords[6*3] = {0.0,1.0,1.0, -0.5,1.0,1.0, 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 654 PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 655 656 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm);CHKERRQ(ierr); 657 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 658 } 659 break; 660 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 661 } 662 break; 663 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 664 } 665 PetscFunctionReturn(0); 666 } 667 668 static PetscErrorCode CustomView(DM dm, PetscViewer v) 669 { 670 DMPlexInterpolatedFlag interpolated; 671 PetscBool distributed; 672 PetscErrorCode ierr; 673 674 PetscFunctionBegin; 675 ierr = DMPlexIsDistributed(dm, &distributed);CHKERRQ(ierr); 676 ierr = DMPlexIsInterpolatedCollective(dm, &interpolated);CHKERRQ(ierr); 677 ierr = PetscViewerASCIIPrintf(v, "DMPlexIsDistributed: %s\n", PetscBools[distributed]);CHKERRQ(ierr); 678 ierr = PetscViewerASCIIPrintf(v, "DMPlexIsInterpolatedCollective: %s\n", DMPlexInterpolatedFlags[interpolated]);CHKERRQ(ierr); 679 PetscFunctionReturn(0); 680 } 681 682 static PetscErrorCode CreateMeshFromFile(MPI_Comm comm, AppCtx *user, DM *dm, DM *serialDM) 683 { 684 const char *filename = user->filename; 685 PetscBool testHeavy = user->testHeavy; 686 PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; 687 PetscBool distributed = PETSC_FALSE; 688 PetscErrorCode ierr; 689 690 PetscFunctionBegin; 691 *serialDM = NULL; 692 if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_FALSE);CHKERRQ(ierr);} 693 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 694 ierr = DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, dm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 695 ierr = PetscLogStagePop();CHKERRQ(ierr); 696 if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_TRUE);CHKERRQ(ierr);} 697 ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr); 698 ierr = PetscPrintf(comm, "DMPlexCreateFromFile produced %s mesh.\n", distributed ? "distributed" : "serial");CHKERRQ(ierr); 699 if (testHeavy && distributed) { 700 ierr = PetscOptionsSetValue(NULL, "-dm_plex_hdf5_force_sequential", NULL);CHKERRQ(ierr); 701 ierr = DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, serialDM);CHKERRQ(ierr); 702 ierr = DMPlexIsDistributed(*serialDM, &distributed);CHKERRQ(ierr); 703 PetscCheckFalse(distributed,comm, PETSC_ERR_PLIB, "unable to create a serial DM from file"); 704 } 705 ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 710 { 711 PetscPartitioner part; 712 PortableBoundary boundary = NULL; 713 DM serialDM = NULL; 714 PetscBool cellSimplex = user->cellSimplex; 715 PetscBool useGenerator = user->useGenerator; 716 PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; 717 PetscBool interpSerial = user->interpolate == AFTER_CREATE ? PETSC_TRUE : PETSC_FALSE; 718 PetscBool interpParallel = user->interpolate == AFTER_DISTRIBUTE ? PETSC_TRUE : PETSC_FALSE; 719 PetscBool testHeavy = user->testHeavy; 720 PetscMPIInt rank; 721 PetscErrorCode ierr; 722 723 PetscFunctionBegin; 724 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 725 if (user->filename[0]) { 726 ierr = CreateMeshFromFile(comm, user, dm, &serialDM);CHKERRQ(ierr); 727 } else if (useGenerator) { 728 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 729 ierr = DMPlexCreateBoxMesh(comm, user->dim, cellSimplex, user->faces, NULL, NULL, NULL, interpCreate, dm);CHKERRQ(ierr); 730 ierr = PetscLogStagePop();CHKERRQ(ierr); 731 } else { 732 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 733 switch (user->dim) { 734 case 1: 735 ierr = CreateMesh_1D(comm, interpCreate, user, dm);CHKERRQ(ierr); 736 break; 737 case 2: 738 if (cellSimplex) { 739 ierr = CreateSimplex_2D(comm, interpCreate, user, dm);CHKERRQ(ierr); 740 } else { 741 ierr = CreateQuad_2D(comm, interpCreate, user, dm);CHKERRQ(ierr); 742 } 743 break; 744 case 3: 745 if (cellSimplex) { 746 ierr = CreateSimplex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr); 747 } else { 748 ierr = CreateHex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr); 749 } 750 break; 751 default: 752 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", user->dim); 753 } 754 ierr = PetscLogStagePop();CHKERRQ(ierr); 755 } 756 PetscCheckFalse(user->ncoords % user->dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "length of coordinates array %D must be divisable by spatial dimension %D", user->ncoords, user->dim); 757 ierr = PetscObjectSetName((PetscObject) *dm, "Original Mesh");CHKERRQ(ierr); 758 ierr = DMViewFromOptions(*dm, NULL, "-orig_dm_view");CHKERRQ(ierr); 759 760 if (interpSerial) { 761 DM idm; 762 763 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);} 764 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 765 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 766 ierr = PetscLogStagePop();CHKERRQ(ierr); 767 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);} 768 ierr = DMDestroy(dm);CHKERRQ(ierr); 769 *dm = idm; 770 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh");CHKERRQ(ierr); 771 ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr); 772 } 773 774 /* Set partitioner options */ 775 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 776 if (part) { 777 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE);CHKERRQ(ierr); 778 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 779 } 780 781 if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 782 if (testHeavy) { 783 PetscBool distributed; 784 785 ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr); 786 if (!serialDM && !distributed) { 787 serialDM = *dm; 788 ierr = PetscObjectReference((PetscObject)*dm);CHKERRQ(ierr); 789 } 790 if (serialDM) { 791 ierr = DMPlexGetExpandedBoundary_Private(serialDM, &boundary);CHKERRQ(ierr); 792 } 793 if (boundary) { 794 /* check DM which has been created in parallel and already interpolated */ 795 ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr); 796 } 797 /* Orient interface because it could be deliberately skipped above. It is idempotent. */ 798 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 799 } 800 if (user->distribute) { 801 DM pdm = NULL; 802 803 /* Redistribute mesh over processes using that partitioner */ 804 ierr = PetscLogStagePush(stage[1]);CHKERRQ(ierr); 805 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 806 ierr = PetscLogStagePop();CHKERRQ(ierr); 807 if (pdm) { 808 ierr = DMDestroy(dm);CHKERRQ(ierr); 809 *dm = pdm; 810 ierr = PetscObjectSetName((PetscObject) *dm, "Redistributed Mesh");CHKERRQ(ierr); 811 ierr = DMViewFromOptions(*dm, NULL, "-dist_dm_view");CHKERRQ(ierr); 812 } 813 814 if (interpParallel) { 815 DM idm; 816 817 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);} 818 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 819 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 820 ierr = PetscLogStagePop();CHKERRQ(ierr); 821 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);} 822 ierr = DMDestroy(dm);CHKERRQ(ierr); 823 *dm = idm; 824 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Redistributed Mesh");CHKERRQ(ierr); 825 ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr); 826 } 827 } 828 if (testHeavy) { 829 if (boundary) { 830 ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr); 831 } 832 /* Orient interface because it could be deliberately skipped above. It is idempotent. */ 833 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 834 } 835 836 ierr = PetscObjectSetName((PetscObject) *dm, "Parallel Mesh");CHKERRQ(ierr); 837 ierr = DMPlexDistributeSetDefault(*dm, PETSC_FALSE);CHKERRQ(ierr); 838 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 839 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 840 841 if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 842 ierr = DMDestroy(&serialDM);CHKERRQ(ierr); 843 ierr = PortableBoundaryDestroy(&boundary);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 #define ps2d(number) ((double) PetscRealPart(number)) 848 static inline PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, const PetscScalar coords[], PetscReal tol) 849 { 850 PetscErrorCode ierr; 851 852 PetscFunctionBegin; 853 PetscCheckFalse(dim > 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3"); 854 if (tol >= 1e-3) { 855 switch (dim) { 856 case 1: ierr = PetscSNPrintf(buf,len,"(%12.3f)",ps2d(coords[0]));CHKERRQ(ierr); 857 case 2: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]));CHKERRQ(ierr); 858 case 3: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]));CHKERRQ(ierr); 859 } 860 } else { 861 switch (dim) { 862 case 1: ierr = PetscSNPrintf(buf,len,"(%12.6f)",ps2d(coords[0]));CHKERRQ(ierr); 863 case 2: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]));CHKERRQ(ierr); 864 case 3: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]));CHKERRQ(ierr); 865 } 866 } 867 PetscFunctionReturn(0); 868 } 869 870 static PetscErrorCode ViewVerticesFromCoords(DM dm, Vec coordsVec, PetscReal tol, PetscViewer viewer) 871 { 872 PetscInt dim, i, npoints; 873 IS pointsIS; 874 const PetscInt *points; 875 const PetscScalar *coords; 876 char coordstr[128]; 877 MPI_Comm comm; 878 PetscMPIInt rank; 879 PetscErrorCode ierr; 880 881 PetscFunctionBegin; 882 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 883 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 884 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 885 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 886 ierr = DMPlexFindVertices(dm, coordsVec, tol, &pointsIS);CHKERRQ(ierr); 887 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 888 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 889 ierr = VecGetArrayRead(coordsVec, &coords);CHKERRQ(ierr); 890 for (i=0; i < npoints; i++) { 891 ierr = coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol);CHKERRQ(ierr); 892 if (rank == 0 && i) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, "-----\n");CHKERRQ(ierr);} 893 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i]);CHKERRQ(ierr); 894 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 895 } 896 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 897 ierr = VecRestoreArrayRead(coordsVec, &coords);CHKERRQ(ierr); 898 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 899 ierr = ISDestroy(&pointsIS);CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user) 904 { 905 IS is; 906 PetscSection *sects; 907 IS *iss; 908 PetscInt d,depth; 909 PetscMPIInt rank; 910 PetscErrorCode ierr; 911 PetscViewer viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer; 912 913 PetscFunctionBegin; 914 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 915 if (user->testExpandPointsEmpty && rank == 0) { 916 ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 917 } else { 918 ierr = ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 919 } 920 ierr = DMPlexGetConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 921 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 922 ierr = PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank);CHKERRQ(ierr); 923 for (d=depth-1; d>=0; d--) { 924 IS checkIS; 925 PetscBool flg; 926 927 ierr = PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d);CHKERRQ(ierr); 928 ierr = PetscSectionView(sects[d], sviewer);CHKERRQ(ierr); 929 ierr = ISView(iss[d], sviewer);CHKERRQ(ierr); 930 /* check reverse operation */ 931 if (d < depth-1) { 932 ierr = DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS);CHKERRQ(ierr); 933 ierr = ISEqualUnsorted(checkIS, iss[d+1], &flg);CHKERRQ(ierr); 934 PetscCheckFalse(!flg,PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS"); 935 ierr = ISDestroy(&checkIS);CHKERRQ(ierr); 936 } 937 } 938 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 939 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 940 ierr = DMPlexRestoreConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 941 ierr = ISDestroy(&is);CHKERRQ(ierr); 942 PetscFunctionReturn(0); 943 } 944 945 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis) 946 { 947 PetscInt n,n1,ncone,numCoveredPoints,o,p,q,start,end; 948 const PetscInt *coveredPoints; 949 const PetscInt *arr, *cone; 950 PetscInt *newarr; 951 PetscErrorCode ierr; 952 953 PetscFunctionBegin; 954 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 955 ierr = PetscSectionGetStorageSize(section, &n1);CHKERRQ(ierr); 956 ierr = PetscSectionGetChart(section, &start, &end);CHKERRQ(ierr); 957 PetscCheckFalse(n != n1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size", n, n1); 958 ierr = ISGetIndices(is, &arr);CHKERRQ(ierr); 959 ierr = PetscMalloc1(end-start, &newarr);CHKERRQ(ierr); 960 for (q=start; q<end; q++) { 961 ierr = PetscSectionGetDof(section, q, &ncone);CHKERRQ(ierr); 962 ierr = PetscSectionGetOffset(section, q, &o);CHKERRQ(ierr); 963 cone = &arr[o]; 964 if (ncone == 1) { 965 numCoveredPoints = 1; 966 p = cone[0]; 967 } else { 968 PetscInt i; 969 p = PETSC_MAX_INT; 970 for (i=0; i<ncone; i++) if (cone[i] < 0) {p = -1; break;} 971 if (p >= 0) { 972 ierr = DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 973 PetscCheckFalse(numCoveredPoints > 1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q); 974 if (numCoveredPoints) p = coveredPoints[0]; 975 else p = -2; 976 ierr = DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 977 } 978 } 979 newarr[q-start] = p; 980 } 981 ierr = ISRestoreIndices(is, &arr);CHKERRQ(ierr); 982 ierr = ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis);CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 986 static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is) 987 { 988 PetscInt d; 989 IS is,newis; 990 PetscErrorCode ierr; 991 992 PetscFunctionBegin; 993 is = boundary_expanded_is; 994 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 995 for (d = 0; d < depth-1; ++d) { 996 ierr = DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis);CHKERRQ(ierr); 997 ierr = ISDestroy(&is);CHKERRQ(ierr); 998 is = newis; 999 } 1000 *boundary_is = is; 1001 PetscFunctionReturn(0); 1002 } 1003 1004 #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; } 1005 1006 static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm) 1007 { 1008 PetscErrorCode ierr; 1009 PetscViewer viewer; 1010 PetscBool flg; 1011 static PetscBool incall = PETSC_FALSE; 1012 PetscViewerFormat format; 1013 1014 PetscFunctionBegin; 1015 if (incall) PetscFunctionReturn(0); 1016 incall = PETSC_TRUE; 1017 ierr = PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg);CHKERRQI(incall,ierr); 1018 if (flg) { 1019 ierr = PetscViewerPushFormat(viewer,format);CHKERRQI(incall,ierr); 1020 ierr = DMLabelView(label, viewer);CHKERRQI(incall,ierr); 1021 ierr = PetscViewerPopFormat(viewer);CHKERRQI(incall,ierr); 1022 ierr = PetscViewerDestroy(&viewer);CHKERRQI(incall,ierr); 1023 } 1024 incall = PETSC_FALSE; 1025 PetscFunctionReturn(0); 1026 } 1027 1028 /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */ 1029 static inline PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is) 1030 { 1031 IS tmpis; 1032 PetscErrorCode ierr; 1033 1034 PetscFunctionBegin; 1035 ierr = DMLabelGetStratumIS(label, value, &tmpis);CHKERRQ(ierr); 1036 if (!tmpis) {ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis);CHKERRQ(ierr);} 1037 ierr = ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is);CHKERRQ(ierr); 1038 ierr = ISDestroy(&tmpis);CHKERRQ(ierr); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 /* currently only for simple PetscSection without fields or constraints */ 1043 static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout) 1044 { 1045 PetscSection sec; 1046 PetscInt chart[2], p; 1047 PetscInt *dofarr; 1048 PetscMPIInt rank; 1049 PetscErrorCode ierr; 1050 1051 PetscFunctionBegin; 1052 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1053 if (rank == rootrank) { 1054 ierr = PetscSectionGetChart(sec0, &chart[0], &chart[1]);CHKERRQ(ierr); 1055 } 1056 ierr = MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1057 ierr = PetscMalloc1(chart[1]-chart[0], &dofarr);CHKERRQ(ierr); 1058 if (rank == rootrank) { 1059 for (p = chart[0]; p < chart[1]; p++) { 1060 ierr = PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]]);CHKERRQ(ierr); 1061 } 1062 } 1063 ierr = MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1064 ierr = PetscSectionCreate(comm, &sec);CHKERRQ(ierr); 1065 ierr = PetscSectionSetChart(sec, chart[0], chart[1]);CHKERRQ(ierr); 1066 for (p = chart[0]; p < chart[1]; p++) { 1067 ierr = PetscSectionSetDof(sec, p, dofarr[p-chart[0]]);CHKERRQ(ierr); 1068 } 1069 ierr = PetscSectionSetUp(sec);CHKERRQ(ierr); 1070 ierr = PetscFree(dofarr);CHKERRQ(ierr); 1071 *secout = sec; 1072 PetscFunctionReturn(0); 1073 } 1074 1075 static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is) 1076 { 1077 IS faces_expanded_is; 1078 PetscErrorCode ierr; 1079 1080 PetscFunctionBegin; 1081 ierr = DMPlexFindVertices(ipdm, bnd->coordinates, 0.0, &faces_expanded_is);CHKERRQ(ierr); 1082 ierr = DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is);CHKERRQ(ierr); 1083 ierr = ISDestroy(&faces_expanded_is);CHKERRQ(ierr); 1084 PetscFunctionReturn(0); 1085 } 1086 1087 /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */ 1088 static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable) 1089 { 1090 PetscOptions options = NULL; 1091 const char *prefix = NULL; 1092 const char opt[] = "-dm_plex_interpolate_orient_interfaces"; 1093 char prefix_opt[512]; 1094 PetscBool flg, set; 1095 static PetscBool wasSetTrue = PETSC_FALSE; 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 if (dm) { 1100 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);CHKERRQ(ierr); 1101 options = ((PetscObject)dm)->options; 1102 } 1103 ierr = PetscStrcpy(prefix_opt, "-");CHKERRQ(ierr); 1104 ierr = PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt));CHKERRQ(ierr); 1105 ierr = PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt));CHKERRQ(ierr); 1106 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1107 if (!enable) { 1108 if (set && flg) wasSetTrue = PETSC_TRUE; 1109 ierr = PetscOptionsSetValue(options, prefix_opt, "0");CHKERRQ(ierr); 1110 } else if (set && !flg) { 1111 if (wasSetTrue) { 1112 ierr = PetscOptionsSetValue(options, prefix_opt, "1");CHKERRQ(ierr); 1113 } else { 1114 /* default is PETSC_TRUE */ 1115 ierr = PetscOptionsClearValue(options, prefix_opt);CHKERRQ(ierr); 1116 } 1117 wasSetTrue = PETSC_FALSE; 1118 } 1119 if (PetscDefined(USE_DEBUG)) { 1120 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1121 PetscCheckFalse(set && flg != enable,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect"); 1122 } 1123 PetscFunctionReturn(0); 1124 } 1125 1126 /* get coordinate description of the whole-domain boundary */ 1127 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary) 1128 { 1129 PortableBoundary bnd0, bnd; 1130 MPI_Comm comm; 1131 DM idm; 1132 DMLabel label; 1133 PetscInt d; 1134 const char boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary"; 1135 IS boundary_is; 1136 IS *boundary_expanded_iss; 1137 PetscMPIInt rootrank = 0; 1138 PetscMPIInt rank, size; 1139 PetscInt value = 1; 1140 DMPlexInterpolatedFlag intp; 1141 PetscBool flg; 1142 PetscErrorCode ierr; 1143 1144 PetscFunctionBegin; 1145 ierr = PetscNew(&bnd);CHKERRQ(ierr); 1146 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1147 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1148 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1149 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 1150 PetscCheckFalse(flg,comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed"); 1151 1152 /* interpolate serial DM if not yet interpolated */ 1153 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1154 if (intp == DMPLEX_INTERPOLATED_FULL) { 1155 idm = dm; 1156 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 1157 } else { 1158 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 1159 ierr = DMViewFromOptions(idm, NULL, "-idm_view");CHKERRQ(ierr); 1160 } 1161 1162 /* mark whole-domain boundary of the serial DM */ 1163 ierr = DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label);CHKERRQ(ierr); 1164 ierr = DMAddLabel(idm, label);CHKERRQ(ierr); 1165 ierr = DMPlexMarkBoundaryFaces(idm, value, label);CHKERRQ(ierr); 1166 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm);CHKERRQ(ierr); 1167 ierr = DMLabelGetStratumIS(label, value, &boundary_is);CHKERRQ(ierr); 1168 1169 /* translate to coordinates */ 1170 ierr = PetscNew(&bnd0);CHKERRQ(ierr); 1171 ierr = DMGetCoordinatesLocalSetUp(idm);CHKERRQ(ierr); 1172 if (rank == rootrank) { 1173 ierr = DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1174 ierr = DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates);CHKERRQ(ierr); 1175 /* self-check */ 1176 { 1177 IS is0; 1178 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0);CHKERRQ(ierr); 1179 ierr = ISEqual(is0, boundary_is, &flg);CHKERRQ(ierr); 1180 PetscCheckFalse(!flg,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS"); 1181 ierr = ISDestroy(&is0);CHKERRQ(ierr); 1182 } 1183 } else { 1184 ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates);CHKERRQ(ierr); 1185 } 1186 1187 { 1188 Vec tmp; 1189 VecScatter sc; 1190 IS xis; 1191 PetscInt n; 1192 1193 /* just convert seq vectors to mpi vector */ 1194 ierr = VecGetLocalSize(bnd0->coordinates, &n);CHKERRQ(ierr); 1195 ierr = MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1196 if (rank == rootrank) { 1197 ierr = VecCreateMPI(comm, n, n, &tmp);CHKERRQ(ierr); 1198 } else { 1199 ierr = VecCreateMPI(comm, 0, n, &tmp);CHKERRQ(ierr); 1200 } 1201 ierr = VecCopy(bnd0->coordinates, tmp);CHKERRQ(ierr); 1202 ierr = VecDestroy(&bnd0->coordinates);CHKERRQ(ierr); 1203 bnd0->coordinates = tmp; 1204 1205 /* replicate coordinates from root rank to all ranks */ 1206 ierr = VecCreateMPI(comm, n, n*size, &bnd->coordinates);CHKERRQ(ierr); 1207 ierr = ISCreateStride(comm, n, 0, 1, &xis);CHKERRQ(ierr); 1208 ierr = VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc);CHKERRQ(ierr); 1209 ierr = VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1210 ierr = VecScatterEnd( sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1211 ierr = VecScatterDestroy(&sc);CHKERRQ(ierr); 1212 ierr = ISDestroy(&xis);CHKERRQ(ierr); 1213 } 1214 bnd->depth = bnd0->depth; 1215 ierr = MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1216 ierr = PetscMalloc1(bnd->depth, &bnd->sections);CHKERRQ(ierr); 1217 for (d=0; d<bnd->depth; d++) { 1218 ierr = PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]);CHKERRQ(ierr); 1219 } 1220 1221 if (rank == rootrank) { 1222 ierr = DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1223 } 1224 ierr = PortableBoundaryDestroy(&bnd0);CHKERRQ(ierr); 1225 ierr = DMRemoveLabelBySelf(idm, &label, PETSC_TRUE);CHKERRQ(ierr); 1226 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1227 ierr = ISDestroy(&boundary_is);CHKERRQ(ierr); 1228 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1229 *boundary = bnd; 1230 PetscFunctionReturn(0); 1231 } 1232 1233 /* get faces of inter-partition interface */ 1234 static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is) 1235 { 1236 MPI_Comm comm; 1237 DMLabel label; 1238 IS part_boundary_faces_is; 1239 const char partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary"; 1240 PetscInt value = 1; 1241 DMPlexInterpolatedFlag intp; 1242 PetscErrorCode ierr; 1243 1244 PetscFunctionBegin; 1245 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1246 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1247 PetscCheckFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1248 1249 /* get ipdm partition boundary (partBoundary) */ 1250 ierr = DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label);CHKERRQ(ierr); 1251 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1252 ierr = DMPlexMarkBoundaryFaces(ipdm, value, label);CHKERRQ(ierr); 1253 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm);CHKERRQ(ierr); 1254 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is);CHKERRQ(ierr); 1255 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1256 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1257 1258 /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */ 1259 ierr = ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is);CHKERRQ(ierr); 1260 ierr = ISDestroy(&part_boundary_faces_is);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 /* compute inter-partition interface including edges and vertices */ 1265 static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is) 1266 { 1267 DMLabel label; 1268 PetscInt value = 1; 1269 const char interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface"; 1270 DMPlexInterpolatedFlag intp; 1271 MPI_Comm comm; 1272 PetscErrorCode ierr; 1273 1274 PetscFunctionBegin; 1275 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1276 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1277 PetscCheckFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1278 1279 ierr = DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label);CHKERRQ(ierr); 1280 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1281 ierr = DMLabelSetStratumIS(label, value, interface_faces_is);CHKERRQ(ierr); 1282 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm);CHKERRQ(ierr); 1283 ierr = DMPlexLabelComplete(ipdm, label);CHKERRQ(ierr); 1284 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm);CHKERRQ(ierr); 1285 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is);CHKERRQ(ierr); 1286 ierr = PetscObjectSetName((PetscObject)*interface_is, "interface_is");CHKERRQ(ierr); 1287 ierr = ISViewFromOptions(*interface_is, NULL, "-interface_is_view");CHKERRQ(ierr); 1288 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1289 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1290 PetscFunctionReturn(0); 1291 } 1292 1293 static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is) 1294 { 1295 PetscInt n; 1296 const PetscInt *arr; 1297 PetscErrorCode ierr; 1298 1299 PetscFunctionBegin; 1300 ierr = PetscSFGetGraph(sf, NULL, &n, &arr, NULL);CHKERRQ(ierr); 1301 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is);CHKERRQ(ierr); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is) 1306 { 1307 PetscInt n; 1308 const PetscInt *rootdegree; 1309 PetscInt *arr; 1310 PetscErrorCode ierr; 1311 1312 PetscFunctionBegin; 1313 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1314 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 1315 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 1316 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr);CHKERRQ(ierr); 1317 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 1318 PetscFunctionReturn(0); 1319 } 1320 1321 static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is) 1322 { 1323 IS pointSF_out_is, pointSF_in_is; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 ierr = PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is);CHKERRQ(ierr); 1328 ierr = PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is);CHKERRQ(ierr); 1329 ierr = ISExpand(pointSF_out_is, pointSF_in_is, is);CHKERRQ(ierr); 1330 ierr = ISDestroy(&pointSF_out_is);CHKERRQ(ierr); 1331 ierr = ISDestroy(&pointSF_in_is);CHKERRQ(ierr); 1332 PetscFunctionReturn(0); 1333 } 1334 1335 #define CHKERRMY(ierr) do {PetscCheckFalse(ierr,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!");} while (0) 1336 1337 static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v) 1338 { 1339 DMLabel label; 1340 PetscSection coordsSection; 1341 Vec coordsVec; 1342 PetscScalar *coordsScalar; 1343 PetscInt coneSize, depth, dim, i, p, npoints; 1344 const PetscInt *points; 1345 PetscErrorCode ierr; 1346 1347 PetscFunctionBegin; 1348 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1349 ierr = DMGetCoordinateSection(dm, &coordsSection);CHKERRQ(ierr); 1350 ierr = DMGetCoordinatesLocal(dm, &coordsVec);CHKERRQ(ierr); 1351 ierr = VecGetArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1352 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 1353 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 1354 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 1355 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 1356 for (i=0; i<npoints; i++) { 1357 p = points[i]; 1358 ierr = DMLabelGetValue(label, p, &depth);CHKERRQ(ierr); 1359 if (!depth) { 1360 PetscInt n, o; 1361 char coordstr[128]; 1362 1363 ierr = PetscSectionGetDof(coordsSection, p, &n);CHKERRQ(ierr); 1364 ierr = PetscSectionGetOffset(coordsSection, p, &o);CHKERRQ(ierr); 1365 ierr = coord2str(coordstr, sizeof(coordstr), n, &coordsScalar[o], 1.0);CHKERRQ(ierr); 1366 ierr = PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr);CHKERRQ(ierr); 1367 } else { 1368 char entityType[16]; 1369 1370 switch (depth) { 1371 case 1: ierr = PetscStrcpy(entityType, "edge");CHKERRQ(ierr); break; 1372 case 2: ierr = PetscStrcpy(entityType, "face");CHKERRQ(ierr); break; 1373 case 3: ierr = PetscStrcpy(entityType, "cell");CHKERRQ(ierr); break; 1374 default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3"); 1375 } 1376 if (depth == dim && dim < 3) { 1377 ierr = PetscStrlcat(entityType, " (cell)", sizeof(entityType));CHKERRQ(ierr); 1378 } 1379 ierr = PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p);CHKERRQ(ierr); 1380 } 1381 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1382 if (coneSize) { 1383 const PetscInt *cone; 1384 IS coneIS; 1385 1386 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1387 ierr = ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS);CHKERRQ(ierr); 1388 ierr = ViewPointsWithType_Internal(dm, coneIS, v);CHKERRQ(ierr); 1389 ierr = ISDestroy(&coneIS);CHKERRQ(ierr); 1390 } 1391 } 1392 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 1393 ierr = VecRestoreArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1394 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 1395 PetscFunctionReturn(0); 1396 } 1397 1398 static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v) 1399 { 1400 PetscBool flg; 1401 PetscInt npoints; 1402 PetscMPIInt rank; 1403 PetscErrorCode ierr; 1404 1405 PetscFunctionBegin; 1406 ierr = PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg);CHKERRQ(ierr); 1407 PetscCheckFalse(!flg,PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer"); 1408 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRMPI(ierr); 1409 ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr); 1410 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 1411 if (npoints) { 1412 ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank);CHKERRQ(ierr); 1413 ierr = ViewPointsWithType_Internal(dm, points, v);CHKERRQ(ierr); 1414 } 1415 ierr = PetscViewerFlush(v);CHKERRQ(ierr); 1416 ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr); 1417 PetscFunctionReturn(0); 1418 } 1419 1420 static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is) 1421 { 1422 PetscSF pointsf; 1423 IS pointsf_is; 1424 PetscBool flg; 1425 MPI_Comm comm; 1426 PetscMPIInt size; 1427 PetscErrorCode ierr; 1428 1429 PetscFunctionBegin; 1430 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1431 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1432 ierr = DMGetPointSF(ipdm, &pointsf);CHKERRQ(ierr); 1433 if (pointsf) { 1434 PetscInt nroots; 1435 ierr = PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1436 if (nroots < 0) pointsf = NULL; /* uninitialized SF */ 1437 } 1438 if (!pointsf) { 1439 PetscInt N=0; 1440 if (interface_is) {ierr = ISGetSize(interface_is, &N);CHKERRQ(ierr);} 1441 PetscCheckFalse(N,comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL"); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 /* get PointSF points as IS pointsf_is */ 1446 ierr = PointSFGetInterfacePoints_Private(pointsf, &pointsf_is);CHKERRQ(ierr); 1447 1448 /* compare pointsf_is with interface_is */ 1449 ierr = ISEqual(interface_is, pointsf_is, &flg);CHKERRQ(ierr); 1450 ierr = MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm);CHKERRMPI(ierr); 1451 if (!flg) { 1452 IS pointsf_extra_is, pointsf_missing_is; 1453 PetscViewer errv = PETSC_VIEWER_STDERR_(comm); 1454 ierr = ISDifference(interface_is, pointsf_is, &pointsf_missing_is);CHKERRMY(ierr); 1455 ierr = ISDifference(pointsf_is, interface_is, &pointsf_extra_is);CHKERRMY(ierr); 1456 ierr = PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n");CHKERRMY(ierr); 1457 ierr = ViewPointsWithType(ipdm, pointsf_missing_is, errv);CHKERRMY(ierr); 1458 ierr = PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n");CHKERRMY(ierr); 1459 ierr = ViewPointsWithType(ipdm, pointsf_extra_is, errv);CHKERRMY(ierr); 1460 ierr = ISDestroy(&pointsf_extra_is);CHKERRMY(ierr); 1461 ierr = ISDestroy(&pointsf_missing_is);CHKERRMY(ierr); 1462 SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above."); 1463 } 1464 ierr = ISDestroy(&pointsf_is);CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 /* remove faces & edges from label, leave just vertices */ 1469 static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points) 1470 { 1471 PetscInt vStart, vEnd; 1472 MPI_Comm comm; 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1477 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1478 ierr = ISGeneralFilter(points, vStart, vEnd);CHKERRQ(ierr); 1479 PetscFunctionReturn(0); 1480 } 1481 1482 /* 1483 DMPlexCheckPointSFHeavy - Thoroughly test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points. 1484 1485 Collective 1486 1487 Input Parameters: 1488 . dm - The DMPlex object 1489 1490 Notes: 1491 The input DMPlex must be serial (one partition has all points, the other partitions have no points). 1492 This is a heavy test which involves DMPlexInterpolate() if the input DM is not interpolated yet, and depends on having a representation of the whole-domain boundary (PortableBoundary), which can be obtained only by DMPlexGetExpandedBoundary_Private() (which involves DMPlexInterpolate() of a sequential DM). 1493 This is mainly intended for debugging/testing purposes. 1494 1495 Algorithm: 1496 1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces() 1497 2. boundary faces are translated into vertices using DMPlexGetConeRecursive() and these are translated into coordinates - this description (aka PortableBoundary) is completely independent of partitioning and point numbering 1498 3. the mesh is distributed or loaded in parallel 1499 4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices() 1500 5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces() 1501 6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary 1502 7. check that interface covered by PointSF (union of inward and outward points) is equal to the partition interface for each rank, otherwise print the difference and throw an error 1503 1504 Level: developer 1505 1506 .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces() 1507 */ 1508 static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd) 1509 { 1510 DM ipdm=NULL; 1511 IS boundary_faces_is, interface_faces_is, interface_is; 1512 DMPlexInterpolatedFlag intp; 1513 MPI_Comm comm; 1514 PetscErrorCode ierr; 1515 1516 PetscFunctionBegin; 1517 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1518 1519 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1520 if (intp == DMPLEX_INTERPOLATED_FULL) { 1521 ipdm = dm; 1522 } else { 1523 /* create temporary interpolated DM if input DM is not interpolated */ 1524 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_FALSE);CHKERRQ(ierr); 1525 ierr = DMPlexInterpolate(dm, &ipdm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */ 1526 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_TRUE);CHKERRQ(ierr); 1527 } 1528 ierr = DMViewFromOptions(ipdm, NULL, "-ipdm_view");CHKERRQ(ierr); 1529 1530 /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */ 1531 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is);CHKERRQ(ierr); 1532 /* get inter-partition interface faces (interface_faces_is)*/ 1533 ierr = DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is);CHKERRQ(ierr); 1534 /* compute inter-partition interface including edges and vertices (interface_is) */ 1535 ierr = DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is);CHKERRQ(ierr); 1536 /* destroy immediate ISs */ 1537 ierr = ISDestroy(&boundary_faces_is);CHKERRQ(ierr); 1538 ierr = ISDestroy(&interface_faces_is);CHKERRQ(ierr); 1539 1540 /* for uninterpolated case, keep just vertices in interface */ 1541 if (!intp) { 1542 ierr = DMPlexISFilterVertices_Private(ipdm, interface_is);CHKERRQ(ierr); 1543 ierr = DMDestroy(&ipdm);CHKERRQ(ierr); 1544 } 1545 1546 /* compare PointSF with the boundary reconstructed from coordinates */ 1547 ierr = DMPlexComparePointSFWithInterface_Private(dm, interface_is);CHKERRQ(ierr); 1548 ierr = PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n");CHKERRQ(ierr); 1549 ierr = ISDestroy(&interface_is);CHKERRQ(ierr); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 int main(int argc, char **argv) 1554 { 1555 DM dm; 1556 AppCtx user; 1557 PetscErrorCode ierr; 1558 1559 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 1560 ierr = PetscLogStageRegister("create",&stage[0]);CHKERRQ(ierr); 1561 ierr = PetscLogStageRegister("distribute",&stage[1]);CHKERRQ(ierr); 1562 ierr = PetscLogStageRegister("interpolate",&stage[2]);CHKERRQ(ierr); 1563 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1564 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1565 if (user.nPointsToExpand) { 1566 ierr = TestExpandPoints(dm, &user);CHKERRQ(ierr); 1567 } 1568 if (user.ncoords) { 1569 Vec coords; 1570 1571 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, user.ncoords, user.ncoords, user.coords, &coords);CHKERRQ(ierr); 1572 ierr = ViewVerticesFromCoords(dm, coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1573 ierr = VecDestroy(&coords);CHKERRQ(ierr); 1574 } 1575 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1576 ierr = PetscFinalize(); 1577 return ierr; 1578 } 1579 1580 /*TEST 1581 1582 testset: 1583 nsize: 2 1584 args: -dm_view ascii::ascii_info_detail 1585 args: -dm_plex_check_all 1586 test: 1587 suffix: 1_tri_dist0 1588 args: -distribute 0 -interpolate {{none create}separate output} 1589 test: 1590 suffix: 1_tri_dist1 1591 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1592 test: 1593 suffix: 1_quad_dist0 1594 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1595 test: 1596 suffix: 1_quad_dist1 1597 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1598 test: 1599 suffix: 1_1d_dist1 1600 args: -dim 1 -distribute 1 1601 1602 testset: 1603 nsize: 3 1604 args: -testnum 1 -interpolate create 1605 args: -dm_plex_check_all 1606 test: 1607 suffix: 2 1608 args: -dm_view ascii::ascii_info_detail 1609 test: 1610 suffix: 2a 1611 args: -dm_plex_check_cones_conform_on_interfaces_verbose 1612 test: 1613 suffix: 2b 1614 args: -test_expand_points 0,1,2,5,6 1615 test: 1616 suffix: 2c 1617 args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty 1618 1619 testset: 1620 # the same as 1% for 3D 1621 nsize: 2 1622 args: -dim 3 -dm_view ascii::ascii_info_detail 1623 args: -dm_plex_check_all 1624 test: 1625 suffix: 4_tet_dist0 1626 args: -distribute 0 -interpolate {{none create}separate output} 1627 test: 1628 suffix: 4_tet_dist1 1629 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1630 test: 1631 suffix: 4_hex_dist0 1632 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1633 test: 1634 suffix: 4_hex_dist1 1635 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1636 1637 test: 1638 # the same as 4_tet_dist0 but test different initial orientations 1639 suffix: 4_tet_test_orient 1640 nsize: 2 1641 args: -dim 3 -distribute 0 1642 args: -dm_plex_check_all 1643 args: -rotate_interface_0 {{0 1 2 11 12 13}} 1644 args: -rotate_interface_1 {{0 1 2 11 12 13}} 1645 1646 testset: 1647 requires: exodusii 1648 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo 1649 args: -dm_view ascii::ascii_info_detail 1650 args: -dm_plex_check_all 1651 args: -custom_view 1652 test: 1653 suffix: 5_seq 1654 nsize: 1 1655 args: -distribute 0 -interpolate {{none create}separate output} 1656 test: 1657 # Detail viewing in a non-distributed mesh is broken because the DMLabelView() is collective, but the label is not shared 1658 suffix: 5_dist0 1659 nsize: 2 1660 args: -distribute 0 -interpolate {{none create}separate output} -dm_view 1661 test: 1662 suffix: 5_dist1 1663 nsize: 2 1664 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1665 1666 testset: 1667 nsize: {{1 2 4}} 1668 args: -use_generator 1669 args: -dm_plex_check_all 1670 args: -distribute -interpolate none 1671 test: 1672 suffix: 6_tri 1673 requires: triangle 1674 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1675 test: 1676 suffix: 6_quad 1677 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1678 test: 1679 suffix: 6_tet 1680 requires: ctetgen 1681 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1682 test: 1683 suffix: 6_hex 1684 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1685 testset: 1686 nsize: {{1 2 4}} 1687 args: -use_generator 1688 args: -dm_plex_check_all 1689 args: -distribute -interpolate create 1690 test: 1691 suffix: 6_int_tri 1692 requires: triangle 1693 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1694 test: 1695 suffix: 6_int_quad 1696 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1697 test: 1698 suffix: 6_int_tet 1699 requires: ctetgen 1700 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1701 test: 1702 suffix: 6_int_hex 1703 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1704 testset: 1705 nsize: {{2 4}} 1706 args: -use_generator 1707 args: -dm_plex_check_all 1708 args: -distribute -interpolate after_distribute 1709 test: 1710 suffix: 6_parint_tri 1711 requires: triangle 1712 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1713 test: 1714 suffix: 6_parint_quad 1715 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1716 test: 1717 suffix: 6_parint_tet 1718 requires: ctetgen 1719 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1720 test: 1721 suffix: 6_parint_hex 1722 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1723 1724 testset: # 7 EXODUS 1725 requires: exodusii 1726 args: -dm_plex_check_all 1727 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo 1728 args: -distribute 1729 test: # seq load, simple partitioner 1730 suffix: 7_exo 1731 nsize: {{1 2 4 5}} 1732 args: -interpolate none 1733 test: # seq load, seq interpolation, simple partitioner 1734 suffix: 7_exo_int_simple 1735 nsize: {{1 2 4 5}} 1736 args: -interpolate create 1737 test: # seq load, seq interpolation, metis partitioner 1738 suffix: 7_exo_int_metis 1739 requires: parmetis 1740 nsize: {{2 4 5}} 1741 args: -interpolate create 1742 args: -petscpartitioner_type parmetis 1743 test: # seq load, simple partitioner, par interpolation 1744 suffix: 7_exo_simple_int 1745 nsize: {{2 4 5}} 1746 args: -interpolate after_distribute 1747 test: # seq load, metis partitioner, par interpolation 1748 suffix: 7_exo_metis_int 1749 requires: parmetis 1750 nsize: {{2 4 5}} 1751 args: -interpolate after_distribute 1752 args: -petscpartitioner_type parmetis 1753 1754 testset: # 7 HDF5 SEQUANTIAL LOAD 1755 requires: hdf5 !complex 1756 args: -dm_plex_check_all 1757 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1758 args: -dm_plex_hdf5_force_sequential 1759 args: -distribute 1760 test: # seq load, simple partitioner 1761 suffix: 7_seq_hdf5_simple 1762 nsize: {{1 2 4 5}} 1763 args: -interpolate none 1764 test: # seq load, seq interpolation, simple partitioner 1765 suffix: 7_seq_hdf5_int_simple 1766 nsize: {{1 2 4 5}} 1767 args: -interpolate after_create 1768 test: # seq load, seq interpolation, metis partitioner 1769 nsize: {{2 4 5}} 1770 suffix: 7_seq_hdf5_int_metis 1771 requires: parmetis 1772 args: -interpolate after_create 1773 args: -petscpartitioner_type parmetis 1774 test: # seq load, simple partitioner, par interpolation 1775 suffix: 7_seq_hdf5_simple_int 1776 nsize: {{2 4 5}} 1777 args: -interpolate after_distribute 1778 test: # seq load, metis partitioner, par interpolation 1779 nsize: {{2 4 5}} 1780 suffix: 7_seq_hdf5_metis_int 1781 requires: parmetis 1782 args: -interpolate after_distribute 1783 args: -petscpartitioner_type parmetis 1784 1785 testset: # 7 HDF5 PARALLEL LOAD 1786 requires: hdf5 !complex 1787 nsize: {{2 4 5}} 1788 args: -dm_plex_check_all 1789 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1790 test: # par load 1791 suffix: 7_par_hdf5 1792 args: -interpolate none 1793 test: # par load, par interpolation 1794 suffix: 7_par_hdf5_int 1795 args: -interpolate after_create 1796 test: # par load, parmetis repartitioner 1797 TODO: Parallel partitioning of uninterpolated meshes not supported 1798 suffix: 7_par_hdf5_parmetis 1799 requires: parmetis 1800 args: -distribute -petscpartitioner_type parmetis 1801 args: -interpolate none 1802 test: # par load, par interpolation, parmetis repartitioner 1803 suffix: 7_par_hdf5_int_parmetis 1804 requires: parmetis 1805 args: -distribute -petscpartitioner_type parmetis 1806 args: -interpolate after_create 1807 test: # par load, parmetis partitioner, par interpolation 1808 TODO: Parallel partitioning of uninterpolated meshes not supported 1809 suffix: 7_par_hdf5_parmetis_int 1810 requires: parmetis 1811 args: -distribute -petscpartitioner_type parmetis 1812 args: -interpolate after_distribute 1813 1814 test: 1815 suffix: 7_hdf5_hierarch 1816 requires: hdf5 ptscotch !complex 1817 nsize: {{2 3 4}separate output} 1818 args: -distribute 1819 args: -interpolate after_create 1820 args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info 1821 args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2 1822 args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch 1823 1824 test: 1825 suffix: 8 1826 requires: hdf5 !complex 1827 nsize: 4 1828 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1829 args: -distribute 0 -interpolate after_create 1830 args: -view_vertices_from_coords 0.,1.,0.,-0.5,1.,0.,0.583,-0.644,0.,-2.,-2.,-2. -view_vertices_from_coords_tol 1e-3 1831 args: -dm_plex_check_all 1832 args: -custom_view 1833 1834 testset: # 9 HDF5 SEQUANTIAL LOAD 1835 requires: hdf5 !complex datafilespath 1836 args: -dm_plex_check_all 1837 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /cells -dm_plex_hdf5_geometry_path /coordinates 1838 args: -dm_plex_hdf5_force_sequential 1839 args: -distribute 1840 test: # seq load, simple partitioner 1841 suffix: 9_seq_hdf5_simple 1842 nsize: {{1 2 4 5}} 1843 args: -interpolate none 1844 test: # seq load, seq interpolation, simple partitioner 1845 suffix: 9_seq_hdf5_int_simple 1846 nsize: {{1 2 4 5}} 1847 args: -interpolate after_create 1848 test: # seq load, seq interpolation, metis partitioner 1849 nsize: {{2 4 5}} 1850 suffix: 9_seq_hdf5_int_metis 1851 requires: parmetis 1852 args: -interpolate after_create 1853 args: -petscpartitioner_type parmetis 1854 test: # seq load, simple partitioner, par interpolation 1855 suffix: 9_seq_hdf5_simple_int 1856 nsize: {{2 4 5}} 1857 args: -interpolate after_distribute 1858 test: # seq load, simple partitioner, par interpolation 1859 # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy(). 1860 # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken. 1861 # We can then provide an intentionally broken mesh instead. 1862 TODO: This test is broken because PointSF is fixed. 1863 suffix: 9_seq_hdf5_simple_int_err 1864 nsize: 4 1865 args: -interpolate after_distribute 1866 filter: sed -e "/PETSC ERROR/,$$d" 1867 test: # seq load, metis partitioner, par interpolation 1868 nsize: {{2 4 5}} 1869 suffix: 9_seq_hdf5_metis_int 1870 requires: parmetis 1871 args: -interpolate after_distribute 1872 args: -petscpartitioner_type parmetis 1873 1874 testset: # 9 HDF5 PARALLEL LOAD 1875 requires: hdf5 !complex datafilespath 1876 nsize: {{2 4 5}} 1877 args: -dm_plex_check_all 1878 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /cells -dm_plex_hdf5_geometry_path /coordinates 1879 test: # par load 1880 suffix: 9_par_hdf5 1881 args: -interpolate none 1882 test: # par load, par interpolation 1883 suffix: 9_par_hdf5_int 1884 args: -interpolate after_create 1885 test: # par load, parmetis repartitioner 1886 TODO: Parallel partitioning of uninterpolated meshes not supported 1887 suffix: 9_par_hdf5_parmetis 1888 requires: parmetis 1889 args: -distribute -petscpartitioner_type parmetis 1890 args: -interpolate none 1891 test: # par load, par interpolation, parmetis repartitioner 1892 suffix: 9_par_hdf5_int_parmetis 1893 requires: parmetis 1894 args: -distribute -petscpartitioner_type parmetis 1895 args: -interpolate after_create 1896 test: # par load, parmetis partitioner, par interpolation 1897 TODO: Parallel partitioning of uninterpolated meshes not supported 1898 suffix: 9_par_hdf5_parmetis_int 1899 requires: parmetis 1900 args: -distribute -petscpartitioner_type parmetis 1901 args: -interpolate after_distribute 1902 1903 testset: # 10 HDF5 PARALLEL LOAD 1904 requires: hdf5 !complex datafilespath 1905 nsize: {{2 4 7}} 1906 args: -dm_plex_check_all 1907 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined2.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /topo -dm_plex_hdf5_geometry_path /geom 1908 test: # par load, par interpolation 1909 suffix: 10_par_hdf5_int 1910 args: -interpolate after_create 1911 TEST*/ 1912