Actual source code: ex16.c

  1: static char help[] = "Tests for creation of submeshes\n\n";

  3: #include <petscdmplex.h>

  5: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
  6: {
  7:   PetscFunctionBeginUser;
  8:   PetscCall(DMCreate(comm, dm));
  9:   PetscCall(DMSetType(*dm, DMPLEX));
 10:   PetscCall(DMSetFromOptions(*dm));
 11:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

 15: // Label half of the cells
 16: static PetscErrorCode CreateHalfCellsLabel(DM dm, PetscBool lower, DMLabel *label)
 17: {
 18:   PetscInt cStart, cEnd, cStartSub, cEndSub;

 20:   PetscFunctionBeginUser;
 21:   PetscCall(DMCreateLabel(dm, "cells"));
 22:   PetscCall(DMGetLabel(dm, "cells", label));
 23:   PetscCall(DMLabelClearStratum(*label, 1));
 24:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 25:   if (lower) {
 26:     cStartSub = cStart;
 27:     cEndSub   = cEnd / 2;
 28:   } else {
 29:     cStartSub = cEnd / 2;
 30:     cEndSub   = cEnd;
 31:   }
 32:   for (PetscInt c = cStartSub; c < cEndSub; ++c) PetscCall(DMLabelSetValue(*label, c, 1));
 33:   PetscCall(DMPlexLabelComplete(dm, *label));
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: // Label everything on the right half of the domain
 38: static PetscErrorCode CreateHalfDomainLabel(DM dm, PetscBool lower, PetscReal height, DMLabel *label)
 39: {
 40:   PetscReal centroid[3];
 41:   PetscInt  cStart, cEnd, cdim;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
 45:   PetscCall(DMCreateLabel(dm, "cells"));
 46:   PetscCall(DMGetLabel(dm, "cells", label));
 47:   PetscCall(DMLabelClearStratum(*label, 1));
 48:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 49:   PetscCall(DMGetCoordinateDim(dm, &cdim));
 50:   for (PetscInt c = cStart; c < cEnd; ++c) {
 51:     PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
 52:     if (height > 0.0 && PetscAbsReal(centroid[1] - height) > PETSC_SMALL) continue;
 53:     if (lower) {
 54:       if (centroid[0] < 0.5) PetscCall(DMLabelSetValue(*label, c, 1));
 55:     } else {
 56:       if (centroid[0] > 0.5) PetscCall(DMLabelSetValue(*label, c, 1));
 57:     }
 58:   }
 59:   PetscCall(DMPlexLabelComplete(dm, *label));
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: // Create a line of faces at a given x value
 64: static PetscErrorCode CreateLineLabel(DM dm, PetscReal x, DMLabel *label)
 65: {
 66:   PetscReal centroid[3];
 67:   PetscInt  fStart, fEnd;

 69:   PetscFunctionBeginUser;
 70:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
 71:   PetscCall(DMCreateLabel(dm, "faces"));
 72:   PetscCall(DMGetLabel(dm, "faces", label));
 73:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
 74:   for (PetscInt f = fStart; f < fEnd; ++f) {
 75:     PetscCall(DMPlexComputeCellGeometryFVM(dm, f, NULL, centroid, NULL));
 76:     if (PetscAbsReal(centroid[0] - x) < PETSC_SMALL) PetscCall(DMLabelSetValue(*label, f, 1));
 77:   }
 78:   PetscCall(DMPlexLabelComplete(dm, *label));
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: static PetscErrorCode CreateVolumeSubmesh(DM dm, PetscBool domain, PetscBool lower, PetscReal height, DM *subdm)
 83: {
 84:   DMLabel label, map;

 86:   PetscFunctionBegin;
 87:   if (domain) PetscCall(CreateHalfDomainLabel(dm, lower, height, &label));
 88:   else PetscCall(CreateHalfCellsLabel(dm, lower, &label));
 89:   PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm));
 90:   PetscCall(PetscObjectSetName((PetscObject)*subdm, "Submesh"));
 91:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subdm, "sub_"));
 92:   PetscCall(DMViewFromOptions(*subdm, NULL, "-dm_view"));
 93:   PetscCall(DMPlexGetSubpointMap(*subdm, &map));
 94:   PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view"));
 95:   PetscFunctionReturn(PETSC_SUCCESS);
 96: }

 98: static PetscErrorCode TestBoundaryField(DM dm)
 99: {
100:   DM       subdm;
101:   DMLabel  label, map;
102:   PetscFV  fvm;
103:   Vec      gv;
104:   PetscInt cdim;

106:   PetscFunctionBeginUser;
107:   PetscCall(CreateLineLabel(dm, 0.5, &label));
108:   PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdm));
109:   PetscCall(PetscObjectSetName((PetscObject)subdm, "Submesh"));
110:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subdm, "sub_"));
111:   PetscCall(DMViewFromOptions(subdm, NULL, "-dm_view"));
112:   PetscCall(DMPlexGetSubpointMap(subdm, &map));
113:   PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view"));

115:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)subdm), &fvm));
116:   PetscCall(PetscObjectSetName((PetscObject)fvm, "testField"));
117:   PetscCall(PetscFVSetNumComponents(fvm, 1));
118:   PetscCall(DMGetCoordinateDim(subdm, &cdim));
119:   PetscCall(PetscFVSetSpatialDimension(fvm, cdim));
120:   PetscCall(PetscFVSetFromOptions(fvm));

122:   PetscCall(DMAddField(subdm, NULL, (PetscObject)fvm));
123:   PetscCall(PetscFVDestroy(&fvm));
124:   PetscCall(DMCreateDS(subdm));

126:   PetscCall(DMCreateGlobalVector(subdm, &gv));
127:   PetscCall(PetscObjectSetName((PetscObject)gv, "potential"));
128:   PetscCall(VecSet(gv, 1.));
129:   PetscCall(VecViewFromOptions(gv, NULL, "-vec_view"));
130:   PetscCall(VecDestroy(&gv));
131:   PetscCall(DMDestroy(&subdm));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: int main(int argc, char **argv)
136: {
137:   DM        dm, subdm;
138:   PetscReal height = -1.0;
139:   PetscBool volume = PETSC_TRUE, domain = PETSC_FALSE;

141:   PetscFunctionBeginUser;
142:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
143:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &height, NULL));
144:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-volume", &volume, NULL));
145:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL));

147:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
148:   if (volume) {
149:     PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_TRUE, height, &subdm));
150:     PetscCall(DMSetFromOptions(subdm));
151:     PetscCall(DMDestroy(&subdm));
152:     PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_FALSE, height, &subdm));
153:     PetscCall(DMSetFromOptions(subdm));
154:     PetscCall(DMDestroy(&subdm));
155:   } else {
156:     PetscCall(TestBoundaryField(dm));
157:   }
158:   PetscCall(DMDestroy(&dm));
159:   PetscCall(PetscFinalize());
160:   return 0;
161: }

163: /*TEST

165:   test:
166:     suffix: 0
167:     requires: triangle
168:     args: -dm_coord_space 0 -sub_dm_plex_check_all \
169:           -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view

171:   test:
172:     suffix: 0_vtk
173:     requires: triangle
174:     args: -dm_coord_space 0 -sub_dm_plex_check_all \
175:           -dm_view vtk: -sub_dm_view vtk: -map_view

177:   # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes
178:   testset:
179:     nsize: 3
180:     requires: parmetis
181:     args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \
182:           -sub_dm_plex_check_all -dm_view -sub_dm_view

184:     test:
185:       suffix: 1

187:     test:
188:       suffix: 2
189:       args: -domain

191:   # This set tests that global numberings can be made when some strata are missing on a process
192:   testset:
193:     nsize: 3
194:     requires: hdf5
195:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
196:           -sub_dm_plex_check_all -sub_dm_view hdf5:subdm.h5

198:     test:
199:       suffix: 3
200:       args: -domain -height 0.625

202:   # This set tests that global numberings can be made when some strata are missing on a process
203:   testset:
204:     nsize: 3
205:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
206:           -sub_dm_plex_check_all -sub_dm_view {{vtk:subdm.vtk: vtk:subdm.vtu :subdm.txt :subdm_d.txt:ascii_info_detail}}

208:     test:
209:       suffix: 3_vtk
210:       args: -domain -height 0.625

212:   # This test checks whether filter can extract a lower-dimensional manifold and output a field on it
213:   testset:
214:     requires: hdf5
215:     args: -volume 0 -dm_plex_simplex 0 -sub_dm_view hdf5:subdm.h5 -vec_view hdf5:subdm.h5::append

217:     test:
218:       suffix: surface_2d
219:       args: -dm_plex_box_faces 5,5

221:     test:
222:       suffix: surface_3d
223:       args: -dm_plex_box_faces 5,5,5

225:     # This test checks that if cells are present in the SF, the dm is marked as having an overlap
226:     test:
227:       suffix: surface_2d_2
228:       nsize: 3
229:       args: -dm_plex_box_faces 5,5 -domain -height 0.625

231: TEST*/