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