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: {
14: options->adapt = PETSC_FALSE;
16: PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
17: PetscOptionsBool("-adapt", "Flag for adaptation of the surface mesh", "ex10.c", options->adapt, &options->adapt, NULL);
18: PetscOptionsEnd();
19: return 0;
20: }
22: static PetscErrorCode CreateDomainLabel(DM dm)
23: {
24: DMLabel label;
25: PetscInt cStart, cEnd, c;
28: DMCreateLabel(dm, "Cell Sets");
29: DMGetLabel(dm, "Cell Sets", &label);
30: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
31: for (c = cStart; c < cEnd; ++c) {
32: PetscReal centroid[3], volume, x, y;
34: DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL);
35: x = centroid[0]; 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)) {DMLabelSetValue(label, c, 1);continue;}
38: /* River channel is (0.1,0.45)--(1.0,0.55) */
39: if ((x >= 0.1 && x <= 1.0) && (y >= 0.45 && y <= 0.55)) {DMLabelSetValue(label, c, 2);continue;}
40: }
41: return 0;
42: }
44: static PetscErrorCode AdaptMesh(DM *dm, AppCtx *ctx)
45: {
46: DM dmCur = *dm;
47: DMLabel label;
48: IS valueIS, vIS;
49: PetscBool hasLabel;
50: const PetscInt *values;
51: PetscReal *volConst; /* Volume constraints for each label value */
52: PetscReal ratio;
53: PetscInt dim, Nv, v, cStart, cEnd, c;
54: PetscBool adapt = PETSC_TRUE;
57: if (!ctx->adapt) return 0;
58: DMHasLabel(*dm, "Cell Sets", &hasLabel);
59: if (!hasLabel) CreateDomainLabel(*dm);
60: DMGetDimension(*dm, &dim);
61: ratio = PetscPowRealInt(0.5, dim);
62: /* Get volume constraints */
63: DMGetLabel(*dm, "Cell Sets", &label);
64: DMLabelGetValueIS(label, &vIS);
65: ISDuplicate(vIS, &valueIS);
66: ISDestroy(&vIS);
67: /* Sorting ruins the label */
68: ISSort(valueIS);
69: ISGetLocalSize(valueIS, &Nv);
70: ISGetIndices(valueIS, &values);
71: PetscMalloc1(Nv, &volConst);
72: for (v = 0; v < Nv; ++v) {
73: char opt[128];
75: volConst[v] = PETSC_MAX_REAL;
76: PetscSNPrintf(opt, 128, "-volume_constraint_%d", (int) values[v]);
77: PetscOptionsGetReal(NULL, NULL, opt, &volConst[v], NULL);
78: }
79: ISRestoreIndices(valueIS, &values);
80: ISDestroy(&valueIS);
81: /* Adapt mesh iteratively */
82: while (adapt) {
83: DM dmAdapt;
84: DMLabel adaptLabel;
85: PetscInt nAdaptLoc[2], nAdapt[2];
87: adapt = PETSC_FALSE;
88: nAdaptLoc[0] = nAdaptLoc[1] = 0;
89: nAdapt[0] = nAdapt[1] = 0;
90: /* Adaptation is not preserving the domain label */
91: DMHasLabel(dmCur, "Cell Sets", &hasLabel);
92: if (!hasLabel) CreateDomainLabel(dmCur);
93: DMGetLabel(dmCur, "Cell Sets", &label);
94: DMLabelGetValueIS(label, &vIS);
95: ISDuplicate(vIS, &valueIS);
96: ISDestroy(&vIS);
97: /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */
98: ISSort(valueIS);
99: ISGetLocalSize(valueIS, &Nv);
100: ISGetIndices(valueIS, &values);
101: /* Construct adaptation label */
102: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
103: DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd);
104: for (c = cStart; c < cEnd; ++c) {
105: PetscReal volume, centroid[3];
106: PetscInt value, vidx;
108: DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL);
109: DMLabelGetValue(label, c, &value);
110: if (value < 0) continue;
111: PetscFindInt(value, Nv, values, &vidx);
113: if (volume > volConst[vidx]) {DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE); ++nAdaptLoc[0];}
114: if (volume < volConst[vidx]*ratio) {DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN); ++nAdaptLoc[1];}
115: }
116: ISRestoreIndices(valueIS, &values);
117: ISDestroy(&valueIS);
118: MPI_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dmCur));
119: if (nAdapt[0]) {
120: PetscInfo(dmCur, "Adapted mesh, marking %D cells for refinement, and %D cells for coarsening\n", nAdapt[0], nAdapt[1]);
121: DMAdaptLabel(dmCur, adaptLabel, &dmAdapt);
122: DMDestroy(&dmCur);
123: DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view");
124: dmCur = dmAdapt;
125: adapt = PETSC_TRUE;
126: }
127: DMLabelDestroy(&adaptLabel);
128: }
129: PetscFree(volConst);
130: *dm = dmCur;
131: return 0;
132: }
134: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
135: {
136: PetscInt dim;
139: /* Create top surface */
140: DMCreate(comm, dm);
141: DMSetType(*dm, DMPLEX);
142: PetscObjectSetOptionsPrefix((PetscObject) *dm, "init_");
143: DMSetFromOptions(*dm);
144: PetscObjectSetOptionsPrefix((PetscObject) *dm, NULL);
145: /* Adapt surface */
146: AdaptMesh(dm, user);
147: /* Extrude surface to get volume mesh */
148: DMGetDimension(*dm, &dim);
149: DMLocalizeCoordinates(*dm);
150: PetscObjectSetName((PetscObject) *dm, "Mesh");
151: DMSetFromOptions(*dm);
152: DMViewFromOptions(*dm, NULL, "-dm_view");
153: return 0;
154: }
156: int main(int argc, char **argv)
157: {
158: DM dm;
159: AppCtx user;
161: PetscInitialize(&argc, &argv, NULL, help);
162: ProcessOptions(PETSC_COMM_WORLD, &user);
163: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
164: DMDestroy(&dm);
165: PetscFinalize();
166: return 0;
167: }
169: /*TEST
171: test:
172: suffix: 0
173: requires: triangle
174: args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 1,1 -dm_extrude 1 -dm_view
176: test: # Regularly refine the surface before extrusion
177: suffix: 1
178: requires: triangle
179: args: -init_dm_plex_dim 2 -init_dm_refine 2 -dm_extrude 1 -dm_view
181: test: # Parallel run
182: suffix: 2
183: requires: triangle
184: nsize: 5
185: args: -init_dm_plex_dim 2 -init_dm_refine 3 -petscpartitioner_type simple -dm_extrude 3 -dm_view
187: test: # adaptively refine the surface before extrusion
188: suffix: 3
189: requires: triangle
190: 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
192: TEST*/