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;

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

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

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

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

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

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

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

 62:   VecDestroy(&g);
 63:   VecDestroy(&g2);
 64:   DMDestroy(&forest);
 65:   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*/