Actual source code: plexinterpolate.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <../src/sys/utils/hash.h>

  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;

 17:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 18:   DMPlexGetConeSize(dm, p, &coneSize);
 19:   DMPlexGetCone(dm, p, &cone);
 20:   DMPlexGetRawFaces_Internal(dm, dim, coneSize, cone, numFaces, faceSize, faces);
 21:   return(0);
 22: }

 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;

 34:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void *) faces);
 35:   return(0);
 36: }

 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;

 51:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
 52:   if (faces) {DMGetWorkArray(dm, PetscSqr(PetscMax(maxConeSize, maxSupportSize)), PETSC_INT, &facesTmp);}
 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:   return(0);
141: }

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;

155:   DMPlexGetDimension(dm, &cellDim);
156:   /* HACK: I need a better way to determine face dimension, or an alternative to GetFaces() */
157:   DMPlexGetSubpointMap(dm, &subpointMap);
158:   if (subpointMap) ++cellDim;
159:   DMPlexGetDepth(dm, &depth);
160:   ++depth;
161:   ++cellDepth;
162:   cellDim -= depth - cellDepth;
163:   PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
164:   for (d = depth-1; d >= faceDepth; --d) {
165:     DMPlexGetDepthStratum(dm, d, &pStart[d+1], &pEnd[d+1]);
166:   }
167:   DMPlexGetDepthStratum(dm, -1, NULL, &pStart[faceDepth]);
168:   pEnd[faceDepth] = pStart[faceDepth];
169:   for (d = faceDepth-1; d >= 0; --d) {
170:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
171:   }
172:   if (pEnd[cellDepth] > pStart[cellDepth]) {DMPlexGetFaces_Internal(dm, cellDim, pStart[cellDepth], NULL, &faceSizeAll, NULL);}
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:   PetscHashIJKLCreate(&faceTable);
175:   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
176:     const PetscInt *cellFaces;
177:     PetscInt        numCellFaces, faceSize, cf;

179:     DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
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;

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:         PetscSortInt(faceSize, (PetscInt *) &key);
194:       }
195:       PetscHashIJKLPut(faceTable, key, &missing, &iter);
196:       if (missing) {PetscHashIJKLSet(faceTable, iter, face++);}
197:     }
198:     DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
199:   }
200:   pEnd[faceDepth] = face;
201:   PetscHashIJKLDestroy(&faceTable);
202:   /* Count new points */
203:   for (d = 0; d <= depth; ++d) {
204:     numPoints += pEnd[d]-pStart[d];
205:   }
206:   DMPlexSetChart(idm, 0, numPoints);
207:   /* Set cone sizes */
208:   for (d = 0; d <= depth; ++d) {
209:     PetscInt coneSize, p;

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:         DMPlexSetConeSize(idm, p, faceSizeAll);
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:         DMPlexGetFaces_Internal(dm, cellDim, p, &coneSize, NULL, NULL);
220:         DMPlexSetConeSize(idm, p, coneSize);
221:       }
222:     } else {
223:       for (p = pStart[d]; p < pEnd[d]; ++p) {
224:         DMPlexGetConeSize(dm, p, &coneSize);
225:         DMPlexSetConeSize(idm, p, coneSize);
226:       }
227:     }
228:   }
229:   DMSetUp(idm);
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:   PetscHashIJKLCreate(&faceTable);
233:   for (d = depth; d > cellDepth; --d) {
234:     const PetscInt *cone;
235:     PetscInt        p;

237:     for (p = pStart[d]; p < pEnd[d]; ++p) {
238:       DMPlexGetCone(dm, p, &cone);
239:       DMPlexSetCone(idm, p, cone);
240:       DMPlexGetConeOrientation(dm, p, &cone);
241:       DMPlexSetConeOrientation(idm, p, cone);
242:     }
243:   }
244:   for (c = pStart[cellDepth], face = pStart[faceDepth]; c < pEnd[cellDepth]; ++c) {
245:     const PetscInt *cellFaces;
246:     PetscInt        numCellFaces, faceSize, cf;

248:     DMPlexGetFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
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;

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:         PetscSortInt(faceSize, (PetscInt *) &key);
263:       }
264:       PetscHashIJKLPut(faceTable, key, &missing, &iter);
265:       if (missing) {
266:         DMPlexSetCone(idm, face, cellFace);
267:         PetscHashIJKLSet(faceTable, iter, face);
268:         DMPlexInsertCone(idm, c, cf, face++);
269:       } else {
270:         const PetscInt *cone;
271:         PetscInt        coneSize, ornt, i, j, f;

273:         PetscHashIJKLGet(faceTable, iter, &f);
274:         DMPlexInsertCone(idm, c, cf, f);
275:         /* Orient face: Do not allow reverse orientation at the first vertex */
276:         DMPlexGetConeSize(idm, f, &coneSize);
277:         DMPlexGetCone(idm, f, &cone);
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:         DMPlexInsertConeOrientation(idm, c, cf, ornt);
295:       }
296:     }
297:     DMPlexRestoreFaces_Internal(dm, cellDim, c, &numCellFaces, &faceSize, &cellFaces);
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:   PetscFree2(pStart,pEnd);
301:   PetscHashIJKLDestroy(&faceTable);
302:   PetscFree2(pStart,pEnd);
303:   DMPlexSymmetrize(idm);
304:   DMPlexStratify(idm);
305:   return(0);
306: }

