Actual source code: ex10.c

  1: static char help[] = "TDycore Mesh Examples\n\n";

  3: #include <petscdmplex.h>

  5: typedef struct {
  6:   PetscBool adapt;   /* Flag for adaptation of the surface mesh */
  7:   PetscBool extrude; /* Flag for extrusion of the suraace mesh */
  8: } AppCtx;

 10: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 11: {

 15:   options->adapt   = PETSC_FALSE;
 16:   options->extrude = PETSC_TRUE;

 18:   PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
 19:   PetscOptionsBool("-adapt", "Flag for adaptation of the surface mesh", "ex10.c", options->adapt, &options->adapt, NULL);
 20:   PetscOptionsBool("-extrude", "Flag for extrusion of the surface mesh", "ex10.c", options->extrude, &options->extrude, NULL);
 21:   PetscOptionsEnd();
 22:   return(0);
 23: }

 25: static PetscErrorCode CreateDomainLabel(DM dm)
 26: {
 27:   DMLabel        label;
 28:   PetscInt       cStart, cEnd, c;

 32:   DMCreateLabel(dm, "Cell Sets");
 33:   DMGetLabel(dm, "Cell Sets", &label);
 34:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 35:   for (c = cStart; c < cEnd; ++c) {
 36:     PetscReal centroid[3], volume, x, y;

 38:     DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL);
 39:     x = centroid[0]; y = centroid[1];
 40:     /* Headwaters are (0.0,0.25)--(0.1,0.75) */
 41:     if ((x >= 0.0 && x <  0.1) && (y >= 0.25 && y <= 0.75)) {DMLabelSetValue(label, c, 1);continue;}
 42:     /* River channel is (0.1,0.45)--(1.0,0.55) */
 43:     if ((x >= 0.1 && x <= 1.0) && (y >= 0.45 && y <= 0.55)) {DMLabelSetValue(label, c, 2);continue;}
 44:   }
 45:   return(0);
 46: }

 48: static PetscErrorCode AdaptMesh(DM *dm, AppCtx *ctx)
 49: {
 50:   DM              dmCur = *dm;
 51:   DMLabel         label;
 52:   IS              valueIS, vIS;
 53:   PetscBool       hasLabel;
 54:   const PetscInt *values;
 55:   PetscReal      *volConst; /* Volume constraints for each label value */
 56:   PetscReal       ratio;
 57:   PetscInt        dim, Nv, v, cStart, cEnd, c;
 58:   PetscBool       adapt = PETSC_TRUE;
 59:   PetscErrorCode  ierr;

 62:   if (!ctx->adapt) return(0);
 63:   DMHasLabel(*dm, "Cell Sets", &hasLabel);
 64:   if (!hasLabel) {CreateDomainLabel(*dm);}
 65:   DMGetDimension(*dm, &dim);
 66:   ratio = PetscPowRealInt(0.5, dim);
 67:   /* Get volume constraints */
 68:   DMGetLabel(*dm, "Cell Sets", &label);
 69:   DMLabelGetValueIS(label, &vIS);
 70:   ISDuplicate(vIS, &valueIS);
 71:   ISDestroy(&vIS);
 72:   /* Sorting ruins the label */
 73:   ISSort(valueIS);
 74:   ISGetLocalSize(valueIS, &Nv);
 75:   ISGetIndices(valueIS, &values);
 76:   PetscMalloc1(Nv, &volConst);
 77:   for (v = 0; v < Nv; ++v) {
 78:     char opt[128];

 80:     volConst[v] = PETSC_MAX_REAL;
 81:     PetscSNPrintf(opt, 128, "-volume_constraint_%d", (int) values[v]);
 82:     PetscOptionsGetReal(NULL, NULL, opt, &volConst[v], NULL);
 83:   }
 84:   ISRestoreIndices(valueIS, &values);
 85:   ISDestroy(&valueIS);
 86:   /* Adapt mesh iteratively */
 87:   while (adapt) {
 88:     DM       dmAdapt;
 89:     DMLabel  adaptLabel;
 90:     PetscInt nAdaptLoc[2], nAdapt[2];

 92:     adapt = PETSC_FALSE;
 93:     nAdaptLoc[0] = nAdaptLoc[1] = 0;
 94:     nAdapt[0]    = nAdapt[1]    = 0;
 95:     /* Adaptation is not preserving the domain label */
 96:     DMHasLabel(dmCur, "Cell Sets", &hasLabel);
 97:     if (!hasLabel) {CreateDomainLabel(dmCur);}
 98:     DMGetLabel(dmCur, "Cell Sets", &label);
 99:     DMLabelGetValueIS(label, &vIS);
100:     ISDuplicate(vIS, &valueIS);
101:     ISDestroy(&vIS);
102:     /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */
103:     ISSort(valueIS);
104:     ISGetLocalSize(valueIS, &Nv);
105:     ISGetIndices(valueIS, &values);
106:     /* Construct adaptation label */
107:     DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
108:     DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd);
109:     for (c = cStart; c < cEnd; ++c) {
110:       PetscReal volume, centroid[3];
111:       PetscInt  value, vidx;

113:       DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL);
114:       DMLabelGetValue(label, c, &value);
115:       if (value < 0) continue;
116:       PetscFindInt(value, Nv, values, &vidx);
117:       if (vidx < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D for cell %D does not exist in label", value, c);
118:       if (volume > volConst[vidx])        {DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);  ++nAdaptLoc[0];}
119:       if (volume < volConst[vidx]*ratio) {DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN); ++nAdaptLoc[1];}
120:     }
121:     ISRestoreIndices(valueIS, &values);
122:     ISDestroy(&valueIS);
123:     MPI_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dmCur));
124:     if (nAdapt[0]) {
125:       PetscInfo2(dmCur, "Adapted mesh, marking %D cells for refinement, and %D cells for coarsening\n", nAdapt[0], nAdapt[1]);
126:       DMAdaptLabel(dmCur, adaptLabel, &dmAdapt);
127:       DMDestroy(&dmCur);
128:       DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view");
129:       dmCur = dmAdapt;
130:       adapt = PETSC_TRUE;
131:     }
132:     DMLabelDestroy(&adaptLabel);
133:   }
134:   PetscFree(volConst);
135:   *dm = dmCur;
136:   return(0);
137: }

