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: `DM`
19: @*/
20: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **Lstart, const PetscReal **L)
21: {
22: PetscFunctionBegin;
24: if (maxCell) *maxCell = dm->maxCell;
25: if (Lstart) *Lstart = dm->Lstart;
26: if (L) *L = dm->L;
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: /*@
31: DMSetPeriodicity - Set the description of mesh periodicity
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
39: Level: developer
41: .seealso: `DM`, `DMGetPeriodicity()`
42: @*/
43: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal Lstart[], const PetscReal L[])
44: {
45: PetscInt dim, d;
47: PetscFunctionBegin;
49: if (maxCell) PetscAssertPointer(maxCell, 2);
50: if (Lstart) PetscAssertPointer(Lstart, 3);
51: if (L) PetscAssertPointer(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(PETSC_SUCCESS);
76: }
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.
81: Input Parameters:
82: + dm - The `DM`
83: . in - The input coordinate point (dim numbers)
84: - endpoint - Include the endpoint L_i
86: Output Parameter:
87: . out - The localized coordinate point (dim numbers)
89: Level: developer
91: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
92: @*/
93: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
94: {
95: PetscInt dim, d;
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) out[d] = in[d] - dm->L[d] * PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]);
112: }
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*
118: 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.
120: Input Parameters:
121: + dm - The `DM`
122: . dim - The spatial dimension
123: . anchor - The anchor point, the input point can be no more than maxCell away from it
124: - in - The input coordinate point (dim numbers)
126: Output Parameter:
127: . out - The localized coordinate point (dim numbers)
129: Level: developer
131: Note:
132: 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: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
135: */
136: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
137: {
138: PetscInt d;
140: PetscFunctionBegin;
141: if (!dm->maxCell) {
142: for (d = 0; d < dim; ++d) out[d] = in[d];
143: } else {
144: for (d = 0; d < dim; ++d) {
145: if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
146: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
147: } else {
148: out[d] = in[d];
149: }
150: }
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
156: {
157: PetscInt d;
159: PetscFunctionBegin;
160: if (!dm->maxCell) {
161: for (d = 0; d < dim; ++d) out[d] = in[d];
162: } else {
163: for (d = 0; d < dim; ++d) {
164: if ((dm->L[d] > 0.0) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
165: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
166: } else {
167: out[d] = in[d];
168: }
169: }
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*
175: 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: Input Parameters:
178: + dm - The `DM`
179: . dim - The spatial dimension
180: . anchor - The anchor point, the input point can be no more than maxCell away from it
181: . in - The input coordinate delta (dim numbers)
182: - out - The input coordinate point (dim numbers)
184: Output Parameter:
185: . out - The localized coordinate in + out
187: Level: developer
189: Note:
190: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually one of the vertices on a containing cell
192: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()`
193: */
194: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
195: {
196: PetscInt d;
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];
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];
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(PETSC_SUCCESS);
217: }
219: /*@
220: DMGetCoordinatesLocalizedLocal - Check if the `DM` coordinates have been localized for cells on this process
222: Not Collective
224: Input Parameter:
225: . dm - The `DM`
227: Output Parameter:
228: . areLocalized - `PETSC_TRUE` if localized
230: Level: developer
232: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()`
233: @*/
234: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized)
235: {
236: PetscFunctionBegin;
238: PetscAssertPointer(areLocalized, 2);
239: *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: DMGetCoordinatesLocalized - Check if the `DM` coordinates have been localized for cells
246: Collective
248: Input Parameter:
249: . dm - The `DM`
251: Output Parameter:
252: . areLocalized - `PETSC_TRUE` if localized
254: Level: developer
256: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()`
257: @*/
258: PetscErrorCode DMGetCoordinatesLocalized(DM dm, PetscBool *areLocalized)
259: {
260: PetscBool localized;
262: PetscFunctionBegin;
264: PetscAssertPointer(areLocalized, 2);
265: PetscCall(DMGetCoordinatesLocalizedLocal(dm, &localized));
266: PetscCall(MPIU_Allreduce(&localized, areLocalized, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: DMGetSparseLocalize - Check if the `DM` coordinates should be localized only for cells near the periodic boundary.
273: Not collective
275: Input Parameter:
276: . dm - The `DM`
278: Output Parameter:
279: . sparse - `PETSC_TRUE` if ony cells near the periodic boundary are localized
281: Level: intermediate
283: .seealso: `DMSetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
284: @*/
285: PetscErrorCode DMGetSparseLocalize(DM dm, PetscBool *sparse)
286: {
287: PetscFunctionBegin;
289: PetscAssertPointer(sparse, 2);
290: *sparse = dm->sparseLocalize;
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: DMSetSparseLocalize - Set the flag indicating that `DM` coordinates should be localized only for cells near the periodic boundary.
297: Logically collective
299: Input Parameters:
300: + dm - The `DM`
301: - sparse - `PETSC_TRUE` if ony cells near the periodic boundary are localized
303: Level: intermediate
305: .seealso: `DMGetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
306: @*/
307: PetscErrorCode DMSetSparseLocalize(DM dm, PetscBool sparse)
308: {
309: PetscFunctionBegin;
312: dm->sparseLocalize = sparse;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@
317: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
319: Collective
321: Input Parameter:
322: . dm - The `DM`
324: Level: developer
326: .seealso: `DM`, `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()`
327: @*/
328: PetscErrorCode DMLocalizeCoordinates(DM dm)
329: {
330: DM cdm, cdgdm, cplex, plex;
331: PetscSection cs, csDG;
332: Vec coordinates, cVec;
333: PetscScalar *coordsDG, *anchor, *localized;
334: const PetscReal *Lstart, *L;
335: PetscInt Nc, vStart, vEnd, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, bs, coordSize;
336: PetscBool isLocalized, sparseLocalize, useDG = PETSC_FALSE, useDGGlobal;
337: PetscInt maxHeight = 0, h;
338: PetscInt *pStart = NULL, *pEnd = NULL;
339: MPI_Comm comm;
341: PetscFunctionBegin;
343: PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
344: PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
345: /* Cannot automatically localize without L and maxCell right now */
346: if (!L) PetscFunctionReturn(PETSC_SUCCESS);
347: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
348: PetscCall(DMGetCoordinatesLocalized(dm, &isLocalized));
349: if (isLocalized) PetscFunctionReturn(PETSC_SUCCESS);
351: PetscCall(DMGetCoordinateDM(dm, &cdm));
352: PetscCall(DMConvert(dm, DMPLEX, &plex));
353: PetscCall(DMConvert(cdm, DMPLEX, &cplex));
354: if (cplex) {
355: PetscCall(DMPlexGetDepthStratum(cplex, 0, &vStart, &vEnd));
356: PetscCall(DMPlexGetMaxProjectionHeight(cplex, &maxHeight));
357: PetscCall(DMGetWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
358: pEnd = &pStart[maxHeight + 1];
359: newStart = vStart;
360: newEnd = vEnd;
361: for (h = 0; h <= maxHeight; h++) {
362: PetscCall(DMPlexGetHeightStratum(cplex, h, &pStart[h], &pEnd[h]));
363: newStart = PetscMin(newStart, pStart[h]);
364: newEnd = PetscMax(newEnd, pEnd[h]);
365: }
366: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
367: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
368: PetscCheck(coordinates, comm, PETSC_ERR_SUP, "Missing local coordinate vector");
369: PetscCall(DMGetCoordinateSection(dm, &cs));
370: PetscCall(VecGetBlockSize(coordinates, &bs));
371: PetscCall(PetscSectionGetChart(cs, &sStart, &sEnd));
373: PetscCall(PetscSectionCreate(comm, &csDG));
374: PetscCall(PetscSectionSetNumFields(csDG, 1));
375: PetscCall(PetscSectionGetFieldComponents(cs, 0, &Nc));
376: PetscCall(PetscSectionSetFieldComponents(csDG, 0, Nc));
377: PetscCall(PetscSectionSetChart(csDG, newStart, newEnd));
378: PetscCheck(bs == Nc, comm, PETSC_ERR_ARG_INCOMP, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " number of components", bs, Nc);
380: PetscCall(DMGetWorkArray(dm, 2 * Nc, MPIU_SCALAR, &anchor));
381: localized = &anchor[Nc];
382: for (h = 0; h <= maxHeight; h++) {
383: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
385: for (c = cStart; c < cEnd; ++c) {
386: PetscScalar *cellCoords = NULL;
387: DMPolytopeType ct;
388: PetscInt dof, d, p;
390: PetscCall(DMPlexGetCellType(plex, c, &ct));
391: if (ct == DM_POLYTOPE_FV_GHOST) continue;
392: PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
393: 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);
394: for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[d];
395: for (p = 0; p < dof / Nc; ++p) {
396: PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], localized));
397: for (d = 0; d < Nc; ++d)
398: if (cellCoords[p * Nc + d] != localized[d]) break;
399: if (d < Nc) break;
400: }
401: if (p < dof / Nc) useDG = PETSC_TRUE;
402: if (p < dof / Nc || !sparseLocalize) {
403: PetscCall(PetscSectionSetDof(csDG, c, dof));
404: PetscCall(PetscSectionSetFieldDof(csDG, c, 0, dof));
405: }
406: PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
407: }
408: }
409: PetscCall(MPIU_Allreduce(&useDG, &useDGGlobal, 1, MPIU_BOOL, MPI_LOR, comm));
410: if (!useDGGlobal) goto end;
412: PetscCall(PetscSectionSetUp(csDG));
413: PetscCall(PetscSectionGetStorageSize(csDG, &coordSize));
414: PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
415: PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
416: PetscCall(VecSetBlockSize(cVec, bs));
417: PetscCall(VecSetSizes(cVec, coordSize, PETSC_DETERMINE));
418: PetscCall(VecSetType(cVec, VECSTANDARD));
419: PetscCall(VecGetArray(cVec, &coordsDG));
420: for (h = 0; h <= maxHeight; h++) {
421: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
423: for (c = cStart; c < cEnd; ++c) {
424: PetscScalar *cellCoords = NULL;
425: PetscInt p = 0, q, dof, cdof, d, offDG;
427: PetscCall(PetscSectionGetDof(csDG, c, &cdof));
428: if (!cdof) continue;
429: PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
430: PetscCall(PetscSectionGetOffset(csDG, c, &offDG));
431: // TODO The coordinates are set in closure order, which might not be the tensor order
432: for (q = 0; q < dof / Nc; ++q) {
433: // Select a trial anchor
434: for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[q * Nc + d];
435: for (p = 0; p < dof / Nc; ++p) {
436: PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], &coordsDG[offDG + p * Nc]));
437: // We need the cell to fit into the torus [lower, lower+L)
438: for (d = 0; d < Nc; ++d)
439: 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;
440: if (d < Nc) break;
441: }
442: if (p == dof / Nc) break;
443: }
444: 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 + " : "");
445: PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
446: }
447: }
448: PetscCall(VecRestoreArray(cVec, &coordsDG));
449: PetscCall(DMClone(cdm, &cdgdm));
450: PetscCall(DMSetCellCoordinateDM(dm, cdgdm));
451: PetscCall(DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG));
452: PetscCall(DMSetCellCoordinatesLocal(dm, cVec));
453: PetscCall(VecDestroy(&cVec));
454: // Convert the discretization
455: {
456: PetscFE fe, dgfe;
457: PetscSpace P;
458: PetscDualSpace Q, dgQ;
459: PetscQuadrature q, fq;
460: PetscClassId id;
462: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
463: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
464: if (id == PETSCFE_CLASSID) {
465: PetscCall(PetscFEGetBasisSpace(fe, &P));
466: PetscCall(PetscObjectReference((PetscObject)P));
467: PetscCall(PetscFEGetDualSpace(fe, &Q));
468: PetscCall(PetscDualSpaceDuplicate(Q, &dgQ));
469: PetscCall(PetscDualSpaceLagrangeSetContinuity(dgQ, PETSC_FALSE));
470: PetscCall(PetscDualSpaceSetUp(dgQ));
471: PetscCall(PetscFEGetQuadrature(fe, &q));
472: PetscCall(PetscObjectReference((PetscObject)q));
473: PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
474: PetscCall(PetscObjectReference((PetscObject)fq));
475: PetscCall(PetscFECreateFromSpaces(P, dgQ, q, fq, &dgfe));
476: PetscCall(DMSetField(cdgdm, 0, NULL, (PetscObject)dgfe));
477: PetscCall(PetscFEDestroy(&dgfe));
478: PetscCall(DMCreateDS(cdgdm));
479: }
480: }
481: PetscCall(DMDestroy(&cdgdm));
483: end:
484: PetscCall(DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor));
485: PetscCall(DMRestoreWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
486: PetscCall(PetscSectionDestroy(&csDG));
487: PetscCall(DMDestroy(&plex));
488: PetscCall(DMDestroy(&cplex));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }