Actual source code: ex1.c

  1: static char help[] = "Test HDF5 input and output.\n\
  2: This exposed a bug with sharing discretizations.\n\n\n";

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

  8: int main(int argc, char **argv)
  9: {
 10:   DM           base, forest, plex;
 11:   Vec          g, g2;
 12:   PetscSection s;
 13:   PetscViewer  viewer;
 14:   PetscReal    diff;
 15:   PetscInt     min_refine = 2, overlap = 0;
 16:   PetscInt     vStart, vEnd, v;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 21:   PetscCall(DMCreate(PETSC_COMM_WORLD, &base));
 22:   PetscCall(DMSetType(base, DMPLEX));
 23:   PetscCall(DMSetFromOptions(base));

 25:   PetscCall(DMCreate(PETSC_COMM_WORLD, &forest));
 26:   PetscCall(DMSetType(forest, DMP4EST));
 27:   PetscCall(DMForestSetBaseDM(forest, base));
 28:   PetscCall(DMForestSetInitialRefinement(forest, min_refine));
 29:   PetscCall(DMForestSetPartitionOverlap(forest, overlap));
 30:   PetscCall(DMSetUp(forest));
 31:   PetscCall(DMDestroy(&base));
 32:   PetscCall(DMViewFromOptions(forest, NULL, "-dm_view"));

 34:   PetscCall(DMConvert(forest, DMPLEX, &plex));
 35:   PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
 36:   PetscCall(DMDestroy(&plex));
 37:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)forest), &s));
 38:   PetscCall(PetscSectionSetChart(s, vStart, vEnd));
 39:   for (v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(s, v, 1));
 40:   PetscCall(PetscSectionSetUp(s));
 41:   PetscCall(DMSetLocalSection(forest, s));
 42:   PetscCall(PetscSectionDestroy(&s));

 44:   PetscCall(DMCreateGlobalVector(forest, &g));
 45:   PetscCall(PetscObjectSetName((PetscObject)g, "g"));
 46:   PetscCall(VecSet(g, 1.0));
 47:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_WRITE, &viewer));
 48:   PetscCall(VecView(g, viewer));
 49:   PetscCall(PetscViewerDestroy(&viewer));

 51:   PetscCall(DMCreateGlobalVector(forest, &g2));
 52:   PetscCall(PetscObjectSetName((PetscObject)g2, "g"));
 53:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_READ, &viewer));
 54:   PetscCall(VecLoad(g2, viewer));
 55:   PetscCall(PetscViewerDestroy(&viewer));

 57:   /*  Check if the data is the same*/
 58:   PetscCall(VecAXPY(g2, -1.0, g));
 59:   PetscCall(VecNorm(g2, NORM_INFINITY, &diff));
 60:   if (diff > PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Check failed: %g\n", (double)diff));

 62:   PetscCall(VecDestroy(&g));
 63:   PetscCall(VecDestroy(&g2));
 64:   PetscCall(DMDestroy(&forest));
 65:   PetscCall(PetscFinalize());
 66:   return 0;
 67: }

 69: /*TEST

 71:   build:
 72:     requires: hdf5 p4est

 74:   test:
 75:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3

 77: TEST*/