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