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