310: /*@C
311:   DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc.

313:   Collective on DM

315:   Input Parameters:
316: + dm - The DMPlex object with only cells and vertices
317: - dmInt - If NULL a new DM is created, otherwise the interpolated DM is put into the given DM

319:   Output Parameter:
320: . dmInt - The complete DMPlex object

322:   Level: intermediate

324: .keywords: mesh
325: .seealso: DMPlexUninterpolate(), DMPlexCreateFromCellList()
326: @*/
327: PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt)
328: {
329:   DM             idm, odm = dm;
330:   PetscInt       depth, dim, d;

334:   PetscLogEventBegin(DMPLEX_Interpolate,dm,0,0,0);
335:   DMPlexGetDepth(dm, &depth);
336:   DMPlexGetDimension(dm, &dim);
337:   if (dim <= 1) {
338:     PetscObjectReference((PetscObject) dm);
339:     idm  = dm;
340:   }
341:   for (d = 1; d < dim; ++d) {
342:     /* Create interpolated mesh */
343:     if ((d == dim-1) && *dmInt) {idm  = *dmInt;}
344:     else                        {DMCreate(PetscObjectComm((PetscObject)dm), &idm);}
345:     DMSetType(idm, DMPLEX);
346:     DMPlexSetDimension(idm, dim);
347:     if (depth > 0) {DMPlexInterpolateFaces_Internal(odm, 1, idm);}
348:     if (odm != dm) {DMDestroy(&odm);}
349:     odm  = idm;
350:   }
351:   *dmInt = idm;
352:   PetscLogEventEnd(DMPLEX_Interpolate,dm,0,0,0);
353:   return(0);
354: }

358: /*@
359:   DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices

361:   Collective on DM

363:   Input Parameter:
364: . dmA - The DMPlex object with initial coordinates

366:   Output Parameter:
367: . dmB - The DMPlex object with copied coordinates

369:   Level: intermediate

371:   Note: This is typically used when adding pieces other than vertices to a mesh

373: .keywords: mesh
374: .seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
375: @*/
376: PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
377: {
378:   Vec            coordinatesA, coordinatesB;
379:   PetscSection   coordSectionA, coordSectionB;
380:   PetscScalar   *coordsA, *coordsB;
381:   PetscInt       spaceDim, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d;

385:   if (dmA == dmB) return(0);
386:   DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA);
387:   DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB);
388:   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);
389:   DMGetCoordinateSection(dmA, &coordSectionA);
390:   DMGetCoordinateSection(dmB, &coordSectionB);
391:   PetscSectionSetNumFields(coordSectionB, 1);
392:   PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim);
393:   PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim);
394:   PetscSectionSetChart(coordSectionB, vStartB, vEndB);
395:   for (v = vStartB; v < vEndB; ++v) {
396:     PetscSectionSetDof(coordSectionB, v, spaceDim);
397:     PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim);
398:   }
399:   PetscSectionSetUp(coordSectionB);
400:   PetscSectionGetStorageSize(coordSectionB, &coordSizeB);
401:   DMGetCoordinatesLocal(dmA, &coordinatesA);
402:   VecCreate(PetscObjectComm((PetscObject) dmB), &coordinatesB);
403:   PetscObjectSetName((PetscObject) coordinatesB, "coordinates");
404:   VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE);
405:   VecSetType(coordinatesB,VECSTANDARD);
406:   VecGetArray(coordinatesA, &coordsA);
407:   VecGetArray(coordinatesB, &coordsB);
408:   for (v = 0; v < vEndB-vStartB; ++v) {
409:     for (d = 0; d < spaceDim; ++d) {
410:       coordsB[v*spaceDim+d] = coordsA[v*spaceDim+d];
411:     }
412:   }
413:   VecRestoreArray(coordinatesA, &coordsA);
414:   VecRestoreArray(coordinatesB, &coordsB);
415:   DMSetCoordinatesLocal(dmB, coordinatesB);
416:   VecDestroy(&coordinatesB);
417:   return(0);
418: }

