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:   Not collective

 10:   Input Parameter:
 11: . dm - The `DM` object

 13:   Output Parameters:
 14: + maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
 15: . Lstart  - If we assume the mesh is a torus, this is the start of each coordinate, or `NULL` for 0.0
 16: - L       - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0

 18:   Level: developer

 20: .seealso: `DM`
 21: @*/
 22: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal *maxCell[], const PetscReal *Lstart[], const PetscReal *L[])
 23: {
 24:   PetscFunctionBegin;
 26:   if (maxCell) *maxCell = dm->maxCell;
 27:   if (Lstart) *Lstart = dm->Lstart;
 28:   if (L) *L = dm->L;
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: /*@
 33:   DMSetPeriodicity - Set the description of mesh periodicity

 35:   Logically Collective

 37:   Input Parameters:
 38: + dm      - The `DM` object
 39: . 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.
 40: . Lstart  - If we assume the mesh is a torus, this is the start of each coordinate, or `NULL` for 0.0
 41: - L       - If we assume the mesh is a torus, this is the length of each coordinate, otherwise it is < 0.0

 43:   Level: developer

 45: .seealso: `DM`, `DMGetPeriodicity()`
 46: @*/
 47: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal Lstart[], const PetscReal L[])
 48: {
 49:   PetscInt dim, d;

 51:   PetscFunctionBegin;
 53:   if (maxCell) PetscAssertPointer(maxCell, 2);
 54:   if (Lstart) PetscAssertPointer(Lstart, 3);
 55:   if (L) PetscAssertPointer(L, 4);
 56:   PetscCall(DMGetDimension(dm, &dim));
 57:   if (maxCell) {
 58:     if (!dm->maxCell) PetscCall(PetscMalloc1(dim, &dm->maxCell));
 59:     for (d = 0; d < dim; ++d) dm->maxCell[d] = maxCell[d];
 60:   } else { /* remove maxCell information to disable automatic computation of localized vertices */
 61:     PetscCall(PetscFree(dm->maxCell));
 62:     dm->maxCell = NULL;
 63:   }
 64:   if (Lstart) {
 65:     if (!dm->Lstart) PetscCall(PetscMalloc1(dim, &dm->Lstart));
 66:     for (d = 0; d < dim; ++d) dm->Lstart[d] = Lstart[d];
 67:   } else { /* remove L information to disable automatic computation of localized vertices */
 68:     PetscCall(PetscFree(dm->Lstart));
 69:     dm->Lstart = NULL;
 70:   }
 71:   if (L) {
 72:     if (!dm->L) PetscCall(PetscMalloc1(dim, &dm->L));
 73:     for (d = 0; d < dim; ++d) dm->L[d] = L[d];
 74:   } else { /* remove L information to disable automatic computation of localized vertices */
 75:     PetscCall(PetscFree(dm->L));
 76:     dm->L = NULL;
 77:   }
 78:   PetscCheck((dm->maxCell && dm->L) || (!dm->maxCell && !dm->L), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot set only one of maxCell/L");
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: /*@
 83:   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.

 85:   Input Parameters:
 86: + dm       - The `DM`
 87: . in       - The input coordinate point (dim numbers)
 88: - endpoint - Include the endpoint L_i

 90:   Output Parameter:
 91: . out - The localized coordinate point (dim numbers)

 93:   Level: developer

 95: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
 96: @*/
 97: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
 98: {
 99:   PetscInt dim, d;

101:   PetscFunctionBegin;
102:   PetscCall(DMGetCoordinateDim(dm, &dim));
103:   if (!dm->maxCell) {
104:     for (d = 0; d < dim; ++d) out[d] = in[d];
105:   } else {
106:     if (endpoint) {
107:       for (d = 0; d < dim; ++d) {
108:         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)) {
109:           out[d] = in[d] - dm->L[d] * (PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]) - 1);
110:         } else {
111:           out[d] = in[d] - dm->L[d] * PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]);
112:         }
113:       }
114:     } else {
115:       for (d = 0; d < dim; ++d) out[d] = in[d] - dm->L[d] * PetscFloorReal(PetscRealPart(in[d]) / dm->L[d]);
116:     }
117:   }
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: static PetscErrorCode DMLocalizeCoordinatePerDimension_Internal(DM dm, PetscInt d, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
122: {
123:   PetscFunctionBegin;
124:   if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
125:     out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
126:   } else {
127:     out[d] = in[d];
128:   }
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: /*
133:   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.

135:   Input Parameters:
136: + dm     - The `DM`
137: . dim    - The spatial dimension
138: . anchor - The anchor point, the input point can be no more than maxCell away from it
139: - in     - The input coordinate point (dim numbers)

141:   Output Parameter:
142: . out - The localized coordinate point (dim numbers)

144:   Level: developer

146:   Note:
147:   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

149: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
150: */
151: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
152: {
153:   PetscInt d;

155:   PetscFunctionBegin;
156:   if (!dm->maxCell) {
157:     for (d = 0; d < dim; ++d) out[d] = in[d];
158:   } else {
159:     for (d = 0; d < dim; ++d) PetscCall(DMLocalizeCoordinatePerDimension_Internal(dm, d, anchor, in, out));
160:   }
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
165: {
166:   PetscInt d;

168:   PetscFunctionBegin;
169:   if (!dm->maxCell) {
170:     for (d = 0; d < dim; ++d) out[d] = in[d];
171:   } else {
172:     for (d = 0; d < dim; ++d) {
173:       if ((dm->L[d] > 0.0) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
174:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
175:       } else {
176:         out[d] = in[d];
177:       }
178:     }
179:   }
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*
184:   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.

186:   Input Parameters:
187: + dm     - The `DM`
188: . dim    - The spatial dimension
189: . anchor - The anchor point, the input point can be no more than maxCell away from it
190: . in     - The input coordinate delta (dim numbers)
191: - out    - The input coordinate point (dim numbers)

193:   Output Parameter:
194: . out    - The localized coordinate in + out

196:   Level: developer

198:   Note:
199:   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

201: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()`
202: */
203: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
204: {
205:   PetscInt d;

207:   PetscFunctionBegin;
208:   if (!dm->maxCell) {
209:     for (d = 0; d < dim; ++d) out[d] += in[d];
210:   } else {
211:     for (d = 0; d < dim; ++d) {
212:       const PetscReal maxC = dm->maxCell[d];

214:       if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > maxC)) {
215:         const PetscScalar newCoord = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];

217:         if (PetscAbsScalar(newCoord - anchor[d]) > maxC)
218:           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]));
219:         out[d] += newCoord;
220:       } else {
221:         out[d] += in[d];
222:       }
223:     }
224:   }
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: /*@
229:   DMGetCoordinatesLocalizedLocal - Check if the `DM` coordinates have been localized for cells on this process

231:   Not Collective

233:   Input Parameter:
234: . dm - The `DM`

236:   Output Parameter:
237: . areLocalized - `PETSC_TRUE` if localized

239:   Level: developer

241: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()`
242: @*/
243: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized)
244: {
245:   PetscFunctionBegin;
247:   PetscAssertPointer(areLocalized, 2);
248:   *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE;
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: /*@
253:   DMGetCoordinatesLocalized - Check if the `DM` coordinates have been localized for cells

255:   Collective

257:   Input Parameter:
258: . dm - The `DM`

260:   Output Parameter:
261: . areLocalized - `PETSC_TRUE` if localized

263:   Level: developer

265: .seealso: `DM`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()`
266: @*/
267: PetscErrorCode DMGetCoordinatesLocalized(DM dm, PetscBool *areLocalized)
268: {
269:   PetscBool localized;

271:   PetscFunctionBegin;
273:   PetscAssertPointer(areLocalized, 2);
274:   PetscCall(DMGetCoordinatesLocalizedLocal(dm, &localized));
275:   PetscCallMPI(MPIU_Allreduce(&localized, areLocalized, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
276:   PetscFunctionReturn(PETSC_SUCCESS);
277: }

279: /*@
280:   DMGetSparseLocalize - Check if the `DM` coordinates should be localized only for cells near the periodic boundary.

282:   Not collective

284:   Input Parameter:
285: . dm - The `DM`

287:   Output Parameter:
288: . sparse - `PETSC_TRUE` if only cells near the periodic boundary are localized

290:   Level: intermediate

292: .seealso: `DMSetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
293: @*/
294: PetscErrorCode DMGetSparseLocalize(DM dm, PetscBool *sparse)
295: {
296:   PetscFunctionBegin;
298:   PetscAssertPointer(sparse, 2);
299:   *sparse = dm->sparseLocalize;
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: /*@
304:   DMSetSparseLocalize - Set the flag indicating that `DM` coordinates should be localized only for cells near the periodic boundary.

306:   Collective

308:   Input Parameters:
309: + dm     - The `DM`
310: - sparse - `PETSC_TRUE` if only cells near the periodic boundary are localized

312:   Level: intermediate

314:   Note:
315:   If previous cell coordinates existed with a different sparse localization then these will be destroyed.

317: .seealso: `DMGetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
318: @*/
319: PetscErrorCode DMSetSparseLocalize(DM dm, PetscBool sparse)
320: {
321:   PetscBool prevSparse, areLocalized;

323:   PetscFunctionBegin;
326:   // Blow away previously localized coordinates if they were created with the other sparse localize option
327:   PetscCall(DMGetSparseLocalize(dm, &prevSparse));
328:   PetscCall(DMGetCoordinatesLocalized(dm, &areLocalized));
329:   if (prevSparse != sparse && areLocalized) PetscCall(DMDestroyCoordinates_Internal(&dm->coordinates[1]));
330:   dm->sparseLocalize = sparse;
331:   PetscFunctionReturn(PETSC_SUCCESS);
332: }

334: /*@
335:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

337:   Collective

339:   Input Parameter:
340: . dm - The `DM`

342:   Level: developer

344: .seealso: `DM`, `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()`
345: @*/
346: PetscErrorCode DMLocalizeCoordinates(DM dm)
347: {
348:   DM               cdm, cdgdm, cplex, plex;
349:   PetscSection     cs, csDG;
350:   Vec              coordinates, cVec;
351:   PetscScalar     *coordsDG, *anchor, *localized;
352:   const PetscReal *Lstart, *L, *maxCell;
353:   PetscInt         Nc, vStart, vEnd, sStart, sEnd, newStart = PETSC_INT_MAX, newEnd = PETSC_INT_MIN, bs, coordSize;
354:   PetscBool        isLocalized, sparseLocalize, useDG = PETSC_FALSE, useDGGlobal;
355:   PetscInt         maxHeight = 0, h;
356:   PetscInt        *pStart = NULL, *pEnd = NULL;
357:   MPI_Comm         comm;

359:   PetscFunctionBegin;
361:   PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
362:   PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
363:   /* Cannot automatically localize without L and maxCell right now */
364:   if (!L) PetscFunctionReturn(PETSC_SUCCESS);
365:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
366:   PetscCall(DMGetCoordinatesLocalized(dm, &isLocalized));
367:   if (isLocalized) PetscFunctionReturn(PETSC_SUCCESS);

369:   PetscCall(DMGetCoordinateDM(dm, &cdm));
370:   PetscCall(DMConvert(dm, DMPLEX, &plex));
371:   PetscCall(DMConvert(cdm, DMPLEX, &cplex));
372:   PetscCheck(cplex, comm, PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
373:   PetscCall(DMPlexGetDepthStratum(cplex, 0, &vStart, &vEnd));
374:   PetscCall(DMPlexGetMaxProjectionHeight(cplex, &maxHeight));
375:   PetscCall(DMGetWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
376:   pEnd     = &pStart[maxHeight + 1];
377:   newStart = vStart;
378:   newEnd   = vEnd;
379:   for (h = 0; h <= maxHeight; h++) {
380:     PetscCall(DMPlexGetHeightStratum(cplex, h, &pStart[h], &pEnd[h]));
381:     newStart = PetscMin(newStart, pStart[h]);
382:     newEnd   = PetscMax(newEnd, pEnd[h]);
383:   }
384:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
385:   PetscCheck(coordinates, comm, PETSC_ERR_SUP, "Missing local coordinate vector");
386:   PetscCall(DMGetCoordinateSection(dm, &cs));
387:   PetscCall(VecGetBlockSize(coordinates, &bs));
388:   PetscCall(PetscSectionGetChart(cs, &sStart, &sEnd));

390:   PetscCall(PetscSectionCreate(comm, &csDG));
391:   PetscCall(PetscSectionSetNumFields(csDG, 1));
392:   PetscCall(PetscSectionGetFieldComponents(cs, 0, &Nc));
393:   PetscCall(PetscSectionSetFieldComponents(csDG, 0, Nc));
394:   PetscCall(PetscSectionSetChart(csDG, newStart, newEnd));
395:   PetscCheck(bs == Nc, comm, PETSC_ERR_ARG_INCOMP, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " number of components", bs, Nc);
396:   for (PetscInt d = 0; d < Nc; ++d) {
397:     PetscCheck(L[d] < 0. || (maxCell[d] > 0. && maxCell[d] < 2 * L[d]), comm, PETSC_ERR_ARG_INCOMP, "Periodic length %g > max cell size %g in dimension %" PetscInt_FMT, (double)L[d], (double)maxCell[d], d);
398:   }

400:   PetscCall(DMGetWorkArray(dm, 2 * Nc, MPIU_SCALAR, &anchor));
401:   localized = &anchor[Nc];
402:   for (h = 0; h <= maxHeight; h++) {
403:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

405:     for (c = cStart; c < cEnd; ++c) {
406:       PetscScalar   *cellCoords = NULL;
407:       DMPolytopeType ct;
408:       PetscInt       dof, d, p;

410:       PetscCall(DMPlexGetCellType(plex, c, &ct));
411:       if (ct == DM_POLYTOPE_FV_GHOST) continue;
412:       PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
413:       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);
414:       for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[d];
415:       for (p = 0; p < dof / Nc; ++p) {
416:         PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], localized));
417:         for (d = 0; d < Nc; ++d)
418:           if (cellCoords[p * Nc + d] != localized[d]) break;
419:         if (d < Nc) break;
420:       }
421:       if (p < dof / Nc) useDG = PETSC_TRUE;
422:       if (p < dof / Nc || !sparseLocalize) {
423:         PetscCall(PetscSectionSetDof(csDG, c, dof));
424:         PetscCall(PetscSectionSetFieldDof(csDG, c, 0, dof));
425:       }
426:       PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
427:     }
428:   }
429:   PetscCallMPI(MPIU_Allreduce(&useDG, &useDGGlobal, 1, MPI_C_BOOL, MPI_LOR, comm));
430:   if (!useDGGlobal) goto end;

432:   PetscCall(PetscSectionSetUp(csDG));
433:   PetscCall(PetscSectionGetStorageSize(csDG, &coordSize));
434:   PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
435:   PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
436:   PetscCall(VecSetBlockSize(cVec, bs));
437:   PetscCall(VecSetSizes(cVec, coordSize, PETSC_DETERMINE));
438:   PetscCall(VecSetType(cVec, VECSTANDARD));
439:   PetscCall(VecGetArray(cVec, &coordsDG));
440:   for (h = 0; h <= maxHeight; h++) {
441:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

443:     for (c = cStart; c < cEnd; ++c) {
444:       PetscScalar *cellCoords = NULL;
445:       PetscInt     p          = 0, q, dof, cdof, d, offDG;

447:       PetscCall(PetscSectionGetDof(csDG, c, &cdof));
448:       if (!cdof) continue;
449:       PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
450:       PetscCall(PetscSectionGetOffset(csDG, c, &offDG));
451:       // Select an anchor for each dimension
452:       // TODO The coordinates are set in closure order, which might not be the tensor order
453:       for (d = 0; d < Nc; ++d) {
454:         const PetscReal start = Lstart ? Lstart[d] : 0.;

456:         for (q = 0; q < dof / Nc; ++q) {
457:           anchor[d] = cellCoords[q * Nc + d];
458:           // Map anchor into [Lstart, Lstart+L)
459:           if (L[d] > 0.0) {
460:             if (PetscRealPart(anchor[d]) < start) anchor[d] += L[d];
461:             else if (PetscRealPart(anchor[d]) > start + L[d]) anchor[d] -= L[d];
462:             PetscCheck(PetscRealPart(anchor[d]) >= start && PetscRealPart(anchor[d]) < start + L[d], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Anchor %g out of bounds for dimension %" PetscInt_FMT, (double)PetscRealPart(anchor[d]), d);
463:           }
464:           for (p = 0; p < dof / Nc; ++p) {
465:             PetscCall(DMLocalizeCoordinatePerDimension_Internal(dm, d, anchor, &cellCoords[p * Nc], &coordsDG[offDG + p * Nc]));
466:             // We need the cell to fit into the torus [lower, lower+L)
467:             if (L[d] > 0. && ((PetscRealPart(coordsDG[offDG + p * Nc + d]) < start) || (PetscRealPart(coordsDG[offDG + p * Nc + d]) > start + L[d]))) break;
468:           }
469:           // Suitable anchor found, continue
470:           if (p == dof / Nc) break;
471:         }
472:         // No suitable anchor found, abort
473:         if (q == dof / Nc) break;
474:       }
475:       PetscCheck(d == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " does not fit into the torus %s[0, L]", c, Lstart ? "Lstart + " : "");
476:       PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
477:     }
478:   }
479:   PetscCall(VecRestoreArray(cVec, &coordsDG));
480:   PetscUseTypeMethod(dm, createcellcoordinatedm, &cdgdm);
481:   PetscCall(DMSetCellCoordinateDM(dm, cdgdm));
482:   PetscCall(DMDestroy(&cdgdm));
483:   // Convert the discretization
484:   {
485:     PetscFE      fe;
486:     PetscClassId id;

488:     PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
489:     PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
490:     if (id == PETSCFE_CLASSID) {
491:       PetscSpace P;
492:       PetscInt   degree;

494:       PetscCall(PetscFEGetBasisSpace(fe, &P));
495:       PetscCall(PetscSpaceGetDegree(P, &degree, NULL));
496:       PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_TRUE, PETSC_FALSE));
497:     }
498:   }
499:   PetscCall(DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG));
500:   PetscCall(DMSetCellCoordinatesLocal(dm, cVec));
501:   PetscCall(VecDestroy(&cVec));

503: end:
504:   PetscCall(DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor));
505:   PetscCall(DMRestoreWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
506:   PetscCall(PetscSectionDestroy(&csDG));
507:   PetscCall(DMDestroy(&plex));
508:   PetscCall(DMDestroy(&cplex));
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }