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