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*/