Actual source code: ex56.c
1: #include <petscdmplex.h>
2: #include <petscviewerhdf5.h>
3: #include <petscsf.h>
5: static const char help[] = "Test DMLabel I/O with PETSc native HDF5 mesh format\n\n";
6: static const char EX[] = "ex56.c";
7: typedef struct {
8: MPI_Comm comm;
9: const char *meshname; /* Mesh name */
10: PetscInt num_labels; /* Asserted number of labels in loaded mesh */
11: PetscBool compare; /* Compare the meshes using DMPlexEqual() and DMCompareLabels() */
12: PetscBool compare_labels; /* Compare labels in the meshes using DMCompareLabels() */
13: PetscBool compare_boundary; /* Check label I/O via boundary vertex coordinates */
14: PetscBool compare_pre_post; /* Compare labels loaded before distribution with those loaded after distribution */
15: char outfile[PETSC_MAX_PATH_LEN]; /* Output file */
16: PetscBool use_low_level_functions; /* Use low level functions for viewing and loading */
17: //TODO This is meant as temporary option; can be removed once we have full parallel loading in place
18: PetscBool distribute_after_topo_load; /* Distribute topology right after DMPlexTopologyLoad(), if use_low_level_functions=true */
19: PetscInt verbose;
20: } AppCtx;
22: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
23: {
25: options->comm = comm;
26: options->num_labels = -1;
27: options->compare = PETSC_FALSE;
28: options->compare_labels = PETSC_FALSE;
29: options->compare_boundary = PETSC_FALSE;
30: options->compare_pre_post = PETSC_FALSE;
31: options->outfile[0] = '\0';
32: options->use_low_level_functions = PETSC_FALSE;
33: options->distribute_after_topo_load = PETSC_FALSE;
34: options->verbose = 0;
36: PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
37: PetscOptionsInt("-num_labels", "Asserted number of labels in meshfile; don't count depth and celltype; -1 to deactivate", EX, options->num_labels, &options->num_labels, NULL);
38: PetscOptionsBool("-compare", "Compare the meshes using DMPlexEqual() and DMCompareLabels()", EX, options->compare, &options->compare, NULL);
39: PetscOptionsBool("-compare_labels", "Compare labels in the meshes using DMCompareLabels()", "ex55.c", options->compare_labels, &options->compare_labels, NULL);
40: PetscOptionsBool("-compare_boundary", "Check label I/O via boundary vertex coordinates", "ex55.c", options->compare_boundary, &options->compare_boundary, NULL);
41: PetscOptionsBool("-compare_pre_post", "Compare labels loaded before distribution with those loaded after distribution", "ex55.c", options->compare_pre_post, &options->compare_pre_post, NULL);
42: PetscOptionsString("-outfile", "Output mesh file", EX, options->outfile, options->outfile, sizeof(options->outfile), NULL);
43: PetscOptionsBool("-use_low_level_functions", "Use low level functions for viewing and loading", EX, options->use_low_level_functions, &options->use_low_level_functions, NULL);
44: PetscOptionsBool("-distribute_after_topo_load", "Distribute topology right after DMPlexTopologyLoad(), if use_low_level_functions=true", EX, options->distribute_after_topo_load, &options->distribute_after_topo_load, NULL);
45: PetscOptionsInt("-verbose", "Verbosity level", EX, options->verbose, &options->verbose, NULL);
46: PetscOptionsEnd();
47: return 0;
48: };
50: static PetscErrorCode CreateMesh(AppCtx *options, DM *newdm)
51: {
52: DM dm;
55: DMCreate(options->comm, &dm);
56: DMSetType(dm, DMPLEX);
57: DMSetFromOptions(dm);
58: PetscObjectGetName((PetscObject)dm, &options->meshname);
59: DMViewFromOptions(dm, NULL, "-dm_view");
60: *newdm = dm;
61: return 0;
62: }
64: static PetscErrorCode SaveMesh(AppCtx *options, DM dm)
65: {
66: PetscViewer v;
69: PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), options->outfile, FILE_MODE_WRITE, &v);
70: if (options->use_low_level_functions) {
71: DMPlexTopologyView(dm, v);
72: DMPlexCoordinatesView(dm, v);
73: DMPlexLabelsView(dm, v);
74: } else {
75: DMView(dm, v);
76: }
77: PetscViewerDestroy(&v);
78: return 0;
79: }
81: typedef enum {
82: NONE = 0,
83: PRE_DIST = 1,
84: POST_DIST = 2
85: } AuxObjLoadMode;
87: static PetscErrorCode LoadMeshLowLevel(AppCtx *options, PetscViewer v, PetscBool explicitDistribute, AuxObjLoadMode mode, DM *newdm)
88: {
89: DM dm;
90: PetscSF sfXC;
93: DMCreate(options->comm, &dm);
94: DMSetType(dm, DMPLEX);
95: PetscObjectSetName((PetscObject)dm, options->meshname);
96: DMPlexTopologyLoad(dm, v, &sfXC);
97: if (mode == PRE_DIST) {
98: DMPlexCoordinatesLoad(dm, v, sfXC);
99: DMPlexLabelsLoad(dm, v, sfXC);
100: }
101: if (explicitDistribute) {
102: DM dmdist;
103: PetscSF sfXB = sfXC, sfBC;
105: DMPlexDistribute(dm, 0, &sfBC, &dmdist);
106: if (dmdist) {
107: const char *name;
109: PetscObjectGetName((PetscObject)dm, &name);
110: PetscObjectSetName((PetscObject)dmdist, name);
111: PetscSFCompose(sfXB, sfBC, &sfXC);
112: PetscSFDestroy(&sfXB);
113: PetscSFDestroy(&sfBC);
114: DMDestroy(&dm);
115: dm = dmdist;
116: }
117: }
118: if (mode == POST_DIST) {
119: DMPlexCoordinatesLoad(dm, v, sfXC);
120: DMPlexLabelsLoad(dm, v, sfXC);
121: }
122: PetscSFDestroy(&sfXC);
123: *newdm = dm;
124: return 0;
125: }
127: static PetscErrorCode LoadMesh(AppCtx *options, DM *dmnew)
128: {
129: DM dm;
130: PetscViewer v;
133: PetscViewerHDF5Open(options->comm, options->outfile, FILE_MODE_READ, &v);
134: if (options->use_low_level_functions) {
135: if (options->compare_pre_post) {
136: DM dm0;
138: LoadMeshLowLevel(options, v, PETSC_TRUE, PRE_DIST, &dm0);
139: LoadMeshLowLevel(options, v, PETSC_TRUE, POST_DIST, &dm);
140: DMCompareLabels(dm0, dm, NULL, NULL);
141: DMDestroy(&dm0);
142: } else {
143: LoadMeshLowLevel(options, v, options->distribute_after_topo_load, POST_DIST, &dm);
144: }
145: } else {
146: DMCreate(options->comm, &dm);
147: DMSetType(dm, DMPLEX);
148: PetscObjectSetName((PetscObject)dm, options->meshname);
149: DMLoad(dm, v);
150: }
151: PetscViewerDestroy(&v);
153: DMSetOptionsPrefix(dm, "load_");
154: DMSetFromOptions(dm);
155: DMViewFromOptions(dm, NULL, "-dm_view");
156: *dmnew = dm;
157: return 0;
158: }
160: static PetscErrorCode CompareMeshes(AppCtx *options, DM dm0, DM dm1)
161: {
162: PetscBool flg;
165: if (options->compare) {
166: DMPlexEqual(dm0, dm1, &flg);
168: PetscPrintf(options->comm, "DMs equal\n");
169: }
170: if (options->compare_labels) {
171: DMCompareLabels(dm0, dm1, NULL, NULL);
172: PetscPrintf(options->comm, "DMLabels equal\n");
173: }
174: return 0;
175: }
177: static PetscErrorCode MarkBoundary(DM dm, const char name[], PetscInt value, PetscBool verticesOnly, DMLabel *label)
178: {
179: DMLabel l;
182: DMLabelCreate(PetscObjectComm((PetscObject)dm), name, &l);
183: DMAddLabel(dm, l);
184: DMPlexMarkBoundaryFaces(dm, value, l);
185: DMPlexLabelComplete(dm, l);
186: if (verticesOnly) {
187: IS points;
188: const PetscInt *idx;
189: PetscInt i, n;
191: DMLabelGetStratumIS(l, value, &points);
192: ISGetLocalSize(points, &n);
193: ISGetIndices(points, &idx);
194: for (i = 0; i < n; i++) {
195: const PetscInt p = idx[i];
196: PetscInt d;
198: DMPlexGetPointDepth(dm, p, &d);
199: if (d != 0) DMLabelClearValue(l, p, value);
200: }
201: ISRestoreIndices(points, &idx);
202: ISDestroy(&points);
203: }
204: if (label) *label = l;
205: else DMLabelDestroy(&l);
206: return 0;
207: }
209: static PetscErrorCode VertexCoordinatesToAll(DM dm, IS vertices, Vec *allCoords)
210: {
211: Vec coords, allCoords_;
212: VecScatter sc;
213: MPI_Comm comm;
216: PetscObjectGetComm((PetscObject)dm, &comm);
217: DMGetCoordinatesLocalSetUp(dm);
218: if (vertices) {
219: DMGetCoordinatesLocalTuple(dm, vertices, NULL, &coords);
220: } else {
221: VecCreateSeq(PETSC_COMM_SELF, 0, &coords);
222: }
223: {
224: PetscInt n;
225: Vec mpivec;
227: VecGetLocalSize(coords, &n);
228: VecCreateMPI(comm, n, PETSC_DECIDE, &mpivec);
229: VecCopy(coords, mpivec);
230: VecDestroy(&coords);
231: coords = mpivec;
232: }
234: VecScatterCreateToAll(coords, &sc, &allCoords_);
235: VecScatterBegin(sc, coords, allCoords_, INSERT_VALUES, SCATTER_FORWARD);
236: VecScatterEnd(sc, coords, allCoords_, INSERT_VALUES, SCATTER_FORWARD);
237: VecScatterDestroy(&sc);
238: VecDestroy(&coords);
239: *allCoords = allCoords_;
240: return 0;
241: }
243: static PetscErrorCode DMLabelGetCoordinateRepresentation(DM dm, DMLabel label, PetscInt value, Vec *allCoords)
244: {
245: IS vertices;
248: DMLabelGetStratumIS(label, value, &vertices);
249: VertexCoordinatesToAll(dm, vertices, allCoords);
250: ISDestroy(&vertices);
251: return 0;
252: }
254: static PetscErrorCode DMLabelCompareWithCoordinateRepresentation(DM dm, DMLabel label, PetscInt value, Vec allCoords, PetscInt verbose)
255: {
256: const char *labelName;
257: IS pointsIS;
258: const PetscInt *points;
259: PetscInt i, n;
260: PetscBool fail = PETSC_FALSE;
261: MPI_Comm comm;
262: PetscMPIInt rank;
266: PetscObjectGetComm((PetscObject)dm, &comm);
267: PetscObjectGetName((PetscObject)label, &labelName);
268: MPI_Comm_rank(comm, &rank);
269: {
270: PetscInt pStart, pEnd;
272: DMPlexGetChart(dm, &pStart, &pEnd);
273: DMLabelCreateIndex(label, pStart, pEnd);
274: }
275: DMPlexFindVertices(dm, allCoords, 0.0, &pointsIS);
276: ISGetIndices(pointsIS, &points);
277: ISGetLocalSize(pointsIS, &n);
278: if (verbose > 1) DMLabelView(label, PETSC_VIEWER_STDOUT_(comm));
279: for (i = 0; i < n; i++) {
280: const PetscInt p = points[i];
281: PetscBool has;
282: PetscInt v;
284: if (p < 0) continue;
285: DMLabelHasPoint(label, p, &has);
286: if (!has) {
287: if (verbose) PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Label \"%s\" does not have point %" PetscInt_FMT "\n", rank, labelName, p);
288: fail = PETSC_TRUE;
289: continue;
290: }
291: DMLabelGetValue(label, p, &v);
292: if (v != value) {
293: if (verbose) PetscSynchronizedFPrintf(comm, PETSC_STDERR, "Point %" PetscInt_FMT " has bad value %" PetscInt_FMT " in label \"%s\"", p, v, labelName);
294: fail = PETSC_TRUE;
295: continue;
296: }
297: if (verbose > 1) PetscSynchronizedPrintf(comm, "[%d] OK point %" PetscInt_FMT "\n", rank, p);
298: }
299: PetscSynchronizedFlush(comm, PETSC_STDOUT);
300: PetscSynchronizedFlush(comm, PETSC_STDERR);
301: ISRestoreIndices(pointsIS, &points);
302: ISDestroy(&pointsIS);
303: MPI_Allreduce(MPI_IN_PLACE, &fail, 1, MPIU_BOOL, MPI_LOR, comm);
305: return 0;
306: }
308: static PetscErrorCode CheckNumLabels(DM dm, AppCtx *ctx)
309: {
310: PetscInt actualNum;
311: PetscBool fail = PETSC_FALSE;
312: MPI_Comm comm;
313: PetscMPIInt rank;
316: if (ctx->num_labels < 0) return 0;
317: PetscObjectGetComm((PetscObject)dm, &comm);
318: MPI_Comm_rank(comm, &rank);
319: DMGetNumLabels(dm, &actualNum);
320: if (ctx->num_labels != actualNum) {
321: fail = PETSC_TRUE;
322: if (ctx->verbose) {
323: PetscInt i;
325: PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Asserted number of labels: %" PetscInt_FMT ", actual: %" PetscInt_FMT "\n", rank, ctx->num_labels, actualNum);
326: for (i = 0; i < actualNum; i++) {
327: DMLabel label;
328: const char *name;
330: DMGetLabelByNum(dm, i, &label);
331: PetscObjectGetName((PetscObject)label, &name);
332: PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Label %" PetscInt_FMT " \"%s\"\n", rank, i, name);
333: }
334: PetscSynchronizedFlush(comm, PETSC_STDERR);
335: }
336: }
337: MPI_Allreduce(MPI_IN_PLACE, &fail, 1, MPIU_BOOL, MPI_LOR, comm);
339: return 0;
340: }
342: static inline PetscErrorCode IncrementNumLabels(AppCtx *ctx)
343: {
345: if (ctx->num_labels >= 0) ctx->num_labels++;
346: return 0;
347: }
349: int main(int argc, char **argv)
350: {
351: const char BOUNDARY_NAME[] = "Boundary";
352: const char BOUNDARY_VERTICES_NAME[] = "BoundaryVertices";
353: const PetscInt BOUNDARY_VALUE = 12345;
354: const PetscInt BOUNDARY_VERTICES_VALUE = 6789;
355: DM dm, dmnew;
356: AppCtx user;
357: Vec allCoords = NULL;
360: PetscInitialize(&argc, &argv, NULL, help);
361: ProcessOptions(PETSC_COMM_WORLD, &user);
362: CreateMesh(&user, &dm);
363: MarkBoundary(dm, BOUNDARY_NAME, BOUNDARY_VALUE, PETSC_FALSE, NULL);
364: IncrementNumLabels(&user);
365: if (user.compare_boundary) {
366: DMLabel label;
368: MarkBoundary(dm, BOUNDARY_VERTICES_NAME, BOUNDARY_VERTICES_VALUE, PETSC_TRUE, &label);
369: IncrementNumLabels(&user);
370: DMLabelGetCoordinateRepresentation(dm, label, BOUNDARY_VERTICES_VALUE, &allCoords);
371: DMLabelDestroy(&label);
372: }
373: SaveMesh(&user, dm);
375: LoadMesh(&user, &dmnew);
376: IncrementNumLabels(&user); /* account for depth label */
377: IncrementNumLabels(&user); /* account for celltype label */
378: CheckNumLabels(dm, &user);
379: CompareMeshes(&user, dm, dmnew);
380: if (user.compare_boundary) {
381: DMLabel label;
383: DMGetLabel(dmnew, BOUNDARY_VERTICES_NAME, &label);
384: DMLabelCompareWithCoordinateRepresentation(dmnew, label, BOUNDARY_VERTICES_VALUE, allCoords, user.verbose);
385: }
386: DMDestroy(&dm);
387: DMDestroy(&dmnew);
388: VecDestroy(&allCoords);
389: PetscFinalize();
390: return 0;
391: }
393: //TODO we can -compare once the new parallel topology format is in place
394: /*TEST
395: build:
396: requires: hdf5
398: # load old format, save in new format, reload, distribute
399: testset:
400: suffix: 1
401: requires: !complex datafilespath
402: args: -dm_plex_name plex
403: args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0
404: args: -dm_plex_interpolate
405: args: -load_dm_plex_check_all
406: args: -use_low_level_functions {{0 1}} -compare_boundary
407: args: -num_labels 1
408: args: -outfile ex56_1.h5
409: nsize: {{1 3}}
410: test:
411: suffix: a
412: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
413: test:
414: suffix: b
415: TODO: broken
416: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
417: test:
418: suffix: c
419: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
420: test:
421: suffix: d
422: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
423: test:
424: suffix: e
425: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
426: test:
427: suffix: f
428: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
430: # load old format, save in new format, reload topology, distribute, load geometry and labels
431: testset:
432: suffix: 2
433: requires: !complex datafilespath
434: args: -dm_plex_name plex
435: args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0
436: args: -dm_plex_interpolate
437: args: -load_dm_plex_check_all
438: args: -use_low_level_functions -load_dm_distribute 0 -distribute_after_topo_load -compare_boundary
439: args: -num_labels 1
440: args: -outfile ex56_2.h5
441: nsize: 3
442: test:
443: suffix: a
444: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
445: test:
446: suffix: b
447: TODO: broken
448: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
449: test:
450: suffix: c
451: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
452: test:
453: suffix: d
454: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
455: test:
456: suffix: e
457: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
458: test:
459: suffix: f
460: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
462: # load old format, save in new format, reload topology, distribute, load geometry and labels
463: testset:
464: suffix: 3
465: requires: !complex datafilespath
466: args: -dm_plex_name plex
467: args: -dm_plex_view_hdf5_storage_version 2.0.0
468: args: -dm_plex_interpolate -load_dm_distribute 0
469: args: -use_low_level_functions -compare_pre_post
470: args: -num_labels 1
471: args: -outfile ex56_3.h5
472: nsize: 3
473: test:
474: suffix: a
475: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
476: test:
477: suffix: b
478: TODO: broken
479: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
480: test:
481: suffix: c
482: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
483: test:
484: suffix: d
485: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
486: test:
487: suffix: e
488: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
489: test:
490: suffix: f
491: args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
492: TEST*/