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 { 5 PetscInt dim, c; 6 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; 17 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; 24 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 } 54 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; 64 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; 85 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; 119 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; 143 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 } 160 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; 175 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 } 189 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 } 213 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]; 228 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 } 247 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; 255 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 } 272 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; 278 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]; 288 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]; 296 297 minMaxFlag[0] = minFlag; 298 minMaxFlag[1] = -maxFlag; 299 PetscCall(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 } 324