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 Not collective
9
10 Input Parameter:
11 . dm - The `DM` object
12
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
17
18 Level: developer
19
20 .seealso: `DM`
21 @*/
DMGetPeriodicity(DM dm,const PetscReal * maxCell[],const PetscReal * Lstart[],const PetscReal * L[])22 PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal *maxCell[], const PetscReal *Lstart[], const PetscReal *L[])
23 {
24 PetscFunctionBegin;
25 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
26 if (maxCell) *maxCell = dm->maxCell;
27 if (Lstart) *Lstart = dm->Lstart;
28 if (L) *L = dm->L;
29 PetscFunctionReturn(PETSC_SUCCESS);
30 }
31
32 /*@
33 DMSetPeriodicity - Set the description of mesh periodicity
34
35 Logically Collective
36
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
42
43 Level: developer
44
45 .seealso: `DM`, `DMGetPeriodicity()`
46 @*/
DMSetPeriodicity(DM dm,const PetscReal maxCell[],const PetscReal Lstart[],const PetscReal L[])47 PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal Lstart[], const PetscReal L[])
48 {
49 PetscInt dim, d;
50
51 PetscFunctionBegin;
52 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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 }
81
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.
84
85 Input Parameters:
86 + dm - The `DM`
87 . in - The input coordinate point (dim numbers)
88 - endpoint - Include the endpoint L_i
89
90 Output Parameter:
91 . out - The localized coordinate point (dim numbers)
92
93 Level: developer
94
95 .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
96 @*/
DMLocalizeCoordinate(DM dm,const PetscScalar in[],PetscBool endpoint,PetscScalar out[])97 PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
98 {
99 PetscInt dim, d;
100
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 }
120
121 /*
122 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.
123
124 Input Parameters:
125 + dm - The `DM`
126 . dim - The spatial dimension
127 . anchor - The anchor point, the input point can be no more than maxCell away from it
128 - in - The input coordinate point (dim numbers)
129
130 Output Parameter:
131 . out - The localized coordinate point (dim numbers)
132
133 Level: developer
134
135 Note:
136 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
137
138 .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeAddCoordinate()`
139 */
DMLocalizeCoordinate_Internal(DM dm,PetscInt dim,const PetscScalar anchor[],const PetscScalar in[],PetscScalar out[])140 PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
141 {
142 PetscInt d;
143
144 PetscFunctionBegin;
145 if (!dm->maxCell) {
146 for (d = 0; d < dim; ++d) out[d] = in[d];
147 } else {
148 for (d = 0; d < dim; ++d) {
149 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
150 out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
151 } else {
152 out[d] = in[d];
153 }
154 }
155 }
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
DMLocalizeCoordinateReal_Internal(DM dm,PetscInt dim,const PetscReal anchor[],const PetscReal in[],PetscReal out[])159 PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
160 {
161 PetscInt d;
162
163 PetscFunctionBegin;
164 if (!dm->maxCell) {
165 for (d = 0; d < dim; ++d) out[d] = in[d];
166 } else {
167 for (d = 0; d < dim; ++d) {
168 if ((dm->L[d] > 0.0) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
169 out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
170 } else {
171 out[d] = in[d];
172 }
173 }
174 }
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /*
179 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.
180
181 Input Parameters:
182 + dm - The `DM`
183 . dim - The spatial dimension
184 . anchor - The anchor point, the input point can be no more than maxCell away from it
185 . in - The input coordinate delta (dim numbers)
186 - out - The input coordinate point (dim numbers)
187
188 Output Parameter:
189 . out - The localized coordinate in + out
190
191 Level: developer
192
193 Note:
194 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
195
196 .seealso: `DM`, `DMLocalizeCoordinates()`, `DMLocalizeCoordinate()`
197 */
DMLocalizeAddCoordinate_Internal(DM dm,PetscInt dim,const PetscScalar anchor[],const PetscScalar in[],PetscScalar out[])198 PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
199 {
200 PetscInt d;
201
202 PetscFunctionBegin;
203 if (!dm->maxCell) {
204 for (d = 0; d < dim; ++d) out[d] += in[d];
205 } else {
206 for (d = 0; d < dim; ++d) {
207 const PetscReal maxC = dm->maxCell[d];
208
209 if ((dm->L[d] > 0.0) && (PetscAbsScalar(anchor[d] - in[d]) > maxC)) {
210 const PetscScalar newCoord = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
211
212 if (PetscAbsScalar(newCoord - anchor[d]) > maxC)
213 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]));
214 out[d] += newCoord;
215 } else {
216 out[d] += in[d];
217 }
218 }
219 }
220 PetscFunctionReturn(PETSC_SUCCESS);
221 }
222
223 /*@
224 DMGetCoordinatesLocalizedLocal - Check if the `DM` coordinates have been localized for cells on this process
225
226 Not Collective
227
228 Input Parameter:
229 . dm - The `DM`
230
231 Output Parameter:
232 . areLocalized - `PETSC_TRUE` if localized
233
234 Level: developer
235
236 .seealso: `DM`, `DMLocalizeCoordinates()`, `DMGetCoordinatesLocalized()`, `DMSetPeriodicity()`
237 @*/
DMGetCoordinatesLocalizedLocal(DM dm,PetscBool * areLocalized)238 PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm, PetscBool *areLocalized)
239 {
240 PetscFunctionBegin;
241 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
242 PetscAssertPointer(areLocalized, 2);
243 *areLocalized = dm->coordinates[1].dim < 0 ? PETSC_FALSE : PETSC_TRUE;
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
247 /*@
248 DMGetCoordinatesLocalized - Check if the `DM` coordinates have been localized for cells
249
250 Collective
251
252 Input Parameter:
253 . dm - The `DM`
254
255 Output Parameter:
256 . areLocalized - `PETSC_TRUE` if localized
257
258 Level: developer
259
260 .seealso: `DM`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`, `DMGetCoordinatesLocalizedLocal()`
261 @*/
DMGetCoordinatesLocalized(DM dm,PetscBool * areLocalized)262 PetscErrorCode DMGetCoordinatesLocalized(DM dm, PetscBool *areLocalized)
263 {
264 PetscBool localized;
265
266 PetscFunctionBegin;
267 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
268 PetscAssertPointer(areLocalized, 2);
269 PetscCall(DMGetCoordinatesLocalizedLocal(dm, &localized));
270 PetscCallMPI(MPIU_Allreduce(&localized, areLocalized, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
274 /*@
275 DMGetSparseLocalize - Check if the `DM` coordinates should be localized only for cells near the periodic boundary.
276
277 Not collective
278
279 Input Parameter:
280 . dm - The `DM`
281
282 Output Parameter:
283 . sparse - `PETSC_TRUE` if only cells near the periodic boundary are localized
284
285 Level: intermediate
286
287 .seealso: `DMSetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
288 @*/
DMGetSparseLocalize(DM dm,PetscBool * sparse)289 PetscErrorCode DMGetSparseLocalize(DM dm, PetscBool *sparse)
290 {
291 PetscFunctionBegin;
292 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
293 PetscAssertPointer(sparse, 2);
294 *sparse = dm->sparseLocalize;
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
298 /*@
299 DMSetSparseLocalize - Set the flag indicating that `DM` coordinates should be localized only for cells near the periodic boundary.
300
301 Logically collective
302
303 Input Parameters:
304 + dm - The `DM`
305 - sparse - `PETSC_TRUE` if only cells near the periodic boundary are localized
306
307 Level: intermediate
308
309 .seealso: `DMGetSparseLocalize()`, `DMLocalizeCoordinates()`, `DMSetPeriodicity()`
310 @*/
DMSetSparseLocalize(DM dm,PetscBool sparse)311 PetscErrorCode DMSetSparseLocalize(DM dm, PetscBool sparse)
312 {
313 PetscFunctionBegin;
314 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
315 PetscValidLogicalCollectiveBool(dm, sparse, 2);
316 dm->sparseLocalize = sparse;
317 PetscFunctionReturn(PETSC_SUCCESS);
318 }
319
320 /*@
321 DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
322
323 Collective
324
325 Input Parameter:
326 . dm - The `DM`
327
328 Level: developer
329
330 .seealso: `DM`, `DMSetPeriodicity()`, `DMLocalizeCoordinate()`, `DMLocalizeAddCoordinate()`
331 @*/
DMLocalizeCoordinates(DM dm)332 PetscErrorCode DMLocalizeCoordinates(DM dm)
333 {
334 DM cdm, cdgdm, cplex, plex;
335 PetscSection cs, csDG;
336 Vec coordinates, cVec;
337 PetscScalar *coordsDG, *anchor, *localized;
338 const PetscReal *Lstart, *L;
339 PetscInt Nc, vStart, vEnd, sStart, sEnd, newStart = PETSC_INT_MAX, newEnd = PETSC_INT_MIN, bs, coordSize;
340 PetscBool isLocalized, sparseLocalize, useDG = PETSC_FALSE, useDGGlobal;
341 PetscInt maxHeight = 0, h;
342 PetscInt *pStart = NULL, *pEnd = NULL;
343 MPI_Comm comm;
344
345 PetscFunctionBegin;
346 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
347 PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
348 PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
349 /* Cannot automatically localize without L and maxCell right now */
350 if (!L) PetscFunctionReturn(PETSC_SUCCESS);
351 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
352 PetscCall(DMGetCoordinatesLocalized(dm, &isLocalized));
353 if (isLocalized) PetscFunctionReturn(PETSC_SUCCESS);
354
355 PetscCall(DMGetCoordinateDM(dm, &cdm));
356 PetscCall(DMConvert(dm, DMPLEX, &plex));
357 PetscCall(DMConvert(cdm, DMPLEX, &cplex));
358 PetscCheck(cplex, comm, PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
359 PetscCall(DMPlexGetDepthStratum(cplex, 0, &vStart, &vEnd));
360 PetscCall(DMPlexGetMaxProjectionHeight(cplex, &maxHeight));
361 PetscCall(DMGetWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
362 pEnd = &pStart[maxHeight + 1];
363 newStart = vStart;
364 newEnd = vEnd;
365 for (h = 0; h <= maxHeight; h++) {
366 PetscCall(DMPlexGetHeightStratum(cplex, h, &pStart[h], &pEnd[h]));
367 newStart = PetscMin(newStart, pStart[h]);
368 newEnd = PetscMax(newEnd, pEnd[h]);
369 }
370 PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
371 PetscCheck(coordinates, comm, PETSC_ERR_SUP, "Missing local coordinate vector");
372 PetscCall(DMGetCoordinateSection(dm, &cs));
373 PetscCall(VecGetBlockSize(coordinates, &bs));
374 PetscCall(PetscSectionGetChart(cs, &sStart, &sEnd));
375
376 PetscCall(PetscSectionCreate(comm, &csDG));
377 PetscCall(PetscSectionSetNumFields(csDG, 1));
378 PetscCall(PetscSectionGetFieldComponents(cs, 0, &Nc));
379 PetscCall(PetscSectionSetFieldComponents(csDG, 0, Nc));
380 PetscCall(PetscSectionSetChart(csDG, newStart, newEnd));
381 PetscCheck(bs == Nc, comm, PETSC_ERR_ARG_INCOMP, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " number of components", bs, Nc);
382
383 PetscCall(DMGetWorkArray(dm, 2 * Nc, MPIU_SCALAR, &anchor));
384 localized = &anchor[Nc];
385 for (h = 0; h <= maxHeight; h++) {
386 PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
387
388 for (c = cStart; c < cEnd; ++c) {
389 PetscScalar *cellCoords = NULL;
390 DMPolytopeType ct;
391 PetscInt dof, d, p;
392
393 PetscCall(DMPlexGetCellType(plex, c, &ct));
394 if (ct == DM_POLYTOPE_FV_GHOST) continue;
395 PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
396 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);
397 for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[d];
398 for (p = 0; p < dof / Nc; ++p) {
399 PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], localized));
400 for (d = 0; d < Nc; ++d)
401 if (cellCoords[p * Nc + d] != localized[d]) break;
402 if (d < Nc) break;
403 }
404 if (p < dof / Nc) useDG = PETSC_TRUE;
405 if (p < dof / Nc || !sparseLocalize) {
406 PetscCall(PetscSectionSetDof(csDG, c, dof));
407 PetscCall(PetscSectionSetFieldDof(csDG, c, 0, dof));
408 }
409 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
410 }
411 }
412 PetscCallMPI(MPIU_Allreduce(&useDG, &useDGGlobal, 1, MPI_C_BOOL, MPI_LOR, comm));
413 if (!useDGGlobal) goto end;
414
415 PetscCall(PetscSectionSetUp(csDG));
416 PetscCall(PetscSectionGetStorageSize(csDG, &coordSize));
417 PetscCall(VecCreate(PETSC_COMM_SELF, &cVec));
418 PetscCall(PetscObjectSetName((PetscObject)cVec, "coordinates"));
419 PetscCall(VecSetBlockSize(cVec, bs));
420 PetscCall(VecSetSizes(cVec, coordSize, PETSC_DETERMINE));
421 PetscCall(VecSetType(cVec, VECSTANDARD));
422 PetscCall(VecGetArray(cVec, &coordsDG));
423 for (h = 0; h <= maxHeight; h++) {
424 PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
425
426 for (c = cStart; c < cEnd; ++c) {
427 PetscScalar *cellCoords = NULL;
428 PetscInt p = 0, q, dof, cdof, d, offDG;
429
430 PetscCall(PetscSectionGetDof(csDG, c, &cdof));
431 if (!cdof) continue;
432 PetscCall(DMPlexVecGetClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
433 PetscCall(PetscSectionGetOffset(csDG, c, &offDG));
434 // TODO The coordinates are set in closure order, which might not be the tensor order
435 for (q = 0; q < dof / Nc; ++q) {
436 // Select a trial anchor
437 for (d = 0; d < Nc; ++d) anchor[d] = cellCoords[q * Nc + d];
438 for (p = 0; p < dof / Nc; ++p) {
439 PetscCall(DMLocalizeCoordinate_Internal(dm, Nc, anchor, &cellCoords[p * Nc], &coordsDG[offDG + p * Nc]));
440 // We need the cell to fit into the torus [lower, lower+L)
441 for (d = 0; d < Nc; ++d)
442 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;
443 if (d < Nc) break;
444 }
445 if (p == dof / Nc) break;
446 }
447 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 + " : "");
448 PetscCall(DMPlexVecRestoreClosure(cplex, cs, coordinates, c, &dof, &cellCoords));
449 }
450 }
451 PetscCall(VecRestoreArray(cVec, &coordsDG));
452 PetscUseTypeMethod(dm, createcellcoordinatedm, &cdgdm);
453 PetscCall(DMSetCellCoordinateDM(dm, cdgdm));
454 PetscCall(DMDestroy(&cdgdm));
455 // Convert the discretization
456 {
457 PetscFE fe;
458 PetscClassId id;
459
460 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
461 PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
462 if (id == PETSCFE_CLASSID) {
463 PetscSpace P;
464 PetscInt degree;
465
466 PetscCall(PetscFEGetBasisSpace(fe, &P));
467 PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
468 PetscCall(DMPlexCreateCoordinateSpace(dm, degree, PETSC_TRUE, PETSC_FALSE));
469 }
470 }
471 PetscCall(DMSetCellCoordinateSection(dm, PETSC_DETERMINE, csDG));
472 PetscCall(DMSetCellCoordinatesLocal(dm, cVec));
473 PetscCall(VecDestroy(&cVec));
474
475 end:
476 PetscCall(DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor));
477 PetscCall(DMRestoreWorkArray(dm, 2 * (maxHeight + 1), MPIU_INT, &pStart));
478 PetscCall(PetscSectionDestroy(&csDG));
479 PetscCall(DMDestroy(&plex));
480 PetscCall(DMDestroy(&cplex));
481 PetscFunctionReturn(PETSC_SUCCESS);
482 }
483