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