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