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