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;
18: PetscInitialize(&argc, &argv, NULL, help);
19: PetscOptionsGetBool(NULL, NULL, "-adapt", &adapt, NULL);
20: PetscOptionsGetBool(NULL, NULL, "-user_section", &userSection, NULL);
22: /* Create a base DMPlex mesh */
23: DMCreate(PETSC_COMM_WORLD, &base);
24: DMSetType(base, DMPLEX);
25: DMSetFromOptions(base);
26: DMViewFromOptions(base, NULL, "-dm_view");
28: /* Convert Plex mesh to Forest and destroy base */
29: DMCreate(PETSC_COMM_WORLD, &forest);
30: DMSetType(forest, DMP4EST);
31: DMForestSetBaseDM(forest, base);
32: DMSetUp(forest);
33: DMDestroy(&base);
34: 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: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
43: DMLabelSetValue(adaptLabel, 0, DM_ADAPT_REFINE);
44: DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest);
45: DMForestSetAdaptivityLabel(postforest, adaptLabel);
46: DMLabelDestroy(&adaptLabel);
47: DMSetUp(postforest);
48: 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: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
61: DMForestGetCellChart(forest, &cStart, &cEnd);
62: for (c = cStart; c < cEnd; ++c) DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);
64: DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest);
65: DMForestSetAdaptivityLabel(postforest, adaptLabel);
66: DMLabelDestroy(&adaptLabel);
67: DMSetUp(postforest);
68: DMDestroy(&forest);
69: forest = postforest;
70: }
71: }
72: DMViewFromOptions(forest, NULL, "-dm_view");
74: /* Setup the section*/
75: if (userSection) {
76: DMConvert(forest, DMPLEX, &plex);
77: DMViewFromOptions(plex, NULL, "-dm_view");
78: DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);
79: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Vertices [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", vStart, vEnd);
80: PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
81: PetscSectionCreate(PetscObjectComm((PetscObject)forest), &s);
82: PetscSectionSetNumFields(s, 1);
83: PetscSectionSetChart(s, vStart, vEnd);
84: for (v = vStart; v < vEnd; ++v) {
85: PetscSectionSetDof(s, v, 1);
86: PetscSectionSetFieldDof(s, v, 0, 1);
87: }
88: PetscSectionSetUp(s);
89: DMSetLocalSection(forest, s);
90: PetscObjectViewFromOptions((PetscObject)s, NULL, "-my_section_view");
91: PetscSectionDestroy(&s);
92: DMDestroy(&plex);
93: } else {
94: PetscFE fe;
95: PetscInt dim;
97: DMGetDimension(forest, &dim);
98: PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 1, PETSC_DETERMINE, &fe);
99: DMAddField(forest, NULL, (PetscObject)fe);
100: PetscFEDestroy(&fe);
101: DMCreateDS(forest);
102: }
104: /* Create the global vector*/
105: DMCreateGlobalVector(forest, &g);
106: PetscObjectSetName((PetscObject)g, "g");
107: VecSet(g, 1.0);
109: /* Test global to local*/
110: Vec l;
111: DMCreateLocalVector(forest, &l);
112: VecZeroEntries(l);
113: DMGlobalToLocal(forest, g, INSERT_VALUES, l);
114: VecZeroEntries(g);
115: DMLocalToGlobal(forest, l, INSERT_VALUES, g);
116: VecDestroy(&l);
118: /* Save a vector*/
119: PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_WRITE, &viewer);
120: VecView(g, viewer);
121: PetscViewerDestroy(&viewer);
123: /* Load another vector to load into*/
124: DMCreateGlobalVector(forest, &g2);
125: PetscObjectSetName((PetscObject)g2, "g");
126: VecZeroEntries(g2);
128: /* Load a vector*/
129: PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_READ, &viewer);
130: VecLoad(g2, viewer);
131: PetscViewerDestroy(&viewer);
133: /* Check if the data is the same*/
134: VecAXPY(g2, -1.0, g);
135: VecNorm(g2, NORM_INFINITY, &nrm);
138: VecDestroy(&g);
139: VecDestroy(&g2);
140: DMDestroy(&forest);
141: 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*/