Actual source code: ex57.c

  1: static char help[] = "Tests for ephemeral meshes.\n";

  3: #include <petscdmplex.h>
  4: #include <petscdmplextransform.h>
  5: #include <petscdmlabelephemeral.h>

  7: /*
  8:   Use

 10:     -ref_dm_view -eph_dm_view

 12:   to view the concrete and ephemeral meshes from the first transformation, and

 14:    -ref_dm_sec_view -eph_dm_sec_view

 16:   for the second.
 17: */

 19: // Should remove when I have an API for everything
 20: #include <petsc/private/dmplextransformimpl.h>

 22: typedef struct {
 23:   DMLabel   active;   /* Label for transform */
 24:   PetscBool second;   /* Flag to execute a second transformation */
 25:   PetscBool concrete; /* Flag to use the concrete mesh for the second transformation */
 26: } AppCtx;

 28: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 29: {
 30:   PetscInt  cells[1024], Nc = 1024;
 31:   PetscBool flg;

 33:   PetscFunctionBeginUser;
 34:   options->active   = NULL;
 35:   options->second   = PETSC_FALSE;
 36:   options->concrete = PETSC_FALSE;

 38:   PetscOptionsBegin(comm, "", "Ephemeral Meshing Options", "DMPLEX");
 39:   PetscCall(PetscOptionsIntArray("-cells", "Cells to mark for transformation", "ex57.c", cells, &Nc, &flg));
 40:   if (flg) {
 41:     PetscCall(DMLabelCreate(comm, "active", &options->active));
 42:     for (PetscInt c = 0; c < Nc; ++c) PetscCall(DMLabelSetValue(options->active, cells[c], DM_ADAPT_REFINE));
 43:   }
 44:   PetscCall(PetscOptionsBool("-second", "Use a second transformation", "ex57.c", options->second, &options->second, &flg));
 45:   PetscCall(PetscOptionsBool("-concrete", "Use concrete mesh for the second transformation", "ex57.c", options->concrete, &options->concrete, &flg));
 46:   PetscOptionsEnd();
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
 51: {
 52:   PetscFunctionBeginUser;
 53:   PetscCall(DMCreate(comm, dm));
 54:   PetscCall(DMSetType(*dm, DMPLEX));
 55:   PetscCall(DMSetFromOptions(*dm));
 56:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
 57:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 61: static PetscErrorCode CreateTransform(DM dm, DMLabel active, const char prefix[], DMPlexTransform *tr)
 62: {
 63:   PetscFunctionBeginUser;
 64:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), tr));
 65:   PetscCall(PetscObjectSetName((PetscObject)*tr, "Transform"));
 66:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*tr, prefix));
 67:   PetscCall(DMPlexTransformSetDM(*tr, dm));
 68:   PetscCall(DMPlexTransformSetActive(*tr, active));
 69:   PetscCall(DMPlexTransformSetFromOptions(*tr));
 70:   PetscCall(DMPlexTransformSetUp(*tr));

 72:   PetscCall(DMSetApplicationContext(dm, *tr));
 73:   PetscCall(PetscObjectViewFromOptions((PetscObject)*tr, NULL, "-dm_plex_transform_view"));
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode CreateEphemeralMesh(DMPlexTransform tr, DM *tdm)
 78: {
 79:   PetscFunctionBegin;
 80:   PetscCall(DMPlexCreateEphemeral(tr, tdm));
 81:   PetscCall(PetscObjectSetName((PetscObject)*tdm, "Ephemeral Mesh"));
 82:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*tdm, "eph_"));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: static PetscErrorCode CreateConcreteMesh(DMPlexTransform tr, DM *rdm)
 87: {
 88:   DM cdm, codm, rcodm;

 90:   PetscFunctionBegin;
 91:   PetscCall(DMPlexTransformGetDM(tr, &cdm));
 92:   PetscCall(DMPlexTransformApply(tr, cdm, rdm));
 93:   PetscCall(DMSetCoarsenLevel(*rdm, cdm->leveldown));
 94:   PetscCall(DMSetRefineLevel(*rdm, cdm->levelup + 1));
 95:   PetscCall(DMCopyDisc(cdm, *rdm));
 96:   PetscCall(DMGetCoordinateDM(cdm, &codm));
 97:   PetscCall(DMGetCoordinateDM(*rdm, &rcodm));
 98:   PetscCall(DMCopyDisc(codm, rcodm));
 99:   PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
100:   PetscCall(DMSetCoarseDM(*rdm, cdm));
101:   PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
102:   if (rdm) {
103:     ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)cdm->data)->printFEM;
104:     ((DM_Plex *)(*rdm)->data)->printL2  = ((DM_Plex *)cdm->data)->printL2;
105:   }
106:   PetscCall(PetscObjectSetName((PetscObject)*rdm, "Concrete Mesh"));
107:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*rdm, "ref_"));
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode CompareMeshes(DM dmA, DM dmB, DM dm)
112: {
113:   PetscInt dim, dimB, pStart, pEnd, pStartB, pEndB;
114:   MPI_Comm comm;

116:   PetscFunctionBegin;
117:   PetscCall(PetscObjectGetComm((PetscObject)dmA, &comm));
118:   PetscCall(DMGetDimension(dmA, &dim));
119:   PetscCall(DMGetDimension(dmB, &dimB));
120:   PetscCheck(dim == dimB, comm, PETSC_ERR_ARG_INCOMP, "Dimension from dmA %" PetscInt_FMT " != %" PetscInt_FMT " from dmB", dim, dimB);
121:   PetscCall(DMPlexGetChart(dmA, &pStart, &pEnd));
122:   PetscCall(DMPlexGetChart(dmB, &pStartB, &pEndB));
123:   PetscCheck(pStart == pStartB && pEnd == pEndB, comm, PETSC_ERR_ARG_INCOMP, "Chart from dmA (%" PetscInt_FMT ", %" PetscInt_FMT ") does not match (%" PetscInt_FMT ", %" PetscInt_FMT ") for dmB", pStart, pEnd, pStartB, pEndB);
124:   for (PetscInt p = pStart; p < pEnd; ++p) {
125:     const PetscInt *cone, *ornt, *coneB, *orntB;
126:     PetscInt        coneSize, coneSizeB;

128:     PetscCall(DMPlexGetConeSize(dmA, p, &coneSize));
129:     PetscCall(DMPlexGetConeSize(dmB, p, &coneSizeB));
130:     PetscCheck(coneSize == coneSizeB, comm, PETSC_ERR_ARG_INCOMP, "Cone size for %" PetscInt_FMT " from dmA %" PetscInt_FMT " does not match %" PetscInt_FMT " for dmB", p, coneSize, coneSizeB);
131:     PetscCall(DMPlexGetOrientedCone(dmA, p, &cone, &ornt));
132:     PetscCall(DMPlexGetOrientedCone(dmB, p, &coneB, &orntB));
133:     for (PetscInt c = 0; c < coneSize; ++c) {
134:       PetscCheck(cone[c] == coneB[c] && ornt[c] == orntB[c], comm, PETSC_ERR_ARG_INCOMP, "Cone point %" PetscInt_FMT " for point %" PetscInt_FMT " from dmA (%" PetscInt_FMT ", %" PetscInt_FMT ") does not match (%" PetscInt_FMT ", %" PetscInt_FMT ") for dmB", c, p, cone[c], ornt[c], coneB[c], orntB[c]);
135:     }
136:     PetscCall(DMPlexRestoreOrientedCone(dmA, p, &cone, &ornt));
137:     PetscCall(DMPlexRestoreOrientedCone(dmB, p, &coneB, &orntB));
138:   }
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: int main(int argc, char *argv[])
143: {
144:   DM              dm, tdm, rdm;
145:   DMPlexTransform tr;
146:   AppCtx          user;
147:   MPI_Comm        comm;

149:   PetscFunctionBeginUser;
150:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
151:   comm = PETSC_COMM_WORLD;
152:   PetscCall(ProcessOptions(comm, &user));
153:   PetscCall(CreateMesh(comm, &dm));
154:   PetscCall(CreateTransform(dm, user.active, "first_", &tr));
155:   PetscCall(CreateEphemeralMesh(tr, &tdm));
156:   PetscCall(CreateConcreteMesh(tr, &rdm));
157:   if (user.second) {
158:     DMPlexTransform tr2;
159:     DM              tdm2, rdm2;

161:     PetscCall(DMViewFromOptions(rdm, NULL, "-dm_sec_view"));
162:     PetscCall(DMViewFromOptions(tdm, NULL, "-dm_sec_view"));
163:     if (user.concrete) PetscCall(CreateTransform(rdm, user.active, "second_", &tr2));
164:     else PetscCall(CreateTransform(tdm, user.active, "second_", &tr2));
165:     PetscCall(CreateEphemeralMesh(tr2, &tdm2));
166:     PetscCall(CreateConcreteMesh(tr2, &rdm2));
167:     PetscCall(DMDestroy(&tdm));
168:     PetscCall(DMDestroy(&rdm));
169:     PetscCall(DMPlexTransformDestroy(&tr2));
170:     tdm = tdm2;
171:     rdm = rdm2;
172:   }
173:   PetscCall(DMViewFromOptions(tdm, NULL, "-dm_view"));
174:   PetscCall(DMViewFromOptions(rdm, NULL, "-dm_view"));
175:   PetscCall(CompareMeshes(rdm, tdm, dm));
176:   PetscCall(DMPlexTransformDestroy(&tr));
177:   PetscCall(DMDestroy(&tdm));
178:   PetscCall(DMDestroy(&rdm));
179:   PetscCall(DMDestroy(&dm));
180:   PetscCall(DMLabelDestroy(&user.active));
181:   PetscCall(PetscFinalize());
182:   return 0;
183: }

185: /*TEST

187:   # Tests for regular refinement of whole meshes
188:   test:
189:     suffix: tri
190:     requires: triangle
191:     args: -first_dm_plex_transform_view ::ascii_info_detail

193:   test:
194:     suffix: quad
195:     args: -dm_plex_simplex 0

197:   # Here I am checking that the 'marker' label is correct for the ephemeral mesh
198:   test:
199:     suffix: tet
200:     requires: ctetgen
201:     args: -dm_plex_dim 3 -ref_dm_view -eph_dm_view -eph_dm_plex_view_labels marker

203:   test:
204:     suffix: hex
205:     args: -dm_plex_dim 3 -dm_plex_simplex 0

207:   # Tests for filter patches
208:   testset:
209:     args: -first_dm_plex_transform_type transform_filter -ref_dm_view

211:     test:
212:       suffix: tri_patch
213:       requires: triangle
214:       args: -cells 0,1,2,4

216:     test:
217:       suffix: quad_patch
218:       args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -cells 0,1,3,4

220:   # Tests for refined filter patches
221:   testset:
222:     args: -first_dm_plex_transform_type transform_filter -ref_dm_view -eph_dm_view -second

224:     test:
225:       suffix: tri_patch_ref
226:       requires: triangle
227:       args: -cells 0,1,2,4

229:     test:
230:       suffix: tri_patch_ref_concrete
231:       requires: triangle
232:       args: -cells 0,1,2,4 -concrete -first_dm_plex_transform_view ::ascii_info_detail

234:   # Tests for boundary layer refinement
235:   test:
236:     suffix: quad_bl
237:     args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_extrude 2 -cells 0,2,4,6,8 \
238:           -first_dm_plex_transform_type refine_boundary_layer -first_dm_plex_transform_bl_splits 4 \
239:           -ref_dm_view

241: TEST*/