Actual source code: ex3.c
1: static char help[] = "Tests adaptive refinement using DMForest, and uses HDF5.\n\n";
3: #include <petscdmforest.h>
4: #include <petscdmplex.h>
5: #include <petscviewerhdf5.h>
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;
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));
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"));
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"));
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;
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;
55: for (i = 0; i < 3; ++i) {
56: DM postforest;
57: DMLabel adaptLabel = NULL;
59: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
61: PetscCall(DMForestGetCellChart(forest, &cStart, &cEnd));
62: for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));
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"));
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;
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: }
104: /* Create the global vector*/
105: PetscCall(DMCreateGlobalVector(forest, &g));
106: PetscCall(PetscObjectSetName((PetscObject)g, "g"));
107: PetscCall(VecSet(g, 1.0));
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));
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));
123: /* Load another vector to load into*/
124: PetscCall(DMCreateGlobalVector(forest, &g2));
125: PetscCall(PetscObjectSetName((PetscObject)g2, "g"));
126: PetscCall(VecZeroEntries(g2));
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));
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);
138: PetscCall(VecDestroy(&g));
139: PetscCall(VecDestroy(&g2));
140: PetscCall(DMDestroy(&forest));
141: PetscCall(PetscFinalize());
142: return 0;
143: }
145: /*TEST
147: build:
148: requires: hdf5 p4est
150: test:
151: suffix: 0
152: nsize: {{1 2 5}}
153: args: -adapt -dm_plex_simplex 0 -dm_plex_box_faces 2,2
155: TEST*/