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 nAdapt[2];
93: adapt = PETSC_FALSE;
94: nAdapt[0] = nAdapt[1] = 0;
95: /* Adaptation is not preserving the domain label */
96: PetscCall(DMHasLabel(dmCur, "Cell Sets", &hasLabel));
97: if (!hasLabel) PetscCall(CreateDomainLabel(dmCur));
98: PetscCall(DMGetLabel(dmCur, "Cell Sets", &label));
99: PetscCall(DMLabelGetValueIS(label, &vIS));
100: PetscCall(ISDuplicate(vIS, &valueIS));
101: PetscCall(ISDestroy(&vIS));
102: /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */
103: PetscCall(ISSort(valueIS));
104: PetscCall(ISGetLocalSize(valueIS, &Nv));
105: PetscCall(ISGetIndices(valueIS, &values));
106: /* Construct adaptation label */
107: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
108: PetscCall(DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd));
109: for (c = cStart; c < cEnd; ++c) {
110: PetscReal volume, centroid[3];
111: PetscInt value, vidx;
113: PetscCall(DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL));
114: PetscCall(DMLabelGetValue(label, c, &value));
115: if (value < 0) continue;
116: PetscCall(PetscFindInt(value, Nv, values, &vidx));
117: PetscCheck(vidx >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " for cell %" PetscInt_FMT " does not exist in label", value, c);
118: if (volume > volConst[vidx]) {
119: PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));
120: ++nAdapt[0];
121: }
122: if (volume < volConst[vidx] * ratio) {
123: PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN));
124: ++nAdapt[1];
125: }
126: }
127: PetscCall(ISRestoreIndices(valueIS, &values));
128: PetscCall(ISDestroy(&valueIS));
129: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dmCur)));
130: if (nAdapt[0]) {
131: PetscCall(PetscInfo(dmCur, "Adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nAdapt[0], nAdapt[1]));
132: PetscCall(DMAdaptLabel(dmCur, adaptLabel, &dmAdapt));
133: PetscCall(DMDestroy(&dmCur));
134: PetscCall(DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view"));
135: dmCur = dmAdapt;
136: adapt = PETSC_TRUE;
137: }
138: PetscCall(DMLabelDestroy(&adaptLabel));
139: }
140: PetscCall(PetscFree(volConst));
141: *dm = dmCur;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
146: {
147: PetscInt dim;
149: PetscFunctionBeginUser;
150: /* Create top surface */
151: PetscCall(DMCreate(comm, dm));
152: PetscCall(DMSetType(*dm, DMPLEX));
153: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "init_"));
154: PetscCall(DMSetFromOptions(*dm));
155: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
156: /* Adapt surface */
157: PetscCall(AdaptMesh(dm, user));
158: /* Extrude surface to get volume mesh */
159: PetscCall(DMGetDimension(*dm, &dim));
160: PetscCall(DMLocalizeCoordinates(*dm));
161: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
162: PetscCall(DMSetFromOptions(*dm));
163: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: int main(int argc, char **argv)
168: {
169: DM dm;
170: AppCtx user;
172: PetscFunctionBeginUser;
173: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
174: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
175: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
176: PetscCall(DMDestroy(&dm));
177: PetscCall(PetscFinalize());
178: return 0;
179: }
181: /*TEST
183: test:
184: suffix: 0
185: requires: triangle
186: args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 1,1 -dm_extrude 1 -dm_view
188: test: # Regularly refine the surface before extrusion
189: suffix: 1
190: requires: triangle
191: args: -init_dm_plex_dim 2 -init_dm_refine 2 -dm_extrude 1 -dm_view
193: test: # Parallel run
194: suffix: 2
195: requires: triangle
196: nsize: 5
197: args: -init_dm_plex_dim 2 -init_dm_refine 3 -petscpartitioner_type simple -dm_extrude 3 -dm_view
199: test: # adaptively refine the surface before extrusion
200: suffix: 3
201: requires: triangle
202: 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
203: output_file: output/empty.out
205: TEST*/