Actual source code: dmperiodicity.c
1: #include <petsc/private/dmimpl.h>
3: #include <petscdmplex.h>
5: /*@C
6: DMGetPeriodicity - Get the description of mesh periodicity
8: Input Parameter:
9: . dm - The DM object
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
16: Level: developer
18: .seealso: `DMGetPeriodicity()`
19: @*/
20: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **Lstart, const PetscReal **L)
21: {
23: if (maxCell) *maxCell = dm->maxCell;
24: if (Lstart) *Lstart = dm->Lstart;
25: if (L) *L = dm->L;
26: return 0;
27: }
29: /*@C
30: DMSetPeriodicity - Set the description of mesh periodicity
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
38: Level: developer
40: .seealso: `DMGetPeriodicity()`
41: @*/
42: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal Lstart[], const PetscReal L[])
43: {
44: PetscInt dim, d;
50: DMGetDimension(dm, &dim);
51: if (maxCell) {
52: if (!dm->maxCell) 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: PetscFree(dm->maxCell);
56: dm->maxCell = NULL;
57: }
58: if (Lstart) {
59: if (!dm->Lstart) 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: PetscFree(dm->Lstart);
63: dm->Lstart = NULL;
64: }
65: if (L) {
66: if (!dm->L) 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: PetscFree(dm->L);
70: dm->L = NULL;
71: }
73: return 0;
74: }
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.
79: Input Parameters:
80: + dm - The DM
81: . in - The input coordinate point (dim numbers)
82: - endpoint - Include the endpoint L_i
84: Output Parameter:
85: . out - The localized coordinate point
87: Level: developer
89: .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
90: @*/
91: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
92: {
93: PetscInt dim, d;
95: 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: return 0;
112: }
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.
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)
123: Output Parameter:
124: . out - The localized coordinate point
126: Level: developer
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
130: .seealso: `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
131: */
132: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
133: {
134: PetscInt d;
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: return 0;
148: }
150: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
151: {
152: PetscInt d;
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: return 0;
166: }
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.
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)
178: Output Parameter:
179: . out - The localized coordinate in + out
181: Level: developer
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
185: .seealso: `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()`
186: */
187: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
188: {
189: PetscInt d;
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];
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];
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: return 0;
209: }
211: /*@
212: DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process
214: Not collective
216: Input Parameter:
217: . dm - The DM
219: Output Parameter:
220: areLocalized - True if localized
222: Level: developer
224: .seealso: `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()`
225: @*/
226: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized)
227: {
230: *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE;
231: return 0;
232: }
234: /*@
235: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
237: Collective on dm
239: Input Parameter:
240: . dm - The DM
242: Output Parameter:
243: areLocalized - True if localized
245: Level: developer
247: .seealso: `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()`
248: @*/
249: PetscErrorCode DMGetCoordinatesLocalized(DM dm, PetscBool *areLocalized)
250: {
251: PetscBool localized;
255: DMGetCoordinatesLocalizedLocal(dm, &localized);
256: MPIU_Allreduce(&localized, areLocalized, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
257: return 0;
258: }
260: /*@
261: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
263: Collective on dm
265: Input Parameter:
266: . dm - The DM
268: Level: developer
270: .seealso: `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()`
271: @*/
272: PetscErrorCode DMLocalizeCoordinates(DM dm)
273: {
274: DM cdm, cdgdm, cplex, plex;
275: PetscSection cs, csDG;
276: Vec coordinates, cVec;
277: PetscScalar *coordsDG, *anchor, *localized;
278: const PetscReal *Lstart, *L;
279: PetscInt Nc, vStart, vEnd, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, bs, coordSize;
280: PetscBool isLocalized, sparseLocalize = dm->sparseLocalize, useDG = PETSC_FALSE, useDGGlobal;
281: PetscInt maxHeight = 0, h;
282: PetscInt *pStart = NULL, *pEnd = NULL;
283: MPI_Comm comm;
286: DMGetPeriodicity(dm, NULL, &Lstart, &L);
287: /* Cannot automatically localize without L and maxCell right now */
288: if (!L) return 0;
289: PetscObjectGetComm((PetscObject)dm, &comm);
290: DMGetCoordinatesLocalized(dm, &isLocalized);
291: if (isLocalized) return 0;
293: DMGetCoordinateDM(dm, &cdm);
294: DMConvert(dm, DMPLEX, &plex);
295: DMConvert(cdm, DMPLEX, &cplex);
296: if (cplex) {
297: DMPlexGetDepthStratum(cplex, 0, &vStart, &vEnd);
298: DMPlexGetMaxProjectionHeight(cplex, &maxHeight);
299: 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: 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: DMGetCoordinatesLocal(dm, &coordinates);
311: DMGetCoordinateSection(dm, &cs);
312: VecGetBlockSize(coordinates, &bs);
313: PetscSectionGetChart(cs, &sStart, &sEnd);
315: PetscSectionCreate(comm, &csDG);
316: PetscSectionSetNumFields(csDG, 1);
317: PetscSectionGetFieldComponents(cs, 0, &Nc);
318: PetscSectionSetFieldComponents(csDG, 0, Nc);
319: PetscSectionSetChart(csDG, newStart, newEnd);
322: 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;
327: for (c = cStart; c < cEnd; ++c) {
328: PetscScalar *cellCoords = NULL;
329: DMPolytopeType ct;
330: PetscInt dof, d, p;
332: DMPlexGetCellType(plex, c, &ct);
333: if (ct == DM_POLYTOPE_FV_GHOST) continue;
334: DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords);
336: for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[d];
337: for (p = 0; p < dof / Nc; ++p) {
338: 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: PetscSectionSetDof(csDG, c, dof);
346: PetscSectionSetFieldDof(csDG, c, 0, dof);
347: }
348: DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords);
349: }
350: }
351: MPI_Allreduce(&useDG, &useDGGlobal, 1, MPIU_BOOL, MPI_LOR, comm);
352: if (!useDGGlobal) goto end;
354: PetscSectionSetUp(csDG);
355: PetscSectionGetStorageSize(csDG, &coordSize);
356: VecCreate(PETSC_COMM_SELF, &cVec);
357: PetscObjectSetName((PetscObject)cVec, "coordinates");
358: VecSetBlockSize(cVec, bs);
359: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
360: VecSetType(cVec, VECSTANDARD);
361: VecGetArray(cVec, &coordsDG);
362: for (h = 0; h <= maxHeight; h++) {
363: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
365: for (c = cStart; c < cEnd; ++c) {
366: PetscScalar *cellCoords = NULL;
367: PetscInt p = 0, q, dof, cdof, d, offDG;
369: PetscSectionGetDof(csDG, c, &cdof);
370: if (!cdof) continue;
371: DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords);
372: 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: 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: }
387: DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords);
388: }
389: }
390: VecRestoreArray(cVec, &coordsDG);
391: DMClone(cdm, &cdgdm);
392: DMSetCellCoordinateDM(dm, cdgdm);
393: DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG);
394: DMSetCellCoordinatesLocal(dm, cVec);
395: 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;
404: DMGetField(cdm, 0, NULL, (PetscObject *)&fe);
405: PetscObjectGetClassId((PetscObject)fe, &id);
406: if (id == PETSCFE_CLASSID) {
407: PetscFEGetBasisSpace(fe, &P);
408: PetscObjectReference((PetscObject)P);
409: PetscFEGetDualSpace(fe, &Q);
410: PetscDualSpaceDuplicate(Q, &dgQ);
411: PetscDualSpaceLagrangeSetContinuity(dgQ, PETSC_FALSE);
412: PetscDualSpaceSetUp(dgQ);
413: PetscFEGetQuadrature(fe, &q);
414: PetscObjectReference((PetscObject)q);
415: PetscFEGetFaceQuadrature(fe, &fq);
416: PetscObjectReference((PetscObject)fq);
417: PetscFECreateFromSpaces(P, dgQ, q, fq, &dgfe);
418: DMSetField(cdgdm, 0, NULL, (PetscObject)dgfe);
419: PetscFEDestroy(&dgfe);
420: DMCreateDS(cdgdm);
421: }
422: }
423: DMDestroy(&cdgdm);
425: end:
426: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
427: DMRestoreWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart);
428: PetscSectionDestroy(&csDG);
429: DMDestroy(&plex);
430: DMDestroy(&cplex);
431: return 0;
432: }