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