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