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: {
12: options->adapt = PETSC_FALSE;
14: PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
15: PetscOptionsBool("-adapt", "Flag for adaptation of the surface mesh", "ex10.c", options->adapt, &options->adapt, NULL);
16: PetscOptionsEnd();
17: return 0;
18: }
20: static PetscErrorCode CreateDomainLabel(DM dm)
21: {
22: DMLabel label;
23: PetscInt cStart, cEnd, c;
26: DMGetCoordinatesLocalSetUp(dm);
27: DMCreateLabel(dm, "Cell Sets");
28: DMGetLabel(dm, "Cell Sets", &label);
29: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
30: for (c = cStart; c < cEnd; ++c) {
31: PetscReal centroid[3], volume, x, y;
33: 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: 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: DMLabelSetValue(label, c, 2);
44: continue;
45: }
46: }
47: return 0;
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;
63: if (!ctx->adapt) return 0;
64: DMHasLabel(*dm, "Cell Sets", &hasLabel);
65: if (!hasLabel) CreateDomainLabel(*dm);
66: DMGetDimension(*dm, &dim);
67: ratio = PetscPowRealInt(0.5, dim);
68: /* Get volume constraints */
69: DMGetLabel(*dm, "Cell Sets", &label);
70: DMLabelGetValueIS(label, &vIS);
71: ISDuplicate(vIS, &valueIS);
72: ISDestroy(&vIS);
73: /* Sorting ruins the label */
74: ISSort(valueIS);
75: ISGetLocalSize(valueIS, &Nv);
76: ISGetIndices(valueIS, &values);
77: PetscMalloc1(Nv, &volConst);
78: for (v = 0; v < Nv; ++v) {
79: char opt[128];
81: volConst[v] = PETSC_MAX_REAL;
82: PetscSNPrintf(opt, 128, "-volume_constraint_%d", (int)values[v]);
83: PetscOptionsGetReal(NULL, NULL, opt, &volConst[v], NULL);
84: }
85: ISRestoreIndices(valueIS, &values);
86: 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: DMHasLabel(dmCur, "Cell Sets", &hasLabel);
98: if (!hasLabel) CreateDomainLabel(dmCur);
99: DMGetLabel(dmCur, "Cell Sets", &label);
100: DMLabelGetValueIS(label, &vIS);
101: ISDuplicate(vIS, &valueIS);
102: ISDestroy(&vIS);
103: /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */
104: ISSort(valueIS);
105: ISGetLocalSize(valueIS, &Nv);
106: ISGetIndices(valueIS, &values);
107: /* Construct adaptation label */
108: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
109: DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd);
110: for (c = cStart; c < cEnd; ++c) {
111: PetscReal volume, centroid[3];
112: PetscInt value, vidx;
114: DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL);
115: DMLabelGetValue(label, c, &value);
116: if (value < 0) continue;
117: PetscFindInt(value, Nv, values, &vidx);
119: if (volume > volConst[vidx]) {
120: DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);
121: ++nAdaptLoc[0];
122: }
123: if (volume < volConst[vidx] * ratio) {
124: DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN);
125: ++nAdaptLoc[1];
126: }
127: }
128: ISRestoreIndices(valueIS, &values);
129: ISDestroy(&valueIS);
130: MPI_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dmCur));
131: if (nAdapt[0]) {
132: PetscInfo(dmCur, "Adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nAdapt[0], nAdapt[1]);
133: DMAdaptLabel(dmCur, adaptLabel, &dmAdapt);
134: DMDestroy(&dmCur);
135: DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view");
136: dmCur = dmAdapt;
137: adapt = PETSC_TRUE;
138: }
139: DMLabelDestroy(&adaptLabel);
140: }
141: PetscFree(volConst);
142: *dm = dmCur;
143: return 0;
144: }
146: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
147: {
148: PetscInt dim;
151: /* Create top surface */
152: DMCreate(comm, dm);
153: DMSetType(*dm, DMPLEX);
154: PetscObjectSetOptionsPrefix((PetscObject)*dm, "init_");
155: DMSetFromOptions(*dm);
156: PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL);
157: /* Adapt surface */
158: AdaptMesh(dm, user);
159: /* Extrude surface to get volume mesh */
160: DMGetDimension(*dm, &dim);
161: DMLocalizeCoordinates(*dm);
162: PetscObjectSetName((PetscObject)*dm, "Mesh");
163: DMSetFromOptions(*dm);
164: DMViewFromOptions(*dm, NULL, "-dm_view");
165: return 0;
166: }
168: int main(int argc, char **argv)
169: {
170: DM dm;
171: AppCtx user;
174: PetscInitialize(&argc, &argv, NULL, help);
175: ProcessOptions(PETSC_COMM_WORLD, &user);
176: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
177: DMDestroy(&dm);
178: 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*/