xref: /petsc/src/dm/impls/plex/plexinterpolate.c (revision 89d949e22d93de971d3ed7d7fdc00f23508da20b) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <../src/sys/utils/hash.h>
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMPlexGetFaces_Internal"
6 /*
7   DMPlexGetFaces_Internal - Gets groups of vertices that correspond to faces for the given cell
8 */
9 PetscErrorCode DMPlexGetFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
10 {
11   const PetscInt *cone = NULL;
12   PetscInt        maxConeSize, maxSupportSize, coneSize;
13   PetscErrorCode  ierr;
14 
15   PetscFunctionBegin;
16   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
17   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
18   ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
19   ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
20   ierr = DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);CHKERRQ(ierr);
21   PetscFunctionReturn(0);
22 }
23 
24 #undef __FUNCT__
25 #define __FUNCT__ "DMPlexRestoreFaces_Internal"
26 /*
27   DMPlexRestoreFaces_Internal - Restores the array
28 */
29 PetscErrorCode DMPlexRestoreFaces_Internal(DM dm, PetscInt dim, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
30 {
31   PetscErrorCode  ierr;
32 
33   PetscFunctionBegin;
34   ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, (void *) faces);CHKERRQ(ierr);
35   PetscFunctionReturn(0);
36 }
37 
38 #undef __FUNCT__
39 #define __FUNCT__ "DMPlexGetRawFaces_Internal"
40 /*
41   DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone
42 */
43 PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, PetscInt dim, PetscInt coneSize, const PetscInt cone[], PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
44 {
45   PetscInt       *facesTmp;
46   PetscInt        maxConeSize, maxSupportSize;
47   PetscErrorCode  ierr;
48 
49   PetscFunctionBegin;
50   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
51   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
52   if (faces) {ierr = DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);CHKERRQ(ierr);}
53   switch (dim) {
54   case 1:
55     switch (coneSize) {
56     case 2:
57       if (faces) {
58         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
59         *faces = facesTmp;
60       }
61       if (numFaces) *numFaces         = 2;
62       if (faceSize) *faceSize         = 1;
63       break;
64     default:
65       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
66     }
67     break;
68   case 2:
69     switch (coneSize) {
70     case 3:
71       if (faces) {
72         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
73         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
74         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
75         *faces = facesTmp;
76       }
77       if (numFaces) *numFaces         = 3;
78       if (faceSize) *faceSize         = 2;
79       break;
80     case 4:
81       /* Vertices follow right hand rule */
82       if (faces) {
83         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
84         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
85         facesTmp[4] = cone[2]; facesTmp[5] = cone[3];
86         facesTmp[6] = cone[3]; facesTmp[7] = cone[0];
87         *faces = facesTmp;
88       }
89       if (numFaces) *numFaces         = 4;
90       if (faceSize) *faceSize         = 2;
91       break;
92     default:
93       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
94     }
95     break;
96   case 3:
97     switch (coneSize) {
98     case 3:
99       if (faces) {
100         facesTmp[0] = cone[0]; facesTmp[1] = cone[1];
101         facesTmp[2] = cone[1]; facesTmp[3] = cone[2];
102         facesTmp[4] = cone[2]; facesTmp[5] = cone[0];
103         *faces = facesTmp;
104       }
105       if (numFaces) *numFaces         = 3;
106       if (faceSize) *faceSize         = 2;
107       break;
108     case 4:
109       /* Vertices of first face follow right hand rule and normal points away from last vertex */
110       if (faces) {
111         facesTmp[0] = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2];
112         facesTmp[3] = cone[0]; facesTmp[4]  = cone[3]; facesTmp[5]  = cone[1];
113         facesTmp[6] = cone[0]; facesTmp[7]  = cone[2]; facesTmp[8]  = cone[3];
114         facesTmp[9] = cone[2]; facesTmp[10] = cone[1]; facesTmp[11] = cone[3];
115         *faces = facesTmp;
116       }
117       if (numFaces) *numFaces         = 4;
118       if (faceSize) *faceSize         = 3;
119       break;
120     case 8:
121       if (faces) {
122         facesTmp[0]  = cone[0]; facesTmp[1]  = cone[1]; facesTmp[2]  = cone[2]; facesTmp[3]  = cone[3]; /* Bottom */
123         facesTmp[4]  = cone[4]; facesTmp[5]  = cone[5]; facesTmp[6]  = cone[6]; facesTmp[7]  = cone[7]; /* Top */
124         facesTmp[8]  = cone[0]; facesTmp[9]  = cone[3]; facesTmp[10] = cone[5]; facesTmp[11] = cone[4]; /* Front */
125         facesTmp[12] = cone[2]; facesTmp[13] = cone[1]; facesTmp[14] = cone[7]; facesTmp[15] = cone[6]; /* Back */
126         facesTmp[16] = cone[3]; facesTmp[17] = cone[2]; facesTmp[18] = cone[6]; facesTmp[19] = cone[5]; /* Right */
127         facesTmp[20] = cone[0]; facesTmp[21] = cone[4]; facesTmp[22] = cone[7]; facesTmp[23] = cone[1]; /* Left */
128         *faces = facesTmp;
129       }
130       if (numFaces) *numFaces         = 6;
131       if (faceSize) *faceSize         = 4;
132       break;
133     default:
134       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %D", coneSize, dim);
135     }
136     break;
137   default:
138     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
139   }
140   PetscFunctionReturn(0);
141 }
142 
143 #undef __FUNCT__
144 #define __FUNCT__ "DMPlexInterpolateFaces_Internal"
145 /* This interpolates faces for cells at some stratum */
146 static PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm)
147 {
148   DMLabel        subpointMap;
149   PetscHashIJKL  faceTable;
150   PetscInt      *pStart, *pEnd;
151   PetscInt       cellDim, depth, faceDepth = cellDepth, numPoints = 0, faceSizeAll = 0, face, c, d;
152   PetscErrorCode ierr;
153 
154   PetscFunctionBegin;
155   ierr = DMPlexGetDimension(dm, &cellDim);CHKERRQ(ierr);
156   /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
157   ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr);
158   if (subpointMap) ++cellDim;
159   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
160   ++depth;
161   ++cellDepth;
162   cellDim -= depth - cellDepth;
163   ierr = PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);CHKERRQ(ierr);
164   for (d = depth-1; d >= faceDepth; --d) {
165     ierr = DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);CHKERRQ(ierr);
166   }
167   ierr = DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);CHKERRQ(ierr);
168   pEnd[faceDepth] = pStart[faceDepth];
169   for (d = faceDepth-1; d >= 0; --d) {
170     ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);
171   }
172   if (pEnd[cellDepth] > pStart[cellDepth]) {ierr = DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);CHKERRQ(ierr);}
173   if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
174   ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr);
175   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
176     const PetscInt *cellFaces;
177     PetscInt        numCellFaces, faceSize, cf;
178 
179     ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr);
180     if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
181     for (cf = 0; cf < numCellFaces; ++cf) {
182       const PetscInt   *cellFace = &cellFaces[cf*faceSize];
183       PetscHashIJKLKey  key;
184       PetscHashIJKLIter missing, iter;
185 
186       if (faceSize == 2) {
187         key.i = PetscMin(cellFace[0], cellFace[1]);
188         key.j = PetscMax(cellFace[0], cellFace[1]);
189         key.k = 0;
190         key.l = 0;
191       } else {
192         key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
193         ierr = PetscSortInt(faceSize, (PetscInt *) &key);
194       }
195       ierr = PetscHashIJKLPut(faceTable, key, &missing, &iter);CHKERRQ(ierr);
196       if (missing) {ierr = PetscHashIJKLSet(faceTable, iter, face++);CHKERRQ(ierr);}
197     }
198     ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr);
199   }
200   pEnd[faceDepth] = face;
201   ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr);
202   /* Count new points */
203   for (d = 0; d <= depth; ++d) {
204     numPoints += pEnd[d]-pStart[d];
205   }
206   ierr = DMPlexSetChart(idm, 0, numPoints);CHKERRQ(ierr);
207   /* Set cone sizes */
208   for (d = 0; d <= depth; ++d) {
209     PetscInt coneSize, p;
210 
211     if (d == faceDepth) {
212       for (p = pStart[d]; p < pEnd[d]; ++p) {
213         /* I see no way to do this if we admit faces of different shapes */
214         ierr = DMPlexSetConeSize(idm, p, faceSizeAll);CHKERRQ(ierr);
215       }
216     } else if (d == cellDepth) {
217       for (p = pStart[d]; p < pEnd[d]; ++p) {
218         /* Number of cell faces may be different from number of cell vertices*/
219         ierr = DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);CHKERRQ(ierr);
220         ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr);
221       }
222     } else {
223       for (p = pStart[d]; p < pEnd[d]; ++p) {
224         ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
225         ierr = DMPlexSetConeSize(idm, p, coneSize);CHKERRQ(ierr);
226       }
227     }
228   }
229   ierr = DMSetUp(idm);CHKERRQ(ierr);
230   /* Get face cones from subsets of cell vertices */
231   if (faceSizeAll > 4) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Do not support interpolation of meshes with faces of %D vertices", faceSizeAll);
232   ierr = PetscHashIJKLCreate(&faceTable);CHKERRQ(ierr);
233   for (d = depth; d > cellDepth; --d) {
234     const PetscInt *cone;
235     PetscInt        p;
236 
237     for (p = pStart[d]; p < pEnd[d]; ++p) {
238       ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
239       ierr = DMPlexSetCone(idm, p, cone);CHKERRQ(ierr);
240       ierr = DMPlexGetConeOrientation(dm, p, &cone);CHKERRQ(ierr);
241       ierr = DMPlexSetConeOrientation(idm, p, cone);CHKERRQ(ierr);
242     }
243   }
244   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
245     const PetscInt *cellFaces;
246     PetscInt        numCellFaces, faceSize, cf;
247 
248     ierr = DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr);
249     if (faceSize != faceSizeAll) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent face for cell %D of size %D != %D", c, faceSize, faceSizeAll);
250     for (cf = 0; cf < numCellFaces; ++cf) {
251       const PetscInt  *cellFace = &cellFaces[cf*faceSize];
252       PetscHashIJKLKey key;
253       PetscHashIJKLIter missing, iter;
254 
255       if (faceSize == 2) {
256         key.i = PetscMin(cellFace[0], cellFace[1]);
257         key.j = PetscMax(cellFace[0], cellFace[1]);
258         key.k = 0;
259         key.l = 0;
260       } else {
261         key.i = cellFace[0]; key.j = cellFace[1]; key.k = cellFace[2]; key.l = faceSize > 3 ? cellFace[3] : 0;
262         ierr = PetscSortInt(faceSize, (PetscInt *) &key);
263       }
264       ierr = PetscHashIJKLPut(faceTable, key, &missing, &iter);CHKERRQ(ierr);
265       if (missing) {
266         ierr = DMPlexSetCone(idm, face, cellFace);CHKERRQ(ierr);
267         ierr = PetscHashIJKLSet(faceTable, iter, face);CHKERRQ(ierr);
268         ierr = DMPlexInsertCone(idm, c, cf, face++);CHKERRQ(ierr);
269       } else {
270         const PetscInt *cone;
271         PetscInt        coneSize, ornt, i, j, f;
272 
273         ierr = PetscHashIJKLGet(faceTable, iter, &f);CHKERRQ(ierr);
274         ierr = DMPlexInsertCone(idm, c, cf, f);CHKERRQ(ierr);
275         /* Orient face: Do not allow reverse orientation at the first vertex */
276         ierr = DMPlexGetConeSize(idm, f, &coneSize);CHKERRQ(ierr);
277         ierr = DMPlexGetCone(idm, f, &cone);CHKERRQ(ierr);
278         if (coneSize != faceSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %D for face %D should be %D", coneSize, f, faceSize);
279         /* - First find the initial vertex */
280         for (i = 0; i < faceSize; ++i) if (cellFace[0] == cone[i]) break;
281         /* - Try forward comparison */
282         for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+j)%faceSize]) break;
283         if (j == faceSize) {
284           if ((faceSize == 2) && (i == 1)) ornt = -2;
285           else                             ornt = i;
286         } else {
287           /* - Try backward comparison */
288           for (j = 0; j < faceSize; ++j) if (cellFace[j] != cone[(i+faceSize-j)%faceSize]) break;
289           if (j == faceSize) {
290             if (i == 0) ornt = -faceSize;
291             else        ornt = -i;
292           } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine face orientation");
293         }
294         ierr = DMPlexInsertConeOrientation(idm, c, cf, ornt);CHKERRQ(ierr);
295       }
296     }
297     ierr = DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);CHKERRQ(ierr);
298   }
299   if (face != pEnd[faceDepth]) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid number of faces %D should be %D", face-pStart[faceDepth], pEnd[faceDepth]-pStart[faceDepth]);
300   ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr);
301   ierr = PetscHashIJKLDestroy(&faceTable);CHKERRQ(ierr);
302   ierr = PetscFree2(pStart,pEnd);CHKERRQ(ierr);
303   ierr = DMPlexSymmetrize(idm);CHKERRQ(ierr);
304   ierr = DMPlexStratify(idm);CHKERRQ(ierr);
305   PetscFunctionReturn(0);
306 }
307 
308 #undef __FUNCT__
309 #define __FUNCT__ "DMPlexInterpolate"
310 /*@
311   DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.
312 
313   Collective on DM
314 
315   Input Parameter:
316 . dm - The DMPlex object with only cells and vertices
317 
318   Output Parameter:
319 . dmInt - The complete DMPlex object
320 
321   Level: intermediate
322 
323 .keywords: mesh
324 .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList()
325 @*/
326 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
327 {
328   DM             idm, odm = dm;
329   PetscInt       depth, dim, d;
330   PetscErrorCode ierr;
331 
332   PetscFunctionBegin;
333   ierr = PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr);
334   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
335   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
336   if (dim <= 1) {
337     ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
338     idm  = dm;
339   }
340   for (d = 1; d < dim; ++d) {
341     /* Create interpolated mesh */
342     ierr = DMCreate(PetscObjectComm((PetscObject)dm), &idm);CHKERRQ(ierr);
343     ierr = DMSetType(idm, DMPLEX);CHKERRQ(ierr);
344     ierr = DMPlexSetDimension(idm, dim);CHKERRQ(ierr);
345     if (depth > 0) {ierr = DMPlexInterpolateFaces_Internal(odm, 1, idm);CHKERRQ(ierr);}
346     if (odm != dm) {ierr = DMDestroy(&odm);CHKERRQ(ierr);}
347     odm  = idm;
348   }
349   *dmInt = idm;
350   ierr = PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);CHKERRQ(ierr);
351   PetscFunctionReturn(0);
352 }
353 
354 #undef __FUNCT__
355 #define __FUNCT__ "DMPlexCopyCoordinates"
356 /*@
357   DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices
358 
359   Collective on DM
360 
361   Input Parameter:
362 . dmA - The DMPlex object with initial coordinates
363 
364   Output Parameter:
365 . dmB - The DMPlex object with copied coordinates
366 
367   Level: intermediate
368 
369   Note: This is typically used when adding pieces other than vertices to a mesh
370 
371 .keywords: mesh
372 .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
373 @*/
374 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
375 {
376   Vec            coordinatesA, coordinatesB;
377   PetscSection   coordSectionA, coordSectionB;
378   PetscScalar   *coordsA, *coordsB;
379   PetscInt       spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;
380   PetscErrorCode ierr;
381 
382   PetscFunctionBegin;
383   if (dmA == dmB) PetscFunctionReturn(0);
384   ierr = DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);CHKERRQ(ierr);
385   ierr = DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);CHKERRQ(ierr);
386   if ((vEndA-vStartA) != (vEndB-vStartB)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %d != %d in the second DM", vEndA-vStartA, vEndB-vStartB);
387   ierr = DMGetCoordinateSection(dmA, &coordSectionA);CHKERRQ(ierr);
388   ierr = DMGetCoordinateSection(dmB, &coordSectionB);CHKERRQ(ierr);
389   ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr);
390   ierr = PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);CHKERRQ(ierr);
391   ierr = PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);CHKERRQ(ierr);
392   ierr = PetscSectionSetChart(coordSectionB, vStartB, vEndB);CHKERRQ(ierr);
393   for (v = vStartB; v < vEndB; ++v) {
394     ierr = PetscSectionSetDof(coordSectionB, v, spaceDim);CHKERRQ(ierr);
395     ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);CHKERRQ(ierr);
396   }
397   ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr);
398   ierr = PetscSectionGetStorageSize(coordSectionB, &coordSizeB);CHKERRQ(ierr);
399   ierr = DMGetCoordinatesLocal(dmA, &coordinatesA);CHKERRQ(ierr);
400   ierr = VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);CHKERRQ(ierr);
401   ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr);
402   ierr = VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);CHKERRQ(ierr);
403   ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr);
404   ierr = VecGetArray(coordinatesA, &coordsA);CHKERRQ(ierr);
405   ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr);
406   for (v = 0; v < vEndB-vStartB; ++v) {
407     for (d = 0; d < spaceDim; ++d) {
408       coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
409     }
410   }
411   ierr = VecRestoreArray(coordinatesA, &coordsA);CHKERRQ(ierr);
412   ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr);
413   ierr = DMSetCoordinatesLocal(dmB, coordinatesB);CHKERRQ(ierr);
414   ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr);
415   PetscFunctionReturn(0);
416 }
417 
418 #undef __FUNCT__
419 #define __FUNCT__ "DMPlexCopyLabels"
420 /*@
421   DMPlexCopyLabels - Copy labels from one mesh to another with a superset of the points
422 
423   Collective on DM
424 
425   Input Parameter:
426 . dmA - The DMPlex object with initial labels
427 
428   Output Parameter:
429 . dmB - The DMPlex object with copied labels
430 
431   Level: intermediate
432 
433   Note: This is typically used when interpolating or otherwise adding to a mesh
434 
435 .keywords: mesh
436 .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
437 @*/
438 PetscErrorCode DMPlexCopyLabels(DM dmA, DM dmB)
439 {
440   PetscInt       numLabels, l;
441   PetscErrorCode ierr;
442 
443   PetscFunctionBegin;
444   if (dmA == dmB) PetscFunctionReturn(0);
445   ierr = DMPlexGetNumLabels(dmA, &numLabels);CHKERRQ(ierr);
446   for (l = 0; l < numLabels; ++l) {
447     DMLabel     label, labelNew;
448     const char *name;
449     PetscBool   flg;
450 
451     ierr = DMPlexGetLabelName(dmA, l, &name);CHKERRQ(ierr);
452     ierr = PetscStrcmp(name, "depth", &flg);CHKERRQ(ierr);
453     if (flg) continue;
454     ierr = DMPlexGetLabel(dmA, name, &label);CHKERRQ(ierr);
455     ierr = DMLabelDuplicate(label, &labelNew);CHKERRQ(ierr);
456     ierr = DMPlexAddLabel(dmB, labelNew);CHKERRQ(ierr);
457   }
458   PetscFunctionReturn(0);
459 }
460 
461 #undef __FUNCT__
462 #define __FUNCT__ "DMPlexUninterpolate"
463 /*@
464   DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh
465 
466   Collective on DM
467 
468   Input Parameter:
469 . dm - The complete DMPlex object
470 
471   Output Parameter:
472 . dmUnint - The DMPlex object with only cells and vertices
473 
474   Level: intermediate
475 
476 .keywords: mesh
477 .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList()
478 @*/
479 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
480 {
481   DM             udm;
482   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;
483   PetscErrorCode ierr;
484 
485   PetscFunctionBegin;
486   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
487   if (dim <= 1) {
488     ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
489     *dmUnint = dm;
490     PetscFunctionReturn(0);
491   }
492   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
493   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
494   ierr = DMCreate(PetscObjectComm((PetscObject) dm), &udm);CHKERRQ(ierr);
495   ierr = DMSetType(udm, DMPLEX);CHKERRQ(ierr);
496   ierr = DMPlexSetDimension(udm, dim);CHKERRQ(ierr);
497   ierr = DMPlexSetChart(udm, cStart, vEnd);CHKERRQ(ierr);
498   for (c = cStart; c < cEnd; ++c) {
499     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
500 
501     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
502     for (cl = 0; cl < closureSize*2; cl += 2) {
503       const PetscInt p = closure[cl];
504 
505       if ((p >= vStart) && (p < vEnd)) ++coneSize;
506     }
507     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
508     ierr = DMPlexSetConeSize(udm, c, coneSize);CHKERRQ(ierr);
509     maxConeSize = PetscMax(maxConeSize, coneSize);
510   }
511   ierr = DMSetUp(udm);CHKERRQ(ierr);
512   ierr = PetscMalloc1(maxConeSize, &cone);CHKERRQ(ierr);
513   for (c = cStart; c < cEnd; ++c) {
514     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
515 
516     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
517     for (cl = 0; cl < closureSize*2; cl += 2) {
518       const PetscInt p = closure[cl];
519 
520       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
521     }
522     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
523     ierr = DMPlexSetCone(udm, c, cone);CHKERRQ(ierr);
524   }
525   ierr = PetscFree(cone);CHKERRQ(ierr);
526   ierr = DMPlexSymmetrize(udm);CHKERRQ(ierr);
527   ierr = DMPlexStratify(udm);CHKERRQ(ierr);
528   *dmUnint = udm;
529   PetscFunctionReturn(0);
530 }
531