1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
DMPlexLabelToVolumeConstraint(DM dm,DMLabel adaptLabel,PetscInt cStart,PetscInt cEnd,PetscReal refRatio,PetscReal maxVolumes[])3 static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
4 {
5 PetscInt dim, c;
6
7 PetscFunctionBegin;
8 PetscCall(DMGetDimension(dm, &dim));
9 refRatio = refRatio == (PetscReal)PETSC_DEFAULT ? (PetscReal)(1 << dim) : refRatio;
10 for (c = cStart; c < cEnd; c++) {
11 PetscReal vol;
12 PetscInt closureSize = 0, cl;
13 PetscInt *closure = NULL;
14 PetscBool anyRefine = PETSC_FALSE;
15 PetscBool anyCoarsen = PETSC_FALSE;
16 PetscBool anyKeep = PETSC_FALSE;
17
18 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL));
19 maxVolumes[c - cStart] = vol;
20 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
21 for (cl = 0; cl < closureSize * 2; cl += 2) {
22 const PetscInt point = closure[cl];
23 PetscInt refFlag;
24
25 PetscCall(DMLabelGetValue(adaptLabel, point, &refFlag));
26 switch (refFlag) {
27 case DM_ADAPT_REFINE:
28 anyRefine = PETSC_TRUE;
29 break;
30 case DM_ADAPT_COARSEN:
31 anyCoarsen = PETSC_TRUE;
32 break;
33 case DM_ADAPT_KEEP:
34 anyKeep = PETSC_TRUE;
35 break;
36 case DM_ADAPT_DETERMINE:
37 break;
38 default:
39 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "DMPlex does not support refinement flag %" PetscInt_FMT, refFlag);
40 }
41 if (anyRefine) break;
42 }
43 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
44 if (anyRefine) {
45 maxVolumes[c - cStart] = vol / refRatio;
46 } else if (anyKeep) {
47 maxVolumes[c - cStart] = vol;
48 } else if (anyCoarsen) {
49 maxVolumes[c - cStart] = vol * refRatio;
50 }
51 }
52 PetscFunctionReturn(PETSC_SUCCESS);
53 }
54
DMPlexLabelToMetricConstraint(DM dm,DMLabel adaptLabel,PetscInt cStart,PetscInt cEnd,PetscInt vStart,PetscInt vEnd,PetscReal refRatio,Vec * metricVec)55 static PetscErrorCode DMPlexLabelToMetricConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscInt vStart, PetscInt vEnd, PetscReal refRatio, Vec *metricVec)
56 {
57 DM udm, coordDM;
58 PetscSection coordSection;
59 Vec coordinates, mb, mx;
60 Mat A;
61 PetscScalar *metric, *eqns;
62 const PetscReal coarseRatio = refRatio == (PetscReal)PETSC_DEFAULT ? PetscSqr(0.5) : 1 / refRatio;
63 PetscInt dim, Nv, Neq, c, v;
64
65 PetscFunctionBegin;
66 PetscCall(DMPlexUninterpolate(dm, &udm));
67 PetscCall(DMGetDimension(dm, &dim));
68 PetscCall(DMGetCoordinateDM(dm, &coordDM));
69 PetscCall(DMGetLocalSection(coordDM, &coordSection));
70 PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
71 Nv = vEnd - vStart;
72 PetscCall(VecCreateSeq(PETSC_COMM_SELF, Nv * PetscSqr(dim), metricVec));
73 PetscCall(VecGetArray(*metricVec, &metric));
74 Neq = (dim * (dim + 1)) / 2;
75 PetscCall(PetscMalloc1(PetscSqr(Neq), &eqns));
76 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A));
77 PetscCall(MatCreateVecs(A, &mx, &mb));
78 PetscCall(VecSet(mb, 1.0));
79 for (c = cStart; c < cEnd; ++c) {
80 const PetscScalar *sol;
81 PetscScalar *cellCoords = NULL;
82 PetscReal e[3], vol;
83 const PetscInt *cone;
84 PetscInt coneSize, cl, i, j, d, r;
85
86 PetscCall(DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords));
87 /* Only works for simplices */
88 for (i = 0, r = 0; i < dim + 1; ++i) {
89 for (j = 0; j < i; ++j, ++r) {
90 for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i * dim + d] - cellCoords[j * dim + d]);
91 /* FORTRAN ORDERING */
92 switch (dim) {
93 case 2:
94 eqns[0 * Neq + r] = PetscSqr(e[0]);
95 eqns[1 * Neq + r] = 2.0 * e[0] * e[1];
96 eqns[2 * Neq + r] = PetscSqr(e[1]);
97 break;
98 case 3:
99 eqns[0 * Neq + r] = PetscSqr(e[0]);
100 eqns[1 * Neq + r] = 2.0 * e[0] * e[1];
101 eqns[2 * Neq + r] = 2.0 * e[0] * e[2];
102 eqns[3 * Neq + r] = PetscSqr(e[1]);
103 eqns[4 * Neq + r] = 2.0 * e[1] * e[2];
104 eqns[5 * Neq + r] = PetscSqr(e[2]);
105 break;
106 }
107 }
108 }
109 PetscCall(MatSetUnfactored(A));
110 PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords));
111 PetscCall(MatLUFactor(A, NULL, NULL, NULL));
112 PetscCall(MatSolve(A, mb, mx));
113 PetscCall(VecGetArrayRead(mx, &sol));
114 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL));
115 PetscCall(DMPlexGetCone(udm, c, &cone));
116 PetscCall(DMPlexGetConeSize(udm, c, &coneSize));
117 for (cl = 0; cl < coneSize; ++cl) {
118 const PetscInt v = cone[cl] - vStart;
119
120 if (dim == 2) {
121 metric[v * 4 + 0] += vol * coarseRatio * sol[0];
122 metric[v * 4 + 1] += vol * coarseRatio * sol[1];
123 metric[v * 4 + 2] += vol * coarseRatio * sol[1];
124 metric[v * 4 + 3] += vol * coarseRatio * sol[2];
125 } else {
126 metric[v * 9 + 0] += vol * coarseRatio * sol[0];
127 metric[v * 9 + 1] += vol * coarseRatio * sol[1];
128 metric[v * 9 + 3] += vol * coarseRatio * sol[1];
129 metric[v * 9 + 2] += vol * coarseRatio * sol[2];
130 metric[v * 9 + 6] += vol * coarseRatio * sol[2];
131 metric[v * 9 + 4] += vol * coarseRatio * sol[3];
132 metric[v * 9 + 5] += vol * coarseRatio * sol[4];
133 metric[v * 9 + 7] += vol * coarseRatio * sol[4];
134 metric[v * 9 + 8] += vol * coarseRatio * sol[5];
135 }
136 }
137 PetscCall(VecRestoreArrayRead(mx, &sol));
138 }
139 for (v = 0; v < Nv; ++v) {
140 const PetscInt *support;
141 PetscInt supportSize, s;
142 PetscReal vol, totVol = 0.0;
143
144 PetscCall(DMPlexGetSupport(udm, v + vStart, &support));
145 PetscCall(DMPlexGetSupportSize(udm, v + vStart, &supportSize));
146 for (s = 0; s < supportSize; ++s) {
147 PetscCall(DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL));
148 totVol += vol;
149 }
150 for (s = 0; s < PetscSqr(dim); ++s) metric[v * PetscSqr(dim) + s] /= totVol;
151 }
152 PetscCall(PetscFree(eqns));
153 PetscCall(VecRestoreArray(*metricVec, &metric));
154 PetscCall(VecDestroy(&mx));
155 PetscCall(VecDestroy(&mb));
156 PetscCall(MatDestroy(&A));
157 PetscCall(DMDestroy(&udm));
158 PetscFunctionReturn(PETSC_SUCCESS);
159 }
160
161 /*
162 Contains the list of registered DMPlexGenerators routines
163 */
DMPlexRefine_Internal(DM dm,PETSC_UNUSED Vec metric,DMLabel adaptLabel,PETSC_UNUSED DMLabel rgLabel,DM * dmRefined)164 PetscErrorCode DMPlexRefine_Internal(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *dmRefined)
165 {
166 DMGeneratorFunctionList fl;
167 PetscErrorCode (*refine)(DM, PetscReal *, DM *);
168 PetscErrorCode (*adapt)(DM, Vec, DMLabel, DMLabel, DM *);
169 PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *);
170 char genname[PETSC_MAX_PATH_LEN], *name = NULL;
171 PetscReal refinementLimit;
172 PetscReal *maxVolumes;
173 PetscInt dim, cStart, cEnd, c;
174 PetscBool flg, flg2, localized;
175
176 PetscFunctionBegin;
177 PetscCall(DMGetCoordinatesLocalized(dm, &localized));
178 PetscCall(DMPlexGetRefinementLimit(dm, &refinementLimit));
179 PetscCall(DMPlexGetRefinementFunction(dm, &refinementFunc));
180 if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) PetscFunctionReturn(PETSC_SUCCESS);
181 PetscCall(DMGetDimension(dm, &dim));
182 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
183 PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_adaptor", genname, sizeof(genname), &flg));
184 if (flg) name = genname;
185 else {
186 PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_generator", genname, sizeof(genname), &flg2));
187 if (flg2) name = genname;
188 }
189
190 fl = DMGenerateList;
191 if (name) {
192 while (fl) {
193 PetscCall(PetscStrcmp(fl->name, name, &flg));
194 if (flg) {
195 refine = fl->refine;
196 adapt = fl->adapt;
197 goto gotit;
198 }
199 fl = fl->next;
200 }
201 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid refiner %s not registered", name);
202 } else {
203 while (fl) {
204 if (fl->dim < 0 || dim - 1 == fl->dim) {
205 refine = fl->refine;
206 adapt = fl->adapt;
207 goto gotit;
208 }
209 fl = fl->next;
210 }
211 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No grid refiner of dimension %" PetscInt_FMT " registered", dim);
212 }
213
214 gotit:
215 switch (dim) {
216 case 1:
217 case 2:
218 case 3:
219 if (adapt) {
220 PetscCall((*adapt)(dm, NULL, adaptLabel, NULL, dmRefined));
221 } else {
222 PetscCall(PetscMalloc1(cEnd - cStart, &maxVolumes));
223 if (adaptLabel) {
224 PetscCall(DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes));
225 } else if (refinementFunc) {
226 for (c = cStart; c < cEnd; ++c) {
227 PetscReal vol, centroid[3];
228
229 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL));
230 PetscCall((*refinementFunc)(centroid, &maxVolumes[c - cStart]));
231 }
232 } else {
233 for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = refinementLimit;
234 }
235 PetscCall((*refine)(dm, maxVolumes, dmRefined));
236 PetscCall(PetscFree(maxVolumes));
237 }
238 break;
239 default:
240 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %" PetscInt_FMT " is not supported.", dim);
241 }
242 PetscCall(DMCopyDisc(dm, *dmRefined));
243 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmRefined));
244 if (localized) PetscCall(DMLocalizeCoordinates(*dmRefined));
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
DMPlexCoarsen_Internal(DM dm,PETSC_UNUSED Vec metric,DMLabel adaptLabel,PETSC_UNUSED DMLabel rgLabel,DM * dmCoarsened)248 PetscErrorCode DMPlexCoarsen_Internal(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *dmCoarsened)
249 {
250 Vec metricVec;
251 PetscInt cStart, cEnd, vStart, vEnd;
252 DMLabel bdLabel = NULL;
253 char bdLabelName[PETSC_MAX_PATH_LEN], rgLabelName[PETSC_MAX_PATH_LEN];
254 PetscBool localized, flg;
255
256 PetscFunctionBegin;
257 PetscCall(DMGetCoordinatesLocalized(dm, &localized));
258 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
259 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
260 PetscCall(DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec));
261 PetscCall(PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg));
262 if (flg) PetscCall(DMGetLabel(dm, bdLabelName, &bdLabel));
263 PetscCall(PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_rg_label", rgLabelName, sizeof(rgLabelName), &flg));
264 if (flg) PetscCall(DMGetLabel(dm, rgLabelName, &rgLabel));
265 PetscCall(DMAdaptMetric(dm, metricVec, bdLabel, rgLabel, dmCoarsened));
266 PetscCall(VecDestroy(&metricVec));
267 PetscCall(DMCopyDisc(dm, *dmCoarsened));
268 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmCoarsened));
269 if (localized) PetscCall(DMLocalizeCoordinates(*dmCoarsened));
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
DMAdaptLabel_Plex(DM dm,PETSC_UNUSED Vec metric,DMLabel adaptLabel,PETSC_UNUSED DMLabel rgLabel,DM * dmAdapted)273 PetscErrorCode DMAdaptLabel_Plex(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *dmAdapted)
274 {
275 IS flagIS;
276 const PetscInt *flags;
277 PetscInt defFlag, minFlag, maxFlag, numFlags, f;
278
279 PetscFunctionBegin;
280 PetscCall(DMLabelGetDefaultValue(adaptLabel, &defFlag));
281 minFlag = defFlag;
282 maxFlag = defFlag;
283 PetscCall(DMLabelGetValueIS(adaptLabel, &flagIS));
284 PetscCall(ISGetLocalSize(flagIS, &numFlags));
285 PetscCall(ISGetIndices(flagIS, &flags));
286 for (f = 0; f < numFlags; ++f) {
287 const PetscInt flag = flags[f];
288
289 minFlag = PetscMin(minFlag, flag);
290 maxFlag = PetscMax(maxFlag, flag);
291 }
292 PetscCall(ISRestoreIndices(flagIS, &flags));
293 PetscCall(ISDestroy(&flagIS));
294 {
295 PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
296
297 minMaxFlag[0] = minFlag;
298 minMaxFlag[1] = -maxFlag;
299 PetscCallMPI(MPIU_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
300 minFlag = minMaxFlagGlobal[0];
301 maxFlag = -minMaxFlagGlobal[1];
302 }
303 if (minFlag == maxFlag) {
304 switch (minFlag) {
305 case DM_ADAPT_DETERMINE:
306 *dmAdapted = NULL;
307 break;
308 case DM_ADAPT_REFINE:
309 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE));
310 PetscCall(DMRefine(dm, MPI_COMM_NULL, dmAdapted));
311 break;
312 case DM_ADAPT_COARSEN:
313 PetscCall(DMCoarsen(dm, MPI_COMM_NULL, dmAdapted));
314 break;
315 default:
316 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "DMPlex does not support refinement flag %" PetscInt_FMT, minFlag);
317 }
318 } else {
319 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_FALSE));
320 PetscCall(DMPlexRefine_Internal(dm, NULL, adaptLabel, NULL, dmAdapted));
321 }
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324