Actual source code: plexadapt.c
1: #include <petsc/private/dmpleximpl.h>
3: static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
4: {
5: PetscInt dim, c;
7: PetscFunctionBegin;
8: PetscCall(DMGetDimension(dm, &dim));
9: refRatio = refRatio == (PetscReal)PETSC_DEFAULT ? (PetscReal)((PetscInt)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;
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;
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: }
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;
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;
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;
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;
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: }
161: /*
162: Contains the list of registered DMPlexGenerators routines
163: */
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;
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: }
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: }
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];
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: }
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;
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: }
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;
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];
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];
297: minMaxFlag[0] = minFlag;
298: minMaxFlag[1] = -maxFlag;
299: PetscCall(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: }