139: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
140: {
141:   PetscInt       dim;

145:   /* Create top surface */
146:   DMCreate(comm, dm);
147:   DMSetType(*dm, DMPLEX);
148:   PetscObjectSetOptionsPrefix((PetscObject) *dm, "init_");
149:   DMSetFromOptions(*dm);
150:   /* Adapt surface */
151:   AdaptMesh(dm, user);
152:   /* Extrude surface to get volume mesh */
153:   DMGetDimension(*dm, &dim);
154:   if (dim < 3 && user->extrude) {
155:     DM edm;

157:     PetscObjectSetOptionsPrefix((PetscObject) *dm, "srf_");
158:     DMSetFromOptions(*dm);
159:     DMViewFromOptions(*dm, NULL, "-dm_view");
160:     PetscObjectSetOptionsPrefix((PetscObject) *dm, NULL);
161:     DMPlexExtrude(*dm, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_TRUE, NULL, PETSC_TRUE, &edm);
162:     DMDestroy(dm);
163:     *dm  = edm;
164:   }
165:   DMLocalizeCoordinates(*dm);
166:   PetscObjectSetName((PetscObject) *dm, "Mesh");
167:   DMSetFromOptions(*dm);
168:   DMViewFromOptions(*dm, NULL, "-dm_view");
169:   return(0);
170: }

172: int main(int argc, char **argv)
173: {
174:   DM             dm;
175:   AppCtx         user;

178:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
179:   ProcessOptions(PETSC_COMM_WORLD, &user);
180:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
181:   DMDestroy(&dm);
182:   PetscFinalize();
183:   return ierr;
184: }

186: /*TEST

188:   test:
189:     suffix: 0
190:     requires: triangle
191:     args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 1,1 -dm_view

193:   test: # Regularly refine the surface before extrusion
194:     suffix: 1
195:     requires: triangle
196:     args: -init_dm_plex_dim 2 -srf_dm_refine 2 -dm_view

198:   test: # Parallel run
199:     suffix: 2
200:     requires: triangle
201:     nsize: 5
202:     args: -init_dm_plex_dim 2 -srf_dm_refine 3 -petscpartitioner_type simple -dm_distribute -dm_plex_extrude_layers 3 -dm_view

204:   test: # adaptively refine the surface before extrusion
205:     suffix: 3
206:     requires: triangle
207:     args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 5,5 -adapt -volume_constraint_1 0.01 -volume_constraint_2 0.000625 -dm_plex_extrude_layers 10

209: TEST*/