1 static char help[] = "Test HDF5 input and output.\n\
2 This exposed a bug with sharing discretizations.\n\n\n";
3
4 #include <petscdmforest.h>
5 #include <petscdmplex.h>
6 #include <petscviewerhdf5.h>
7
main(int argc,char ** argv)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;
17
18 PetscFunctionBeginUser;
19 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
20
21 PetscCall(DMCreate(PETSC_COMM_WORLD, &base));
22 PetscCall(DMSetType(base, DMPLEX));
23 PetscCall(DMSetFromOptions(base));
24
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"));
33
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));
43
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));
50
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));
56
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));
61
62 PetscCall(VecDestroy(&g));
63 PetscCall(VecDestroy(&g2));
64 PetscCall(DMDestroy(&forest));
65 PetscCall(PetscFinalize());
66 return 0;
67 }
68
69 /*TEST
70
71 build:
72 requires: hdf5 p4est
73
74 test:
75 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3
76 output_file: output/empty.out
77
78 TEST*/
79