422: /*@
423:   DMPlexCopyLabels - Copy labels from one mesh to another with a superset of the points

425:   Collective on DM

427:   Input Parameter:
428: . dmA - The DMPlex object with initial labels

430:   Output Parameter:
431: . dmB - The DMPlex object with copied labels

433:   Level: intermediate

435:   Note: This is typically used when interpolating or otherwise adding to a mesh

437: .keywords: mesh
438: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
439: @*/
440: PetscErrorCode DMPlexCopyLabels(DM dmA, DM dmB)
441: {
442:   PetscInt       numLabels, l;

446:   if (dmA == dmB) return(0);
447:   DMPlexGetNumLabels(dmA, &numLabels);
448:   for (l = 0; l < numLabels; ++l) {
449:     DMLabel     label, labelNew;
450:     const char *name;
451:     PetscBool   flg;

453:     DMPlexGetLabelName(dmA, l, &name);
454:     PetscStrcmp(name, "depth", &flg);
455:     if (flg) continue;
456:     DMPlexGetLabel(dmA, name, &label);
457:     DMLabelDuplicate(label, &labelNew);
458:     DMPlexAddLabel(dmB, labelNew);
459:   }
460:   return(0);
461: }

465: /*@
466:   DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh

468:   Collective on DM

470:   Input Parameter:
471: . dm - The complete DMPlex object

473:   Output Parameter:
474: . dmUnint - The DMPlex object with only cells and vertices

476:   Level: intermediate

478: .keywords: mesh
479: .seealso: DMPlexInterpolate(), DMPlexCreateFromCellList()
480: @*/
481: PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint)
482: {
483:   DM             udm;
484:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone;

488:   DMPlexGetDimension(dm, &dim);
489:   if (dim <= 1) {
490:     PetscObjectReference((PetscObject) dm);
491:     *dmUnint = dm;
492:     return(0);
493:   }
494:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
495:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
496:   DMCreate(PetscObjectComm((PetscObject) dm), &udm);
497:   DMSetType(udm, DMPLEX);
498:   DMPlexSetDimension(udm, dim);
499:   DMPlexSetChart(udm, cStart, vEnd);
500:   for (c = cStart; c < cEnd; ++c) {
501:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

503:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
504:     for (cl = 0; cl < closureSize*2; cl += 2) {
505:       const PetscInt p = closure[cl];

507:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
508:     }
509:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
510:     DMPlexSetConeSize(udm, c, coneSize);
511:     maxConeSize = PetscMax(maxConeSize, coneSize);
512:   }
513:   DMSetUp(udm);
514:   PetscMalloc1(maxConeSize, &cone);
515:   for (c = cStart; c < cEnd; ++c) {
516:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

518:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
519:     for (cl = 0; cl < closureSize*2; cl += 2) {
520:       const PetscInt p = closure[cl];

522:       if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p;
523:     }
524:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
525:     DMPlexSetCone(udm, c, cone);
526:   }
527:   PetscFree(cone);
528:   DMPlexSymmetrize(udm);
529:   DMPlexStratify(udm);
530:   *dmUnint = udm;
531:   return(0);
532: }