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