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*/