1 static char help[] = "Tests adaptive refinement using DMForest, and uses HDF5.\n\n";
2
3 #include <petscdmforest.h>
4 #include <petscdmplex.h>
5 #include <petscviewerhdf5.h>
6
main(int argc,char ** argv)7 int main(int argc, char **argv)
8 {
9 DM base, forest, plex;
10 PetscSection s;
11 PetscViewer viewer;
12 Vec g = NULL, g2 = NULL;
13 PetscReal nrm;
14 PetscBool adapt = PETSC_FALSE, userSection = PETSC_FALSE;
15 PetscInt vStart, vEnd, v, i;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
19 PetscCall(PetscOptionsGetBool(NULL, NULL, "-adapt", &adapt, NULL));
20 PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_section", &userSection, NULL));
21
22 /* Create a base DMPlex mesh */
23 PetscCall(DMCreate(PETSC_COMM_WORLD, &base));
24 PetscCall(DMSetType(base, DMPLEX));
25 PetscCall(DMSetFromOptions(base));
26 PetscCall(DMViewFromOptions(base, NULL, "-dm_view"));
27
28 /* Convert Plex mesh to Forest and destroy base */
29 PetscCall(DMCreate(PETSC_COMM_WORLD, &forest));
30 PetscCall(DMSetType(forest, DMP4EST));
31 PetscCall(DMForestSetBaseDM(forest, base));
32 PetscCall(DMSetUp(forest));
33 PetscCall(DMDestroy(&base));
34 PetscCall(DMViewFromOptions(forest, NULL, "-dm_view"));
35
36 if (adapt) {
37 /* Adaptively refine the cell 0 of the mesh */
38 for (i = 0; i < 3; ++i) {
39 DM postforest;
40 DMLabel adaptLabel = NULL;
41
42 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
43 PetscCall(DMLabelSetValue(adaptLabel, 0, DM_ADAPT_REFINE));
44 PetscCall(DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest));
45 PetscCall(DMForestSetAdaptivityLabel(postforest, adaptLabel));
46 PetscCall(DMLabelDestroy(&adaptLabel));
47 PetscCall(DMSetUp(postforest));
48 PetscCall(DMDestroy(&forest));
49 forest = postforest;
50 }
51 } else {
52 /* Adaptively refine all cells of the mesh */
53 PetscInt cStart, cEnd, c;
54
55 for (i = 0; i < 3; ++i) {
56 DM postforest;
57 DMLabel adaptLabel = NULL;
58
59 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
60
61 PetscCall(DMForestGetCellChart(forest, &cStart, &cEnd));
62 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));
63
64 PetscCall(DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest));
65 PetscCall(DMForestSetAdaptivityLabel(postforest, adaptLabel));
66 PetscCall(DMLabelDestroy(&adaptLabel));
67 PetscCall(DMSetUp(postforest));
68 PetscCall(DMDestroy(&forest));
69 forest = postforest;
70 }
71 }
72 PetscCall(DMViewFromOptions(forest, NULL, "-dm_view"));
73
74 /* Setup the section*/
75 if (userSection) {
76 PetscCall(DMConvert(forest, DMPLEX, &plex));
77 PetscCall(DMViewFromOptions(plex, NULL, "-dm_view"));
78 PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd));
79 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Vertices [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", vStart, vEnd));
80 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
81 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)forest), &s));
82 PetscCall(PetscSectionSetNumFields(s, 1));
83 PetscCall(PetscSectionSetChart(s, vStart, vEnd));
84 for (v = vStart; v < vEnd; ++v) {
85 PetscCall(PetscSectionSetDof(s, v, 1));
86 PetscCall(PetscSectionSetFieldDof(s, v, 0, 1));
87 }
88 PetscCall(PetscSectionSetUp(s));
89 PetscCall(DMSetLocalSection(forest, s));
90 PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-my_section_view"));
91 PetscCall(PetscSectionDestroy(&s));
92 PetscCall(DMDestroy(&plex));
93 } else {
94 PetscFE fe;
95 PetscInt dim;
96
97 PetscCall(DMGetDimension(forest, &dim));
98 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 1, PETSC_DETERMINE, &fe));
99 PetscCall(DMAddField(forest, NULL, (PetscObject)fe));
100 PetscCall(PetscFEDestroy(&fe));
101 PetscCall(DMCreateDS(forest));
102 }
103
104 /* Create the global vector*/
105 PetscCall(DMCreateGlobalVector(forest, &g));
106 PetscCall(PetscObjectSetName((PetscObject)g, "g"));
107 PetscCall(VecSet(g, 1.0));
108
109 /* Test global to local*/
110 Vec l;
111 PetscCall(DMCreateLocalVector(forest, &l));
112 PetscCall(VecZeroEntries(l));
113 PetscCall(DMGlobalToLocal(forest, g, INSERT_VALUES, l));
114 PetscCall(VecZeroEntries(g));
115 PetscCall(DMLocalToGlobal(forest, l, INSERT_VALUES, g));
116 PetscCall(VecDestroy(&l));
117
118 /* Save a vector*/
119 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_WRITE, &viewer));
120 PetscCall(VecView(g, viewer));
121 PetscCall(PetscViewerDestroy(&viewer));
122
123 /* Load another vector to load into*/
124 PetscCall(DMCreateGlobalVector(forest, &g2));
125 PetscCall(PetscObjectSetName((PetscObject)g2, "g"));
126 PetscCall(VecZeroEntries(g2));
127
128 /* Load a vector*/
129 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_READ, &viewer));
130 PetscCall(VecLoad(g2, viewer));
131 PetscCall(PetscViewerDestroy(&viewer));
132
133 /* Check if the data is the same*/
134 PetscCall(VecAXPY(g2, -1.0, g));
135 PetscCall(VecNorm(g2, NORM_INFINITY, &nrm));
136 PetscCheck(PetscAbsReal(nrm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Invalid difference norm %g", (double)nrm);
137
138 PetscCall(VecDestroy(&g));
139 PetscCall(VecDestroy(&g2));
140 PetscCall(DMDestroy(&forest));
141 PetscCall(PetscFinalize());
142 return 0;
143 }
144
145 /*TEST
146
147 build:
148 requires: hdf5 p4est
149
150 test:
151 suffix: 0
152 nsize: {{1 2 5}}
153 args: -adapt -dm_plex_simplex 0 -dm_plex_box_faces 2,2
154 output_file: output/empty.out
155
156 TEST*/
157