1 static char help[] = "Tests for parallel mesh loading and parallel topological interpolation\n\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 /* List of test meshes 5 6 Network 7 ------- 8 Test 0 (2 ranks): 9 10 network=0: 11 --------- 12 cell 0 cell 1 cell 2 nCells-1 (edge) 13 0 ------ 1 ------ 2 ------ 3 -- -- v -- -- nCells (vertex) 14 15 vertex distribution: 16 rank 0: 0 1 17 rank 1: 2 3 ... nCells 18 cell(edge) distribution: 19 rank 0: 0 1 20 rank 1: 2 ... nCells-1 21 22 network=1: 23 --------- 24 v2 25 ^ 26 | 27 cell 2 28 | 29 v0 --cell 0--> v3--cell 1--> v1 30 31 vertex distribution: 32 rank 0: 0 1 3 33 rank 1: 2 34 cell(edge) distribution: 35 rank 0: 0 1 36 rank 1: 2 37 38 example: 39 mpiexec -n 2 ./ex18 -distribute 1 -dim 1 -orig_dm_view -dist_dm_view -dist_dm_view -petscpartitioner_type parmetis -ncells 50 40 41 Triangle 42 -------- 43 Test 0 (2 ranks): 44 Two triangles sharing a face 45 46 2 47 / | \ 48 / | \ 49 / | \ 50 0 0 | 1 3 51 \ | / 52 \ | / 53 \ | / 54 1 55 56 vertex distribution: 57 rank 0: 0 1 58 rank 1: 2 3 59 cell distribution: 60 rank 0: 0 61 rank 1: 1 62 63 Test 1 (3 ranks): 64 Four triangles partitioned across 3 ranks 65 66 0 _______ 3 67 | \ / | 68 | \ 1 / | 69 | \ / | 70 | 0 2 2 | 71 | / \ | 72 | / 3 \ | 73 | / \ | 74 1 ------- 4 75 76 vertex distribution: 77 rank 0: 0 1 78 rank 1: 2 3 79 rank 2: 4 80 cell distribution: 81 rank 0: 0 82 rank 1: 1 83 rank 2: 2 3 84 85 Test 2 (3 ranks): 86 Four triangles partitioned across 3 ranks 87 88 1 _______ 3 89 | \ / | 90 | \ 1 / | 91 | \ / | 92 | 0 0 2 | 93 | / \ | 94 | / 3 \ | 95 | / \ | 96 2 ------- 4 97 98 vertex distribution: 99 rank 0: 0 1 100 rank 1: 2 3 101 rank 2: 4 102 cell distribution: 103 rank 0: 0 104 rank 1: 1 105 rank 2: 2 3 106 107 Tetrahedron 108 ----------- 109 Test 0: 110 Two tets sharing a face 111 112 cell 3 _______ cell 113 0 / | \ \ 1 114 / | \ \ 115 / | \ \ 116 0----|----4-----2 117 \ | / / 118 \ | / / 119 \ | / / 120 1------- 121 y 122 | x 123 |/ 124 *----z 125 126 vertex distribution: 127 rank 0: 0 1 128 rank 1: 2 3 4 129 cell distribution: 130 rank 0: 0 131 rank 1: 1 132 133 Quadrilateral 134 ------------- 135 Test 0 (2 ranks): 136 Two quads sharing a face 137 138 3-------2-------5 139 | | | 140 | 0 | 1 | 141 | | | 142 0-------1-------4 143 144 vertex distribution: 145 rank 0: 0 1 2 146 rank 1: 3 4 5 147 cell distribution: 148 rank 0: 0 149 rank 1: 1 150 151 TODO Test 1: 152 A quad and a triangle sharing a face 153 154 5-------4 155 | | \ 156 | 0 | \ 157 | | 1 \ 158 2-------3----6 159 160 Hexahedron 161 ---------- 162 Test 0 (2 ranks): 163 Two hexes sharing a face 164 165 cell 7-------------6-------------11 cell 166 0 /| /| /| 1 167 / | F1 / | F7 / | 168 / | / | / | 169 4-------------5-------------10 | 170 | | F4 | | F10 | | 171 | | | | | | 172 |F5 | |F3 | |F9 | 173 | | F2 | | F8 | | 174 | 3---------|---2---------|---9 175 | / | / | / 176 | / F0 | / F6 | / 177 |/ |/ |/ 178 0-------------1-------------8 179 180 vertex distribution: 181 rank 0: 0 1 2 3 4 5 182 rank 1: 6 7 8 9 10 11 183 cell distribution: 184 rank 0: 0 185 rank 1: 1 186 187 */ 188 189 typedef enum {NONE, CREATE, AFTER_CREATE, AFTER_DISTRIBUTE} InterpType; 190 191 typedef struct { 192 PetscInt debug; /* The debugging level */ 193 PetscInt testNum; /* Indicates the mesh to create */ 194 PetscInt dim; /* The topological mesh dimension */ 195 PetscBool cellSimplex; /* Use simplices or hexes */ 196 PetscBool distribute; /* Distribute the mesh */ 197 InterpType interpolate; /* Interpolate the mesh before or after DMPlexDistribute() */ 198 PetscBool useGenerator; /* Construct mesh with a mesh generator */ 199 PetscBool testOrientIF; /* Test for different original interface orientations */ 200 PetscBool testHeavy; /* Run the heavy PointSF test */ 201 PetscBool customView; /* Show results of DMPlexIsInterpolated() etc. */ 202 PetscInt ornt[2]; /* Orientation of interface on rank 0 and rank 1 */ 203 PetscInt faces[3]; /* Number of faces per dimension for generator */ 204 PetscReal coords[128]; 205 PetscReal coordsTol; 206 PetscInt ncoords; 207 PetscInt pointsToExpand[128]; 208 PetscInt nPointsToExpand; 209 PetscBool testExpandPointsEmpty; 210 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 211 } AppCtx; 212 213 struct _n_PortableBoundary { 214 Vec coordinates; 215 PetscInt depth; 216 PetscSection *sections; 217 }; 218 typedef struct _n_PortableBoundary * PortableBoundary; 219 220 #if defined(PETSC_USE_LOG) 221 static PetscLogStage stage[3]; 222 #endif 223 224 static PetscErrorCode DMPlexCheckPointSFHeavy(DM, PortableBoundary); 225 static PetscErrorCode DMPlexSetOrientInterface_Private(DM,PetscBool); 226 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM, PortableBoundary *); 227 static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM, IS, PetscSection, IS *); 228 229 static PetscErrorCode PortableBoundaryDestroy(PortableBoundary *bnd) 230 { 231 PetscInt d; 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 if (!*bnd) PetscFunctionReturn(0); 236 ierr = VecDestroy(&(*bnd)->coordinates);CHKERRQ(ierr); 237 for (d=0; d < (*bnd)->depth; d++) { 238 ierr = PetscSectionDestroy(&(*bnd)->sections[d]);CHKERRQ(ierr); 239 } 240 ierr = PetscFree((*bnd)->sections);CHKERRQ(ierr); 241 ierr = PetscFree(*bnd);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 246 { 247 const char *interpTypes[4] = {"none", "create", "after_create", "after_distribute"}; 248 PetscInt interp=NONE, dim; 249 PetscBool flg1, flg2; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 options->debug = 0; 254 options->testNum = 0; 255 options->dim = 2; 256 options->cellSimplex = PETSC_TRUE; 257 options->distribute = PETSC_FALSE; 258 options->interpolate = NONE; 259 options->useGenerator = PETSC_FALSE; 260 options->testOrientIF = PETSC_FALSE; 261 options->testHeavy = PETSC_TRUE; 262 options->customView = PETSC_FALSE; 263 options->testExpandPointsEmpty = PETSC_FALSE; 264 options->ornt[0] = 0; 265 options->ornt[1] = 0; 266 options->faces[0] = 2; 267 options->faces[1] = 2; 268 options->faces[2] = 2; 269 options->filename[0] = '\0'; 270 options->coordsTol = PETSC_DEFAULT; 271 272 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 273 ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 274 ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 275 ierr = PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 276 ierr = PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL);CHKERRQ(ierr); 277 ierr = PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL);CHKERRQ(ierr); 278 options->interpolate = (InterpType) interp; 279 if (!options->distribute && options->interpolate == AFTER_DISTRIBUTE) SETERRQ(comm, PETSC_ERR_SUP, "-interpolate after_distribute needs -distribute 1"); 280 ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr); 281 options->ncoords = 128; 282 ierr = PetscOptionsRealArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL);CHKERRQ(ierr); 283 ierr = PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL);CHKERRQ(ierr); 284 options->nPointsToExpand = 128; 285 ierr = PetscOptionsIntArray("-test_expand_points", "Expand given array of DAG point using DMPlexGetConeRecursive() and print results", "ex18.c", options->pointsToExpand, &options->nPointsToExpand, NULL);CHKERRQ(ierr); 286 if (options->nPointsToExpand) { 287 ierr = PetscOptionsBool("-test_expand_points_empty", "For -test_expand_points, rank 0 will have empty input array", "ex18.c", options->testExpandPointsEmpty, &options->testExpandPointsEmpty, NULL);CHKERRQ(ierr); 288 } 289 ierr = PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL);CHKERRQ(ierr); 290 ierr = PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL);CHKERRQ(ierr); 291 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3);CHKERRQ(ierr); 292 dim = 3; 293 ierr = PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2);CHKERRQ(ierr); 294 if (flg2) { 295 if (flg1 && dim != options->dim) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "specified -dim %D is not equal to length %D of -faces (note that -dim can be omitted)", options->dim, dim); 296 options->dim = dim; 297 } 298 ierr = PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 299 ierr = PetscOptionsBoundedInt("-rotate_interface_0", "Rotation (relative orientation) of interface on rank 0; implies -interpolate create -distribute 0", "ex18.c", options->ornt[0], &options->ornt[0], &options->testOrientIF,0);CHKERRQ(ierr); 300 ierr = PetscOptionsBoundedInt("-rotate_interface_1", "Rotation (relative orientation) of interface on rank 1; implies -interpolate create -distribute 0", "ex18.c", options->ornt[1], &options->ornt[1], &flg2,0);CHKERRQ(ierr); 301 if (flg2 != options->testOrientIF) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set"); 302 if (options->testOrientIF) { 303 PetscInt i; 304 for (i=0; i<2; i++) { 305 if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10); /* 11 12 13 become -1 -2 -3 */ 306 } 307 options->filename[0] = 0; 308 options->useGenerator = PETSC_FALSE; 309 options->dim = 3; 310 options->cellSimplex = PETSC_TRUE; 311 options->interpolate = CREATE; 312 options->distribute = PETSC_FALSE; 313 } 314 ierr = PetscOptionsEnd();CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) 319 { 320 PetscInt testNum = user->testNum; 321 PetscMPIInt rank,size; 322 PetscErrorCode ierr; 323 PetscInt numCorners=2,i; 324 PetscInt numCells,numVertices,network; 325 PetscInt *cells; 326 PetscReal *coords; 327 328 PetscFunctionBegin; 329 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 330 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 331 if (size > 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum); 332 333 numCells = 3; 334 ierr = PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL);CHKERRQ(ierr); 335 if (numCells < 3) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells); 336 337 if (size == 1) { 338 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);CHKERRMPI(ierr); 415 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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);CHKERRMPI(ierr); 536 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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 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 = DMPlexOrientPoint(*dm, ifp[rank], user->ornt[rank]);CHKERRQ(ierr); 575 ierr = DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view");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);CHKERRMPI(ierr); 591 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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);CHKERRMPI(ierr); 634 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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 == 0 && 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);CHKERRMPI(ierr); 908 if (user->testExpandPointsEmpty && rank == 0) { 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);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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));CHKERRQ(ierr); 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);CHKERRMPI(ierr); 1173 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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"); 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"); 1435 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(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 # Detail viewing in a non-distributed mesh is broken because the DMLabelView() is collective, but the label is not shared 1681 suffix: 5_dist0 1682 nsize: 2 1683 args: -distribute 0 -interpolate {{none create}separate output} -dm_view 1684 test: 1685 suffix: 5_dist1 1686 nsize: 2 1687 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1688 1689 testset: 1690 nsize: {{1 2 4}} 1691 args: -use_generator 1692 args: -dm_plex_check_all 1693 args: -distribute -interpolate none 1694 test: 1695 suffix: 6_tri 1696 requires: triangle 1697 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1698 test: 1699 suffix: 6_quad 1700 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1701 test: 1702 suffix: 6_tet 1703 requires: ctetgen 1704 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1705 test: 1706 suffix: 6_hex 1707 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1708 testset: 1709 nsize: {{1 2 4}} 1710 args: -use_generator 1711 args: -dm_plex_check_all 1712 args: -distribute -interpolate create 1713 test: 1714 suffix: 6_int_tri 1715 requires: triangle 1716 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1717 test: 1718 suffix: 6_int_quad 1719 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1720 test: 1721 suffix: 6_int_tet 1722 requires: ctetgen 1723 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1724 test: 1725 suffix: 6_int_hex 1726 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1727 testset: 1728 nsize: {{2 4}} 1729 args: -use_generator 1730 args: -dm_plex_check_all 1731 args: -distribute -interpolate after_distribute 1732 test: 1733 suffix: 6_parint_tri 1734 requires: triangle 1735 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1736 test: 1737 suffix: 6_parint_quad 1738 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1739 test: 1740 suffix: 6_parint_tet 1741 requires: ctetgen 1742 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1743 test: 1744 suffix: 6_parint_hex 1745 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1746 1747 testset: # 7 EXODUS 1748 requires: exodusii 1749 args: -dm_plex_check_all 1750 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo 1751 args: -distribute 1752 test: # seq load, simple partitioner 1753 suffix: 7_exo 1754 nsize: {{1 2 4 5}} 1755 args: -interpolate none 1756 test: # seq load, seq interpolation, simple partitioner 1757 suffix: 7_exo_int_simple 1758 nsize: {{1 2 4 5}} 1759 args: -interpolate create 1760 test: # seq load, seq interpolation, metis partitioner 1761 suffix: 7_exo_int_metis 1762 requires: parmetis 1763 nsize: {{2 4 5}} 1764 args: -interpolate create 1765 args: -petscpartitioner_type parmetis 1766 test: # seq load, simple partitioner, par interpolation 1767 suffix: 7_exo_simple_int 1768 nsize: {{2 4 5}} 1769 args: -interpolate after_distribute 1770 test: # seq load, metis partitioner, par interpolation 1771 suffix: 7_exo_metis_int 1772 requires: parmetis 1773 nsize: {{2 4 5}} 1774 args: -interpolate after_distribute 1775 args: -petscpartitioner_type parmetis 1776 1777 testset: # 7 HDF5 SEQUANTIAL LOAD 1778 requires: hdf5 !complex 1779 args: -dm_plex_check_all 1780 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1781 args: -dm_plex_hdf5_force_sequential 1782 args: -distribute 1783 test: # seq load, simple partitioner 1784 suffix: 7_seq_hdf5_simple 1785 nsize: {{1 2 4 5}} 1786 args: -interpolate none 1787 test: # seq load, seq interpolation, simple partitioner 1788 suffix: 7_seq_hdf5_int_simple 1789 nsize: {{1 2 4 5}} 1790 args: -interpolate after_create 1791 test: # seq load, seq interpolation, metis partitioner 1792 nsize: {{2 4 5}} 1793 suffix: 7_seq_hdf5_int_metis 1794 requires: parmetis 1795 args: -interpolate after_create 1796 args: -petscpartitioner_type parmetis 1797 test: # seq load, simple partitioner, par interpolation 1798 suffix: 7_seq_hdf5_simple_int 1799 nsize: {{2 4 5}} 1800 args: -interpolate after_distribute 1801 test: # seq load, metis partitioner, par interpolation 1802 nsize: {{2 4 5}} 1803 suffix: 7_seq_hdf5_metis_int 1804 requires: parmetis 1805 args: -interpolate after_distribute 1806 args: -petscpartitioner_type parmetis 1807 1808 testset: # 7 HDF5 PARALLEL LOAD 1809 requires: hdf5 !complex 1810 nsize: {{2 4 5}} 1811 args: -dm_plex_check_all 1812 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1813 test: # par load 1814 suffix: 7_par_hdf5 1815 args: -interpolate none 1816 test: # par load, par interpolation 1817 suffix: 7_par_hdf5_int 1818 args: -interpolate after_create 1819 test: # par load, parmetis repartitioner 1820 TODO: Parallel partitioning of uninterpolated meshes not supported 1821 suffix: 7_par_hdf5_parmetis 1822 requires: parmetis 1823 args: -distribute -petscpartitioner_type parmetis 1824 args: -interpolate none 1825 test: # par load, par interpolation, parmetis repartitioner 1826 suffix: 7_par_hdf5_int_parmetis 1827 requires: parmetis 1828 args: -distribute -petscpartitioner_type parmetis 1829 args: -interpolate after_create 1830 test: # par load, parmetis partitioner, par interpolation 1831 TODO: Parallel partitioning of uninterpolated meshes not supported 1832 suffix: 7_par_hdf5_parmetis_int 1833 requires: parmetis 1834 args: -distribute -petscpartitioner_type parmetis 1835 args: -interpolate after_distribute 1836 1837 test: 1838 suffix: 7_hdf5_hierarch 1839 requires: hdf5 ptscotch !complex 1840 nsize: {{2 3 4}separate output} 1841 args: -distribute 1842 args: -interpolate after_create 1843 args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info 1844 args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2 1845 args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch 1846 1847 test: 1848 suffix: 8 1849 requires: hdf5 !complex 1850 nsize: 4 1851 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1852 args: -distribute 0 -interpolate after_create 1853 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 1854 args: -dm_plex_check_all 1855 args: -custom_view 1856 1857 testset: # 9 HDF5 SEQUANTIAL LOAD 1858 requires: hdf5 !complex datafilespath 1859 args: -dm_plex_check_all 1860 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 1861 args: -dm_plex_hdf5_force_sequential 1862 args: -distribute 1863 test: # seq load, simple partitioner 1864 suffix: 9_seq_hdf5_simple 1865 nsize: {{1 2 4 5}} 1866 args: -interpolate none 1867 test: # seq load, seq interpolation, simple partitioner 1868 suffix: 9_seq_hdf5_int_simple 1869 nsize: {{1 2 4 5}} 1870 args: -interpolate after_create 1871 test: # seq load, seq interpolation, metis partitioner 1872 nsize: {{2 4 5}} 1873 suffix: 9_seq_hdf5_int_metis 1874 requires: parmetis 1875 args: -interpolate after_create 1876 args: -petscpartitioner_type parmetis 1877 test: # seq load, simple partitioner, par interpolation 1878 suffix: 9_seq_hdf5_simple_int 1879 nsize: {{2 4 5}} 1880 args: -interpolate after_distribute 1881 test: # seq load, simple partitioner, par interpolation 1882 # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy(). 1883 # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken. 1884 # We can then provide an intentionally broken mesh instead. 1885 TODO: This test is broken because PointSF is fixed. 1886 suffix: 9_seq_hdf5_simple_int_err 1887 nsize: 4 1888 args: -interpolate after_distribute 1889 filter: sed -e "/PETSC ERROR/,$$d" 1890 test: # seq load, metis partitioner, par interpolation 1891 nsize: {{2 4 5}} 1892 suffix: 9_seq_hdf5_metis_int 1893 requires: parmetis 1894 args: -interpolate after_distribute 1895 args: -petscpartitioner_type parmetis 1896 1897 testset: # 9 HDF5 PARALLEL LOAD 1898 requires: hdf5 !complex datafilespath 1899 nsize: {{2 4 5}} 1900 args: -dm_plex_check_all 1901 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 1902 test: # par load 1903 suffix: 9_par_hdf5 1904 args: -interpolate none 1905 test: # par load, par interpolation 1906 suffix: 9_par_hdf5_int 1907 args: -interpolate after_create 1908 test: # par load, parmetis repartitioner 1909 TODO: Parallel partitioning of uninterpolated meshes not supported 1910 suffix: 9_par_hdf5_parmetis 1911 requires: parmetis 1912 args: -distribute -petscpartitioner_type parmetis 1913 args: -interpolate none 1914 test: # par load, par interpolation, parmetis repartitioner 1915 suffix: 9_par_hdf5_int_parmetis 1916 requires: parmetis 1917 args: -distribute -petscpartitioner_type parmetis 1918 args: -interpolate after_create 1919 test: # par load, parmetis partitioner, par interpolation 1920 TODO: Parallel partitioning of uninterpolated meshes not supported 1921 suffix: 9_par_hdf5_parmetis_int 1922 requires: parmetis 1923 args: -distribute -petscpartitioner_type parmetis 1924 args: -interpolate after_distribute 1925 1926 testset: # 10 HDF5 PARALLEL LOAD 1927 requires: hdf5 !complex datafilespath 1928 nsize: {{2 4 7}} 1929 args: -dm_plex_check_all 1930 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 1931 test: # par load, par interpolation 1932 suffix: 10_par_hdf5_int 1933 args: -interpolate after_create 1934 TEST*/ 1935