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: }