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