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: }