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