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: {
  8:   DMCreate(comm, dm);
  9:   DMSetType(*dm, DMPLEX);
 10:   DMSetFromOptions(*dm);
 11:   DMViewFromOptions(*dm, NULL, "-dm_view");
 12:   return 0;
 13: }

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

 21:   DMCreateLabel(dm, "cells");
 22:   DMGetLabel(dm, "cells", label);
 23:   DMLabelClearStratum(*label, 1);
 24:   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) DMLabelSetValue(*label, c, 1);
 33:   DMPlexLabelComplete(dm, *label);
 34:   return 0;
 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;

 44:   DMGetCoordinatesLocalSetUp(dm);
 45:   DMCreateLabel(dm, "cells");
 46:   DMGetLabel(dm, "cells", label);
 47:   DMLabelClearStratum(*label, 1);
 48:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 49:   DMGetCoordinateDim(dm, &cdim);
 50:   for (PetscInt c = cStart; c < cEnd; ++c) {
 51:     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) DMLabelSetValue(*label, c, 1);
 55:     } else {
 56:       if (centroid[0] > 0.5) DMLabelSetValue(*label, c, 1);
 57:     }
 58:   }
 59:   DMPlexLabelComplete(dm, *label);
 60:   return 0;
 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;

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

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

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

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

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

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

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

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

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

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

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

162: /*TEST

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

170:   # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes
171:   testset:
172:     nsize: 3
173:     requires: parmetis
174:     args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \
175:           -sub_dm_plex_check_all -dm_view -sub_dm_view

177:     test:
178:       suffix: 1

180:     test:
181:       suffix: 2
182:       args: -domain

184:   # This set tests that global numberings can be made when some strata are missing on a process
185:   testset:
186:     nsize: 3
187:     requires: hdf5
188:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \
189:           -sub_dm_plex_check_all -sub_dm_view hdf5:subdm.h5

191:     test:
192:       suffix: 3
193:       args: -domain -height 0.625

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

200:     test:
201:       suffix: surface_2d
202:       args: -dm_plex_box_faces 5,5

204:     test:
205:       suffix: surface_3d
206:       args: -dm_plex_box_faces 5,5,5

208:     # This test checks that if cells are present in the SF, the dm is marked as having an overlap
209:     test:
210:       suffix: surface_2d_2
211:       nsize: 3
212:       args: -dm_plex_box_faces 5,5 -domain -height 0.625

214: TEST*/