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