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, PETSC_MAX_PATH_LEN, 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 int *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 double *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 = DMPlexCreateFromCellList(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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 int 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 = DMPlexCreateFromCellListParallel(comm, user->dim, numCells, numVertices, 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 defined(PETSC_USE_DEBUG) 1145 { 1146 ierr = PetscOptionsGetBool(options, prefix, opt, &flg, &set);CHKERRQ(ierr); 1147 if (PetscUnlikely(set && flg != enable)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect"); 1148 } 1149 #endif 1150 PetscFunctionReturn(0); 1151 } 1152 1153 /* get coordinate description of the whole-domain boundary */ 1154 static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary) 1155 { 1156 PortableBoundary bnd0, bnd; 1157 MPI_Comm comm; 1158 DM idm; 1159 DMLabel label; 1160 PetscInt d; 1161 const char boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary"; 1162 IS boundary_is; 1163 IS *boundary_expanded_iss; 1164 PetscMPIInt rootrank = 0; 1165 PetscMPIInt rank, size; 1166 PetscInt value = 1; 1167 DMPlexInterpolatedFlag intp; 1168 PetscBool flg; 1169 PetscErrorCode ierr; 1170 1171 PetscFunctionBegin; 1172 ierr = PetscNew(&bnd);CHKERRQ(ierr); 1173 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1174 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1175 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1176 ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr); 1177 if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed"); 1178 1179 /* interpolate serial DM if not yet interpolated */ 1180 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1181 if (intp == DMPLEX_INTERPOLATED_FULL) { 1182 idm = dm; 1183 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 1184 } else { 1185 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 1186 ierr = DMViewFromOptions(idm, NULL, "-idm_view");CHKERRQ(ierr); 1187 } 1188 1189 /* mark whole-domain boundary of the serial DM */ 1190 ierr = DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label);CHKERRQ(ierr); 1191 ierr = DMAddLabel(idm, label);CHKERRQ(ierr); 1192 ierr = DMPlexMarkBoundaryFaces(idm, value, label);CHKERRQ(ierr); 1193 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm);CHKERRQ(ierr); 1194 ierr = DMLabelGetStratumIS(label, value, &boundary_is);CHKERRQ(ierr); 1195 1196 /* translate to coordinates */ 1197 ierr = PetscNew(&bnd0);CHKERRQ(ierr); 1198 ierr = DMGetCoordinatesLocalSetUp(idm);CHKERRQ(ierr); 1199 if (rank == rootrank) { 1200 ierr = DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1201 ierr = DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates);CHKERRQ(ierr); 1202 /* self-check */ 1203 { 1204 IS is0; 1205 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0);CHKERRQ(ierr); 1206 ierr = ISEqual(is0, boundary_is, &flg);CHKERRQ(ierr); 1207 if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS"); 1208 ierr = ISDestroy(&is0);CHKERRQ(ierr); 1209 } 1210 } else { 1211 ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates);CHKERRQ(ierr); 1212 } 1213 1214 { 1215 Vec tmp; 1216 VecScatter sc; 1217 IS xis; 1218 PetscInt n; 1219 1220 /* just convert seq vectors to mpi vector */ 1221 ierr = VecGetLocalSize(bnd0->coordinates, &n);CHKERRQ(ierr); 1222 ierr = MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1223 if (rank == rootrank) { 1224 ierr = VecCreateMPI(comm, n, n, &tmp);CHKERRQ(ierr); 1225 } else { 1226 ierr = VecCreateMPI(comm, 0, n, &tmp);CHKERRQ(ierr); 1227 } 1228 ierr = VecCopy(bnd0->coordinates, tmp);CHKERRQ(ierr); 1229 ierr = VecDestroy(&bnd0->coordinates);CHKERRQ(ierr); 1230 bnd0->coordinates = tmp; 1231 1232 /* replicate coordinates from root rank to all ranks */ 1233 ierr = VecCreateMPI(comm, n, n*size, &bnd->coordinates);CHKERRQ(ierr); 1234 ierr = ISCreateStride(comm, n, 0, 1, &xis);CHKERRQ(ierr); 1235 ierr = VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc);CHKERRQ(ierr); 1236 ierr = VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1237 ierr = VecScatterEnd( sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1238 ierr = VecScatterDestroy(&sc);CHKERRQ(ierr); 1239 ierr = ISDestroy(&xis);CHKERRQ(ierr); 1240 } 1241 bnd->depth = bnd0->depth; 1242 ierr = MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm);CHKERRQ(ierr); 1243 ierr = PetscMalloc1(bnd->depth, &bnd->sections);CHKERRQ(ierr); 1244 for (d=0; d<bnd->depth; d++) { 1245 ierr = PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d]);CHKERRQ(ierr); 1246 } 1247 1248 if (rank == rootrank) { 1249 ierr = DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections);CHKERRQ(ierr); 1250 } 1251 ierr = PortableBoundaryDestroy(&bnd0);CHKERRQ(ierr); 1252 ierr = DMRemoveLabelBySelf(idm, &label, PETSC_TRUE);CHKERRQ(ierr); 1253 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1254 ierr = ISDestroy(&boundary_is);CHKERRQ(ierr); 1255 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1256 *boundary = bnd; 1257 PetscFunctionReturn(0); 1258 } 1259 1260 /* get faces of inter-partition interface */ 1261 static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is) 1262 { 1263 MPI_Comm comm; 1264 DMLabel label; 1265 IS part_boundary_faces_is; 1266 const char partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary"; 1267 PetscInt value = 1; 1268 DMPlexInterpolatedFlag intp; 1269 PetscErrorCode ierr; 1270 1271 PetscFunctionBegin; 1272 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1273 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1274 if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1275 1276 /* get ipdm partition boundary (partBoundary) */ 1277 ierr = DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label);CHKERRQ(ierr); 1278 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1279 ierr = DMPlexMarkBoundaryFaces(ipdm, value, label);CHKERRQ(ierr); 1280 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm);CHKERRQ(ierr); 1281 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is);CHKERRQ(ierr); 1282 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1283 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1284 1285 /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */ 1286 ierr = ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is);CHKERRQ(ierr); 1287 ierr = ISDestroy(&part_boundary_faces_is);CHKERRQ(ierr); 1288 PetscFunctionReturn(0); 1289 } 1290 1291 /* compute inter-partition interface including edges and vertices */ 1292 static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is) 1293 { 1294 DMLabel label; 1295 PetscInt value = 1; 1296 const char interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface"; 1297 DMPlexInterpolatedFlag intp; 1298 MPI_Comm comm; 1299 PetscErrorCode ierr; 1300 1301 PetscFunctionBegin; 1302 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1303 ierr = DMPlexIsInterpolatedCollective(ipdm, &intp);CHKERRQ(ierr); 1304 if (intp != DMPLEX_INTERPOLATED_FULL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); 1305 1306 ierr = DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label);CHKERRQ(ierr); 1307 ierr = DMAddLabel(ipdm, label);CHKERRQ(ierr); 1308 ierr = DMLabelSetStratumIS(label, value, interface_faces_is);CHKERRQ(ierr); 1309 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm);CHKERRQ(ierr); 1310 ierr = DMPlexLabelComplete(ipdm, label);CHKERRQ(ierr); 1311 ierr = DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm);CHKERRQ(ierr); 1312 ierr = DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is);CHKERRQ(ierr); 1313 ierr = PetscObjectSetName((PetscObject)*interface_is, "interface_is");CHKERRQ(ierr); 1314 ierr = ISViewFromOptions(*interface_is, NULL, "-interface_is_view");CHKERRQ(ierr); 1315 ierr = DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE);CHKERRQ(ierr); 1316 ierr = DMLabelDestroy(&label);CHKERRQ(ierr); 1317 PetscFunctionReturn(0); 1318 } 1319 1320 static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is) 1321 { 1322 PetscInt n; 1323 const PetscInt *arr; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 ierr = PetscSFGetGraph(sf, NULL, &n, &arr, NULL);CHKERRQ(ierr); 1328 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is);CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is) 1333 { 1334 PetscInt n; 1335 const PetscInt *rootdegree; 1336 PetscInt *arr; 1337 PetscErrorCode ierr; 1338 1339 PetscFunctionBegin; 1340 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 1341 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 1342 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 1343 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr);CHKERRQ(ierr); 1344 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 1345 PetscFunctionReturn(0); 1346 } 1347 1348 static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is) 1349 { 1350 IS pointSF_out_is, pointSF_in_is; 1351 PetscErrorCode ierr; 1352 1353 PetscFunctionBegin; 1354 ierr = PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is);CHKERRQ(ierr); 1355 ierr = PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is);CHKERRQ(ierr); 1356 ierr = ISExpand(pointSF_out_is, pointSF_in_is, is);CHKERRQ(ierr); 1357 ierr = ISDestroy(&pointSF_out_is);CHKERRQ(ierr); 1358 ierr = ISDestroy(&pointSF_in_is);CHKERRQ(ierr); 1359 PetscFunctionReturn(0); 1360 } 1361 1362 #define MYCHKERRQ(ierr) do {if (PetscUnlikely(ierr)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!");} while (0) 1363 1364 static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v) 1365 { 1366 DMLabel label; 1367 PetscSection coordsSection; 1368 Vec coordsVec; 1369 PetscScalar *coordsScalar; 1370 PetscInt coneSize, depth, dim, i, p, npoints; 1371 const PetscInt *points; 1372 PetscErrorCode ierr; 1373 1374 PetscFunctionBegin; 1375 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1376 ierr = DMGetCoordinateSection(dm, &coordsSection);CHKERRQ(ierr); 1377 ierr = DMGetCoordinatesLocal(dm, &coordsVec);CHKERRQ(ierr); 1378 ierr = VecGetArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1379 ierr = ISGetLocalSize(pointsIS, &npoints);CHKERRQ(ierr); 1380 ierr = ISGetIndices(pointsIS, &points);CHKERRQ(ierr); 1381 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 1382 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 1383 for (i=0; i<npoints; i++) { 1384 p = points[i]; 1385 ierr = DMLabelGetValue(label, p, &depth);CHKERRQ(ierr); 1386 if (!depth) { 1387 PetscInt c, n, o; 1388 PetscReal coords[3]; 1389 char coordstr[128]; 1390 1391 ierr = PetscSectionGetDof(coordsSection, p, &n);CHKERRQ(ierr); 1392 ierr = PetscSectionGetOffset(coordsSection, p, &o);CHKERRQ(ierr); 1393 for (c=0; c<n; c++) coords[c] = PetscRealPart(coordsScalar[o+c]); 1394 ierr = coord2str(coordstr, sizeof(coordstr), n, coords, 1.0);CHKERRQ(ierr); 1395 ierr = PetscViewerASCIISynchronizedPrintf(v, "vertex %D w/ coordinates %s\n", p, coordstr);CHKERRQ(ierr); 1396 } else { 1397 char entityType[16]; 1398 1399 switch (depth) { 1400 case 1: ierr = PetscStrcpy(entityType, "edge");CHKERRQ(ierr); break; 1401 case 2: ierr = PetscStrcpy(entityType, "face");CHKERRQ(ierr); break; 1402 case 3: ierr = PetscStrcpy(entityType, "cell");CHKERRQ(ierr); break; 1403 default: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for depth <= 3");CHKERRQ(ierr); 1404 } 1405 if (depth == dim && dim < 3) { 1406 ierr = PetscStrlcat(entityType, " (cell)", sizeof(entityType));CHKERRQ(ierr); 1407 } 1408 ierr = PetscViewerASCIISynchronizedPrintf(v, "%s %D\n", entityType, p);CHKERRQ(ierr); 1409 } 1410 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1411 if (coneSize) { 1412 const PetscInt *cone; 1413 IS coneIS; 1414 1415 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1416 ierr = ISCreateGeneral(PETSC_COMM_SELF, coneSize, cone, PETSC_USE_POINTER, &coneIS);CHKERRQ(ierr); 1417 ierr = ViewPointsWithType_Internal(dm, coneIS, v);CHKERRQ(ierr); 1418 ierr = ISDestroy(&coneIS);CHKERRQ(ierr); 1419 } 1420 } 1421 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 1422 ierr = VecRestoreArray(coordsVec, &coordsScalar);CHKERRQ(ierr); 1423 ierr = ISRestoreIndices(pointsIS, &points);CHKERRQ(ierr); 1424 PetscFunctionReturn(0); 1425 } 1426 1427 static PetscErrorCode ViewPointsWithType(DM dm, IS points, PetscViewer v) 1428 { 1429 PetscBool flg; 1430 PetscInt npoints; 1431 PetscMPIInt rank; 1432 PetscErrorCode ierr; 1433 1434 PetscFunctionBegin; 1435 ierr = PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &flg);CHKERRQ(ierr); 1436 if (!flg) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Only for ASCII viewer");CHKERRQ(ierr); 1437 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRQ(ierr); 1438 ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr); 1439 ierr = ISGetLocalSize(points, &npoints);CHKERRQ(ierr); 1440 if (npoints) { 1441 ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank);CHKERRQ(ierr); 1442 ierr = ViewPointsWithType_Internal(dm, points, v);CHKERRQ(ierr); 1443 } 1444 ierr = PetscViewerFlush(v);CHKERRQ(ierr); 1445 ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr); 1446 PetscFunctionReturn(0); 1447 } 1448 1449 static PetscErrorCode DMPlexComparePointSFWithInterface_Private(DM ipdm, IS interface_is) 1450 { 1451 PetscSF pointsf; 1452 IS pointsf_is; 1453 PetscBool flg; 1454 MPI_Comm comm; 1455 PetscMPIInt size; 1456 PetscErrorCode ierr; 1457 1458 PetscFunctionBegin; 1459 ierr = PetscObjectGetComm((PetscObject)ipdm, &comm);CHKERRQ(ierr); 1460 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1461 ierr = DMGetPointSF(ipdm, &pointsf);CHKERRQ(ierr); 1462 if (pointsf) { 1463 PetscInt nroots; 1464 ierr = PetscSFGetGraph(pointsf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 1465 if (nroots < 0) pointsf = NULL; /* uninitialized SF */ 1466 } 1467 if (!pointsf) { 1468 PetscInt N=0; 1469 if (interface_is) {ierr = ISGetSize(interface_is, &N);CHKERRQ(ierr);} 1470 if (N) SETERRQ(comm, PETSC_ERR_PLIB, "interface_is should be NULL or empty for PointSF being NULL"); 1471 PetscFunctionReturn(0); 1472 } 1473 1474 /* get PointSF points as IS pointsf_is */ 1475 ierr = PointSFGetInterfacePoints_Private(pointsf, &pointsf_is);CHKERRQ(ierr); 1476 1477 /* compare pointsf_is with interface_is */ 1478 ierr = ISEqual(interface_is, pointsf_is, &flg);CHKERRQ(ierr); 1479 ierr = MPI_Allreduce(MPI_IN_PLACE,&flg,1,MPIU_BOOL,MPI_LAND,comm);CHKERRQ(ierr); 1480 if (!flg) { 1481 IS pointsf_extra_is, pointsf_missing_is; 1482 PetscViewer errv = PETSC_VIEWER_STDERR_(comm); 1483 ierr = ISDifference(interface_is, pointsf_is, &pointsf_missing_is);MYCHKERRQ(ierr); 1484 ierr = ISDifference(pointsf_is, interface_is, &pointsf_extra_is);MYCHKERRQ(ierr); 1485 ierr = PetscViewerASCIIPrintf(errv, "Points missing in PointSF:\n");MYCHKERRQ(ierr); 1486 ierr = ViewPointsWithType(ipdm, pointsf_missing_is, errv);MYCHKERRQ(ierr); 1487 ierr = PetscViewerASCIIPrintf(errv, "Extra points in PointSF:\n");MYCHKERRQ(ierr); 1488 ierr = ViewPointsWithType(ipdm, pointsf_extra_is, errv);MYCHKERRQ(ierr); 1489 ierr = ISDestroy(&pointsf_extra_is);MYCHKERRQ(ierr); 1490 ierr = ISDestroy(&pointsf_missing_is);MYCHKERRQ(ierr); 1491 SETERRQ(comm, PETSC_ERR_PLIB, "PointSF is wrong! See details above."); 1492 } 1493 ierr = ISDestroy(&pointsf_is);CHKERRQ(ierr); 1494 PetscFunctionReturn(0); 1495 } 1496 1497 /* remove faces & edges from label, leave just vertices */ 1498 static PetscErrorCode DMPlexISFilterVertices_Private(DM dm, IS points) 1499 { 1500 PetscInt vStart, vEnd; 1501 MPI_Comm comm; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1506 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1507 ierr = ISGeneralFilter(points, vStart, vEnd);CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /* 1512 DMPlexCheckPointSFHeavy - Thorough test that the PointSF after parallel DMPlexInterpolate() includes exactly all interface points. 1513 1514 Collective 1515 1516 Input Parameters: 1517 . dm - The DMPlex object 1518 1519 Notes: 1520 The input DMPlex must be serial (one partition has all points, the other partitions have no points). 1521 This is a heavy test which involves DMPlexInterpolate() if the input DM is not interpolated yet, and depends on having a representation of the whole-domain boundary (PortableBoundary), which can be obtained only by DMPlexGetExpandedBoundary_Private() (which involves DMPlexInterpolate() of a sequential DM). 1522 This is mainly intended for debugging/testing purposes. 1523 1524 Algorithm: 1525 1. boundary faces of the serial version of the whole mesh are found using DMPlexMarkBoundaryFaces() 1526 2. boundary faces are translated into vertices using DMPlexGetConeRecursive() and these are translated into coordinates - this description (aka PortableBoundary) is completely independent of partitioning and point numbering 1527 3. the mesh is distributed or loaded in parallel 1528 4. boundary faces of the distributed mesh are reconstructed from PortableBoundary using DMPlexFindVertices() 1529 5. partition boundary faces of the parallel mesh are found using DMPlexMarkBoundaryFaces() 1530 6. partition interfaces are computed as set difference of partition boundary faces minus the reconstructed boundary 1531 7. check that interface covered by PointSF (union of inward and outward points) is equal to the partition interface for each rank, otherwise print the difference and throw an error 1532 1533 Level: developer 1534 1535 .seealso: DMGetPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces() 1536 */ 1537 static PetscErrorCode DMPlexCheckPointSFHeavy(DM dm, PortableBoundary bnd) 1538 { 1539 DM ipdm=NULL; 1540 IS boundary_faces_is, interface_faces_is, interface_is; 1541 DMPlexInterpolatedFlag intp; 1542 MPI_Comm comm; 1543 PetscErrorCode ierr; 1544 1545 PetscFunctionBegin; 1546 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 1547 1548 ierr = DMPlexIsInterpolatedCollective(dm, &intp);CHKERRQ(ierr); 1549 if (intp == DMPLEX_INTERPOLATED_FULL) { 1550 ipdm = dm; 1551 } else { 1552 /* create temporary interpolated DM if input DM is not interpolated */ 1553 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_FALSE);CHKERRQ(ierr); 1554 ierr = DMPlexInterpolate(dm, &ipdm);CHKERRQ(ierr); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexComparePointSFWithInterface_Private() below */ 1555 ierr = DMPlexSetOrientInterface_Private(dm, PETSC_TRUE);CHKERRQ(ierr); 1556 } 1557 ierr = DMViewFromOptions(ipdm, NULL, "-ipdm_view");CHKERRQ(ierr); 1558 1559 /* recover ipdm whole-domain boundary faces from the expanded vertices coordinates */ 1560 ierr = DMPlexExpandedVerticesCoordinatesToFaces_Private(ipdm, bnd, &boundary_faces_is);CHKERRQ(ierr); 1561 /* get inter-partition interface faces (interface_faces_is)*/ 1562 ierr = DMPlexGetInterfaceFaces_Private(ipdm, boundary_faces_is, &interface_faces_is);CHKERRQ(ierr); 1563 /* compute inter-partition interface including edges and vertices (interface_is) */ 1564 ierr = DMPlexComputeCompleteInterface_Private(ipdm, interface_faces_is, &interface_is);CHKERRQ(ierr); 1565 /* destroy immediate ISs */ 1566 ierr = ISDestroy(&boundary_faces_is);CHKERRQ(ierr); 1567 ierr = ISDestroy(&interface_faces_is);CHKERRQ(ierr); 1568 1569 /* for uninterpolated case, keep just vertices in interface */ 1570 if (!intp) { 1571 ierr = DMPlexISFilterVertices_Private(ipdm, interface_is);CHKERRQ(ierr); 1572 ierr = DMDestroy(&ipdm);CHKERRQ(ierr); 1573 } 1574 1575 /* compare PointSF with the boundary reconstructed from coordinates */ 1576 ierr = DMPlexComparePointSFWithInterface_Private(dm, interface_is);CHKERRQ(ierr); 1577 ierr = PetscPrintf(comm, "DMPlexCheckPointSFHeavy PASSED\n");CHKERRQ(ierr); 1578 ierr = ISDestroy(&interface_is);CHKERRQ(ierr); 1579 PetscFunctionReturn(0); 1580 } 1581 1582 int main(int argc, char **argv) 1583 { 1584 DM dm; 1585 AppCtx user; 1586 PetscErrorCode ierr; 1587 1588 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 1589 ierr = PetscLogStageRegister("create",&stage[0]);CHKERRQ(ierr); 1590 ierr = PetscLogStageRegister("distribute",&stage[1]);CHKERRQ(ierr); 1591 ierr = PetscLogStageRegister("interpolate",&stage[2]);CHKERRQ(ierr); 1592 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1593 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1594 if (user.nPointsToExpand) { 1595 ierr = TestExpandPoints(dm, &user);CHKERRQ(ierr); 1596 } 1597 if (user.ncoords) { 1598 ierr = ViewVerticesFromCoords(dm, user.ncoords/user.dim, user.coords, user.coordsTol, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1599 } 1600 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1601 ierr = PetscFinalize(); 1602 return ierr; 1603 } 1604 1605 /*TEST 1606 1607 testset: 1608 nsize: 2 1609 args: -dm_view ascii::ascii_info_detail 1610 args: -dm_plex_check_all 1611 test: 1612 suffix: 1_tri_dist0 1613 args: -distribute 0 -interpolate {{none create}separate output} 1614 test: 1615 suffix: 1_tri_dist1 1616 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1617 test: 1618 suffix: 1_quad_dist0 1619 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1620 test: 1621 suffix: 1_quad_dist1 1622 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1623 test: 1624 suffix: 1_1d_dist1 1625 args: -dim 1 -distribute 1 1626 1627 testset: 1628 nsize: 3 1629 args: -testnum 1 -interpolate create 1630 args: -dm_plex_check_all 1631 test: 1632 suffix: 2 1633 args: -dm_view ascii::ascii_info_detail 1634 test: 1635 suffix: 2a 1636 args: -dm_plex_check_cones_conform_on_interfaces_verbose 1637 test: 1638 suffix: 2b 1639 args: -test_expand_points 0,1,2,5,6 1640 test: 1641 suffix: 2c 1642 args: -test_expand_points 0,1,2,5,6 -test_expand_points_empty 1643 1644 testset: 1645 # the same as 1% for 3D 1646 nsize: 2 1647 args: -dim 3 -dm_view ascii::ascii_info_detail 1648 args: -dm_plex_check_all 1649 test: 1650 suffix: 4_tet_dist0 1651 args: -distribute 0 -interpolate {{none create}separate output} 1652 test: 1653 suffix: 4_tet_dist1 1654 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1655 test: 1656 suffix: 4_hex_dist0 1657 args: -cell_simplex 0 -distribute 0 -interpolate {{none create}separate output} 1658 test: 1659 suffix: 4_hex_dist1 1660 args: -cell_simplex 0 -distribute 1 -interpolate {{none create after_distribute}separate output} 1661 1662 test: 1663 # the same as 4_tet_dist0 but test different initial orientations 1664 suffix: 4_tet_test_orient 1665 nsize: 2 1666 args: -dim 3 -distribute 0 1667 args: -dm_plex_check_all 1668 args: -rotate_interface_0 {{0 1 2 11 12 13}} 1669 args: -rotate_interface_1 {{0 1 2 11 12 13}} 1670 1671 testset: 1672 requires: exodusii 1673 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/TwoQuads.exo 1674 args: -dm_view ascii::ascii_info_detail 1675 args: -dm_plex_check_all 1676 args: -custom_view 1677 test: 1678 suffix: 5_seq 1679 nsize: 1 1680 args: -distribute 0 -interpolate {{none create}separate output} 1681 test: 1682 suffix: 5_dist0 1683 nsize: 2 1684 args: -distribute 0 -interpolate {{none create}separate output} 1685 test: 1686 suffix: 5_dist1 1687 nsize: 2 1688 args: -distribute 1 -interpolate {{none create after_distribute}separate output} 1689 1690 testset: 1691 nsize: {{1 2 4}} 1692 args: -use_generator 1693 args: -dm_plex_check_all 1694 args: -distribute -interpolate none 1695 test: 1696 suffix: 6_tri 1697 requires: triangle 1698 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1699 test: 1700 suffix: 6_quad 1701 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1702 test: 1703 suffix: 6_tet 1704 requires: ctetgen 1705 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1706 test: 1707 suffix: 6_hex 1708 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1709 testset: 1710 nsize: {{1 2 4}} 1711 args: -use_generator 1712 args: -dm_plex_check_all 1713 args: -distribute -interpolate create 1714 test: 1715 suffix: 6_int_tri 1716 requires: triangle 1717 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1718 test: 1719 suffix: 6_int_quad 1720 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1721 test: 1722 suffix: 6_int_tet 1723 requires: ctetgen 1724 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1725 test: 1726 suffix: 6_int_hex 1727 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1728 testset: 1729 nsize: {{2 4}} 1730 args: -use_generator 1731 args: -dm_plex_check_all 1732 args: -distribute -interpolate after_distribute 1733 test: 1734 suffix: 6_parint_tri 1735 requires: triangle 1736 args: -faces {{2,2 1,3 7,4}} -cell_simplex 1 -dm_plex_generator triangle 1737 test: 1738 suffix: 6_parint_quad 1739 args: -faces {{2,2 1,3 7,4}} -cell_simplex 0 1740 test: 1741 suffix: 6_parint_tet 1742 requires: ctetgen 1743 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 1 -dm_plex_generator ctetgen 1744 test: 1745 suffix: 6_parint_hex 1746 args: -faces {{2,2,2 1,3,5 3,4,7}} -cell_simplex 0 1747 1748 testset: # 7 EXODUS 1749 requires: exodusii 1750 args: -dm_plex_check_all 1751 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo 1752 args: -distribute 1753 test: # seq load, simple partitioner 1754 suffix: 7_exo 1755 nsize: {{1 2 4 5}} 1756 args: -interpolate none 1757 test: # seq load, seq interpolation, simple partitioner 1758 suffix: 7_exo_int_simple 1759 nsize: {{1 2 4 5}} 1760 args: -interpolate create 1761 test: # seq load, seq interpolation, metis partitioner 1762 suffix: 7_exo_int_metis 1763 requires: parmetis 1764 nsize: {{2 4 5}} 1765 args: -interpolate create 1766 args: -petscpartitioner_type parmetis 1767 test: # seq load, simple partitioner, par interpolation 1768 suffix: 7_exo_simple_int 1769 nsize: {{2 4 5}} 1770 args: -interpolate after_distribute 1771 test: # seq load, metis partitioner, par interpolation 1772 suffix: 7_exo_metis_int 1773 requires: parmetis 1774 nsize: {{2 4 5}} 1775 args: -interpolate after_distribute 1776 args: -petscpartitioner_type parmetis 1777 1778 testset: # 7 HDF5 SEQUANTIAL LOAD 1779 requires: hdf5 !complex 1780 args: -dm_plex_check_all 1781 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1782 args: -dm_plex_hdf5_force_sequential 1783 args: -distribute 1784 test: # seq load, simple partitioner 1785 suffix: 7_seq_hdf5_simple 1786 nsize: {{1 2 4 5}} 1787 args: -interpolate none 1788 test: # seq load, seq interpolation, simple partitioner 1789 suffix: 7_seq_hdf5_int_simple 1790 nsize: {{1 2 4 5}} 1791 args: -interpolate after_create 1792 test: # seq load, seq interpolation, metis partitioner 1793 nsize: {{2 4 5}} 1794 suffix: 7_seq_hdf5_int_metis 1795 requires: parmetis 1796 args: -interpolate after_create 1797 args: -petscpartitioner_type parmetis 1798 test: # seq load, simple partitioner, par interpolation 1799 suffix: 7_seq_hdf5_simple_int 1800 nsize: {{2 4 5}} 1801 args: -interpolate after_distribute 1802 test: # seq load, metis partitioner, par interpolation 1803 nsize: {{2 4 5}} 1804 suffix: 7_seq_hdf5_metis_int 1805 requires: parmetis 1806 args: -interpolate after_distribute 1807 args: -petscpartitioner_type parmetis 1808 1809 testset: # 7 HDF5 PARALLEL LOAD 1810 requires: hdf5 !complex 1811 nsize: {{2 4 5}} 1812 args: -dm_plex_check_all 1813 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1814 test: # par load 1815 suffix: 7_par_hdf5 1816 args: -interpolate none 1817 test: # par load, par interpolation 1818 suffix: 7_par_hdf5_int 1819 args: -interpolate after_create 1820 test: # par load, parmetis repartitioner 1821 TODO: Parallel partitioning of uninterpolated meshes not supported 1822 suffix: 7_par_hdf5_parmetis 1823 requires: parmetis 1824 args: -distribute -petscpartitioner_type parmetis 1825 args: -interpolate none 1826 test: # par load, par interpolation, parmetis repartitioner 1827 suffix: 7_par_hdf5_int_parmetis 1828 requires: parmetis 1829 args: -distribute -petscpartitioner_type parmetis 1830 args: -interpolate after_create 1831 test: # par load, parmetis partitioner, par interpolation 1832 TODO: Parallel partitioning of uninterpolated meshes not supported 1833 suffix: 7_par_hdf5_parmetis_int 1834 requires: parmetis 1835 args: -distribute -petscpartitioner_type parmetis 1836 args: -interpolate after_distribute 1837 1838 test: 1839 suffix: 7_hdf5_hierarch 1840 requires: hdf5 ptscotch !complex 1841 nsize: {{2 3 4}separate output} 1842 args: -distribute 1843 args: -interpolate after_create 1844 args: -petscpartitioner_type matpartitioning -petscpartitioner_view ::ascii_info 1845 args: -mat_partitioning_type hierarch -mat_partitioning_hierarchical_nfineparts 2 1846 args: -mat_partitioning_hierarchical_coarseparttype ptscotch -mat_partitioning_hierarchical_fineparttype ptscotch 1847 1848 test: 1849 suffix: 8 1850 requires: hdf5 !complex 1851 nsize: 4 1852 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -dm_plex_create_from_hdf5_xdmf 1853 args: -distribute 0 -interpolate after_create 1854 args: -view_vertices_from_coords 0.,1.,0.,-0.5,1.,0.,0.583,-0.644,0.,-2.,-2.,-2. -view_vertices_from_coords_tol 1e-3 1855 args: -dm_plex_check_all 1856 args: -custom_view 1857 1858 testset: # 9 HDF5 SEQUANTIAL LOAD 1859 requires: hdf5 !complex datafilespath 1860 args: -dm_plex_check_all 1861 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /cells -dm_plex_hdf5_geometry_path /coordinates 1862 args: -dm_plex_hdf5_force_sequential 1863 args: -distribute 1864 test: # seq load, simple partitioner 1865 suffix: 9_seq_hdf5_simple 1866 nsize: {{1 2 4 5}} 1867 args: -interpolate none 1868 test: # seq load, seq interpolation, simple partitioner 1869 suffix: 9_seq_hdf5_int_simple 1870 nsize: {{1 2 4 5}} 1871 args: -interpolate after_create 1872 test: # seq load, seq interpolation, metis partitioner 1873 nsize: {{2 4 5}} 1874 suffix: 9_seq_hdf5_int_metis 1875 requires: parmetis 1876 args: -interpolate after_create 1877 args: -petscpartitioner_type parmetis 1878 test: # seq load, simple partitioner, par interpolation 1879 suffix: 9_seq_hdf5_simple_int 1880 nsize: {{2 4 5}} 1881 args: -interpolate after_distribute 1882 test: # seq load, simple partitioner, par interpolation 1883 # This is like 9_seq_hdf5_simple_int but testing error output of DMPlexCheckPointSFHeavy(). 1884 # Once 9_seq_hdf5_simple_int gets fixed, this one gets broken. 1885 # We can then provide an intentionally broken mesh instead. 1886 TODO: This test is broken because PointSF is fixed. 1887 suffix: 9_seq_hdf5_simple_int_err 1888 nsize: 4 1889 args: -interpolate after_distribute 1890 filter: sed -e "/PETSC ERROR/,$$d" 1891 test: # seq load, metis partitioner, par interpolation 1892 nsize: {{2 4 5}} 1893 suffix: 9_seq_hdf5_metis_int 1894 requires: parmetis 1895 args: -interpolate after_distribute 1896 args: -petscpartitioner_type parmetis 1897 1898 testset: # 9 HDF5 PARALLEL LOAD 1899 requires: hdf5 !complex datafilespath 1900 nsize: {{2 4 5}} 1901 args: -dm_plex_check_all 1902 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /cells -dm_plex_hdf5_geometry_path /coordinates 1903 test: # par load 1904 suffix: 9_par_hdf5 1905 args: -interpolate none 1906 test: # par load, par interpolation 1907 suffix: 9_par_hdf5_int 1908 args: -interpolate after_create 1909 test: # par load, parmetis repartitioner 1910 TODO: Parallel partitioning of uninterpolated meshes not supported 1911 suffix: 9_par_hdf5_parmetis 1912 requires: parmetis 1913 args: -distribute -petscpartitioner_type parmetis 1914 args: -interpolate none 1915 test: # par load, par interpolation, parmetis repartitioner 1916 suffix: 9_par_hdf5_int_parmetis 1917 requires: parmetis 1918 args: -distribute -petscpartitioner_type parmetis 1919 args: -interpolate after_create 1920 test: # par load, parmetis partitioner, par interpolation 1921 TODO: Parallel partitioning of uninterpolated meshes not supported 1922 suffix: 9_par_hdf5_parmetis_int 1923 requires: parmetis 1924 args: -distribute -petscpartitioner_type parmetis 1925 args: -interpolate after_distribute 1926 1927 testset: # 10 HDF5 PARALLEL LOAD 1928 requires: hdf5 !complex datafilespath 1929 nsize: {{2 4 7}} 1930 args: -dm_plex_check_all 1931 args: -filename ${DATAFILESPATH}/meshes/cube-hexahedra-refined2.h5 -dm_plex_create_from_hdf5_xdmf -dm_plex_hdf5_topology_path /topo -dm_plex_hdf5_geometry_path /geom 1932 test: # par load, par interpolation 1933 suffix: 10_par_hdf5_int 1934 args: -interpolate after_create 1935 TEST*/ 1936