1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 2 3 #include <petscdmplex.h> 4 5 /*@C 6 DMGetPeriodicity - Get the description of mesh periodicity 7 8 Input Parameter: 9 . dm - The DM object 10 11 Output Parameters: 12 + maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates 13 . Lstart - If we assume the mesh is a torus, this is the start of each coordinate, or NULL for 0.0 14 - L - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0 15 16 Level: developer 17 18 .seealso: `DMGetPeriodicity()` 19 @*/ 20 PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **Lstart, const PetscReal **L) { 21 PetscFunctionBegin; 22 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 23 if (maxCell) *maxCell = dm->maxCell; 24 if (Lstart) *Lstart = dm->Lstart; 25 if (L) *L = dm->L; 26 PetscFunctionReturn(0); 27 } 28 29 /*@C 30 DMSetPeriodicity - Set the description of mesh periodicity 31 32 Input Parameters: 33 + dm - The DM object 34 . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates. Pass NULL to remove such information. 35 . Lstart - If we assume the mesh is a torus, this is the start of each coordinate, or NULL for 0.0 36 - L - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0 37 38 Level: developer 39 40 .seealso: `DMGetPeriodicity()` 41 @*/ 42 PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal Lstart[], const PetscReal L[]) { 43 PetscInt dim, d; 44 45 PetscFunctionBegin; 46 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 47 if (maxCell) { PetscValidRealPointer(maxCell, 2); } 48 if (Lstart) { PetscValidRealPointer(Lstart, 3); } 49 if (L) { PetscValidRealPointer(L, 4); } 50 PetscCall(DMGetDimension(dm, &dim)); 51 if (maxCell) { 52 if (!dm->maxCell) PetscCall(PetscMalloc1(dim, &dm->maxCell)); 53 for (d = 0; d < dim; ++d) dm->maxCell[d] = maxCell[d]; 54 } else { /* remove maxCell information to disable automatic computation of localized vertices */ 55 PetscCall(PetscFree(dm->maxCell)); 56 dm->maxCell = NULL; 57 } 58 if (Lstart) { 59 if (!dm->Lstart) PetscCall(PetscMalloc1(dim, &dm->Lstart)); 60 for (d = 0; d < dim; ++d) dm->Lstart[d] = Lstart[d]; 61 } else { /* remove L information to disable automatic computation of localized vertices */ 62 PetscCall(PetscFree(dm->Lstart)); 63 dm->Lstart = NULL; 64 } 65 if (L) { 66 if (!dm->L) PetscCall(PetscMalloc1(dim, &dm->L)); 67 for (d = 0; d < dim; ++d) dm->L[d] = L[d]; 68 } else { /* remove L information to disable automatic computation of localized vertices */ 69 PetscCall(PetscFree(dm->L)); 70 dm->L = NULL; 71 } 72 PetscCheck((dm->maxCell && dm->L) || (!dm->maxCell && !dm->L), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot set only one of maxCell/L"); 73 PetscFunctionReturn(0); 74 } 75 76 /*@ 77 DMLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension. 78 79 Input Parameters: 80 + dm - The DM 81 . in - The input coordinate point (dim numbers) 82 - endpoint - Include the endpoint L_i 83 84 Output Parameter: 85 . out - The localized coordinate point 86 87 Level: developer 88 89 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()` 90 @*/ 91 PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[]) { 92 PetscInt dim, d; 93 94 PetscFunctionBegin; 95 PetscCall(DMGetCoordinateDim(dm, &dim)); 96 if (!dm->maxCell) { 97 for (d = 0; d < dim; ++d) out[d] = in[d]; 98 } else { 99 if (endpoint) { 100 for (d = 0; d < dim; ++d) { 101 if ((PetscAbsReal(PetscRealPart(in[d]) / dm->L[d] - PetscFloorReal(PetscRealPart(in[d]) / dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d]) / dm->L[d] > PETSC_SMALL)) { 102 out[d] = in[d] - dm->L[d] * (PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]) - 1); 103 } else { 104 out[d] = in[d] - dm->L[d] * PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]); 105 } 106 } 107 } else { 108 for (d = 0; d < dim; ++d) { out[d] = in[d] - dm->L[d] * PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]); } 109 } 110 } 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 DMLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer. 116 117 Input Parameters: 118 + dm - The DM 119 . dim - The spatial dimension 120 . anchor - The anchor point, the input point can be no more than maxCell away from it 121 - in - The input coordinate point (dim numbers) 122 123 Output Parameter: 124 . out - The localized coordinate point 125 126 Level: developer 127 128 Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell 129 130 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()` 131 */ 132 PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) { 133 PetscInt d; 134 135 PetscFunctionBegin; 136 if (!dm->maxCell) { 137 for (d = 0; d < dim; ++d) out[d] = in[d]; 138 } else { 139 for (d = 0; d < dim; ++d) { 140 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) { 141 out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 142 } else { 143 out[d] = in[d]; 144 } 145 } 146 } 147 PetscFunctionReturn(0); 148 } 149 150 PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[]) { 151 PetscInt d; 152 153 PetscFunctionBegin; 154 if (!dm->maxCell) { 155 for (d = 0; d < dim; ++d) out[d] = in[d]; 156 } else { 157 for (d = 0; d < dim; ++d) { 158 if ((dm->L[d] > 0.0) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) { 159 out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d]; 160 } else { 161 out[d] = in[d]; 162 } 163 } 164 } 165 PetscFunctionReturn(0); 166 } 167 168 /* 169 DMLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer. 170 171 Input Parameters: 172 + dm - The DM 173 . dim - The spatial dimension 174 . anchor - The anchor point, the input point can be no more than maxCell away from it 175 . in - The input coordinate delta (dim numbers) 176 - out - The input coordinate point (dim numbers) 177 178 Output Parameter: 179 . out - The localized coordinate in + out 180 181 Level: developer 182 183 Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell 184 185 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()` 186 */ 187 PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) { 188 PetscInt d; 189 190 PetscFunctionBegin; 191 if (!dm->maxCell) { 192 for (d = 0; d < dim; ++d) out[d] += in[d]; 193 } else { 194 for (d = 0; d < dim; ++d) { 195 const PetscReal maxC = dm->maxCell[d]; 196 197 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > maxC)) { 198 const PetscScalar newCoord = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 199 200 if (PetscAbsScalar(newCoord - anchor[d]) > maxC) 201 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT "-Coordinate %g more than %g away from anchor %g", d, (double)PetscRealPart(in[d]), (double)maxC, (double)PetscRealPart(anchor[d])); 202 out[d] += newCoord; 203 } else { 204 out[d] += in[d]; 205 } 206 } 207 } 208 PetscFunctionReturn(0); 209 } 210 211 /*@ 212 DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process 213 214 Not collective 215 216 Input Parameter: 217 . dm - The DM 218 219 Output Parameter: 220 areLocalized - True if localized 221 222 Level: developer 223 224 .seealso: `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()` 225 @*/ 226 PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized) { 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 229 PetscValidBoolPointer(areLocalized, 2); 230 *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE; 231 PetscFunctionReturn(0); 232 } 233 234 /*@ 235 DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells 236 237 Collective on dm 238 239 Input Parameter: 240 . dm - The DM 241 242 Output Parameter: 243 areLocalized - True if localized 244 245 Level: developer 246 247 .seealso: `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()` 248 @*/ 249 PetscErrorCode DMGetCoordinatesLocalized(DM dm, PetscBool *areLocalized) { 250 PetscBool localized; 251 252 PetscFunctionBegin; 253 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 254 PetscValidBoolPointer(areLocalized, 2); 255 PetscCall(DMGetCoordinatesLocalizedLocal(dm, &localized)); 256 PetscCall(MPIU_Allreduce(&localized, areLocalized, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm))); 257 PetscFunctionReturn(0); 258 } 259 260 /*@ 261 DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces 262 263 Collective on dm 264 265 Input Parameter: 266 . dm - The DM 267 268 Level: developer 269 270 .seealso: `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()` 271 @*/ 272 PetscErrorCode DMLocalizeCoordinates(DM dm) { 273 DM cdm, cdgdm, cplex, plex; 274 PetscSection cs, csDG; 275 Vec coordinates, cVec; 276 PetscScalar *coordsDG, *anchor, *localized; 277 const PetscReal *Lstart, *L; 278 PetscInt Nc, vStart, vEnd, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, bs, coordSize; 279 PetscBool isLocalized, sparseLocalize = dm->sparseLocalize, useDG = PETSC_FALSE, useDGGlobal; 280 PetscInt maxHeight = 0, h; 281 PetscInt *pStart = NULL, *pEnd = NULL; 282 MPI_Comm comm; 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 286 PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L)); 287 /* Cannot automatically localize without L and maxCell right now */ 288 if (!L) PetscFunctionReturn(0); 289 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 290 PetscCall(DMGetCoordinatesLocalized(dm, &isLocalized)); 291 if (isLocalized) PetscFunctionReturn(0); 292 293 PetscCall(DMGetCoordinateDM(dm, &cdm)); 294 PetscCall(DMConvert(dm, DMPLEX, &plex)); 295 PetscCall(DMConvert(cdm, DMPLEX, &cplex)); 296 if (cplex) { 297 PetscCall(DMPlexGetDepthStratum(cplex, 0, &vStart, &vEnd)); 298 PetscCall(DMPlexGetMaxProjectionHeight(cplex, &maxHeight)); 299 PetscCall(DMGetWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart)); 300 pEnd = &pStart[maxHeight + 1]; 301 newStart = vStart; 302 newEnd = vEnd; 303 for (h = 0; h <= maxHeight; h++) { 304 PetscCall(DMPlexGetHeightStratum(cplex, h, &pStart[h], &pEnd[h])); 305 newStart = PetscMin(newStart, pStart[h]); 306 newEnd = PetscMax(newEnd, pEnd[h]); 307 } 308 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM"); 309 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 310 PetscCheck(coordinates, comm, PETSC_ERR_SUP, "Missing local coordinate vector"); 311 PetscCall(DMGetCoordinateSection(dm, &cs)); 312 PetscCall(VecGetBlockSize(coordinates, &bs)); 313 PetscCall(PetscSectionGetChart(cs, &sStart, &sEnd)); 314 315 PetscCall(PetscSectionCreate(comm, &csDG)); 316 PetscCall(PetscSectionSetNumFields(csDG, 1)); 317 PetscCall(PetscSectionGetFieldComponents(cs, 0, &Nc)); 318 PetscCall(PetscSectionSetFieldComponents(csDG, 0, Nc)); 319 PetscCall(PetscSectionSetChart(csDG, newStart, newEnd)); 320 PetscCheck(bs == Nc, comm, PETSC_ERR_ARG_INCOMP, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " number of components", bs, Nc); 321 322 PetscCall(DMGetWorkArray(dm, 2 * Nc, MPIU_SCALAR, &anchor)); 323 localized = &anchor[Nc]; 324 for (h = 0; h <= maxHeight; h++) { 325 PetscInt cStart = pStart[h], cEnd = pEnd[h], c; 326 327 for (c = cStart; c < cEnd; ++c) { 328 PetscScalar *cellCoords = NULL; 329 DMPolytopeType ct; 330 PetscInt dof, d, p; 331 332 PetscCall(DMPlexGetCellType(plex, c, &ct)); 333 if (ct == DM_POLYTOPE_FV_GHOST) continue; 334 PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 335 PetscCheck(!(dof % Nc), comm, PETSC_ERR_ARG_INCOMP, "Coordinate size on cell %" PetscInt_FMT " closure %" PetscInt_FMT " not divisible by %" PetscInt_FMT " number of components", c, dof, Nc); 336 for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[d]; 337 for (p = 0; p < dof / Nc; ++p) { 338 PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], localized)); 339 for (d = 0; d < Nc; ++d) 340 if (cellCoords[p * Nc + d] != localized[d]) break; 341 if (d < Nc) break; 342 } 343 if (p < dof / Nc) useDG = PETSC_TRUE; 344 if (p < dof / Nc || !sparseLocalize) { 345 PetscCall(PetscSectionSetDof(csDG, c, dof)); 346 PetscCall(PetscSectionSetFieldDof(csDG, c, 0, dof)); 347 } 348 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 349 } 350 } 351 PetscCallMPI(MPI_Allreduce(&useDG, &useDGGlobal, 1, MPIU_BOOL, MPI_LOR, comm)); 352 if (!useDGGlobal) goto end; 353 354 PetscCall(PetscSectionSetUp(csDG)); 355 PetscCall(PetscSectionGetStorageSize(csDG, &coordSize)); 356 PetscCall(VecCreate(PETSC_COMM_SELF, &cVec)); 357 PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates")); 358 PetscCall(VecSetBlockSize(cVec, bs)); 359 PetscCall(VecSetSizes(cVec, coordSize, PETSC_DETERMINE)); 360 PetscCall(VecSetType(cVec, VECSTANDARD)); 361 PetscCall(VecGetArray(cVec, &coordsDG)); 362 for (h = 0; h <= maxHeight; h++) { 363 PetscInt cStart = pStart[h], cEnd = pEnd[h], c; 364 365 for (c = cStart; c < cEnd; ++c) { 366 PetscScalar *cellCoords = NULL; 367 PetscInt p = 0, q, dof, cdof, d, offDG; 368 369 PetscCall(PetscSectionGetDof(csDG, c, &cdof)); 370 if (!cdof) continue; 371 PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 372 PetscCall(PetscSectionGetOffset(csDG, c, &offDG)); 373 // TODO The coordinates are set in closure order, which might not be the tensor order 374 for (q = 0; q < dof / Nc; ++q) { 375 // Select a trial anchor 376 for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[q * Nc + d]; 377 for (p = 0; p < dof / Nc; ++p) { 378 PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], &coordsDG[offDG + p * Nc])); 379 // We need the cell to fit into the torus [lower, lower+L) 380 for (d = 0; d < Nc; ++d) 381 if (L[d] > 0. && ((PetscRealPart(coordsDG[offDG + p * Nc + d]) < (Lstart ? Lstart[d] : 0.)) || (PetscRealPart(coordsDG[offDG + p * Nc + d]) > (Lstart ? Lstart[d] : 0.) + L[d]))) break; 382 if (d < Nc) break; 383 } 384 if (p == dof / Nc) break; 385 } 386 PetscCheck(p == dof / Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " does not fit into the torus %s[0, L]", c, Lstart ? "Lstart + " : ""); 387 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 388 } 389 } 390 PetscCall(VecRestoreArray(cVec, &coordsDG)); 391 PetscCall(DMClone(cdm, &cdgdm)); 392 PetscCall(DMSetCellCoordinateDM(dm, cdgdm)); 393 PetscCall(DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG)); 394 PetscCall(DMSetCellCoordinatesLocal(dm, cVec)); 395 PetscCall(VecDestroy(&cVec)); 396 // Convert the discretization 397 { 398 PetscFE fe, dgfe; 399 PetscSpace P; 400 PetscDualSpace Q, dgQ; 401 PetscQuadrature q, fq; 402 PetscClassId id; 403 404 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe)); 405 PetscCall(PetscObjectGetClassId((PetscObject)fe, &id)); 406 if (id == PETSCFE_CLASSID) { 407 PetscCall(PetscFEGetBasisSpace(fe, &P)); 408 PetscCall(PetscObjectReference((PetscObject)P)); 409 PetscCall(PetscFEGetDualSpace(fe, &Q)); 410 PetscCall(PetscDualSpaceDuplicate(Q, &dgQ)); 411 PetscCall(PetscDualSpaceLagrangeSetContinuity(dgQ, PETSC_FALSE)); 412 PetscCall(PetscDualSpaceSetUp(dgQ)); 413 PetscCall(PetscFEGetQuadrature(fe, &q)); 414 PetscCall(PetscObjectReference((PetscObject)q)); 415 PetscCall(PetscFEGetFaceQuadrature(fe, &fq)); 416 PetscCall(PetscObjectReference((PetscObject)fq)); 417 PetscCall(PetscFECreateFromSpaces(P, dgQ, q, fq, &dgfe)); 418 PetscCall(DMSetField(cdgdm, 0, NULL, (PetscObject)dgfe)); 419 PetscCall(PetscFEDestroy(&dgfe)); 420 PetscCall(DMCreateDS(cdgdm)); 421 } 422 } 423 PetscCall(DMDestroy(&cdgdm)); 424 425 end: 426 PetscCall(DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor)); 427 PetscCall(DMRestoreWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart)); 428 PetscCall(PetscSectionDestroy(&csDG)); 429 PetscCall(DMDestroy(&plex)); 430 PetscCall(DMDestroy(&cplex)); 431 PetscFunctionReturn(0); 432 } 433