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 static PetscLogStage stage[3]; 221 222 static PetscErrorCode DMPlexCheckPointSFHeavy(DM, PortableBoundary); 223 static PetscErrorCode DMPlexSetOrientInterface_Private(DM,PetscBool); 224 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM, PortableBoundary *); 225 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM, IS, PetscSection, IS *); 226 227 static PetscErrorCode PortableBoundaryDestroy(PortableBoundary *bnd) 228 { 229 PetscInt d; 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 if (!*bnd) PetscFunctionReturn(0); 234 ierr = VecDestroy(&(*bnd)->coordinates);CHKERRQ(ierr); 235 for (d=0; d < (*bnd)->depth; d++) { 236 ierr = PetscSectionDestroy(&(*bnd)->sections[d]);CHKERRQ(ierr); 237 } 238 ierr = PetscFree((*bnd)->sections);CHKERRQ(ierr); 239 ierr = PetscFree(*bnd);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 244 { 245 const char *interpTypes[4] = {"none", "create", "after_create", "after_distribute"}; 246 PetscInt interp=NONE, dim; 247 PetscBool flg1, flg2; 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 options->debug = 0; 252 options->testNum = 0; 253 options->dim = 2; 254 options->cellSimplex = PETSC_TRUE; 255 options->distribute = PETSC_FALSE; 256 options->interpolate = NONE; 257 options->useGenerator = PETSC_FALSE; 258 options->testOrientIF = PETSC_FALSE; 259 options->testHeavy = PETSC_TRUE; 260 options->customView = PETSC_FALSE; 261 options->testExpandPointsEmpty = PETSC_FALSE; 262 options->ornt[0] = 0; 263 options->ornt[1] = 0; 264 options->faces[0] = 2; 265 options->faces[1] = 2; 266 options->faces[2] = 2; 267 options->filename[0] = '\0'; 268 options->coordsTol = PETSC_DEFAULT; 269 270 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 271 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 272 ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 273 ierr = PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 274 ierr = PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL);CHKERRQ(ierr); 275 ierr = PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL);CHKERRQ(ierr); 276 options->interpolate = (InterpType) interp; 277 if (!options->distribute && options->interpolate == AFTER_DISTRIBUTE) SETERRQ(comm, PETSC_ERR_SUP, "-interpolate after_distribute needs -distribute 1"); 278 ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr); 279 options->ncoords = 128; 280 ierr = PetscOptionsRealArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL);CHKERRQ(ierr); 281 ierr = PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL);CHKERRQ(ierr); 282 options->nPointsToExpand = 128; 283 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); 284 if (options->nPointsToExpand) { 285 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); 286 } 287 ierr = PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL);CHKERRQ(ierr); 288 ierr = PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL);CHKERRQ(ierr); 289 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3);CHKERRQ(ierr); 290 dim = 3; 291 ierr = PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2);CHKERRQ(ierr); 292 if (flg2) { 293 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); 294 options->dim = dim; 295 } 296 ierr = PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 297 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); 298 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); 299 if (flg2 != options->testOrientIF) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set"); 300 if (options->testOrientIF) { 301 PetscInt i; 302 for (i=0; i<2; i++) { 303 if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10); /* 11 12 13 become -1 -2 -3 */ 304 } 305 options->filename[0] = 0; 306 options->useGenerator = PETSC_FALSE; 307 options->dim = 3; 308 options->cellSimplex = PETSC_TRUE; 309 options->interpolate = CREATE; 310 options->distribute = PETSC_FALSE; 311 } 312 ierr = PetscOptionsEnd(); 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 317 { 318 PetscInt testNum = user->testNum; 319 PetscMPIInt rank,size; 320 PetscErrorCode ierr; 321 PetscInt numCorners=2,i; 322 PetscInt numCells,numVertices,network; 323 PetscInt *cells; 324 PetscReal *coords; 325 326 PetscFunctionBegin; 327 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 328 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 329 if (size > 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum); 330 331 numCells = 3; 332 ierr = PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL);CHKERRQ(ierr); 333 if (numCells < 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells); 334 335 if (size == 1) { 336 PetscReal *dcoords; 337 numVertices = numCells + 1; 338 ierr = PetscMalloc2(2*numCells,&cells,2*numVertices,&dcoords);CHKERRQ(ierr); 339 for (i=0; i<numCells; i++) { 340 cells[2*i] = i; cells[2*i+1] = i + 1; 341 dcoords[2*i] = i; dcoords[2*i+1] = i + 1; 342 } 343 344 ierr = DMPlexCreateFromCellListPetsc(comm, user->dim, numCells, numVertices, numCorners, PETSC_FALSE, cells, user->dim, dcoords, dm);CHKERRQ(ierr); 345 ierr = PetscFree2(cells,coords);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 network = 0; 350 ierr = PetscOptionsGetInt(NULL, NULL, "-network_case", &network, NULL);CHKERRQ(ierr); 351 if (network == 0) { 352 switch (rank) { 353 case 0: 354 { 355 numCells = 2; 356 numVertices = numCells; 357 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 358 cells[0] = 0; cells[1] = 1; 359 cells[2] = 1; cells[3] = 2; 360 coords[0] = 0.; coords[1] = 1.; 361 coords[2] = 1.; coords[3] = 2.; 362 } 363 break; 364 case 1: 365 { 366 numCells -= 2; 367 numVertices = numCells + 1; 368 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 369 for (i=0; i<numCells; i++) { 370 cells[2*i] = 2+i; cells[2*i+1] = 2 + i + 1; 371 coords[2*i] = 2+i; coords[2*i+1] = 2 + i + 1; 372 } 373 } 374 break; 375 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 376 } 377 } else { /* network_case = 1 */ 378 /* ----------------------- */ 379 switch (rank) { 380 case 0: 381 { 382 numCells = 2; 383 numVertices = 3; 384 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 385 cells[0] = 0; cells[1] = 3; 386 cells[2] = 3; cells[3] = 1; 387 } 388 break; 389 case 1: 390 { 391 numCells = 1; 392 numVertices = 1; 393 ierr = PetscMalloc2(2*numCells,&cells,2*numCells,&coords);CHKERRQ(ierr); 394 cells[0] = 3; cells[1] = 2; 395 } 396 break; 397 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 398 } 399 } 400 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, PETSC_FALSE, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 401 ierr = PetscFree2(cells,coords);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 406 { 407 PetscInt testNum = user->testNum, p; 408 PetscMPIInt rank, size; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 413 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 414 switch (testNum) { 415 case 0: 416 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 417 switch (rank) { 418 case 0: 419 { 420 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 421 const PetscInt cells[3] = {0, 1, 2}; 422 PetscReal coords[4] = {-0.5, 0.5, 0.0, 0.0}; 423 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 424 425 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 426 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 427 } 428 break; 429 case 1: 430 { 431 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 432 const PetscInt cells[3] = {1, 3, 2}; 433 PetscReal coords[4] = {0.0, 1.0, 0.5, 0.5}; 434 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 435 436 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 437 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 438 } 439 break; 440 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 441 } 442 break; 443 case 1: 444 if (size != 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); 445 switch (rank) { 446 case 0: 447 { 448 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 449 const PetscInt cells[3] = {0, 1, 2}; 450 PetscReal coords[4] = {0.0, 1.0, 0.0, 0.0}; 451 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 452 453 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 454 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 455 } 456 break; 457 case 1: 458 { 459 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 460 const PetscInt cells[3] = {0, 2, 3}; 461 PetscReal coords[4] = {0.5, 0.5, 1.0, 1.0}; 462 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 463 464 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 465 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 466 } 467 break; 468 case 2: 469 { 470 const PetscInt numCells = 2, numVertices = 1, numCorners = 3; 471 const PetscInt cells[6] = {2, 4, 3, 2, 1, 4}; 472 PetscReal coords[2] = {1.0, 0.0}; 473 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 474 475 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 476 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 477 } 478 break; 479 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 480 } 481 break; 482 case 2: 483 if (size != 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); 484 switch (rank) { 485 case 0: 486 { 487 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 488 const PetscInt cells[3] = {1, 2, 0}; 489 PetscReal coords[4] = {0.5, 0.5, 0.0, 1.0}; 490 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 491 492 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 493 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 494 } 495 break; 496 case 1: 497 { 498 const PetscInt numCells = 1, numVertices = 2, numCorners = 3; 499 const PetscInt cells[3] = {1, 0, 3}; 500 PetscReal coords[4] = {0.0, 0.0, 1.0, 1.0}; 501 PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; 502 503 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 504 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 505 } 506 break; 507 case 2: 508 { 509 const PetscInt numCells = 2, numVertices = 1, numCorners = 3; 510 const PetscInt cells[6] = {0, 4, 3, 0, 2, 4}; 511 PetscReal coords[2] = {1.0, 0.0}; 512 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 513 514 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 515 for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 516 } 517 break; 518 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 519 } 520 break; 521 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 522 } 523 PetscFunctionReturn(0); 524 } 525 526 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 527 { 528 PetscInt testNum = user->testNum, p; 529 PetscMPIInt rank, size; 530 PetscErrorCode ierr; 531 532 PetscFunctionBegin; 533 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 534 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 535 switch (testNum) { 536 case 0: 537 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 538 switch (rank) { 539 case 0: 540 { 541 const PetscInt numCells = 1, numVertices = 2, numCorners = 4; 542 const PetscInt cells[4] = {0, 2, 1, 3}; 543 PetscReal coords[6] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0}; 544 PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; 545 546 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 547 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 548 } 549 break; 550 case 1: 551 { 552 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 553 const PetscInt cells[4] = {1, 2, 4, 3}; 554 PetscReal coords[9] = {1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; 555 PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; 556 557 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 558 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 559 } 560 break; 561 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 562 } 563 break; 564 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 565 } 566 if (user->testOrientIF) { 567 PetscInt start; 568 PetscBool reverse; 569 PetscInt ifp[] = {8, 6}; 570 571 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh before orientation");CHKERRQ(ierr); 572 ierr = DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view");CHKERRQ(ierr); 573 /* rotate interface face ifp[rank] by given orientation ornt[rank] */ 574 ierr = DMPlexFixFaceOrientations_Translate_Private(user->ornt[rank], &start, &reverse);CHKERRQ(ierr); 575 ierr = DMPlexOrientCell_Internal(*dm, ifp[rank], start, reverse);CHKERRQ(ierr); 576 ierr = DMPlexCheckFaces(*dm, 0);CHKERRQ(ierr); 577 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 578 ierr = PetscPrintf(comm, "Orientation test PASSED\n");CHKERRQ(ierr); 579 } 580 PetscFunctionReturn(0); 581 } 582 583 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 584 { 585 PetscInt testNum = user->testNum, p; 586 PetscMPIInt rank, size; 587 PetscErrorCode ierr; 588 589 PetscFunctionBegin; 590 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 591 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 592 switch (testNum) { 593 case 0: 594 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 595 switch (rank) { 596 case 0: 597 { 598 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 599 const PetscInt cells[4] = {0, 1, 2, 3}; 600 PetscReal coords[6] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0}; 601 PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; 602 603 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 604 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 605 } 606 break; 607 case 1: 608 { 609 const PetscInt numCells = 1, numVertices = 3, numCorners = 4; 610 const PetscInt cells[4] = {1, 4, 5, 2}; 611 PetscReal coords[6] = {-0.5, 1.0, 0.5, 0.0, 0.5, 1.0}; 612 PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; 613 614 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 615 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 616 } 617 break; 618 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 619 } 620 break; 621 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 622 } 623 PetscFunctionReturn(0); 624 } 625 626 static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 627 { 628 PetscInt testNum = user->testNum, p; 629 PetscMPIInt rank, size; 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 634 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 635 switch (testNum) { 636 case 0: 637 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); 638 switch (rank) { 639 case 0: 640 { 641 const PetscInt numCells = 1, numVertices = 6, numCorners = 8; 642 const PetscInt cells[8] = {0, 3, 2, 1, 4, 5, 6, 7}; 643 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}; 644 PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 645 646 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 647 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 648 } 649 break; 650 case 1: 651 { 652 const PetscInt numCells = 1, numVertices = 6, numCorners = 8; 653 const PetscInt cells[8] = {1, 2, 9, 8, 5, 10, 11, 6}; 654 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}; 655 PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 656 657 ierr = DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, dm);CHKERRQ(ierr); 658 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 659 } 660 break; 661 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); 662 } 663 break; 664 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); 665 } 666 PetscFunctionReturn(0); 667 } 668 669 static PetscErrorCode CustomView(DM dm, PetscViewer v) 670 { 671 DMPlexInterpolatedFlag interpolated; 672 PetscBool distributed; 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 ierr = DMPlexIsDistributed(dm, &distributed);CHKERRQ(ierr); 677 ierr = DMPlexIsInterpolatedCollective(dm, &interpolated);CHKERRQ(ierr); 678 ierr = PetscViewerASCIIPrintf(v, "DMPlexIsDistributed: %s\n", PetscBools[distributed]);CHKERRQ(ierr); 679 ierr = PetscViewerASCIIPrintf(v, "DMPlexIsInterpolatedCollective: %s\n", DMPlexInterpolatedFlags[interpolated]);CHKERRQ(ierr); 680 PetscFunctionReturn(0); 681 } 682 683 static PetscErrorCode CreateMeshFromFile(MPI_Comm comm, AppCtx *user, DM *dm, DM *serialDM) 684 { 685 const char *filename = user->filename; 686 PetscBool testHeavy = user->testHeavy; 687 PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; 688 PetscBool distributed = PETSC_FALSE; 689 PetscErrorCode ierr; 690 691 PetscFunctionBegin; 692 *serialDM = NULL; 693 if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_FALSE);CHKERRQ(ierr);} 694 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 695 ierr = DMPlexCreateFromFile(comm, filename, interpCreate, dm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 696 ierr = PetscLogStagePop();CHKERRQ(ierr); 697 if (testHeavy && interpCreate) {ierr = DMPlexSetOrientInterface_Private(NULL, PETSC_TRUE);CHKERRQ(ierr);} 698 ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr); 699 ierr = PetscPrintf(comm, "DMPlexCreateFromFile produced %s mesh.\n", distributed ? "distributed" : "serial");CHKERRQ(ierr); 700 if (testHeavy && distributed) { 701 ierr = PetscOptionsSetValue(NULL, "-dm_plex_hdf5_force_sequential", NULL);CHKERRQ(ierr); 702 ierr = DMPlexCreateFromFile(comm, filename, interpCreate, serialDM);CHKERRQ(ierr); 703 ierr = DMPlexIsDistributed(*serialDM, &distributed);CHKERRQ(ierr); 704 if (distributed) SETERRQ(comm, PETSC_ERR_PLIB, "unable to create a serial DM from file"); 705 } 706 ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 710 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 711 { 712 PetscPartitioner part; 713 PortableBoundary boundary = NULL; 714 DM serialDM = NULL; 715 PetscBool cellSimplex = user->cellSimplex; 716 PetscBool useGenerator = user->useGenerator; 717 PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; 718 PetscBool interpSerial = user->interpolate == AFTER_CREATE ? PETSC_TRUE : PETSC_FALSE; 719 PetscBool interpParallel = user->interpolate == AFTER_DISTRIBUTE ? PETSC_TRUE : PETSC_FALSE; 720 PetscBool testHeavy = user->testHeavy; 721 PetscMPIInt rank; 722 PetscErrorCode ierr; 723 724 PetscFunctionBegin; 725 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 726 if (user->filename[0]) { 727 ierr = CreateMeshFromFile(comm, user, dm, &serialDM);CHKERRQ(ierr); 728 } else if (useGenerator) { 729 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 730 ierr = DMPlexCreateBoxMesh(comm, user->dim, cellSimplex, user->faces, NULL, NULL, NULL, interpCreate, dm);CHKERRQ(ierr); 731 ierr = PetscLogStagePop();CHKERRQ(ierr); 732 } else { 733 ierr = PetscLogStagePush(stage[0]);CHKERRQ(ierr); 734 switch (user->dim) { 735 case 1: 736 ierr = CreateMesh_1D(comm, interpCreate, user, dm);CHKERRQ(ierr); 737 break; 738 case 2: 739 if (cellSimplex) { 740 ierr = CreateSimplex_2D(comm, interpCreate, user, dm);CHKERRQ(ierr); 741 } else { 742 ierr = CreateQuad_2D(comm, interpCreate, user, dm);CHKERRQ(ierr); 743 } 744 break; 745 case 3: 746 if (cellSimplex) { 747 ierr = CreateSimplex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr); 748 } else { 749 ierr = CreateHex_3D(comm, interpCreate, user, dm);CHKERRQ(ierr); 750 } 751 break; 752 default: 753 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", user->dim); 754 } 755 ierr = PetscLogStagePop();CHKERRQ(ierr); 756 } 757 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); 758 ierr = PetscObjectSetName((PetscObject) *dm, "Original Mesh");CHKERRQ(ierr); 759 ierr = DMViewFromOptions(*dm, NULL, "-orig_dm_view");CHKERRQ(ierr); 760 761 if (interpSerial) { 762 DM idm; 763 764 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);} 765 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 766 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 767 ierr = PetscLogStagePop();CHKERRQ(ierr); 768 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);} 769 ierr = DMDestroy(dm);CHKERRQ(ierr); 770 *dm = idm; 771 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh");CHKERRQ(ierr); 772 ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr); 773 } 774 775 /* Set partitioner options */ 776 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 777 if (part) { 778 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE);CHKERRQ(ierr); 779 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 780 } 781 782 if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 783 if (testHeavy) { 784 PetscBool distributed; 785 786 ierr = DMPlexIsDistributed(*dm, &distributed);CHKERRQ(ierr); 787 if (!serialDM && !distributed) { 788 serialDM = *dm; 789 ierr = PetscObjectReference((PetscObject)*dm);CHKERRQ(ierr); 790 } 791 if (serialDM) { 792 ierr = DMPlexGetExpandedBoundary_Private(serialDM, &boundary);CHKERRQ(ierr); 793 } 794 if (boundary) { 795 /* check DM which has been created in parallel and already interpolated */ 796 ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr); 797 } 798 /* Orient interface because it could be deliberately skipped above. It is idempotent. */ 799 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 800 } 801 if (user->distribute) { 802 DM pdm = NULL; 803 804 /* Redistribute mesh over processes using that partitioner */ 805 ierr = PetscLogStagePush(stage[1]);CHKERRQ(ierr); 806 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 807 ierr = PetscLogStagePop();CHKERRQ(ierr); 808 if (pdm) { 809 ierr = DMDestroy(dm);CHKERRQ(ierr); 810 *dm = pdm; 811 ierr = PetscObjectSetName((PetscObject) *dm, "Redistributed Mesh");CHKERRQ(ierr); 812 ierr = DMViewFromOptions(*dm, NULL, "-dist_dm_view");CHKERRQ(ierr); 813 } 814 815 if (interpParallel) { 816 DM idm; 817 818 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE);CHKERRQ(ierr);} 819 ierr = PetscLogStagePush(stage[2]);CHKERRQ(ierr); 820 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ 821 ierr = PetscLogStagePop();CHKERRQ(ierr); 822 if (testHeavy) {ierr = DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE);CHKERRQ(ierr);} 823 ierr = DMDestroy(dm);CHKERRQ(ierr); 824 *dm = idm; 825 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Redistributed Mesh");CHKERRQ(ierr); 826 ierr = DMViewFromOptions(*dm, NULL, "-intp_dm_view");CHKERRQ(ierr); 827 } 828 } 829 if (testHeavy) { 830 if (boundary) { 831 ierr = DMPlexCheckPointSFHeavy(*dm, boundary);CHKERRQ(ierr); 832 } 833 /* Orient interface because it could be deliberately skipped above. It is idempotent. */ 834 ierr = DMPlexOrientInterface_Internal(*dm);CHKERRQ(ierr); 835 } 836 837 ierr = PetscObjectSetName((PetscObject) *dm, "Parallel Mesh");CHKERRQ(ierr); 838 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 839 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 840 841 if (user->customView) {ierr = CustomView(*dm, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 842 ierr = DMDestroy(&serialDM);CHKERRQ(ierr); 843 ierr = PortableBoundaryDestroy(&boundary);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 PETSC_STATIC_INLINE PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, PetscReal *coords, PetscReal tol) 848 { 849 PetscErrorCode ierr; 850 851 PetscFunctionBegin; 852 if (dim > 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3"); 853 if (tol >= 1e-3) { 854 switch (dim) { 855 case 1: ierr = PetscSNPrintf(buf,len,"(%12.3f)",(double)coords[0]);CHKERRQ(ierr); 856 case 2: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",(double)coords[0],(double)coords[1]);CHKERRQ(ierr); 857 case 3: ierr = PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",(double)coords[0],(double)coords[1],(double)coords[2]);CHKERRQ(ierr); 858 } 859 } else { 860 switch (dim) { 861 case 1: ierr = PetscSNPrintf(buf,len,"(%12.6f)",(double)coords[0]);CHKERRQ(ierr); 862 case 2: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",(double)coords[0],(double)coords[1]);CHKERRQ(ierr); 863 case 3: ierr = PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",(double)coords[0],(double)coords[1],(double)coords[2]);CHKERRQ(ierr); 864 } 865 } 866 PetscFunctionReturn(0); 867 } 868 869 static PetscErrorCode ViewVerticesFromCoords(DM dm, PetscInt npoints, PetscReal coords[], PetscReal tol, PetscViewer viewer) 870 { 871 PetscInt dim, i; 872 PetscInt *points; 873 char coordstr[128]; 874 MPI_Comm comm; 875 PetscMPIInt rank; 876 PetscErrorCode ierr; 877 878 PetscFunctionBegin; 879 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 880 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 881 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 882 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 883 ierr = PetscMalloc1(npoints, &points);CHKERRQ(ierr); 884 ierr = DMPlexFindVertices(dm, npoints, coords, tol, points);CHKERRQ(ierr); 885 for (i=0; i < npoints; i++) { 886 ierr = coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol);CHKERRQ(ierr); 887 if (!rank && i) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, "-----\n");CHKERRQ(ierr);} 888 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i]);CHKERRQ(ierr); 889 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 890 } 891 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 892 ierr = PetscFree(points);CHKERRQ(ierr); 893 PetscFunctionReturn(0); 894 } 895 896 static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user) 897 { 898 IS is; 899 PetscSection *sects; 900 IS *iss; 901 PetscInt d,depth; 902 PetscMPIInt rank; 903 PetscErrorCode ierr; 904 PetscViewer viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer; 905 906 PetscFunctionBegin; 907 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 908 if (user->testExpandPointsEmpty && !rank) { 909 ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 910 } else { 911 ierr = ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is);CHKERRQ(ierr); 912 } 913 ierr = DMPlexGetConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 914 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 915 ierr = PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank);CHKERRQ(ierr); 916 for (d=depth-1; d>=0; d--) { 917 IS checkIS; 918 PetscBool flg; 919 920 ierr = PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d);CHKERRQ(ierr); 921 ierr = PetscSectionView(sects[d], sviewer);CHKERRQ(ierr); 922 ierr = ISView(iss[d], sviewer);CHKERRQ(ierr); 923 /* check reverse operation */ 924 if (d < depth-1) { 925 ierr = DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS);CHKERRQ(ierr); 926 ierr = ISEqualUnsorted(checkIS, iss[d+1], &flg);CHKERRQ(ierr); 927 if (!flg) SETERRQ(PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS"); 928 ierr = ISDestroy(&checkIS);CHKERRQ(ierr); 929 } 930 } 931 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 932 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 933 ierr = DMPlexRestoreConeRecursive(dm, is, &depth, &iss, §s);CHKERRQ(ierr); 934 ierr = ISDestroy(&is);CHKERRQ(ierr); 935 PetscFunctionReturn(0); 936 } 937 938 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis) 939 { 940 PetscInt n,n1,ncone,numCoveredPoints,o,p,q,start,end; 941 const PetscInt *coveredPoints; 942 const PetscInt *arr, *cone; 943 PetscInt *newarr; 944 PetscErrorCode ierr; 945 946 PetscFunctionBegin; 947 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 948 ierr = PetscSectionGetStorageSize(section, &n1);CHKERRQ(ierr); 949 ierr = PetscSectionGetChart(section, &start, &end);CHKERRQ(ierr); 950 if (n != n1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size\n", n, n1); 951 ierr = ISGetIndices(is, &arr);CHKERRQ(ierr); 952 ierr = PetscMalloc1(end-start, &newarr);CHKERRQ(ierr); 953 for (q=start; q<end; q++) { 954 ierr = PetscSectionGetDof(section, q, &ncone);CHKERRQ(ierr); 955 ierr = PetscSectionGetOffset(section, q, &o);CHKERRQ(ierr); 956 cone = &arr[o]; 957 if (ncone == 1) { 958 numCoveredPoints = 1; 959 p = cone[0]; 960 } else { 961 PetscInt i; 962 p = PETSC_MAX_INT; 963 for (i=0; i<ncone; i++) if (cone[i] < 0) {p = -1; break;} 964 if (p >= 0) { 965 ierr = DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 966 if (numCoveredPoints > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q); 967 if (numCoveredPoints) p = coveredPoints[0]; 968 else p = -2; 969 ierr = DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 970 } 971 } 972 newarr[q-start] = p; 973 } 974 ierr = ISRestoreIndices(is, &arr);CHKERRQ(ierr); 975 ierr = ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis);CHKERRQ(ierr); 976 PetscFunctionReturn(0); 977 } 978 979 static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is) 980 { 981 PetscInt d; 982 IS is,newis; 983 PetscErrorCode ierr; 984 985 PetscFunctionBegin; 986 is = boundary_expanded_is; 987 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 988 for (d = 0; d < depth-1; ++d) { 989 ierr = DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis);CHKERRQ(ierr); 990 ierr = ISDestroy(&is);CHKERRQ(ierr); 991 is = newis; 992 } 993 *boundary_is = is; 994 PetscFunctionReturn(0); 995 } 996 997 #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; } 998 999 static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm) 1000 { 1001 PetscErrorCode ierr; 1002 PetscViewer viewer; 1003 PetscBool flg; 1004 static PetscBool incall = PETSC_FALSE; 1005 PetscViewerFormat format; 1006 1007 PetscFunctionBegin; 1008 if (incall) PetscFunctionReturn(0); 1009 incall = PETSC_TRUE; 1010 ierr = PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg);CHKERRQI(incall,ierr); 1011 if (flg) { 1012 ierr = PetscViewerPushFormat(viewer,format);CHKERRQI(incall,ierr); 1013 ierr = DMLabelView(label, viewer);CHKERRQI(incall,ierr); 1014 ierr = PetscViewerPopFormat(viewer);CHKERRQI(incall,ierr); 1015 ierr = PetscViewerDestroy(&viewer);CHKERRQI(incall,ierr); 1016 } 1017 incall = PETSC_FALSE; 1018 PetscFunctionReturn(0); 1019 } 1020 1021 /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */ 1022 PETSC_STATIC_INLINE PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is) 1023 { 1024 IS tmpis; 1025 PetscErrorCode ierr; 1026 1027 PetscFunctionBegin; 1028 ierr = DMLabelGetStratumIS(label, value, &tmpis);CHKERRQ(ierr); 1029 if (!tmpis) {ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis);CHKERRQ(ierr);} 1030 ierr = ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is);CHKERRQ(ierr); 1031 ierr = ISDestroy(&tmpis);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 /* currently only for simple PetscSection without fields or constraints */ 1036 static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout) 1037 { 1038 PetscSection sec; 1039 PetscInt chart[2], p; 1040 PetscInt *dofarr; 1041 PetscMPIInt rank; 1042 PetscErrorCode ierr; 1043 1044 PetscFunctionBegin; 1045 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1046 if (rank == rootrank) { 1047 ierr = PetscSectionGetChart(sec0, &chart[0], &chart[1]);CHKERRQ(ierr); 1048 } 1049 ierr = MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1050 ierr = PetscMalloc1(chart[1]-chart[0], &dofarr);CHKERRQ(ierr); 1051 if (rank == rootrank) { 1052 for (p = chart[0]; p < chart[1]; p++) { 1053 ierr = PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]]);CHKERRQ(ierr); 1054 } 1055 } 1056 ierr = MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1057 ierr = PetscSectionCreate(comm, &sec);CHKERRQ(ierr); 1058 ierr = PetscSectionSetChart(sec, chart[0], chart[1]);CHKERRQ(ierr); 1059 for (p = chart[0]; p < chart[1]; p++) { 1060 ierr = PetscSectionSetDof(sec, p, dofarr[p-chart[0]]);CHKERRQ(ierr); 1061 } 1062 ierr = PetscSectionSetUp(sec);CHKERRQ(ierr); 1063 ierr = PetscFree(dofarr);CHKERRQ(ierr); 1064 *secout = sec; 1065 PetscFunctionReturn(0); 1066 } 1067 1068 static PetscErrorCode VecToPetscReal_Private(Vec vec, PetscReal *rvals[]) 1069 { 1070 PetscInt n; 1071 const PetscScalar *svals; 1072 PetscErrorCode ierr; 1073 1074 PetscFunctionBegin; 1075 ierr = VecGetLocalSize(vec, &n);CHKERRQ(ierr); 1076 ierr = VecGetArrayRead(vec, &svals);CHKERRQ(ierr); 1077 ierr = PetscMalloc1(n, rvals);CHKERRQ(ierr); 1078 #if defined(PETSC_USE_COMPLEX) 1079 { 1080 PetscInt i; 1081 for (i=0; i<n; i++) (*rvals)[i] = PetscRealPart(svals[i]); 1082 } 1083 #else 1084 ierr = PetscMemcpy(*rvals, svals, n*sizeof(PetscReal)); 1085 #endif 1086 ierr = VecRestoreArrayRead(vec, &svals);CHKERRQ(ierr); 1087 PetscFunctionReturn(0); 1088 } 1089 1090 static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is) 1091 { 1092 PetscInt dim, ncoords, npoints; 1093 PetscReal *rcoords; 1094 PetscInt *points; 1095 IS faces_expanded_is; 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 ierr = DMGetDimension(ipdm, &dim);CHKERRQ(ierr); 1100 ierr = VecGetLocalSize(bnd->coordinates, &ncoords);CHKERRQ(ierr); 1101 ierr = VecToPetscReal_Private(bnd->coordinates, &rcoords);CHKERRQ(ierr); 1102 npoints = ncoords / dim; 1103 ierr = PetscMalloc1(npoints, &points);CHKERRQ(ierr); 1104 ierr = DMPlexFindVertices(ipdm, npoints, rcoords, 0.0, points);CHKERRQ(ierr); 1105 ierr = ISCreateGeneral(PETSC_COMM_SELF, npoints, points, PETSC_OWN_POINTER, &faces_expanded_is);CHKERRQ(ierr); 1106 ierr = DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is);CHKERRQ(ierr); 1107 ierr = PetscFree(rcoords);CHKERRQ(ierr); 1108 ierr = ISDestroy(&faces_expanded_is);CHKERRQ(ierr); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */ 1113 static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable) 1114 { 1115 PetscOptions options = NULL; 1116 const char *prefix = NULL; 1117 const char opt[] = "-dm_plex_interpolate_orient_interfaces"; 1118 char prefix_opt[512]; 1119 PetscBool flg, set; 1120 static PetscBool wasSetTrue = PETSC_FALSE; 1121 PetscErrorCode ierr; 1122 1123 PetscFunctionBegin; 1124 if (dm) { 1125 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);CHKERRQ(ierr); 1126 options = ((PetscObject)dm)->options; 1127 } 1128 ierr = PetscStrcpy(prefix_opt, "-");CHKERRQ(ierr); 1129 ierr = PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt));CHKERRQ(ierr); 1130 ierr = PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt));CHKERRQ(ierr); 1131 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1132 if (!enable) { 1133 if (set && flg) wasSetTrue = PETSC_TRUE; 1134 ierr = PetscOptionsSetValue(options, prefix_opt, "0");CHKERRQ(ierr); 1135 } else if (set && !flg) { 1136 if (wasSetTrue) { 1137 ierr = PetscOptionsSetValue(options, prefix_opt, "1");CHKERRQ(ierr); 1138 } else { 1139 /* default is PETSC_TRUE */ 1140 ierr = PetscOptionsClearValue(options, prefix_opt);CHKERRQ(ierr); 1141 } 1142 wasSetTrue = PETSC_FALSE; 1143 } 1144 if (PetscDefined(USE_DEBUG)) { 1145 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1146 if (PetscUnlikely(set && flg != enable)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect"); 1147 } 1148 PetscFunctionReturn(0); 1149 } 1150 1151 /* get coordinate description of the whole-domain boundary */ 1152 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary) 1153 { 1154 PortableBoundary bnd0, bnd; 1155 MPI_Comm comm; 1156 DM idm; 1157 DMLabel label; 1158 PetscInt d; 1159 const char boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary"; 1160 IS boundary_is; 1161 IS *boundary_expanded_iss; 1162 PetscMPIInt rootrank = 0; 1163 PetscMPIInt rank, size; 1164 PetscInt value = 1; 1165 DMPlexInterpolatedFlag intp; 1166 PetscBool flg; 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = PetscNew(&bnd);CHKERRQ(ierr); 1171 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1172 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1173 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1174 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 1175 if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed"); 1176 1177 /* interpolate serial DM if not yet interpolated */ 1178 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1179 if (intp == DMPLEX_INTERPOLATED_FULL) { 1180 idm = dm; 1181 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 1182 } else { 1183 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 1184 ierr = DMViewFromOptions(idm, NULL, "-idm_view");CHKERRQ(ierr); 1185 } 1186 1187 /* mark whole-domain boundary of the serial DM */ 1188 ierr = DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label);CHKERRQ(ierr); 1189 ierr = DMAddLabel(idm, label);CHKERRQ(ierr); 1190 ierr = DMPlexMarkBoundaryFaces(idm, value, label);CHKERRQ(ierr); 1191 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm);CHKERRQ(ierr); 1192 ierr = DMLabelGetStratumIS(label, value, &boundary_is);CHKERRQ(ierr); 1193 1194 /* translate to coordinates */ 1195 ierr = PetscNew(&bnd0);CHKERRQ(ierr); 1196 ierr = DMGetCoordinatesLocalSetUp(idm);CHKERRQ(ierr); 1197 if (rank == rootrank) { 1198 ierr = DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1199 ierr = DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates);CHKERRQ(ierr); 1200 /* self-check */ 1201 { 1202 IS is0; 1203 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0);CHKERRQ(ierr); 1204 ierr = ISEqual(is0, boundary_is, &flg);CHKERRQ(ierr); 1205 if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS"); 1206 ierr = ISDestroy(&is0);CHKERRQ(ierr); 1207 } 1208 } else { 1209 ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates);CHKERRQ(ierr); 1210 } 1211 1212 { 1213 Vec tmp; 1214 VecScatter sc; 1215 IS xis; 1216 PetscInt n; 1217 1218 /* just convert seq vectors to mpi vector */ 1219 ierr = VecGetLocalSize(bnd0->coordinates, &n);CHKERRQ(ierr); 1220 ierr = MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1221 if (rank == rootrank) { 1222 ierr = VecCreateMPI(comm, n, n, &tmp);CHKERRQ(ierr); 1223 } else { 1224 ierr = VecCreateMPI(comm, 0, n, &tmp);CHKERRQ(ierr); 1225 } 1226 ierr = VecCopy(bnd0->coordinates, tmp);CHKERRQ(ierr); 1227 ierr = VecDestroy(&bnd0->coordinates);CHKERRQ(ierr); 1228 bnd0->coordinates = tmp; 1229 1230 /* replicate coordinates from root rank to all ranks */ 1231 ierr = VecCreateMPI(comm, n, n*size, &bnd->coordinates);CHKERRQ(ierr); 1232 ierr = ISCreateStride(comm, n, 0, 1, &xis);CHKERRQ(ierr); 1233 ierr = VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc);CHKERRQ(ierr); 1234 ierr = VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1235 ierr = VecScatterEnd( sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1236 ierr = VecScatterDestroy(&sc);CHKERRQ(ierr); 1237 ierr = ISDestroy(&xis);CHKERRQ(ierr); 1238 } 1239 bnd->depth = bnd0->depth; 1240 ierr = MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1241 ierr = PetscMalloc1(bnd->depth, &bnd->sections);CHKERRQ(ierr); 1242 for (d=0; d<bnd->depth; d++) { 1243 ierr = PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]);CHKERRQ(ierr); 1244 } 1245 1246 if (rank == rootrank) { 1247 ierr = DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1248 } 1249 ierr = PortableBoundaryDestroy(&bnd0);CHKERRQ(ierr); 1250 ierr = DMRemoveLabelBySelf(idm, &label, PETSC_TRUE);CHKERRQ(ierr); 1251 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1252 ierr = ISDestroy(&boundary_is);CHKERRQ(ierr); 1253 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1254 *boundary = bnd; 1255 PetscFunctionReturn(0); 1256 } 1257 1258 /* get faces of inter-partition interface */ 1259 static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is) 1260 { 1261 MPI_Comm comm; 1262 DMLabel label; 1263 IS part_boundary_faces_is; 1264 const char partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary"; 1265 PetscInt value = 1; 1266 DMPlexInterpolatedFlag intp; 1267 PetscErrorCode ierr; 1268 1269 PetscFunctionBegin; 1270 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1271 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1272 if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1273 1274 /* get ipdm partition boundary (partBoundary) */ 1275 ierr = DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label);CHKERRQ(ierr); 1276 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1277 ierr = DMPlexMarkBoundaryFaces(ipdm, value, label);CHKERRQ(ierr); 1278 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm);CHKERRQ(ierr); 1279 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is);CHKERRQ(ierr); 1280 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1281 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1282 1283 /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */ 1284 ierr = ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is);CHKERRQ(ierr); 1285 ierr = ISDestroy(&part_boundary_faces_is);CHKERRQ(ierr); 1286 PetscFunctionReturn(0); 1287 } 1288 1289 /* compute inter-partition interface including edges and vertices */ 1290 static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is) 1291 { 1292 DMLabel label; 1293 PetscInt value = 1; 1294 const char interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface"; 1295 DMPlexInterpolatedFlag intp; 1296 MPI_Comm comm; 1297 PetscErrorCode ierr; 1298 1299 PetscFunctionBegin; 1300 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1301 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1302 if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1303 1304 ierr = DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label);CHKERRQ(ierr); 1305 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1306 ierr = DMLabelSetStratumIS(label, value, interface_faces_is);CHKERRQ(ierr); 1307 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm);CHKERRQ(ierr); 1308 ierr = DMPlexLabelComplete(ipdm, label);CHKERRQ(ierr); 1309 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm);CHKERRQ(ierr); 1310 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is);CHKERRQ(ierr); 1311 ierr = PetscObjectSetName((PetscObject)*interface_is, "interface_is");CHKERRQ(ierr); 1312 ierr = ISViewFromOptions(*interface_is, NULL, "-interface_is_view");CHKERRQ(ierr); 1313 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1314 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1315 PetscFunctionReturn(0); 1316 } 1317 1318 static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is) 1319 { 1320 PetscInt n; 1321 const PetscInt *arr; 1322 PetscErrorCode ierr; 1323 1324 PetscFunctionBegin; 1325 ierr = PetscSFGetGraph(sf, NULL, &n, &arr, NULL);CHKERRQ(ierr); 1326 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is);CHKERRQ(ierr); 1327 PetscFunctionReturn(0); 1328 } 1329 1330 static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is) 1331 { 1332 PetscInt n; 1333 const PetscInt *rootdegree; 1334 PetscInt *arr; 1335 PetscErrorCode ierr; 1336 1337 PetscFunctionBegin; 1338 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1339 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 1340 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 1341 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr);CHKERRQ(ierr); 1342 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345 1346 static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is) 1347 { 1348 IS pointSF_out_is, pointSF_in_is; 1349 PetscErrorCode ierr; 1350 1351 PetscFunctionBegin; 1352 ierr = PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is);CHKERRQ(ierr); 1353 ierr = PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is);CHKERRQ(ierr); 1354 ierr = ISExpand(pointSF_out_is, pointSF_in_is, is);CHKERRQ(ierr); 1355 ierr = ISDestroy(&pointSF_out_is);CHKERRQ(ierr); 1356 ierr = ISDestroy(&pointSF_in_is);CHKERRQ(ierr); 1357 PetscFunctionReturn(0); 1358 } 1359 1360 #define MYCHKERRQ(ierr) do {if (PetscUnlikely(ierr)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!");} while (0) 1361 1362 static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v) 1363 { 1364 DMLabel label; 1365 PetscSection coordsSection; 1366 Vec coordsVec; 1367 PetscScalar *coordsScalar; 1368 PetscInt coneSize, depth, dim, i, p, npoints; 1369 const PetscInt *points; 1370 PetscErrorCode ierr; 1371 1372 PetscFunctionBegin; 1373 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1374 ierr = DMGetCoordinateSection(dm, &coordsSection);CHKERRQ(ierr); 1375 ierr = DMGetCoordinatesLocal(dm, &coordsVec);CHKERRQ(ierr); 1376 ierr = VecGetArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1377 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 1378 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 1379 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 1380 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 1381 for (i=0; i<npoints; i++) { 1382 p = points[i]; 1383 ierr = DMLabelGetValue(label, p, &depth);CHKERRQ(ierr); 1384 if (!depth) { 1385 PetscInt c, n, o; 1386 PetscReal coords[3]; 1387 char coordstr[128]; 1388 1389 ierr = PetscSectionGetDof(coordsSection, p, &n);CHKERRQ(ierr); 1390 ierr = PetscSectionGetOffset(coordsSection, p, &o);CHKERRQ(ierr); 1391 for (c=0; c<n; c++) coords[c] = PetscRealPart(coordsScalar[o+c]); 1392 ierr = coord2str(coordstr, sizeof(coordstr), n, coords, 1.0);CHKERRQ(ierr); 1393 ierr = PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr);CHKERRQ(ierr); 1394 } else { 1395 char entityType[16]; 1396 1397 switch (depth) { 1398 case 1: ierr = PetscStrcpy(entityType, "edge");CHKERRQ(ierr); break; 1399 case 2: ierr = PetscStrcpy(entityType, "face");CHKERRQ(ierr); break; 1400 case 3: ierr = PetscStrcpy(entityType, "cell");CHKERRQ(ierr); break; 1401 default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3");CHKERRQ(ierr); 1402 } 1403 if (depth == dim && dim < 3) { 1404 ierr = PetscStrlcat(entityType, " (cell)", sizeof(entityType));CHKERRQ(ierr); 1405 } 1406 ierr = PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p);CHKERRQ(ierr); 1407 } 1408 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1409 if (coneSize) { 1410 const PetscInt *cone; 1411 IS coneIS; 1412 1413 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1414 ierr = ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS);CHKERRQ(ierr); 1415 ierr = ViewPointsWithType_Internal(dm, coneIS, v);CHKERRQ(ierr); 1416 ierr = ISDestroy(&coneIS);CHKERRQ(ierr); 1417 } 1418 } 1419 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 1420 ierr = VecRestoreArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1421 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 1422 PetscFunctionReturn(0); 1423 } 1424 1425 static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v) 1426 { 1427 PetscBool flg; 1428 PetscInt npoints; 1429 PetscMPIInt rank; 1430 PetscErrorCode ierr; 1431 1432 PetscFunctionBegin; 1433 ierr = PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg);CHKERRQ(ierr); 1434 if (!flg) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer");CHKERRQ(ierr); 1435 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRQ(ierr); 1436 ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr); 1437 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 1438 if (npoints) { 1439 ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank);CHKERRQ(ierr); 1440 ierr = ViewPointsWithType_Internal(dm, points, v);CHKERRQ(ierr); 1441 } 1442 ierr = PetscViewerFlush(v);CHKERRQ(ierr); 1443 ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is) 1448 { 1449 PetscSF pointsf; 1450 IS pointsf_is; 1451 PetscBool flg; 1452 MPI_Comm comm; 1453 PetscMPIInt size; 1454 PetscErrorCode ierr; 1455 1456 PetscFunctionBegin; 1457 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1458 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1459 ierr = DMGetPointSF(ipdm, &pointsf);CHKERRQ(ierr); 1460 if (pointsf) { 1461 PetscInt nroots; 1462 ierr = PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1463 if (nroots < 0) pointsf = NULL; /* uninitialized SF */ 1464 } 1465 if (!pointsf) { 1466 PetscInt N=0; 1467 if (interface_is) {ierr = ISGetSize(interface_is, &N);CHKERRQ(ierr);} 1468 if (N) SETERRQ(comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL"); 1469 PetscFunctionReturn(0); 1470 } 1471 1472 /* get PointSF points as IS pointsf_is */ 1473 ierr = PointSFGetInterfacePoints_Private(pointsf, &pointsf_is);CHKERRQ(ierr); 1474 1475 /* compare pointsf_is with interface_is */ 1476 ierr = ISEqual(interface_is, pointsf_is, &flg);CHKERRQ(ierr); 1477 ierr = MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm);CHKERRQ(ierr); 1478 if (!flg) { 1479 IS pointsf_extra_is, pointsf_missing_is; 1480 PetscViewer errv = PETSC_VIEWER_STDERR_(comm); 1481 ierr = ISDifference(interface_is, pointsf_is, &pointsf_missing_is);MYCHKERRQ(ierr); 1482 ierr = ISDifference(pointsf_is, interface_is, &pointsf_extra_is);MYCHKERRQ(ierr); 1483 ierr = PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n");MYCHKERRQ(ierr); 1484 ierr = ViewPointsWithType(ipdm, pointsf_missing_is, errv);MYCHKERRQ(ierr); 1485 ierr = PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n");MYCHKERRQ(ierr); 1486 ierr = ViewPointsWithType(ipdm, pointsf_extra_is, errv);MYCHKERRQ(ierr); 1487 ierr = ISDestroy(&pointsf_extra_is);MYCHKERRQ(ierr); 1488 ierr = ISDestroy(&pointsf_missing_is);MYCHKERRQ(ierr); 1489 SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above."); 1490 } 1491 ierr = ISDestroy(&pointsf_is);CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 /* remove faces & edges from label, leave just vertices */ 1496 static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points) 1497 { 1498 PetscInt vStart, vEnd; 1499 MPI_Comm comm; 1500 PetscErrorCode ierr; 1501 1502 PetscFunctionBegin; 1503 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1504 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1505 ierr = ISGeneralFilter(points, vStart, vEnd);CHKERRQ(ierr); 1506 PetscFunctionReturn(0); 1507 } 1508 1509 /* 1510 DMPlexCheckPointSFHeavy - Thorough test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points. 1511 1512 Collective 1513 1514 Input Parameters: 1515 . dm - The DMPlex object 1516 1517 Notes: 1518 The input DMPlex must be serial (one partition has all points, the other partitions have no points). 1519 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). 1520 This is mainly intended for debugging/testing purposes. 1521 1522 Algorithm: 1523 1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces() 1524 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 1525 3. the mesh is distributed or loaded in parallel 1526 4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices() 1527 5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces() 1528 6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary 1529 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 1530 1531 Level: developer 1532 1533 .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces() 1534 */ 1535 static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd) 1536 { 1537 DM ipdm=NULL; 1538 IS boundary_faces_is, interface_faces_is, interface_is; 1539 DMPlexInterpolatedFlag intp; 1540 MPI_Comm comm; 1541 PetscErrorCode ierr; 1542 1543 PetscFunctionBegin; 1544 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1545 1546 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1547 if (intp == DMPLEX_INTERPOLATED_FULL) { 1548 ipdm = dm; 1549 } else { 1550 /* create temporary interpolated DM if input DM is not interpolated */ 1551 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_FALSE);CHKERRQ(ierr); 1552 ierr = DMPlexInterpolate(dm, &ipdm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */ 1553 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_TRUE);CHKERRQ(ierr); 1554 } 1555 ierr = DMViewFromOptions(ipdm, NULL, "-ipdm_view");CHKERRQ(ierr); 1556 1557 /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */ 1558 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is);CHKERRQ(ierr); 1559 /* get inter-partition interface faces (interface_faces_is)*/ 1560 ierr = DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is);CHKERRQ(ierr); 1561 /* compute inter-partition interface including edges and vertices (interface_is) */ 1562 ierr = DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is);CHKERRQ(ierr); 1563 /* destroy immediate ISs */ 1564 ierr = ISDestroy(&boundary_faces_is);CHKERRQ(ierr); 1565 ierr = ISDestroy(&interface_faces_is);CHKERRQ(ierr); 1566 1567 /* for uninterpolated case, keep just vertices in interface */ 1568 if (!intp) { 1569 ierr = DMPlexISFilterVertices_Private(ipdm, interface_is);CHKERRQ(ierr); 1570 ierr = DMDestroy(&ipdm);CHKERRQ(ierr); 1571 } 1572 1573 /* compare PointSF with the boundary reconstructed from coordinates */ 1574 ierr = DMPlexComparePointSFWithInterface_Private(dm, interface_is);CHKERRQ(ierr); 1575 ierr = PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n");CHKERRQ(ierr); 1576 ierr = ISDestroy(&interface_is);CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 int main(int argc, char **argv) 1581 { 1582 DM dm; 1583 AppCtx user; 1584 PetscErrorCode ierr; 1585 1586 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 1587 ierr = PetscLogStageRegister("create",&stage[0]);CHKERRQ(ierr); 1588 ierr = PetscLogStageRegister("distribute",&stage[1]);CHKERRQ(ierr); 1589 ierr = PetscLogStageRegister("interpolate",&stage[2]);CHKERRQ(ierr); 1590 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1591 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1592 if (user.nPointsToExpand) { 1593 ierr = TestExpandPoints(dm, &user);CHKERRQ(ierr); 1594 } 1595 if (user.ncoords) { 1596 ierr = ViewVerticesFromCoords(dm, user.ncoords/user.dim, user.coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1597 } 1598 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1599 ierr = PetscFinalize(); 1600 return ierr; 1601 } 1602 1603 /*TEST 1604 1605 testset: 1606 nsize: 2 1607 args: -dm_view ascii::ascii_info_detail 1608 args: -dm_plex_check_all 1609 test: 1610 suffix: 1_tri_dist0 1611 args: -distribute 0 -interpolate {{none create}separate output} 1612 test: 1613 suffix: 1_tri_dist1 1614 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1615 test: 1616 suffix: 1_quad_dist0 1617 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1618 test: 1619 suffix: 1_quad_dist1 1620 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1621 test: 1622 suffix: 1_1d_dist1 1623 args: -dim 1 -distribute 1 1624 1625 testset: 1626 nsize: 3 1627 args: -testnum 1 -interpolate create 1628 args: -dm_plex_check_all 1629 test: 1630 suffix: 2 1631 args: -dm_view ascii::ascii_info_detail 1632 test: 1633 suffix: 2a 1634 args: -dm_plex_check_cones_conform_on_interfaces_verbose 1635 test: 1636 suffix: 2b 1637 args: -test_expand_points 0,1,2,5,6 1638 test: 1639 suffix: 2c 1640 args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty 1641 1642 testset: 1643 # the same as 1% for 3D 1644 nsize: 2 1645 args: -dim 3 -dm_view ascii::ascii_info_detail 1646 args: -dm_plex_check_all 1647 test: 1648 suffix: 4_tet_dist0 1649 args: -distribute 0 -interpolate {{none create}separate output} 1650 test: 1651 suffix: 4_tet_dist1 1652 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1653 test: 1654 suffix: 4_hex_dist0 1655 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1656 test: 1657 suffix: 4_hex_dist1 1658 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1659 1660 test: 1661 # the same as 4_tet_dist0 but test different initial orientations 1662 suffix: 4_tet_test_orient 1663 nsize: 2 1664 args: -dim 3 -distribute 0 1665 args: -dm_plex_check_all 1666 args: -rotate_interface_0 {{0 1 2 11 12 13}} 1667 args: -rotate_interface_1 {{0 1 2 11 12 13}} 1668 1669 testset: 1670 requires: exodusii 1671 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo 1672 args: -dm_view ascii::ascii_info_detail 1673 args: -dm_plex_check_all 1674 args: -custom_view 1675 test: 1676 suffix: 5_seq 1677 nsize: 1 1678 args: -distribute 0 -interpolate {{none create}separate output} 1679 test: 1680 suffix: 5_dist0 1681 nsize: 2 1682 args: -distribute 0 -interpolate {{none create}separate output} 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_plex_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_plex_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_plex_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_plex_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_plex_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_plex_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