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 PetscAssertFalse(!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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 PetscAssertFalse(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 = DMSetFromOptions(*dm);CHKERRQ(ierr); 838 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 839 840 if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 841 ierr = DMDestroy(&serialDM);CHKERRQ(ierr); 842 ierr = PortableBoundaryDestroy(&boundary);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 846 #define ps2d(number) ((double) PetscRealPart(number)) 847 static inline PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, const PetscScalar coords[], PetscReal tol) 848 { 849 PetscErrorCode ierr; 850 851 PetscFunctionBegin; 852 PetscAssertFalse(dim > 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3"); 853 if (tol >= 1e-3) { 854 switch (dim) { 855 case 1: ierr = PetscSNPrintf(buf,len,"(%12.3f)",ps2d(coords[0]));CHKERRQ(ierr); 856 case 2: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]));CHKERRQ(ierr); 857 case 3: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]));CHKERRQ(ierr); 858 } 859 } else { 860 switch (dim) { 861 case 1: ierr = PetscSNPrintf(buf,len,"(%12.6f)",ps2d(coords[0]));CHKERRQ(ierr); 862 case 2: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]));CHKERRQ(ierr); 863 case 3: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]));CHKERRQ(ierr); 864 } 865 } 866 PetscFunctionReturn(0); 867 } 868 869 static PetscErrorCode ViewVerticesFromCoords(DM dm, Vec coordsVec, PetscReal tol, PetscViewer viewer) 870 { 871 PetscInt dim, i, npoints; 872 IS pointsIS; 873 const PetscInt *points; 874 const PetscScalar *coords; 875 char coordstr[128]; 876 MPI_Comm comm; 877 PetscMPIInt rank; 878 PetscErrorCode ierr; 879 880 PetscFunctionBegin; 881 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 882 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 883 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 884 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 885 ierr = DMPlexFindVertices(dm, coordsVec, tol, &pointsIS);CHKERRQ(ierr); 886 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 887 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 888 ierr = VecGetArrayRead(coordsVec, &coords);CHKERRQ(ierr); 889 for (i=0; i < npoints; i++) { 890 ierr = coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol);CHKERRQ(ierr); 891 if (rank == 0 && i) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, "-----\n");CHKERRQ(ierr);} 892 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i]);CHKERRQ(ierr); 893 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 894 } 895 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 896 ierr = VecRestoreArrayRead(coordsVec, &coords);CHKERRQ(ierr); 897 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 898 ierr = ISDestroy(&pointsIS);CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user) 903 { 904 IS is; 905 PetscSection *sects; 906 IS *iss; 907 PetscInt d,depth; 908 PetscMPIInt rank; 909 PetscErrorCode ierr; 910 PetscViewer viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer; 911 912 PetscFunctionBegin; 913 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 914 if (user->testExpandPointsEmpty && rank == 0) { 915 ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 916 } else { 917 ierr = ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 918 } 919 ierr = DMPlexGetConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 920 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 921 ierr = PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank);CHKERRQ(ierr); 922 for (d=depth-1; d>=0; d--) { 923 IS checkIS; 924 PetscBool flg; 925 926 ierr = PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d);CHKERRQ(ierr); 927 ierr = PetscSectionView(sects[d], sviewer);CHKERRQ(ierr); 928 ierr = ISView(iss[d], sviewer);CHKERRQ(ierr); 929 /* check reverse operation */ 930 if (d < depth-1) { 931 ierr = DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS);CHKERRQ(ierr); 932 ierr = ISEqualUnsorted(checkIS, iss[d+1], &flg);CHKERRQ(ierr); 933 PetscAssertFalse(!flg,PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS"); 934 ierr = ISDestroy(&checkIS);CHKERRQ(ierr); 935 } 936 } 937 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 938 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 939 ierr = DMPlexRestoreConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 940 ierr = ISDestroy(&is);CHKERRQ(ierr); 941 PetscFunctionReturn(0); 942 } 943 944 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis) 945 { 946 PetscInt n,n1,ncone,numCoveredPoints,o,p,q,start,end; 947 const PetscInt *coveredPoints; 948 const PetscInt *arr, *cone; 949 PetscInt *newarr; 950 PetscErrorCode ierr; 951 952 PetscFunctionBegin; 953 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 954 ierr = PetscSectionGetStorageSize(section, &n1);CHKERRQ(ierr); 955 ierr = PetscSectionGetChart(section, &start, &end);CHKERRQ(ierr); 956 PetscAssertFalse(n != n1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size", n, n1); 957 ierr = ISGetIndices(is, &arr);CHKERRQ(ierr); 958 ierr = PetscMalloc1(end-start, &newarr);CHKERRQ(ierr); 959 for (q=start; q<end; q++) { 960 ierr = PetscSectionGetDof(section, q, &ncone);CHKERRQ(ierr); 961 ierr = PetscSectionGetOffset(section, q, &o);CHKERRQ(ierr); 962 cone = &arr[o]; 963 if (ncone == 1) { 964 numCoveredPoints = 1; 965 p = cone[0]; 966 } else { 967 PetscInt i; 968 p = PETSC_MAX_INT; 969 for (i=0; i<ncone; i++) if (cone[i] < 0) {p = -1; break;} 970 if (p >= 0) { 971 ierr = DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 972 PetscAssertFalse(numCoveredPoints > 1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q); 973 if (numCoveredPoints) p = coveredPoints[0]; 974 else p = -2; 975 ierr = DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 976 } 977 } 978 newarr[q-start] = p; 979 } 980 ierr = ISRestoreIndices(is, &arr);CHKERRQ(ierr); 981 ierr = ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis);CHKERRQ(ierr); 982 PetscFunctionReturn(0); 983 } 984 985 static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is) 986 { 987 PetscInt d; 988 IS is,newis; 989 PetscErrorCode ierr; 990 991 PetscFunctionBegin; 992 is = boundary_expanded_is; 993 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 994 for (d = 0; d < depth-1; ++d) { 995 ierr = DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis);CHKERRQ(ierr); 996 ierr = ISDestroy(&is);CHKERRQ(ierr); 997 is = newis; 998 } 999 *boundary_is = is; 1000 PetscFunctionReturn(0); 1001 } 1002 1003 #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; } 1004 1005 static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm) 1006 { 1007 PetscErrorCode ierr; 1008 PetscViewer viewer; 1009 PetscBool flg; 1010 static PetscBool incall = PETSC_FALSE; 1011 PetscViewerFormat format; 1012 1013 PetscFunctionBegin; 1014 if (incall) PetscFunctionReturn(0); 1015 incall = PETSC_TRUE; 1016 ierr = PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg);CHKERRQI(incall,ierr); 1017 if (flg) { 1018 ierr = PetscViewerPushFormat(viewer,format);CHKERRQI(incall,ierr); 1019 ierr = DMLabelView(label, viewer);CHKERRQI(incall,ierr); 1020 ierr = PetscViewerPopFormat(viewer);CHKERRQI(incall,ierr); 1021 ierr = PetscViewerDestroy(&viewer);CHKERRQI(incall,ierr); 1022 } 1023 incall = PETSC_FALSE; 1024 PetscFunctionReturn(0); 1025 } 1026 1027 /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */ 1028 static inline PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is) 1029 { 1030 IS tmpis; 1031 PetscErrorCode ierr; 1032 1033 PetscFunctionBegin; 1034 ierr = DMLabelGetStratumIS(label, value, &tmpis);CHKERRQ(ierr); 1035 if (!tmpis) {ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis);CHKERRQ(ierr);} 1036 ierr = ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is);CHKERRQ(ierr); 1037 ierr = ISDestroy(&tmpis);CHKERRQ(ierr); 1038 PetscFunctionReturn(0); 1039 } 1040 1041 /* currently only for simple PetscSection without fields or constraints */ 1042 static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout) 1043 { 1044 PetscSection sec; 1045 PetscInt chart[2], p; 1046 PetscInt *dofarr; 1047 PetscMPIInt rank; 1048 PetscErrorCode ierr; 1049 1050 PetscFunctionBegin; 1051 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1052 if (rank == rootrank) { 1053 ierr = PetscSectionGetChart(sec0, &chart[0], &chart[1]);CHKERRQ(ierr); 1054 } 1055 ierr = MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1056 ierr = PetscMalloc1(chart[1]-chart[0], &dofarr);CHKERRQ(ierr); 1057 if (rank == rootrank) { 1058 for (p = chart[0]; p < chart[1]; p++) { 1059 ierr = PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]]);CHKERRQ(ierr); 1060 } 1061 } 1062 ierr = MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1063 ierr = PetscSectionCreate(comm, &sec);CHKERRQ(ierr); 1064 ierr = PetscSectionSetChart(sec, chart[0], chart[1]);CHKERRQ(ierr); 1065 for (p = chart[0]; p < chart[1]; p++) { 1066 ierr = PetscSectionSetDof(sec, p, dofarr[p-chart[0]]);CHKERRQ(ierr); 1067 } 1068 ierr = PetscSectionSetUp(sec);CHKERRQ(ierr); 1069 ierr = PetscFree(dofarr);CHKERRQ(ierr); 1070 *secout = sec; 1071 PetscFunctionReturn(0); 1072 } 1073 1074 static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is) 1075 { 1076 IS faces_expanded_is; 1077 PetscErrorCode ierr; 1078 1079 PetscFunctionBegin; 1080 ierr = DMPlexFindVertices(ipdm, bnd->coordinates, 0.0, &faces_expanded_is);CHKERRQ(ierr); 1081 ierr = DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is);CHKERRQ(ierr); 1082 ierr = ISDestroy(&faces_expanded_is);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */ 1087 static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable) 1088 { 1089 PetscOptions options = NULL; 1090 const char *prefix = NULL; 1091 const char opt[] = "-dm_plex_interpolate_orient_interfaces"; 1092 char prefix_opt[512]; 1093 PetscBool flg, set; 1094 static PetscBool wasSetTrue = PETSC_FALSE; 1095 PetscErrorCode ierr; 1096 1097 PetscFunctionBegin; 1098 if (dm) { 1099 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);CHKERRQ(ierr); 1100 options = ((PetscObject)dm)->options; 1101 } 1102 ierr = PetscStrcpy(prefix_opt, "-");CHKERRQ(ierr); 1103 ierr = PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt));CHKERRQ(ierr); 1104 ierr = PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt));CHKERRQ(ierr); 1105 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1106 if (!enable) { 1107 if (set && flg) wasSetTrue = PETSC_TRUE; 1108 ierr = PetscOptionsSetValue(options, prefix_opt, "0");CHKERRQ(ierr); 1109 } else if (set && !flg) { 1110 if (wasSetTrue) { 1111 ierr = PetscOptionsSetValue(options, prefix_opt, "1");CHKERRQ(ierr); 1112 } else { 1113 /* default is PETSC_TRUE */ 1114 ierr = PetscOptionsClearValue(options, prefix_opt);CHKERRQ(ierr); 1115 } 1116 wasSetTrue = PETSC_FALSE; 1117 } 1118 if (PetscDefined(USE_DEBUG)) { 1119 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1120 PetscAssertFalse(set && flg != enable,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect"); 1121 } 1122 PetscFunctionReturn(0); 1123 } 1124 1125 /* get coordinate description of the whole-domain boundary */ 1126 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary) 1127 { 1128 PortableBoundary bnd0, bnd; 1129 MPI_Comm comm; 1130 DM idm; 1131 DMLabel label; 1132 PetscInt d; 1133 const char boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary"; 1134 IS boundary_is; 1135 IS *boundary_expanded_iss; 1136 PetscMPIInt rootrank = 0; 1137 PetscMPIInt rank, size; 1138 PetscInt value = 1; 1139 DMPlexInterpolatedFlag intp; 1140 PetscBool flg; 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 ierr = PetscNew(&bnd);CHKERRQ(ierr); 1145 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1146 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1147 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1148 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 1149 PetscAssertFalse(flg,comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed"); 1150 1151 /* interpolate serial DM if not yet interpolated */ 1152 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1153 if (intp == DMPLEX_INTERPOLATED_FULL) { 1154 idm = dm; 1155 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 1156 } else { 1157 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 1158 ierr = DMViewFromOptions(idm, NULL, "-idm_view");CHKERRQ(ierr); 1159 } 1160 1161 /* mark whole-domain boundary of the serial DM */ 1162 ierr = DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label);CHKERRQ(ierr); 1163 ierr = DMAddLabel(idm, label);CHKERRQ(ierr); 1164 ierr = DMPlexMarkBoundaryFaces(idm, value, label);CHKERRQ(ierr); 1165 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm);CHKERRQ(ierr); 1166 ierr = DMLabelGetStratumIS(label, value, &boundary_is);CHKERRQ(ierr); 1167 1168 /* translate to coordinates */ 1169 ierr = PetscNew(&bnd0);CHKERRQ(ierr); 1170 ierr = DMGetCoordinatesLocalSetUp(idm);CHKERRQ(ierr); 1171 if (rank == rootrank) { 1172 ierr = DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1173 ierr = DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates);CHKERRQ(ierr); 1174 /* self-check */ 1175 { 1176 IS is0; 1177 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0);CHKERRQ(ierr); 1178 ierr = ISEqual(is0, boundary_is, &flg);CHKERRQ(ierr); 1179 PetscAssertFalse(!flg,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS"); 1180 ierr = ISDestroy(&is0);CHKERRQ(ierr); 1181 } 1182 } else { 1183 ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates);CHKERRQ(ierr); 1184 } 1185 1186 { 1187 Vec tmp; 1188 VecScatter sc; 1189 IS xis; 1190 PetscInt n; 1191 1192 /* just convert seq vectors to mpi vector */ 1193 ierr = VecGetLocalSize(bnd0->coordinates, &n);CHKERRQ(ierr); 1194 ierr = MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1195 if (rank == rootrank) { 1196 ierr = VecCreateMPI(comm, n, n, &tmp);CHKERRQ(ierr); 1197 } else { 1198 ierr = VecCreateMPI(comm, 0, n, &tmp);CHKERRQ(ierr); 1199 } 1200 ierr = VecCopy(bnd0->coordinates, tmp);CHKERRQ(ierr); 1201 ierr = VecDestroy(&bnd0->coordinates);CHKERRQ(ierr); 1202 bnd0->coordinates = tmp; 1203 1204 /* replicate coordinates from root rank to all ranks */ 1205 ierr = VecCreateMPI(comm, n, n*size, &bnd->coordinates);CHKERRQ(ierr); 1206 ierr = ISCreateStride(comm, n, 0, 1, &xis);CHKERRQ(ierr); 1207 ierr = VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc);CHKERRQ(ierr); 1208 ierr = VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1209 ierr = VecScatterEnd( sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1210 ierr = VecScatterDestroy(&sc);CHKERRQ(ierr); 1211 ierr = ISDestroy(&xis);CHKERRQ(ierr); 1212 } 1213 bnd->depth = bnd0->depth; 1214 ierr = MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm);CHKERRMPI(ierr); 1215 ierr = PetscMalloc1(bnd->depth, &bnd->sections);CHKERRQ(ierr); 1216 for (d=0; d<bnd->depth; d++) { 1217 ierr = PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]);CHKERRQ(ierr); 1218 } 1219 1220 if (rank == rootrank) { 1221 ierr = DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1222 } 1223 ierr = PortableBoundaryDestroy(&bnd0);CHKERRQ(ierr); 1224 ierr = DMRemoveLabelBySelf(idm, &label, PETSC_TRUE);CHKERRQ(ierr); 1225 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1226 ierr = ISDestroy(&boundary_is);CHKERRQ(ierr); 1227 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1228 *boundary = bnd; 1229 PetscFunctionReturn(0); 1230 } 1231 1232 /* get faces of inter-partition interface */ 1233 static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is) 1234 { 1235 MPI_Comm comm; 1236 DMLabel label; 1237 IS part_boundary_faces_is; 1238 const char partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary"; 1239 PetscInt value = 1; 1240 DMPlexInterpolatedFlag intp; 1241 PetscErrorCode ierr; 1242 1243 PetscFunctionBegin; 1244 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1245 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1246 PetscAssertFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1247 1248 /* get ipdm partition boundary (partBoundary) */ 1249 ierr = DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label);CHKERRQ(ierr); 1250 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1251 ierr = DMPlexMarkBoundaryFaces(ipdm, value, label);CHKERRQ(ierr); 1252 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm);CHKERRQ(ierr); 1253 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is);CHKERRQ(ierr); 1254 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1255 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1256 1257 /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */ 1258 ierr = ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is);CHKERRQ(ierr); 1259 ierr = ISDestroy(&part_boundary_faces_is);CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262 1263 /* compute inter-partition interface including edges and vertices */ 1264 static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is) 1265 { 1266 DMLabel label; 1267 PetscInt value = 1; 1268 const char interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface"; 1269 DMPlexInterpolatedFlag intp; 1270 MPI_Comm comm; 1271 PetscErrorCode ierr; 1272 1273 PetscFunctionBegin; 1274 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1275 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1276 PetscAssertFalse(intp != DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1277 1278 ierr = DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label);CHKERRQ(ierr); 1279 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1280 ierr = DMLabelSetStratumIS(label, value, interface_faces_is);CHKERRQ(ierr); 1281 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm);CHKERRQ(ierr); 1282 ierr = DMPlexLabelComplete(ipdm, label);CHKERRQ(ierr); 1283 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm);CHKERRQ(ierr); 1284 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is);CHKERRQ(ierr); 1285 ierr = PetscObjectSetName((PetscObject)*interface_is, "interface_is");CHKERRQ(ierr); 1286 ierr = ISViewFromOptions(*interface_is, NULL, "-interface_is_view");CHKERRQ(ierr); 1287 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1288 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is) 1293 { 1294 PetscInt n; 1295 const PetscInt *arr; 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 ierr = PetscSFGetGraph(sf, NULL, &n, &arr, NULL);CHKERRQ(ierr); 1300 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is);CHKERRQ(ierr); 1301 PetscFunctionReturn(0); 1302 } 1303 1304 static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is) 1305 { 1306 PetscInt n; 1307 const PetscInt *rootdegree; 1308 PetscInt *arr; 1309 PetscErrorCode ierr; 1310 1311 PetscFunctionBegin; 1312 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1313 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 1314 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 1315 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr);CHKERRQ(ierr); 1316 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 1317 PetscFunctionReturn(0); 1318 } 1319 1320 static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is) 1321 { 1322 IS pointSF_out_is, pointSF_in_is; 1323 PetscErrorCode ierr; 1324 1325 PetscFunctionBegin; 1326 ierr = PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is);CHKERRQ(ierr); 1327 ierr = PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is);CHKERRQ(ierr); 1328 ierr = ISExpand(pointSF_out_is, pointSF_in_is, is);CHKERRQ(ierr); 1329 ierr = ISDestroy(&pointSF_out_is);CHKERRQ(ierr); 1330 ierr = ISDestroy(&pointSF_in_is);CHKERRQ(ierr); 1331 PetscFunctionReturn(0); 1332 } 1333 1334 #define CHKERRMY(ierr) do {PetscAssertFalse(ierr,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!");} while (0) 1335 1336 static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v) 1337 { 1338 DMLabel label; 1339 PetscSection coordsSection; 1340 Vec coordsVec; 1341 PetscScalar *coordsScalar; 1342 PetscInt coneSize, depth, dim, i, p, npoints; 1343 const PetscInt *points; 1344 PetscErrorCode ierr; 1345 1346 PetscFunctionBegin; 1347 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1348 ierr = DMGetCoordinateSection(dm, &coordsSection);CHKERRQ(ierr); 1349 ierr = DMGetCoordinatesLocal(dm, &coordsVec);CHKERRQ(ierr); 1350 ierr = VecGetArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1351 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 1352 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 1353 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 1354 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 1355 for (i=0; i<npoints; i++) { 1356 p = points[i]; 1357 ierr = DMLabelGetValue(label, p, &depth);CHKERRQ(ierr); 1358 if (!depth) { 1359 PetscInt n, o; 1360 char coordstr[128]; 1361 1362 ierr = PetscSectionGetDof(coordsSection, p, &n);CHKERRQ(ierr); 1363 ierr = PetscSectionGetOffset(coordsSection, p, &o);CHKERRQ(ierr); 1364 ierr = coord2str(coordstr, sizeof(coordstr), n, &coordsScalar[o], 1.0);CHKERRQ(ierr); 1365 ierr = PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr);CHKERRQ(ierr); 1366 } else { 1367 char entityType[16]; 1368 1369 switch (depth) { 1370 case 1: ierr = PetscStrcpy(entityType, "edge");CHKERRQ(ierr); break; 1371 case 2: ierr = PetscStrcpy(entityType, "face");CHKERRQ(ierr); break; 1372 case 3: ierr = PetscStrcpy(entityType, "cell");CHKERRQ(ierr); break; 1373 default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3"); 1374 } 1375 if (depth == dim && dim < 3) { 1376 ierr = PetscStrlcat(entityType, " (cell)", sizeof(entityType));CHKERRQ(ierr); 1377 } 1378 ierr = PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p);CHKERRQ(ierr); 1379 } 1380 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1381 if (coneSize) { 1382 const PetscInt *cone; 1383 IS coneIS; 1384 1385 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1386 ierr = ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS);CHKERRQ(ierr); 1387 ierr = ViewPointsWithType_Internal(dm, coneIS, v);CHKERRQ(ierr); 1388 ierr = ISDestroy(&coneIS);CHKERRQ(ierr); 1389 } 1390 } 1391 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 1392 ierr = VecRestoreArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1393 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 1394 PetscFunctionReturn(0); 1395 } 1396 1397 static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v) 1398 { 1399 PetscBool flg; 1400 PetscInt npoints; 1401 PetscMPIInt rank; 1402 PetscErrorCode ierr; 1403 1404 PetscFunctionBegin; 1405 ierr = PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg);CHKERRQ(ierr); 1406 PetscAssertFalse(!flg,PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer"); 1407 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRMPI(ierr); 1408 ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr); 1409 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 1410 if (npoints) { 1411 ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank);CHKERRQ(ierr); 1412 ierr = ViewPointsWithType_Internal(dm, points, v);CHKERRQ(ierr); 1413 } 1414 ierr = PetscViewerFlush(v);CHKERRQ(ierr); 1415 ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is) 1420 { 1421 PetscSF pointsf; 1422 IS pointsf_is; 1423 PetscBool flg; 1424 MPI_Comm comm; 1425 PetscMPIInt size; 1426 PetscErrorCode ierr; 1427 1428 PetscFunctionBegin; 1429 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1430 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 1431 ierr = DMGetPointSF(ipdm, &pointsf);CHKERRQ(ierr); 1432 if (pointsf) { 1433 PetscInt nroots; 1434 ierr = PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1435 if (nroots < 0) pointsf = NULL; /* uninitialized SF */ 1436 } 1437 if (!pointsf) { 1438 PetscInt N=0; 1439 if (interface_is) {ierr = ISGetSize(interface_is, &N);CHKERRQ(ierr);} 1440 PetscAssertFalse(N,comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL"); 1441 PetscFunctionReturn(0); 1442 } 1443 1444 /* get PointSF points as IS pointsf_is */ 1445 ierr = PointSFGetInterfacePoints_Private(pointsf, &pointsf_is);CHKERRQ(ierr); 1446 1447 /* compare pointsf_is with interface_is */ 1448 ierr = ISEqual(interface_is, pointsf_is, &flg);CHKERRQ(ierr); 1449 ierr = MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm);CHKERRMPI(ierr); 1450 if (!flg) { 1451 IS pointsf_extra_is, pointsf_missing_is; 1452 PetscViewer errv = PETSC_VIEWER_STDERR_(comm); 1453 ierr = ISDifference(interface_is, pointsf_is, &pointsf_missing_is);CHKERRMY(ierr); 1454 ierr = ISDifference(pointsf_is, interface_is, &pointsf_extra_is);CHKERRMY(ierr); 1455 ierr = PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n");CHKERRMY(ierr); 1456 ierr = ViewPointsWithType(ipdm, pointsf_missing_is, errv);CHKERRMY(ierr); 1457 ierr = PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n");CHKERRMY(ierr); 1458 ierr = ViewPointsWithType(ipdm, pointsf_extra_is, errv);CHKERRMY(ierr); 1459 ierr = ISDestroy(&pointsf_extra_is);CHKERRMY(ierr); 1460 ierr = ISDestroy(&pointsf_missing_is);CHKERRMY(ierr); 1461 SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above."); 1462 } 1463 ierr = ISDestroy(&pointsf_is);CHKERRQ(ierr); 1464 PetscFunctionReturn(0); 1465 } 1466 1467 /* remove faces & edges from label, leave just vertices */ 1468 static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points) 1469 { 1470 PetscInt vStart, vEnd; 1471 MPI_Comm comm; 1472 PetscErrorCode ierr; 1473 1474 PetscFunctionBegin; 1475 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1476 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1477 ierr = ISGeneralFilter(points, vStart, vEnd);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 /* 1482 DMPlexCheckPointSFHeavy - Thoroughly test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points. 1483 1484 Collective 1485 1486 Input Parameters: 1487 . dm - The DMPlex object 1488 1489 Notes: 1490 The input DMPlex must be serial (one partition has all points, the other partitions have no points). 1491 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). 1492 This is mainly intended for debugging/testing purposes. 1493 1494 Algorithm: 1495 1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces() 1496 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 1497 3. the mesh is distributed or loaded in parallel 1498 4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices() 1499 5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces() 1500 6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary 1501 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 1502 1503 Level: developer 1504 1505 .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces() 1506 */ 1507 static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd) 1508 { 1509 DM ipdm=NULL; 1510 IS boundary_faces_is, interface_faces_is, interface_is; 1511 DMPlexInterpolatedFlag intp; 1512 MPI_Comm comm; 1513 PetscErrorCode ierr; 1514 1515 PetscFunctionBegin; 1516 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1517 1518 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1519 if (intp == DMPLEX_INTERPOLATED_FULL) { 1520 ipdm = dm; 1521 } else { 1522 /* create temporary interpolated DM if input DM is not interpolated */ 1523 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_FALSE);CHKERRQ(ierr); 1524 ierr = DMPlexInterpolate(dm, &ipdm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */ 1525 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_TRUE);CHKERRQ(ierr); 1526 } 1527 ierr = DMViewFromOptions(ipdm, NULL, "-ipdm_view");CHKERRQ(ierr); 1528 1529 /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */ 1530 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is);CHKERRQ(ierr); 1531 /* get inter-partition interface faces (interface_faces_is)*/ 1532 ierr = DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is);CHKERRQ(ierr); 1533 /* compute inter-partition interface including edges and vertices (interface_is) */ 1534 ierr = DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is);CHKERRQ(ierr); 1535 /* destroy immediate ISs */ 1536 ierr = ISDestroy(&boundary_faces_is);CHKERRQ(ierr); 1537 ierr = ISDestroy(&interface_faces_is);CHKERRQ(ierr); 1538 1539 /* for uninterpolated case, keep just vertices in interface */ 1540 if (!intp) { 1541 ierr = DMPlexISFilterVertices_Private(ipdm, interface_is);CHKERRQ(ierr); 1542 ierr = DMDestroy(&ipdm);CHKERRQ(ierr); 1543 } 1544 1545 /* compare PointSF with the boundary reconstructed from coordinates */ 1546 ierr = DMPlexComparePointSFWithInterface_Private(dm, interface_is);CHKERRQ(ierr); 1547 ierr = PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n");CHKERRQ(ierr); 1548 ierr = ISDestroy(&interface_is);CHKERRQ(ierr); 1549 PetscFunctionReturn(0); 1550 } 1551 1552 int main(int argc, char **argv) 1553 { 1554 DM dm; 1555 AppCtx user; 1556 PetscErrorCode ierr; 1557 1558 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 1559 ierr = PetscLogStageRegister("create",&stage[0]);CHKERRQ(ierr); 1560 ierr = PetscLogStageRegister("distribute",&stage[1]);CHKERRQ(ierr); 1561 ierr = PetscLogStageRegister("interpolate",&stage[2]);CHKERRQ(ierr); 1562 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1563 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1564 if (user.nPointsToExpand) { 1565 ierr = TestExpandPoints(dm, &user);CHKERRQ(ierr); 1566 } 1567 if (user.ncoords) { 1568 Vec coords; 1569 1570 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, user.ncoords, user.ncoords, user.coords, &coords);CHKERRQ(ierr); 1571 ierr = ViewVerticesFromCoords(dm, coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1572 ierr = VecDestroy(&coords);CHKERRQ(ierr); 1573 } 1574 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1575 ierr = PetscFinalize(); 1576 return ierr; 1577 } 1578 1579 /*TEST 1580 1581 testset: 1582 nsize: 2 1583 args: -dm_view ascii::ascii_info_detail 1584 args: -dm_plex_check_all 1585 test: 1586 suffix: 1_tri_dist0 1587 args: -distribute 0 -interpolate {{none create}separate output} 1588 test: 1589 suffix: 1_tri_dist1 1590 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1591 test: 1592 suffix: 1_quad_dist0 1593 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1594 test: 1595 suffix: 1_quad_dist1 1596 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1597 test: 1598 suffix: 1_1d_dist1 1599 args: -dim 1 -distribute 1 1600 1601 testset: 1602 nsize: 3 1603 args: -testnum 1 -interpolate create 1604 args: -dm_plex_check_all 1605 test: 1606 suffix: 2 1607 args: -dm_view ascii::ascii_info_detail 1608 test: 1609 suffix: 2a 1610 args: -dm_plex_check_cones_conform_on_interfaces_verbose 1611 test: 1612 suffix: 2b 1613 args: -test_expand_points 0,1,2,5,6 1614 test: 1615 suffix: 2c 1616 args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty 1617 1618 testset: 1619 # the same as 1% for 3D 1620 nsize: 2 1621 args: -dim 3 -dm_view ascii::ascii_info_detail 1622 args: -dm_plex_check_all 1623 test: 1624 suffix: 4_tet_dist0 1625 args: -distribute 0 -interpolate {{none create}separate output} 1626 test: 1627 suffix: 4_tet_dist1 1628 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1629 test: 1630 suffix: 4_hex_dist0 1631 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1632 test: 1633 suffix: 4_hex_dist1 1634 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1635 1636 test: 1637 # the same as 4_tet_dist0 but test different initial orientations 1638 suffix: 4_tet_test_orient 1639 nsize: 2 1640 args: -dim 3 -distribute 0 1641 args: -dm_plex_check_all 1642 args: -rotate_interface_0 {{0 1 2 11 12 13}} 1643 args: -rotate_interface_1 {{0 1 2 11 12 13}} 1644 1645 testset: 1646 requires: exodusii 1647 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo 1648 args: -dm_view ascii::ascii_info_detail 1649 args: -dm_plex_check_all 1650 args: -custom_view 1651 test: 1652 suffix: 5_seq 1653 nsize: 1 1654 args: -distribute 0 -interpolate {{none create}separate output} 1655 test: 1656 # Detail viewing in a non-distributed mesh is broken because the DMLabelView() is collective, but the label is not shared 1657 suffix: 5_dist0 1658 nsize: 2 1659 args: -distribute 0 -interpolate {{none create}separate output} -dm_view 1660 test: 1661 suffix: 5_dist1 1662 nsize: 2 1663 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1664 1665 testset: 1666 nsize: {{1 2 4}} 1667 args: -use_generator 1668 args: -dm_plex_check_all 1669 args: -distribute -interpolate none 1670 test: 1671 suffix: 6_tri 1672 requires: triangle 1673 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1674 test: 1675 suffix: 6_quad 1676 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1677 test: 1678 suffix: 6_tet 1679 requires: ctetgen 1680 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1681 test: 1682 suffix: 6_hex 1683 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1684 testset: 1685 nsize: {{1 2 4}} 1686 args: -use_generator 1687 args: -dm_plex_check_all 1688 args: -distribute -interpolate create 1689 test: 1690 suffix: 6_int_tri 1691 requires: triangle 1692 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1693 test: 1694 suffix: 6_int_quad 1695 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1696 test: 1697 suffix: 6_int_tet 1698 requires: ctetgen 1699 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1700 test: 1701 suffix: 6_int_hex 1702 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1703 testset: 1704 nsize: {{2 4}} 1705 args: -use_generator 1706 args: -dm_plex_check_all 1707 args: -distribute -interpolate after_distribute 1708 test: 1709 suffix: 6_parint_tri 1710 requires: triangle 1711 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_generator triangle 1712 test: 1713 suffix: 6_parint_quad 1714 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1715 test: 1716 suffix: 6_parint_tet 1717 requires: ctetgen 1718 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_generator ctetgen 1719 test: 1720 suffix: 6_parint_hex 1721 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1722 1723 testset: # 7 EXODUS 1724 requires: exodusii 1725 args: -dm_plex_check_all 1726 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo 1727 args: -distribute 1728 test: # seq load, simple partitioner 1729 suffix: 7_exo 1730 nsize: {{1 2 4 5}} 1731 args: -interpolate none 1732 test: # seq load, seq interpolation, simple partitioner 1733 suffix: 7_exo_int_simple 1734 nsize: {{1 2 4 5}} 1735 args: -interpolate create 1736 test: # seq load, seq interpolation, metis partitioner 1737 suffix: 7_exo_int_metis 1738 requires: parmetis 1739 nsize: {{2 4 5}} 1740 args: -interpolate create 1741 args: -petscpartitioner_type parmetis 1742 test: # seq load, simple partitioner, par interpolation 1743 suffix: 7_exo_simple_int 1744 nsize: {{2 4 5}} 1745 args: -interpolate after_distribute 1746 test: # seq load, metis partitioner, par interpolation 1747 suffix: 7_exo_metis_int 1748 requires: parmetis 1749 nsize: {{2 4 5}} 1750 args: -interpolate after_distribute 1751 args: -petscpartitioner_type parmetis 1752 1753 testset: # 7 HDF5 SEQUANTIAL LOAD 1754 requires: hdf5 !complex 1755 args: -dm_plex_check_all 1756 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1757 args: -dm_plex_hdf5_force_sequential 1758 args: -distribute 1759 test: # seq load, simple partitioner 1760 suffix: 7_seq_hdf5_simple 1761 nsize: {{1 2 4 5}} 1762 args: -interpolate none 1763 test: # seq load, seq interpolation, simple partitioner 1764 suffix: 7_seq_hdf5_int_simple 1765 nsize: {{1 2 4 5}} 1766 args: -interpolate after_create 1767 test: # seq load, seq interpolation, metis partitioner 1768 nsize: {{2 4 5}} 1769 suffix: 7_seq_hdf5_int_metis 1770 requires: parmetis 1771 args: -interpolate after_create 1772 args: -petscpartitioner_type parmetis 1773 test: # seq load, simple partitioner, par interpolation 1774 suffix: 7_seq_hdf5_simple_int 1775 nsize: {{2 4 5}} 1776 args: -interpolate after_distribute 1777 test: # seq load, metis partitioner, par interpolation 1778 nsize: {{2 4 5}} 1779 suffix: 7_seq_hdf5_metis_int 1780 requires: parmetis 1781 args: -interpolate after_distribute 1782 args: -petscpartitioner_type parmetis 1783 1784 testset: # 7 HDF5 PARALLEL LOAD 1785 requires: hdf5 !complex 1786 nsize: {{2 4 5}} 1787 args: -dm_plex_check_all 1788 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1789 test: # par load 1790 suffix: 7_par_hdf5 1791 args: -interpolate none 1792 test: # par load, par interpolation 1793 suffix: 7_par_hdf5_int 1794 args: -interpolate after_create 1795 test: # par load, parmetis repartitioner 1796 TODO: Parallel partitioning of uninterpolated meshes not supported 1797 suffix: 7_par_hdf5_parmetis 1798 requires: parmetis 1799 args: -distribute -petscpartitioner_type parmetis 1800 args: -interpolate none 1801 test: # par load, par interpolation, parmetis repartitioner 1802 suffix: 7_par_hdf5_int_parmetis 1803 requires: parmetis 1804 args: -distribute -petscpartitioner_type parmetis 1805 args: -interpolate after_create 1806 test: # par load, parmetis partitioner, par interpolation 1807 TODO: Parallel partitioning of uninterpolated meshes not supported 1808 suffix: 7_par_hdf5_parmetis_int 1809 requires: parmetis 1810 args: -distribute -petscpartitioner_type parmetis 1811 args: -interpolate after_distribute 1812 1813 test: 1814 suffix: 7_hdf5_hierarch 1815 requires: hdf5 ptscotch !complex 1816 nsize: {{2 3 4}separate output} 1817 args: -distribute 1818 args: -interpolate after_create 1819 args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info 1820 args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2 1821 args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch 1822 1823 test: 1824 suffix: 8 1825 requires: hdf5 !complex 1826 nsize: 4 1827 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1828 args: -distribute 0 -interpolate after_create 1829 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 1830 args: -dm_plex_check_all 1831 args: -custom_view 1832 1833 testset: # 9 HDF5 SEQUANTIAL LOAD 1834 requires: hdf5 !complex datafilespath 1835 args: -dm_plex_check_all 1836 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 1837 args: -dm_plex_hdf5_force_sequential 1838 args: -distribute 1839 test: # seq load, simple partitioner 1840 suffix: 9_seq_hdf5_simple 1841 nsize: {{1 2 4 5}} 1842 args: -interpolate none 1843 test: # seq load, seq interpolation, simple partitioner 1844 suffix: 9_seq_hdf5_int_simple 1845 nsize: {{1 2 4 5}} 1846 args: -interpolate after_create 1847 test: # seq load, seq interpolation, metis partitioner 1848 nsize: {{2 4 5}} 1849 suffix: 9_seq_hdf5_int_metis 1850 requires: parmetis 1851 args: -interpolate after_create 1852 args: -petscpartitioner_type parmetis 1853 test: # seq load, simple partitioner, par interpolation 1854 suffix: 9_seq_hdf5_simple_int 1855 nsize: {{2 4 5}} 1856 args: -interpolate after_distribute 1857 test: # seq load, simple partitioner, par interpolation 1858 # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy(). 1859 # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken. 1860 # We can then provide an intentionally broken mesh instead. 1861 TODO: This test is broken because PointSF is fixed. 1862 suffix: 9_seq_hdf5_simple_int_err 1863 nsize: 4 1864 args: -interpolate after_distribute 1865 filter: sed -e "/PETSC ERROR/,$$d" 1866 test: # seq load, metis partitioner, par interpolation 1867 nsize: {{2 4 5}} 1868 suffix: 9_seq_hdf5_metis_int 1869 requires: parmetis 1870 args: -interpolate after_distribute 1871 args: -petscpartitioner_type parmetis 1872 1873 testset: # 9 HDF5 PARALLEL LOAD 1874 requires: hdf5 !complex datafilespath 1875 nsize: {{2 4 5}} 1876 args: -dm_plex_check_all 1877 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 1878 test: # par load 1879 suffix: 9_par_hdf5 1880 args: -interpolate none 1881 test: # par load, par interpolation 1882 suffix: 9_par_hdf5_int 1883 args: -interpolate after_create 1884 test: # par load, parmetis repartitioner 1885 TODO: Parallel partitioning of uninterpolated meshes not supported 1886 suffix: 9_par_hdf5_parmetis 1887 requires: parmetis 1888 args: -distribute -petscpartitioner_type parmetis 1889 args: -interpolate none 1890 test: # par load, par interpolation, parmetis repartitioner 1891 suffix: 9_par_hdf5_int_parmetis 1892 requires: parmetis 1893 args: -distribute -petscpartitioner_type parmetis 1894 args: -interpolate after_create 1895 test: # par load, parmetis partitioner, par interpolation 1896 TODO: Parallel partitioning of uninterpolated meshes not supported 1897 suffix: 9_par_hdf5_parmetis_int 1898 requires: parmetis 1899 args: -distribute -petscpartitioner_type parmetis 1900 args: -interpolate after_distribute 1901 1902 testset: # 10 HDF5 PARALLEL LOAD 1903 requires: hdf5 !complex datafilespath 1904 nsize: {{2 4 7}} 1905 args: -dm_plex_check_all 1906 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 1907 test: # par load, par interpolation 1908 suffix: 10_par_hdf5_int 1909 args: -interpolate after_create 1910 TEST*/ 1911