xref: /petsc/src/dm/impls/plex/tests/ex56.c (revision 8112c1cbf372cb53bf7c5aca994a84a6a303db4d)
11d5c9e23SVaclav Hapla #include <petscdmplex.h>
21d5c9e23SVaclav Hapla #include <petscviewerhdf5.h>
31d5c9e23SVaclav Hapla #include <petscsf.h>
41d5c9e23SVaclav Hapla 
5e22c0602SVaclav Hapla static const char help[] = "Test DMLabel I/O with PETSc native HDF5 mesh format\n\n";
61d5c9e23SVaclav Hapla static const char EX[]   = "ex56.c";
71d5c9e23SVaclav Hapla typedef struct {
81d5c9e23SVaclav Hapla   MPI_Comm    comm;
91d5c9e23SVaclav Hapla   const char *meshname;                    /* Mesh name */
105c820da0SVaclav Hapla   PetscInt    num_labels;                  /* Asserted number of labels in loaded mesh */
111d5c9e23SVaclav Hapla   PetscBool   compare;                     /* Compare the meshes using DMPlexEqual() and DMCompareLabels() */
121d5c9e23SVaclav Hapla   PetscBool   compare_labels;              /* Compare labels in the meshes using DMCompareLabels() */
131d5c9e23SVaclav Hapla   PetscBool   compare_boundary;            /* Check label I/O via boundary vertex coordinates */
141d5c9e23SVaclav Hapla   PetscBool   compare_pre_post;            /* Compare labels loaded before distribution with those loaded after distribution */
151d5c9e23SVaclav Hapla   char        outfile[PETSC_MAX_PATH_LEN]; /* Output file */
161d5c9e23SVaclav Hapla   PetscBool   use_low_level_functions;     /* Use low level functions for viewing and loading */
171d5c9e23SVaclav Hapla   //TODO This is meant as temporary option; can be removed once we have full parallel loading in place
181d5c9e23SVaclav Hapla   PetscBool distribute_after_topo_load; /* Distribute topology right after DMPlexTopologyLoad(), if use_low_level_functions=true */
195c820da0SVaclav Hapla   PetscInt  verbose;
201d5c9e23SVaclav Hapla } AppCtx;
211d5c9e23SVaclav Hapla 
ProcessOptions(MPI_Comm comm,AppCtx * options)22d71ae5a4SJacob Faibussowitsch static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
23d71ae5a4SJacob Faibussowitsch {
241d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
251d5c9e23SVaclav Hapla   options->comm                       = comm;
265c820da0SVaclav Hapla   options->num_labels                 = -1;
271d5c9e23SVaclav Hapla   options->compare                    = PETSC_FALSE;
281d5c9e23SVaclav Hapla   options->compare_labels             = PETSC_FALSE;
291d5c9e23SVaclav Hapla   options->compare_boundary           = PETSC_FALSE;
301d5c9e23SVaclav Hapla   options->compare_pre_post           = PETSC_FALSE;
311d5c9e23SVaclav Hapla   options->outfile[0]                 = '\0';
321d5c9e23SVaclav Hapla   options->use_low_level_functions    = PETSC_FALSE;
331d5c9e23SVaclav Hapla   options->distribute_after_topo_load = PETSC_FALSE;
345c820da0SVaclav Hapla   options->verbose                    = 0;
351d5c9e23SVaclav Hapla 
36d0609cedSBarry Smith   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
375c820da0SVaclav Hapla   PetscCall(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));
389566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-compare", "Compare the meshes using DMPlexEqual() and DMCompareLabels()", EX, options->compare, &options->compare, NULL));
399566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-compare_labels", "Compare labels in the meshes using DMCompareLabels()", "ex55.c", options->compare_labels, &options->compare_labels, NULL));
409566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-compare_boundary", "Check label I/O via boundary vertex coordinates", "ex55.c", options->compare_boundary, &options->compare_boundary, NULL));
419566063dSJacob Faibussowitsch   PetscCall(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));
429566063dSJacob Faibussowitsch   PetscCall(PetscOptionsString("-outfile", "Output mesh file", EX, options->outfile, options->outfile, sizeof(options->outfile), NULL));
439566063dSJacob Faibussowitsch   PetscCall(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));
449566063dSJacob Faibussowitsch   PetscCall(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));
455c820da0SVaclav Hapla   PetscCall(PetscOptionsInt("-verbose", "Verbosity level", EX, options->verbose, &options->verbose, NULL));
46d0609cedSBarry Smith   PetscOptionsEnd();
473ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
48f4f49eeaSPierre Jolivet }
491d5c9e23SVaclav Hapla 
CreateMesh(AppCtx * options,DM * newdm)50d71ae5a4SJacob Faibussowitsch static PetscErrorCode CreateMesh(AppCtx *options, DM *newdm)
51d71ae5a4SJacob Faibussowitsch {
521d5c9e23SVaclav Hapla   DM dm;
531d5c9e23SVaclav Hapla 
541d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
559566063dSJacob Faibussowitsch   PetscCall(DMCreate(options->comm, &dm));
569566063dSJacob Faibussowitsch   PetscCall(DMSetType(dm, DMPLEX));
579566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(dm));
589566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetName((PetscObject)dm, &options->meshname));
599566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
601d5c9e23SVaclav Hapla   *newdm = dm;
613ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
621d5c9e23SVaclav Hapla }
631d5c9e23SVaclav Hapla 
SaveMesh(AppCtx * options,DM dm)64d71ae5a4SJacob Faibussowitsch static PetscErrorCode SaveMesh(AppCtx *options, DM dm)
65d71ae5a4SJacob Faibussowitsch {
661d5c9e23SVaclav Hapla   PetscViewer v;
671d5c9e23SVaclav Hapla 
681d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
699566063dSJacob Faibussowitsch   PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), options->outfile, FILE_MODE_WRITE, &v));
701d5c9e23SVaclav Hapla   if (options->use_low_level_functions) {
719566063dSJacob Faibussowitsch     PetscCall(DMPlexTopologyView(dm, v));
729566063dSJacob Faibussowitsch     PetscCall(DMPlexCoordinatesView(dm, v));
739566063dSJacob Faibussowitsch     PetscCall(DMPlexLabelsView(dm, v));
741d5c9e23SVaclav Hapla   } else {
759566063dSJacob Faibussowitsch     PetscCall(DMView(dm, v));
761d5c9e23SVaclav Hapla   }
779566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&v));
783ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
791d5c9e23SVaclav Hapla }
801d5c9e23SVaclav Hapla 
819371c9d4SSatish Balay typedef enum {
829371c9d4SSatish Balay   NONE      = 0,
839371c9d4SSatish Balay   PRE_DIST  = 1,
849371c9d4SSatish Balay   POST_DIST = 2
859371c9d4SSatish Balay } AuxObjLoadMode;
861d5c9e23SVaclav Hapla 
LoadMeshLowLevel(AppCtx * options,PetscViewer v,PetscBool explicitDistribute,AuxObjLoadMode mode,DM * newdm)87d71ae5a4SJacob Faibussowitsch static PetscErrorCode LoadMeshLowLevel(AppCtx *options, PetscViewer v, PetscBool explicitDistribute, AuxObjLoadMode mode, DM *newdm)
88d71ae5a4SJacob Faibussowitsch {
891d5c9e23SVaclav Hapla   DM      dm;
901d5c9e23SVaclav Hapla   PetscSF sfXC;
911d5c9e23SVaclav Hapla 
921d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
939566063dSJacob Faibussowitsch   PetscCall(DMCreate(options->comm, &dm));
949566063dSJacob Faibussowitsch   PetscCall(DMSetType(dm, DMPLEX));
959566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)dm, options->meshname));
969566063dSJacob Faibussowitsch   PetscCall(DMPlexTopologyLoad(dm, v, &sfXC));
971d5c9e23SVaclav Hapla   if (mode == PRE_DIST) {
989566063dSJacob Faibussowitsch     PetscCall(DMPlexCoordinatesLoad(dm, v, sfXC));
999566063dSJacob Faibussowitsch     PetscCall(DMPlexLabelsLoad(dm, v, sfXC));
1001d5c9e23SVaclav Hapla   }
1011d5c9e23SVaclav Hapla   if (explicitDistribute) {
1021d5c9e23SVaclav Hapla     DM      dmdist;
1031d5c9e23SVaclav Hapla     PetscSF sfXB = sfXC, sfBC;
1041d5c9e23SVaclav Hapla 
1059566063dSJacob Faibussowitsch     PetscCall(DMPlexDistribute(dm, 0, &sfBC, &dmdist));
1061d5c9e23SVaclav Hapla     if (dmdist) {
1071d5c9e23SVaclav Hapla       const char *name;
1081d5c9e23SVaclav Hapla 
1099566063dSJacob Faibussowitsch       PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1109566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)dmdist, name));
1119566063dSJacob Faibussowitsch       PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
1129566063dSJacob Faibussowitsch       PetscCall(PetscSFDestroy(&sfXB));
1139566063dSJacob Faibussowitsch       PetscCall(PetscSFDestroy(&sfBC));
1149566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm));
1151d5c9e23SVaclav Hapla       dm = dmdist;
1161d5c9e23SVaclav Hapla     }
1171d5c9e23SVaclav Hapla   }
1181d5c9e23SVaclav Hapla   if (mode == POST_DIST) {
1199566063dSJacob Faibussowitsch     PetscCall(DMPlexLabelsLoad(dm, v, sfXC));
12021027e53SStefano Zampini     PetscCall(DMPlexCoordinatesLoad(dm, v, sfXC));
1211d5c9e23SVaclav Hapla   }
1229566063dSJacob Faibussowitsch   PetscCall(PetscSFDestroy(&sfXC));
1231d5c9e23SVaclav Hapla   *newdm = dm;
1243ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1251d5c9e23SVaclav Hapla }
1261d5c9e23SVaclav Hapla 
LoadMesh(AppCtx * options,DM * dmnew)127d71ae5a4SJacob Faibussowitsch static PetscErrorCode LoadMesh(AppCtx *options, DM *dmnew)
128d71ae5a4SJacob Faibussowitsch {
1291d5c9e23SVaclav Hapla   DM          dm;
1301d5c9e23SVaclav Hapla   PetscViewer v;
1311d5c9e23SVaclav Hapla 
1321d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
1339566063dSJacob Faibussowitsch   PetscCall(PetscViewerHDF5Open(options->comm, options->outfile, FILE_MODE_READ, &v));
1341d5c9e23SVaclav Hapla   if (options->use_low_level_functions) {
1351d5c9e23SVaclav Hapla     if (options->compare_pre_post) {
1361d5c9e23SVaclav Hapla       DM dm0;
1371d5c9e23SVaclav Hapla 
1389566063dSJacob Faibussowitsch       PetscCall(LoadMeshLowLevel(options, v, PETSC_TRUE, PRE_DIST, &dm0));
1399566063dSJacob Faibussowitsch       PetscCall(LoadMeshLowLevel(options, v, PETSC_TRUE, POST_DIST, &dm));
1409566063dSJacob Faibussowitsch       PetscCall(DMCompareLabels(dm0, dm, NULL, NULL));
1419566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm0));
1421d5c9e23SVaclav Hapla     } else {
1439566063dSJacob Faibussowitsch       PetscCall(LoadMeshLowLevel(options, v, options->distribute_after_topo_load, POST_DIST, &dm));
1441d5c9e23SVaclav Hapla     }
1451d5c9e23SVaclav Hapla   } else {
1469566063dSJacob Faibussowitsch     PetscCall(DMCreate(options->comm, &dm));
1479566063dSJacob Faibussowitsch     PetscCall(DMSetType(dm, DMPLEX));
1489566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)dm, options->meshname));
1499566063dSJacob Faibussowitsch     PetscCall(DMLoad(dm, v));
1501d5c9e23SVaclav Hapla   }
1519566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&v));
15207dc484cSMatthew G. Knepley   DMLabel celltypes;
15307dc484cSMatthew G. Knepley   PetscCall(DMPlexGetCellTypeLabel(dm, &celltypes));
1541d5c9e23SVaclav Hapla 
1559566063dSJacob Faibussowitsch   PetscCall(DMSetOptionsPrefix(dm, "load_"));
1569566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(dm));
1579566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1581d5c9e23SVaclav Hapla   *dmnew = dm;
1593ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1601d5c9e23SVaclav Hapla }
1611d5c9e23SVaclav Hapla 
CompareMeshes(AppCtx * options,DM dm0,DM dm1)162d71ae5a4SJacob Faibussowitsch static PetscErrorCode CompareMeshes(AppCtx *options, DM dm0, DM dm1)
163d71ae5a4SJacob Faibussowitsch {
1641d5c9e23SVaclav Hapla   PetscBool flg;
1651d5c9e23SVaclav Hapla 
1661d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
1671d5c9e23SVaclav Hapla   if (options->compare) {
1689566063dSJacob Faibussowitsch     PetscCall(DMPlexEqual(dm0, dm1, &flg));
16928b400f6SJacob Faibussowitsch     PetscCheck(flg, options->comm, PETSC_ERR_ARG_INCOMP, "DMs are not equal");
1709566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(options->comm, "DMs equal\n"));
1711d5c9e23SVaclav Hapla   }
1721d5c9e23SVaclav Hapla   if (options->compare_labels) {
1739566063dSJacob Faibussowitsch     PetscCall(DMCompareLabels(dm0, dm1, NULL, NULL));
1749566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(options->comm, "DMLabels equal\n"));
1751d5c9e23SVaclav Hapla   }
1763ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1771d5c9e23SVaclav Hapla }
1781d5c9e23SVaclav Hapla 
MarkBoundary(DM dm,const char name[],PetscInt value,PetscBool verticesOnly,DMLabel * label)179d71ae5a4SJacob Faibussowitsch static PetscErrorCode MarkBoundary(DM dm, const char name[], PetscInt value, PetscBool verticesOnly, DMLabel *label)
180d71ae5a4SJacob Faibussowitsch {
1811d5c9e23SVaclav Hapla   DMLabel l;
1825c820da0SVaclav Hapla 
1835c820da0SVaclav Hapla   PetscFunctionBeginUser;
1845c820da0SVaclav Hapla   PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)dm), name, &l));
1855c820da0SVaclav Hapla   PetscCall(DMAddLabel(dm, l));
1865c820da0SVaclav Hapla   PetscCall(DMPlexMarkBoundaryFaces(dm, value, l));
1875c820da0SVaclav Hapla   PetscCall(DMPlexLabelComplete(dm, l));
1885c820da0SVaclav Hapla   if (verticesOnly) {
1891d5c9e23SVaclav Hapla     IS              points;
1901d5c9e23SVaclav Hapla     const PetscInt *idx;
191704d4d0eSMatthew G. Knepley     PetscInt        i, n = 0;
1921d5c9e23SVaclav Hapla 
1939566063dSJacob Faibussowitsch     PetscCall(DMLabelGetStratumIS(l, value, &points));
194704d4d0eSMatthew G. Knepley     if (points) {
1959566063dSJacob Faibussowitsch       PetscCall(ISGetLocalSize(points, &n));
1969566063dSJacob Faibussowitsch       PetscCall(ISGetIndices(points, &idx));
197704d4d0eSMatthew G. Knepley     }
1981d5c9e23SVaclav Hapla     for (i = 0; i < n; i++) {
1991d5c9e23SVaclav Hapla       const PetscInt p = idx[i];
2001d5c9e23SVaclav Hapla       PetscInt       d;
2011d5c9e23SVaclav Hapla 
2029566063dSJacob Faibussowitsch       PetscCall(DMPlexGetPointDepth(dm, p, &d));
20348a46eb9SPierre Jolivet       if (d != 0) PetscCall(DMLabelClearValue(l, p, value));
2041d5c9e23SVaclav Hapla     }
205704d4d0eSMatthew G. Knepley     if (points) PetscCall(ISRestoreIndices(points, &idx));
2069566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&points));
2075c820da0SVaclav Hapla   }
2085c820da0SVaclav Hapla   if (label) *label = l;
2095c820da0SVaclav Hapla   else PetscCall(DMLabelDestroy(&l));
2103ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2111d5c9e23SVaclav Hapla }
2121d5c9e23SVaclav Hapla 
VertexCoordinatesToAll(DM dm,IS vertices,Vec * allCoords)213d71ae5a4SJacob Faibussowitsch static PetscErrorCode VertexCoordinatesToAll(DM dm, IS vertices, Vec *allCoords)
214d71ae5a4SJacob Faibussowitsch {
2151d5c9e23SVaclav Hapla   Vec        coords, allCoords_;
2161d5c9e23SVaclav Hapla   VecScatter sc;
2171d5c9e23SVaclav Hapla   MPI_Comm   comm;
2181d5c9e23SVaclav Hapla 
2191d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
2209566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2219566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocalSetUp(dm));
2221d5c9e23SVaclav Hapla   if (vertices) {
2239566063dSJacob Faibussowitsch     PetscCall(DMGetCoordinatesLocalTuple(dm, vertices, NULL, &coords));
2241d5c9e23SVaclav Hapla   } else {
22577433607SBarry Smith     PetscCall(VecCreateFromOptions(PETSC_COMM_SELF, NULL, 1, 0, 0, &coords));
2261d5c9e23SVaclav Hapla   }
2271d5c9e23SVaclav Hapla   {
2281d5c9e23SVaclav Hapla     PetscInt n;
2291d5c9e23SVaclav Hapla     Vec      mpivec;
2301d5c9e23SVaclav Hapla 
2319566063dSJacob Faibussowitsch     PetscCall(VecGetLocalSize(coords, &n));
23277433607SBarry Smith     PetscCall(VecCreateFromOptions(comm, NULL, 1, n, PETSC_DECIDE, &mpivec));
2339566063dSJacob Faibussowitsch     PetscCall(VecCopy(coords, mpivec));
2349566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&coords));
2351d5c9e23SVaclav Hapla     coords = mpivec;
2361d5c9e23SVaclav Hapla   }
2371d5c9e23SVaclav Hapla 
2389566063dSJacob Faibussowitsch   PetscCall(VecScatterCreateToAll(coords, &sc, &allCoords_));
2399566063dSJacob Faibussowitsch   PetscCall(VecScatterBegin(sc, coords, allCoords_, INSERT_VALUES, SCATTER_FORWARD));
2409566063dSJacob Faibussowitsch   PetscCall(VecScatterEnd(sc, coords, allCoords_, INSERT_VALUES, SCATTER_FORWARD));
2419566063dSJacob Faibussowitsch   PetscCall(VecScatterDestroy(&sc));
2429566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&coords));
2431d5c9e23SVaclav Hapla   *allCoords = allCoords_;
2443ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2451d5c9e23SVaclav Hapla }
2461d5c9e23SVaclav Hapla 
DMLabelGetCoordinateRepresentation(DM dm,DMLabel label,PetscInt value,Vec * allCoords)247d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMLabelGetCoordinateRepresentation(DM dm, DMLabel label, PetscInt value, Vec *allCoords)
248d71ae5a4SJacob Faibussowitsch {
2491d5c9e23SVaclav Hapla   IS vertices;
2501d5c9e23SVaclav Hapla 
2511d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
2525c820da0SVaclav Hapla   PetscCall(DMLabelGetStratumIS(label, value, &vertices));
2539566063dSJacob Faibussowitsch   PetscCall(VertexCoordinatesToAll(dm, vertices, allCoords));
2549566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&vertices));
2553ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
2561d5c9e23SVaclav Hapla }
2571d5c9e23SVaclav Hapla 
DMLabelCompareWithCoordinateRepresentation(DM dm,DMLabel label,PetscInt value,Vec allCoords,PetscInt verbose)258d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMLabelCompareWithCoordinateRepresentation(DM dm, DMLabel label, PetscInt value, Vec allCoords, PetscInt verbose)
259d71ae5a4SJacob Faibussowitsch {
2605c820da0SVaclav Hapla   const char     *labelName;
2611d5c9e23SVaclav Hapla   IS              pointsIS;
2621d5c9e23SVaclav Hapla   const PetscInt *points;
2631d5c9e23SVaclav Hapla   PetscInt        i, n;
2641d5c9e23SVaclav Hapla   PetscBool       fail = PETSC_FALSE;
2651d5c9e23SVaclav Hapla   MPI_Comm        comm;
2661d5c9e23SVaclav Hapla   PetscMPIInt     rank;
2671d5c9e23SVaclav Hapla 
2681d5c9e23SVaclav Hapla   PetscFunctionBeginUser;
2695c820da0SVaclav Hapla   PetscCheck(label, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Label not loaded");
2709566063dSJacob Faibussowitsch   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2715c820da0SVaclav Hapla   PetscCall(PetscObjectGetName((PetscObject)label, &labelName));
2729566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2731d5c9e23SVaclav Hapla   {
2741d5c9e23SVaclav Hapla     PetscInt pStart, pEnd;
2751d5c9e23SVaclav Hapla 
2769566063dSJacob Faibussowitsch     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2779566063dSJacob Faibussowitsch     PetscCall(DMLabelCreateIndex(label, pStart, pEnd));
2781d5c9e23SVaclav Hapla   }
2799566063dSJacob Faibussowitsch   PetscCall(DMPlexFindVertices(dm, allCoords, 0.0, &pointsIS));
2809566063dSJacob Faibussowitsch   PetscCall(ISGetIndices(pointsIS, &points));
2819566063dSJacob Faibussowitsch   PetscCall(ISGetLocalSize(pointsIS, &n));
2825c820da0SVaclav Hapla   if (verbose > 1) PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm)));
2831d5c9e23SVaclav Hapla   for (i = 0; i < n; i++) {
2841d5c9e23SVaclav Hapla     const PetscInt p = points[i];
2851d5c9e23SVaclav Hapla     PetscBool      has;
2861d5c9e23SVaclav Hapla     PetscInt       v;
2871d5c9e23SVaclav Hapla 
2881d5c9e23SVaclav Hapla     if (p < 0) continue;
2899566063dSJacob Faibussowitsch     PetscCall(DMLabelHasPoint(label, p, &has));
2901d5c9e23SVaclav Hapla     if (!has) {
29163a3b9bcSJacob Faibussowitsch       if (verbose) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Label \"%s\" does not have point %" PetscInt_FMT "\n", rank, labelName, p));
2921d5c9e23SVaclav Hapla       fail = PETSC_TRUE;
2931d5c9e23SVaclav Hapla       continue;
2941d5c9e23SVaclav Hapla     }
2959566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValue(label, p, &v));
2965c820da0SVaclav Hapla     if (v != value) {
29763a3b9bcSJacob Faibussowitsch       if (verbose) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "Point %" PetscInt_FMT " has bad value %" PetscInt_FMT " in label \"%s\"", p, v, labelName));
2981d5c9e23SVaclav Hapla       fail = PETSC_TRUE;
2991d5c9e23SVaclav Hapla       continue;
3001d5c9e23SVaclav Hapla     }
30163a3b9bcSJacob Faibussowitsch     if (verbose > 1) PetscCall(PetscSynchronizedPrintf(comm, "[%d] OK point %" PetscInt_FMT "\n", rank, p));
3021d5c9e23SVaclav Hapla   }
3039566063dSJacob Faibussowitsch   PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
3049566063dSJacob Faibussowitsch   PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
3059566063dSJacob Faibussowitsch   PetscCall(ISRestoreIndices(pointsIS, &points));
3069566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&pointsIS));
307*5440e5dcSBarry Smith   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &fail, 1, MPI_C_BOOL, MPI_LOR, comm));
3085c820da0SVaclav Hapla   PetscCheck(!fail, comm, PETSC_ERR_PLIB, "Label \"%s\" was not loaded correctly%s", labelName, verbose ? " - see details above" : "");
3093ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3105c820da0SVaclav Hapla }
3115c820da0SVaclav Hapla 
CheckNumLabels(DM dm,AppCtx * ctx)312d71ae5a4SJacob Faibussowitsch static PetscErrorCode CheckNumLabels(DM dm, AppCtx *ctx)
313d71ae5a4SJacob Faibussowitsch {
3145c820da0SVaclav Hapla   PetscInt    actualNum;
3155c820da0SVaclav Hapla   PetscBool   fail = PETSC_FALSE;
3165c820da0SVaclav Hapla   MPI_Comm    comm;
3175c820da0SVaclav Hapla   PetscMPIInt rank;
3185c820da0SVaclav Hapla 
3195c820da0SVaclav Hapla   PetscFunctionBeginUser;
3203ba16761SJacob Faibussowitsch   if (ctx->num_labels < 0) PetscFunctionReturn(PETSC_SUCCESS);
3215c820da0SVaclav Hapla   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3225c820da0SVaclav Hapla   PetscCallMPI(MPI_Comm_rank(comm, &rank));
3235c820da0SVaclav Hapla   PetscCall(DMGetNumLabels(dm, &actualNum));
3245c820da0SVaclav Hapla   if (ctx->num_labels != actualNum) {
3255c820da0SVaclav Hapla     fail = PETSC_TRUE;
3265c820da0SVaclav Hapla     if (ctx->verbose) {
3275c820da0SVaclav Hapla       PetscInt i;
3285c820da0SVaclav Hapla 
32963a3b9bcSJacob Faibussowitsch       PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Asserted number of labels: %" PetscInt_FMT ", actual: %" PetscInt_FMT "\n", rank, ctx->num_labels, actualNum));
3305c820da0SVaclav Hapla       for (i = 0; i < actualNum; i++) {
3315c820da0SVaclav Hapla         DMLabel     label;
3325c820da0SVaclav Hapla         const char *name;
3335c820da0SVaclav Hapla 
3345c820da0SVaclav Hapla         PetscCall(DMGetLabelByNum(dm, i, &label));
3355c820da0SVaclav Hapla         PetscCall(PetscObjectGetName((PetscObject)label, &name));
33663a3b9bcSJacob Faibussowitsch         PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Label %" PetscInt_FMT " \"%s\"\n", rank, i, name));
3375c820da0SVaclav Hapla       }
3385c820da0SVaclav Hapla       PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
3395c820da0SVaclav Hapla     }
3405c820da0SVaclav Hapla   }
341*5440e5dcSBarry Smith   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &fail, 1, MPI_C_BOOL, MPI_LOR, comm));
3425c820da0SVaclav Hapla   PetscCheck(!fail, comm, PETSC_ERR_PLIB, "Wrong number of labels%s", ctx->verbose ? " - see details above" : "");
3433ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3445c820da0SVaclav Hapla }
3455c820da0SVaclav Hapla 
IncrementNumLabels(AppCtx * ctx)346d71ae5a4SJacob Faibussowitsch static inline PetscErrorCode IncrementNumLabels(AppCtx *ctx)
347d71ae5a4SJacob Faibussowitsch {
3485c820da0SVaclav Hapla   PetscFunctionBeginUser;
3495c820da0SVaclav Hapla   if (ctx->num_labels >= 0) ctx->num_labels++;
3503ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3511d5c9e23SVaclav Hapla }
3521d5c9e23SVaclav Hapla 
main(int argc,char ** argv)353d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
354d71ae5a4SJacob Faibussowitsch {
3555c820da0SVaclav Hapla   const char     BOUNDARY_NAME[]          = "Boundary";
3565c820da0SVaclav Hapla   const char     BOUNDARY_VERTICES_NAME[] = "BoundaryVertices";
3575c820da0SVaclav Hapla   const PetscInt BOUNDARY_VALUE           = 12345;
3585c820da0SVaclav Hapla   const PetscInt BOUNDARY_VERTICES_VALUE  = 6789;
3591d5c9e23SVaclav Hapla   DM             dm, dmnew;
3601d5c9e23SVaclav Hapla   AppCtx         user;
3611d5c9e23SVaclav Hapla   Vec            allCoords = NULL;
3621d5c9e23SVaclav Hapla 
363327415f7SBarry Smith   PetscFunctionBeginUser;
3649566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
3659566063dSJacob Faibussowitsch   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
3669566063dSJacob Faibussowitsch   PetscCall(CreateMesh(&user, &dm));
3675c820da0SVaclav Hapla   PetscCall(MarkBoundary(dm, BOUNDARY_NAME, BOUNDARY_VALUE, PETSC_FALSE, NULL));
3685c820da0SVaclav Hapla   PetscCall(IncrementNumLabels(&user));
3691d5c9e23SVaclav Hapla   if (user.compare_boundary) {
3705c820da0SVaclav Hapla     DMLabel label;
3715c820da0SVaclav Hapla 
3725c820da0SVaclav Hapla     PetscCall(MarkBoundary(dm, BOUNDARY_VERTICES_NAME, BOUNDARY_VERTICES_VALUE, PETSC_TRUE, &label));
3735c820da0SVaclav Hapla     PetscCall(IncrementNumLabels(&user));
3745c820da0SVaclav Hapla     PetscCall(DMLabelGetCoordinateRepresentation(dm, label, BOUNDARY_VERTICES_VALUE, &allCoords));
3755c820da0SVaclav Hapla     PetscCall(DMLabelDestroy(&label));
3761d5c9e23SVaclav Hapla   }
3779566063dSJacob Faibussowitsch   PetscCall(SaveMesh(&user, dm));
3785c820da0SVaclav Hapla 
3799566063dSJacob Faibussowitsch   PetscCall(LoadMesh(&user, &dmnew));
3805c820da0SVaclav Hapla   PetscCall(IncrementNumLabels(&user)); /* account for depth label */
3815c820da0SVaclav Hapla   PetscCall(IncrementNumLabels(&user)); /* account for celltype label */
3825c820da0SVaclav Hapla   PetscCall(CheckNumLabels(dm, &user));
3839566063dSJacob Faibussowitsch   PetscCall(CompareMeshes(&user, dm, dmnew));
3841d5c9e23SVaclav Hapla   if (user.compare_boundary) {
3855c820da0SVaclav Hapla     DMLabel label;
3865c820da0SVaclav Hapla 
3875c820da0SVaclav Hapla     PetscCall(DMGetLabel(dmnew, BOUNDARY_VERTICES_NAME, &label));
3885c820da0SVaclav Hapla     PetscCall(DMLabelCompareWithCoordinateRepresentation(dmnew, label, BOUNDARY_VERTICES_VALUE, allCoords, user.verbose));
3891d5c9e23SVaclav Hapla   }
3909566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
3919566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmnew));
3929566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&allCoords));
3939566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
394b122ec5aSJacob Faibussowitsch   return 0;
3951d5c9e23SVaclav Hapla }
3961d5c9e23SVaclav Hapla 
3971d5c9e23SVaclav Hapla //TODO we can -compare once the new parallel topology format is in place
3981d5c9e23SVaclav Hapla /*TEST
3991d5c9e23SVaclav Hapla   build:
4001d5c9e23SVaclav Hapla     requires: hdf5
4011d5c9e23SVaclav Hapla 
4021d5c9e23SVaclav Hapla   # load old format, save in new format, reload, distribute
4031d5c9e23SVaclav Hapla   testset:
4041d5c9e23SVaclav Hapla     suffix: 1
4051d5c9e23SVaclav Hapla     requires: !complex datafilespath
4061d5c9e23SVaclav Hapla     args: -dm_plex_name plex
4071d5c9e23SVaclav Hapla     args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0
408704d4d0eSMatthew G. Knepley     args: -dm_plex_interpolate -petscpartitioner_type simple
4091d5c9e23SVaclav Hapla     args: -load_dm_plex_check_all
410e600fa54SMatthew G. Knepley     args: -use_low_level_functions {{0 1}} -compare_boundary
4115c820da0SVaclav Hapla     args: -num_labels 1
4121d5c9e23SVaclav Hapla     args: -outfile ex56_1.h5
4131d5c9e23SVaclav Hapla     nsize: {{1 3}}
4143886731fSPierre Jolivet     output_file: output/empty.out
4151d5c9e23SVaclav Hapla     test:
4161d5c9e23SVaclav Hapla       suffix: a
4171d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
4181d5c9e23SVaclav Hapla     test:
4191d5c9e23SVaclav Hapla       suffix: b
4201d5c9e23SVaclav Hapla       TODO: broken
4211d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
4221d5c9e23SVaclav Hapla     test:
4231d5c9e23SVaclav Hapla       suffix: c
4241d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
4251d5c9e23SVaclav Hapla     test:
4261d5c9e23SVaclav Hapla       suffix: d
4275c820da0SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
4281d5c9e23SVaclav Hapla     test:
4291d5c9e23SVaclav Hapla       suffix: e
4301d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
4311d5c9e23SVaclav Hapla     test:
4321d5c9e23SVaclav Hapla       suffix: f
4331d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
4341d5c9e23SVaclav Hapla 
4351d5c9e23SVaclav Hapla   # load old format, save in new format, reload topology, distribute, load geometry and labels
4361d5c9e23SVaclav Hapla   testset:
4371d5c9e23SVaclav Hapla     suffix: 2
4381d5c9e23SVaclav Hapla     requires: !complex datafilespath
4391d5c9e23SVaclav Hapla     args: -dm_plex_name plex
4401d5c9e23SVaclav Hapla     args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0
441704d4d0eSMatthew G. Knepley     args: -dm_plex_interpolate -petscpartitioner_type simple
4421d5c9e23SVaclav Hapla     args: -load_dm_plex_check_all
443e600fa54SMatthew G. Knepley     args: -use_low_level_functions -load_dm_distribute 0 -distribute_after_topo_load -compare_boundary
4445c820da0SVaclav Hapla     args: -num_labels 1
4451d5c9e23SVaclav Hapla     args: -outfile ex56_2.h5
4461d5c9e23SVaclav Hapla     nsize: 3
4473886731fSPierre Jolivet     output_file: output/empty.out
4481d5c9e23SVaclav Hapla     test:
4491d5c9e23SVaclav Hapla       suffix: a
4501d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
4511d5c9e23SVaclav Hapla     test:
4521d5c9e23SVaclav Hapla       suffix: b
4531d5c9e23SVaclav Hapla       TODO: broken
4541d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
4551d5c9e23SVaclav Hapla     test:
4561d5c9e23SVaclav Hapla       suffix: c
4571d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
4581d5c9e23SVaclav Hapla     test:
4591d5c9e23SVaclav Hapla       suffix: d
4605c820da0SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
4611d5c9e23SVaclav Hapla     test:
4621d5c9e23SVaclav Hapla       suffix: e
4631d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
4641d5c9e23SVaclav Hapla     test:
4651d5c9e23SVaclav Hapla       suffix: f
4661d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
4671d5c9e23SVaclav Hapla 
4681d5c9e23SVaclav Hapla   # load old format, save in new format, reload topology, distribute, load geometry and labels
4691d5c9e23SVaclav Hapla   testset:
4701d5c9e23SVaclav Hapla     suffix: 3
4711d5c9e23SVaclav Hapla     requires: !complex datafilespath
4721d5c9e23SVaclav Hapla     args: -dm_plex_name plex
4731d5c9e23SVaclav Hapla     args: -dm_plex_view_hdf5_storage_version 2.0.0
474704d4d0eSMatthew G. Knepley     args: -dm_plex_interpolate -load_dm_distribute 0 -petscpartitioner_type simple
4751d5c9e23SVaclav Hapla     args: -use_low_level_functions -compare_pre_post
4765c820da0SVaclav Hapla     args: -num_labels 1
4771d5c9e23SVaclav Hapla     args: -outfile ex56_3.h5
4781d5c9e23SVaclav Hapla     nsize: 3
4793886731fSPierre Jolivet     output_file: output/empty.out
4801d5c9e23SVaclav Hapla     test:
4811d5c9e23SVaclav Hapla       suffix: a
4821d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5
4831d5c9e23SVaclav Hapla     test:
4841d5c9e23SVaclav Hapla       suffix: b
4851d5c9e23SVaclav Hapla       TODO: broken
4861d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5
4871d5c9e23SVaclav Hapla     test:
4881d5c9e23SVaclav Hapla       suffix: c
4891d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5
4901d5c9e23SVaclav Hapla     test:
4911d5c9e23SVaclav Hapla       suffix: d
4925c820da0SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 -num_labels 0
4931d5c9e23SVaclav Hapla     test:
4941d5c9e23SVaclav Hapla       suffix: e
4951d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5
4961d5c9e23SVaclav Hapla     test:
4971d5c9e23SVaclav Hapla       suffix: f
4981d5c9e23SVaclav Hapla       args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5
4991d5c9e23SVaclav Hapla TEST*/
500