1: static char help[] = "Reload a DMPLEX converted from a DMFOREST (p8est) from HDF5 in parallel and check its point SF.\n\
  2: A DMPLEX converted from a DMFOREST orders its strata cells, faces, edges, vertices instead of cells, vertices, faces, edges.\n\n";

  4: #include <petscdmplex.h>
  5: #include <petscdmforest.h>
  6: #include <petscviewerhdf5.h>
  7: #include <petscsf.h>

  9: int main(int argc, char **argv)
 10: {
 11:   DM          base, forest, plex, dm;
 12:   PetscViewer viewer;
 13:   PetscSF     sfXC;
 14:   PetscInt    faces[3] = {4, 4, 4};

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 19:   // Unrefined p8est forest over a box, converted to DMPLEX (strata ordered cells, faces, edges, vertices)
 20:   PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 3, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &base));
 21:   PetscCall(DMCreate(PETSC_COMM_WORLD, &forest));
 22:   PetscCall(DMSetType(forest, DMP8EST));
 23:   PetscCall(DMForestSetBaseDM(forest, base));
 24:   PetscCall(DMSetUp(forest));
 25:   PetscCall(DMConvert(forest, DMPLEX, &plex));
 26:   PetscCall(PetscObjectSetName((PetscObject)plex, "plex"));

 28:   // Save the topology, with storage version 3.1.0 (the default 1.0.0 loader reorders the strata to cells, vertices, faces, edges)
 29:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest_pointsf.h5", FILE_MODE_WRITE, &viewer));
 30:   {
 31:     struct _n_DMPlexStorageVersion version = {3, 1, 0};

 33:     PetscCall(PetscViewerHDF5SetDMPlexStorageVersionWriting(viewer, &version));
 34:   }
 35:   PetscCall(DMPlexTopologyView(plex, viewer));
 36:   PetscCall(PetscViewerDestroy(&viewer));

 38:   // Load it back; the loaded DM has no tree, so DMPlexCheckPointSF() performs its depth check
 39:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest_pointsf.h5", FILE_MODE_READ, &viewer));
 40:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 41:   PetscCall(DMSetType(dm, DMPLEX));
 42:   PetscCall(PetscObjectSetName((PetscObject)dm, "plex"));
 43:   PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXC));
 44:   PetscCall(PetscViewerDestroy(&viewer));

 46:   // Without order-independent depth lookup this fails on 2+ ranks with "Leaf .. has depth 2 but remote (.., ..) depth is 3"
 47:   PetscCall(DMPlexCheckPointSF(dm, NULL, PETSC_FALSE));
 48:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Point SF of the reloaded forest-ordered DMPLEX is consistent\n"));

 50:   PetscCall(PetscSFDestroy(&sfXC));
 51:   PetscCall(DMDestroy(&dm));
 52:   PetscCall(DMDestroy(&plex));
 53:   PetscCall(DMDestroy(&forest));
 54:   PetscCall(DMDestroy(&base));
 55:   PetscCall(PetscFinalize());
 56:   return 0;
 57: }

 59: /*TEST

 61:   build:
 62:     requires: p4est hdf5

 64:   test:
 65:     suffix: 0
 66:     nsize: 2

 68: TEST*/