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: } AppCtx;

  9: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 10: {
 11:   PetscFunctionBeginUser;
 12:   options->adapt = PETSC_FALSE;

 14:   PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
 15:   PetscCall(PetscOptionsBool("-adapt", "Flag for adaptation of the surface mesh", "ex10.c", options->adapt, &options->adapt, NULL));
 16:   PetscOptionsEnd();
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: static PetscErrorCode CreateDomainLabel(DM dm)
 21: {
 22:   DMLabel  label;
 23:   PetscInt cStart, cEnd, c;

 25:   PetscFunctionBeginUser;
 26:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
 27:   PetscCall(DMCreateLabel(dm, "Cell Sets"));
 28:   PetscCall(DMGetLabel(dm, "Cell Sets", &label));
 29:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 30:   for (c = cStart; c < cEnd; ++c) {
 31:     PetscReal centroid[3], volume, x, y;

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

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

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

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

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

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

146: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
147: {
148:   PetscInt dim;

150:   PetscFunctionBeginUser;
151:   /* Create top surface */
152:   PetscCall(DMCreate(comm, dm));
153:   PetscCall(DMSetType(*dm, DMPLEX));
154:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "init_"));
155:   PetscCall(DMSetFromOptions(*dm));
156:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
157:   /* Adapt surface */
158:   PetscCall(AdaptMesh(dm, user));
159:   /* Extrude surface to get volume mesh */
160:   PetscCall(DMGetDimension(*dm, &dim));
161:   PetscCall(DMLocalizeCoordinates(*dm));
162:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
163:   PetscCall(DMSetFromOptions(*dm));
164:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: int main(int argc, char **argv)
169: {
170:   DM     dm;
171:   AppCtx user;

173:   PetscFunctionBeginUser;
174:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
175:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
176:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
177:   PetscCall(DMDestroy(&dm));
178:   PetscCall(PetscFinalize());
179:   return 0;
180: }

182: /*TEST

184:   test:
185:     suffix: 0
186:     requires: triangle
187:     args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 1,1 -dm_extrude 1 -dm_view

189:   test: # Regularly refine the surface before extrusion
190:     suffix: 1
191:     requires: triangle
192:     args: -init_dm_plex_dim 2 -init_dm_refine 2 -dm_extrude 1 -dm_view

194:   test: # Parallel run
195:     suffix: 2
196:     requires: triangle
197:     nsize: 5
198:     args: -init_dm_plex_dim 2 -init_dm_refine 3 -petscpartitioner_type simple -dm_extrude 3 -dm_view

200:   test: # adaptively refine the surface before extrusion
201:     suffix: 3
202:     requires: triangle
203:     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_extrude 10

205: TEST*/