1 static char help[] = "TDycore Mesh Examples\n\n";
2
3 #include <petscdmplex.h>
4
5 typedef struct {
6 PetscBool adapt; /* Flag for adaptation of the surface mesh */
7 } AppCtx;
8
ProcessOptions(MPI_Comm comm,AppCtx * options)9 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
10 {
11 PetscFunctionBeginUser;
12 options->adapt = PETSC_FALSE;
13
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 }
19
CreateDomainLabel(DM dm)20 static PetscErrorCode CreateDomainLabel(DM dm)
21 {
22 DMLabel label;
23 PetscInt cStart, cEnd, c;
24
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;
32
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 }
49
AdaptMesh(DM * dm,AppCtx * ctx)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;
61
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];
80
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];
92
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;
113
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 PetscCallMPI(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 }
145
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)146 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
147 {
148 PetscInt dim;
149
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 }
167
main(int argc,char ** argv)168 int main(int argc, char **argv)
169 {
170 DM dm;
171 AppCtx user;
172
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 }
181
182 /*TEST
183
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
188
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
193
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
199
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
204 output_file: output/empty.out
205
206 TEST*/
207