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 - L - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0 14 15 Level: developer 16 17 .seealso: `DMGetPeriodicity()` 18 @*/ 19 PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L) 20 { 21 PetscFunctionBegin; 22 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 23 if (L) *L = dm->L; 24 if (maxCell) *maxCell = dm->maxCell; 25 PetscFunctionReturn(0); 26 } 27 28 /*@C 29 DMSetPeriodicity - Set the description of mesh periodicity 30 31 Input Parameters: 32 + dm - The DM object 33 . 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. 34 - L - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0 35 36 Level: developer 37 38 .seealso: `DMGetPeriodicity()` 39 @*/ 40 PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[]) 41 { 42 PetscInt dim, d; 43 44 PetscFunctionBegin; 45 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 46 if (maxCell) {PetscValidRealPointer(maxCell,2);} 47 if (L) {PetscValidRealPointer(L,3);} 48 PetscCall(DMGetDimension(dm, &dim)); 49 if (maxCell) { 50 if (!dm->maxCell) PetscCall(PetscMalloc1(dim, &dm->maxCell)); 51 for (d = 0; d < dim; ++d) dm->maxCell[d] = maxCell[d]; 52 } else { /* remove maxCell information to disable automatic computation of localized vertices */ 53 PetscCall(PetscFree(dm->maxCell)); 54 dm->maxCell = NULL; 55 } 56 if (L) { 57 if (!dm->L) PetscCall(PetscMalloc1(dim, &dm->L)); 58 for (d = 0; d < dim; ++d) dm->L[d] = L[d]; 59 } else { /* remove L information to disable automatic computation of localized vertices */ 60 PetscCall(PetscFree(dm->L)); 61 dm->L = NULL; 62 } 63 PetscCheck((dm->maxCell && dm->L) || (!dm->maxCell && !dm->L), PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Cannot set only one of maxCell/L"); 64 PetscFunctionReturn(0); 65 } 66 67 /*@ 68 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. 69 70 Input Parameters: 71 + dm - The DM 72 . in - The input coordinate point (dim numbers) 73 - endpoint - Include the endpoint L_i 74 75 Output Parameter: 76 . out - The localized coordinate point 77 78 Level: developer 79 80 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()` 81 @*/ 82 PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[]) 83 { 84 PetscInt dim, d; 85 86 PetscFunctionBegin; 87 PetscCall(DMGetCoordinateDim(dm, &dim)); 88 if (!dm->maxCell) { 89 for (d = 0; d < dim; ++d) out[d] = in[d]; 90 } else { 91 if (endpoint) { 92 for (d = 0; d < dim; ++d) { 93 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)) { 94 out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1); 95 } else { 96 out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]); 97 } 98 } 99 } else { 100 for (d = 0; d < dim; ++d) { 101 out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]); 102 } 103 } 104 } 105 PetscFunctionReturn(0); 106 } 107 108 /* 109 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. 110 111 Input Parameters: 112 + dm - The DM 113 . dim - The spatial dimension 114 . anchor - The anchor point, the input point can be no more than maxCell away from it 115 - in - The input coordinate point (dim numbers) 116 117 Output Parameter: 118 . out - The localized coordinate point 119 120 Level: developer 121 122 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 123 124 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()` 125 */ 126 PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) 127 { 128 PetscInt d; 129 130 PetscFunctionBegin; 131 if (!dm->maxCell) { 132 for (d = 0; d < dim; ++d) out[d] = in[d]; 133 } else { 134 for (d = 0; d < dim; ++d) { 135 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) { 136 out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 137 } else { 138 out[d] = in[d]; 139 } 140 } 141 } 142 PetscFunctionReturn(0); 143 } 144 145 PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[]) 146 { 147 PetscInt d; 148 149 PetscFunctionBegin; 150 if (!dm->maxCell) { 151 for (d = 0; d < dim; ++d) out[d] = in[d]; 152 } else { 153 for (d = 0; d < dim; ++d) { 154 if ((dm->L[d] > 0.0) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) { 155 out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d]; 156 } else { 157 out[d] = in[d]; 158 } 159 } 160 } 161 PetscFunctionReturn(0); 162 } 163 164 /* 165 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. 166 167 Input Parameters: 168 + dm - The DM 169 . dim - The spatial dimension 170 . anchor - The anchor point, the input point can be no more than maxCell away from it 171 . in - The input coordinate delta (dim numbers) 172 - out - The input coordinate point (dim numbers) 173 174 Output Parameter: 175 . out - The localized coordinate in + out 176 177 Level: developer 178 179 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 180 181 .seealso: `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()` 182 */ 183 PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) 184 { 185 PetscInt d; 186 187 PetscFunctionBegin; 188 if (!dm->maxCell) { 189 for (d = 0; d < dim; ++d) out[d] += in[d]; 190 } else { 191 for (d = 0; d < dim; ++d) { 192 const PetscReal maxC = dm->maxCell[d]; 193 194 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > maxC)) { 195 const PetscScalar newCoord = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 196 197 if (PetscAbsScalar(newCoord - anchor[d]) > maxC) 198 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])); 199 out[d] += newCoord; 200 } else { 201 out[d] += in[d]; 202 } 203 } 204 } 205 PetscFunctionReturn(0); 206 } 207 208 /*@ 209 DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process 210 211 Not collective 212 213 Input Parameter: 214 . dm - The DM 215 216 Output Parameter: 217 areLocalized - True if localized 218 219 Level: developer 220 221 .seealso: `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()` 222 @*/ 223 PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized) 224 { 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 227 PetscValidBoolPointer(areLocalized, 2); 228 *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE; 229 PetscFunctionReturn(0); 230 } 231 232 /*@ 233 DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells 234 235 Collective on dm 236 237 Input Parameter: 238 . dm - The DM 239 240 Output Parameter: 241 areLocalized - True if localized 242 243 Level: developer 244 245 .seealso: `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()` 246 @*/ 247 PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized) 248 { 249 PetscBool localized; 250 251 PetscFunctionBegin; 252 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 253 PetscValidBoolPointer(areLocalized, 2); 254 PetscCall(DMGetCoordinatesLocalizedLocal(dm, &localized)); 255 PetscCall(MPIU_Allreduce(&localized, areLocalized, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm))); 256 PetscFunctionReturn(0); 257 } 258 259 /*@ 260 DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces 261 262 Collective on dm 263 264 Input Parameter: 265 . dm - The DM 266 267 Level: developer 268 269 .seealso: `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()` 270 @*/ 271 PetscErrorCode DMLocalizeCoordinates(DM dm) 272 { 273 DM cdm, cdgdm, cplex, plex; 274 PetscSection cs, csDG; 275 Vec coordinates, cVec; 276 PetscScalar *coordsDG, *anchor, *localized; 277 const PetscReal *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, &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) if (cellCoords[p*Nc + d] != localized[d]) break; 340 if (d < Nc) break; 341 } 342 if (p < dof/Nc) useDG = PETSC_TRUE; 343 if (p < dof/Nc || !sparseLocalize) { 344 PetscCall(PetscSectionSetDof(csDG, c, dof)); 345 PetscCall(PetscSectionSetFieldDof(csDG, c, 0, dof)); 346 } 347 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 348 } 349 } 350 PetscCallMPI(MPI_Allreduce(&useDG, &useDGGlobal, 1, MPIU_BOOL, MPI_LOR, comm)); 351 if (!useDGGlobal) goto end; 352 353 PetscCall(PetscSectionSetUp(csDG)); 354 PetscCall(PetscSectionGetStorageSize(csDG, &coordSize)); 355 PetscCall(VecCreate(PETSC_COMM_SELF, &cVec)); 356 PetscCall(PetscObjectSetName((PetscObject) cVec, "coordinates")); 357 PetscCall(VecSetBlockSize(cVec, bs)); 358 PetscCall(VecSetSizes(cVec, coordSize, PETSC_DETERMINE)); 359 PetscCall(VecSetType(cVec, VECSTANDARD)); 360 PetscCall(VecGetArray(cVec, &coordsDG)); 361 for (h = 0; h <= maxHeight; h++) { 362 PetscInt cStart = pStart[h], cEnd = pEnd[h], c; 363 364 for (c = cStart; c < cEnd; ++c) { 365 PetscScalar *cellCoords = NULL; 366 PetscInt p = 0, q, dof, cdof, d, offDG; 367 368 PetscCall(PetscSectionGetDof(csDG, c, &cdof)); 369 if (!cdof) continue; 370 PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 371 PetscCall(PetscSectionGetOffset(csDG, c, &offDG)); 372 // We need the cell to fit into [0, L] 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 for (d = 0; d < Nc; ++d) 380 if (L[d] > 0. && ((PetscRealPart(coordsDG[offDG + p*Nc + d]) < 0.) || (PetscRealPart(coordsDG[offDG + p*Nc + d]) > L[d]))) break; 381 if (d < Nc) break; 382 } 383 if (p == dof/Nc) break; 384 } 385 PetscCheck(p == dof/Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " does not fit into the torus [0, L]", c); 386 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords)); 387 } 388 } 389 PetscCall(VecRestoreArray(cVec, &coordsDG)); 390 PetscCall(DMClone(cdm, &cdgdm)); 391 PetscCall(DMSetCellCoordinateDM(dm, cdgdm)); 392 PetscCall(DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG)); 393 PetscCall(DMSetCellCoordinatesLocal(dm, cVec)); 394 PetscCall(VecDestroy(&cVec)); 395 // Convert the discretization 396 { 397 PetscFE fe, dgfe; 398 PetscSpace P; 399 PetscDualSpace Q, dgQ; 400 PetscQuadrature q, fq; 401 PetscClassId id; 402 403 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *) &fe)); 404 PetscCall(PetscObjectGetClassId((PetscObject) fe, &id)); 405 if (id == PETSCFE_CLASSID) { 406 PetscCall(PetscFEGetBasisSpace(fe, &P)); 407 PetscCall(PetscObjectReference((PetscObject) P)); 408 PetscCall(PetscFEGetDualSpace(fe, &Q)); 409 PetscCall(PetscDualSpaceDuplicate(Q, &dgQ)); 410 PetscCall(PetscDualSpaceLagrangeSetContinuity(dgQ, PETSC_FALSE)); 411 PetscCall(PetscDualSpaceSetUp(dgQ)); 412 PetscCall(PetscFEGetQuadrature(fe, &q)); 413 PetscCall(PetscObjectReference((PetscObject) q)); 414 PetscCall(PetscFEGetFaceQuadrature(fe, &fq)); 415 PetscCall(PetscObjectReference((PetscObject) fq)); 416 PetscCall(PetscFECreateFromSpaces(P, dgQ, q, fq, &dgfe)); 417 PetscCall(DMSetField(cdgdm, 0, NULL, (PetscObject) dgfe)); 418 PetscCall(PetscFEDestroy(&dgfe)); 419 PetscCall(DMCreateDS(cdgdm)); 420 } 421 } 422 PetscCall(DMDestroy(&cdgdm)); 423 424 end: 425 PetscCall(DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor)); 426 PetscCall(DMRestoreWorkArray(dm, 2*(maxHeight + 1), MPIU_INT, &pStart)); 427 PetscCall(PetscSectionDestroy(&csDG)); 428 PetscCall(DMDestroy(&plex)); 429 PetscCall(DMDestroy(&cplex)); 430 PetscFunctionReturn(0); 431 } 432