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