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