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