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 CHKERRQ(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 CHKERRQ(DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL)); 19 maxVolumes[c - cStart] = vol; 20 CHKERRQ(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 CHKERRQ(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 %D", refFlag); 37 } 38 if (anyRefine) break; 39 } 40 CHKERRQ(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 CHKERRQ(DMPlexUninterpolate(dm, &udm)); 64 CHKERRQ(DMGetDimension(dm, &dim)); 65 CHKERRQ(DMGetCoordinateDM(dm, &coordDM)); 66 CHKERRQ(DMGetLocalSection(coordDM, &coordSection)); 67 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 68 Nv = vEnd - vStart; 69 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec)); 70 CHKERRQ(VecGetArray(*metricVec, &metric)); 71 Neq = (dim*(dim+1))/2; 72 CHKERRQ(PetscMalloc1(PetscSqr(Neq), &eqns)); 73 CHKERRQ(MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A)); 74 CHKERRQ(MatCreateVecs(A, &mx, &mb)); 75 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetUnfactored(A)); 107 CHKERRQ(DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords)); 108 CHKERRQ(MatLUFactor(A, NULL, NULL, NULL)); 109 CHKERRQ(MatSolve(A, mb, mx)); 110 CHKERRQ(VecGetArrayRead(mx, &sol)); 111 CHKERRQ(DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL)); 112 CHKERRQ(DMPlexGetCone(udm, c, &cone)); 113 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexGetSupport(udm, v+vStart, &support)); 142 CHKERRQ(DMPlexGetSupportSize(udm, v+vStart, &supportSize)); 143 for (s = 0; s < supportSize; ++s) {CHKERRQ(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 CHKERRQ(PetscFree(eqns)); 147 CHKERRQ(VecRestoreArray(*metricVec, &metric)); 148 CHKERRQ(VecDestroy(&mx)); 149 CHKERRQ(VecDestroy(&mb)); 150 CHKERRQ(MatDestroy(&A)); 151 CHKERRQ(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 CHKERRQ(DMGetCoordinatesLocalized(dm, &localized)); 172 CHKERRQ(DMPlexGetRefinementLimit(dm, &refinementLimit)); 173 CHKERRQ(DMPlexGetRefinementFunction(dm, &refinementFunc)); 174 if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) PetscFunctionReturn(0); 175 CHKERRQ(DMGetDimension(dm, &dim)); 176 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 177 CHKERRQ(PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_adaptor", genname, sizeof(genname), &flg)); 178 if (flg) name = genname; 179 else { 180 CHKERRQ(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 CHKERRQ(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 %D registered",dim); 206 } 207 208 gotit: 209 switch (dim) { 210 case 1: 211 case 2: 212 case 3: 213 if (adapt) { 214 CHKERRQ((*adapt)(dm, NULL, adaptLabel, NULL, dmRefined)); 215 } else { 216 CHKERRQ(PetscMalloc1(cEnd - cStart, &maxVolumes)); 217 if (adaptLabel) { 218 CHKERRQ(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 CHKERRQ(DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL)); 224 CHKERRQ((*refinementFunc)(centroid, &maxVolumes[c-cStart])); 225 } 226 } else { 227 for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit; 228 } 229 CHKERRQ((*refine)(dm, maxVolumes, dmRefined)); 230 CHKERRQ(PetscFree(maxVolumes)); 231 } 232 break; 233 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %D is not supported.", dim); 234 } 235 CHKERRQ(DMCopyDisc(dm, *dmRefined)); 236 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *dmRefined)); 237 if (localized) CHKERRQ(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 CHKERRQ(DMGetCoordinatesLocalized(dm, &localized)); 251 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 252 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 253 CHKERRQ(DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec)); 254 CHKERRQ(PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg)); 255 if (flg) CHKERRQ(DMGetLabel(dm, bdLabelName, &bdLabel)); 256 CHKERRQ(PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_rg_label", rgLabelName, sizeof(rgLabelName), &flg)); 257 if (flg) CHKERRQ(DMGetLabel(dm, rgLabelName, &rgLabel)); 258 CHKERRQ(DMAdaptMetric(dm, metricVec, bdLabel, rgLabel, dmCoarsened)); 259 CHKERRQ(VecDestroy(&metricVec)); 260 CHKERRQ(DMCopyDisc(dm, *dmCoarsened)); 261 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *dmCoarsened)); 262 if (localized) CHKERRQ(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 CHKERRQ(DMLabelGetDefaultValue(adaptLabel, &defFlag)); 274 minFlag = defFlag; 275 maxFlag = defFlag; 276 CHKERRQ(DMLabelGetValueIS(adaptLabel, &flagIS)); 277 CHKERRQ(ISGetLocalSize(flagIS, &numFlags)); 278 CHKERRQ(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 CHKERRQ(ISRestoreIndices(flagIS, &flags)); 286 CHKERRQ(ISDestroy(&flagIS)); 287 { 288 PetscInt minMaxFlag[2], minMaxFlagGlobal[2]; 289 290 minMaxFlag[0] = minFlag; 291 minMaxFlag[1] = -maxFlag; 292 CHKERRMPI(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 CHKERRQ(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 302 CHKERRQ(DMRefine(dm, MPI_COMM_NULL, dmAdapted));break; 303 case DM_ADAPT_COARSEN: 304 CHKERRQ(DMCoarsen(dm, MPI_COMM_NULL, dmAdapted));break; 305 default: 306 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D", minFlag); 307 } 308 } else { 309 CHKERRQ(DMPlexSetRefinementUniform(dm, PETSC_FALSE)); 310 CHKERRQ(DMPlexRefine_Internal(dm, NULL, adaptLabel, NULL, dmAdapted)); 311 } 312 PetscFunctionReturn(0); 313 } 314