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