Actual source code: plexsubmesh.c

petsc-dev 2014-04-15
Report Typos and Errors
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petscsf.h>

  6: /*@
  7:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

  9:   Not Collective

 11:   Input Parameter:
 12: . dm - The original DM

 14:   Output Parameter:
 15: . label - The DMLabel marking boundary faces with value 1

 17:   Level: developer

 19: .seealso: DMLabelCreate(), DMPlexCreateLabel()
 20: @*/
 21: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
 22: {
 23:   PetscInt       fStart, fEnd, f;

 28:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
 29:   for (f = fStart; f < fEnd; ++f) {
 30:     PetscInt supportSize;

 32:     DMPlexGetSupportSize(dm, f, &supportSize);
 33:     if (supportSize == 1) {
 34:       DMLabelSetValue(label, f, 1);
 35:     }
 36:   }
 37:   return(0);
 38: }

 42: /*@
 43:   DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface

 45:   Input Parameters:
 46: + dm - The DM
 47: - label - A DMLabel marking the surface points

 49:   Output Parameter:
 50: . label - A DMLabel marking all surface points in the transitive closure

 52:   Level: developer

 54: .seealso: DMPlexLabelCohesiveComplete()
 55: @*/
 56: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
 57: {
 58:   IS              valueIS;
 59:   const PetscInt *values;
 60:   PetscInt        numValues, v;
 61:   PetscErrorCode  ierr;

 64:   DMLabelGetNumValues(label, &numValues);
 65:   DMLabelGetValueIS(label, &valueIS);
 66:   ISGetIndices(valueIS, &values);
 67:   for (v = 0; v < numValues; ++v) {
 68:     IS              pointIS;
 69:     const PetscInt *points;
 70:     PetscInt        numPoints, p;

 72:     DMLabelGetStratumSize(label, values[v], &numPoints);
 73:     DMLabelGetStratumIS(label, values[v], &pointIS);
 74:     ISGetIndices(pointIS, &points);
 75:     for (p = 0; p < numPoints; ++p) {
 76:       PetscInt *closure = NULL;
 77:       PetscInt  closureSize, c;

 79:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 80:       for (c = 0; c < closureSize*2; c += 2) {
 81:         DMLabelSetValue(label, closure[c], values[v]);
 82:       }
 83:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 84:     }
 85:     ISRestoreIndices(pointIS, &points);
 86:     ISDestroy(&pointIS);
 87:   }
 88:   ISRestoreIndices(valueIS, &values);
 89:   ISDestroy(&valueIS);
 90:   return(0);
 91: }

 95: /*@
 96:   DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face

 98:   Input Parameters:
 99: + dm - The DM
100: - label - A DMLabel marking the surface points

102:   Output Parameter:
103: . label - A DMLabel incorporating cells

105:   Level: developer

107:   Note: The cells allow FEM boundary conditions to be applied using the cell geometry

109: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
110: @*/
111: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
112: {
113:   IS              valueIS;
114:   const PetscInt *values;
115:   PetscInt        numValues, v, cStart, cEnd;
116:   PetscErrorCode  ierr;

119:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
120:   DMLabelGetNumValues(label, &numValues);
121:   DMLabelGetValueIS(label, &valueIS);
122:   ISGetIndices(valueIS, &values);
123:   for (v = 0; v < numValues; ++v) {
124:     IS              pointIS;
125:     const PetscInt *points;
126:     PetscInt        numPoints, p;

128:     DMLabelGetStratumSize(label, values[v], &numPoints);
129:     DMLabelGetStratumIS(label, values[v], &pointIS);
130:     ISGetIndices(pointIS, &points);
131:     for (p = 0; p < numPoints; ++p) {
132:       PetscInt *closure = NULL;
133:       PetscInt  closureSize, point;

135:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
136:       point = closure[(closureSize-1)*2];
137:       if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]);}
138:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
139:     }
140:     ISRestoreIndices(pointIS, &points);
141:     ISDestroy(&pointIS);
142:   }
143:   ISRestoreIndices(valueIS, &values);
144:   ISDestroy(&valueIS);
145:   return(0);
146: }

150: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthMax[], PetscInt depthEnd[], PetscInt depthShift[])
151: {
152:   if (depth < 0) return p;
153:   /* Normal Cells                 */ if (p < depthMax[depth])                return p;
154:   /* Hybrid Cells+Normal Vertices */ if (p < depthMax[0])                    return p + depthShift[depth];
155:   /* Hybrid Vertices+Normal Faces */ if (depth < 2 || p < depthMax[depth-1]) return p + depthShift[depth] + depthShift[0];
156:   /* Hybrid Faces+Normal Edges    */ if (depth < 3 || p < depthMax[depth-2]) return p + depthShift[depth] + depthShift[0] + depthShift[depth-1];
157:   /* Hybrid Edges                 */                                         return p + depthShift[depth] + depthShift[0] + depthShift[depth-1] + depthShift[depth-2];
158: }

162: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
163: {
164:   PetscInt      *depthMax, *depthEnd;
165:   PetscInt       depth = 0, d, pStart, pEnd, p;

169:   DMPlexGetDepth(dm, &depth);
170:   if (depth < 0) return(0);
171:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
172:   /* Step 1: Expand chart */
173:   DMPlexGetChart(dm, &pStart, &pEnd);
174:   DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
175:   for (d = 0; d <= depth; ++d) {
176:     pEnd += depthShift[d];
177:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
178:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
179:   }
180:   DMPlexSetChart(dmNew, pStart, pEnd);
181:   /* Step 2: Set cone and support sizes */
182:   for (d = 0; d <= depth; ++d) {
183:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
184:     for (p = pStart; p < pEnd; ++p) {
185:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);
186:       PetscInt size;

188:       DMPlexGetConeSize(dm, p, &size);
189:       DMPlexSetConeSize(dmNew, newp, size);
190:       DMPlexGetSupportSize(dm, p, &size);
191:       DMPlexSetSupportSize(dmNew, newp, size);
192:     }
193:   }
194:   PetscFree2(depthMax,depthEnd);
195:   return(0);
196: }

200: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
201: {
202:   PetscInt      *depthEnd, *depthMax, *newpoints;
203:   PetscInt       depth = 0, d, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

207:   DMPlexGetDepth(dm, &depth);
208:   if (depth < 0) return(0);
209:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
210:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
211:   PetscMalloc3(depth+1,&depthEnd,depth+1,&depthMax,PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
212:   DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
213:   for (d = 0; d <= depth; ++d) {
214:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
215:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
216:   }
217:   /* Step 5: Set cones and supports */
218:   DMPlexGetChart(dm, &pStart, &pEnd);
219:   for (p = pStart; p < pEnd; ++p) {
220:     const PetscInt *points = NULL, *orientations = NULL;
221:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);

223:     DMPlexGetConeSize(dm, p, &size);
224:     DMPlexGetCone(dm, p, &points);
225:     DMPlexGetConeOrientation(dm, p, &orientations);
226:     for (i = 0; i < size; ++i) {
227:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
228:     }
229:     DMPlexSetCone(dmNew, newp, newpoints);
230:     DMPlexSetConeOrientation(dmNew, newp, orientations);
231:     DMPlexGetSupportSize(dm, p, &size);
232:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
233:     DMPlexGetSupport(dm, p, &points);
234:     for (i = 0; i < size; ++i) {
235:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
236:     }
237:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
238:     DMPlexSetSupport(dmNew, newp, newpoints);
239:   }
240:   PetscFree3(depthEnd,depthMax,newpoints);
241:   return(0);
242: }

246: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
247: {
248:   PetscSection   coordSection, newCoordSection;
249:   Vec            coordinates, newCoordinates;
250:   PetscScalar   *coords, *newCoords;
251:   PetscInt      *depthEnd, coordSize;
252:   PetscInt       dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v;

256:   DMPlexGetDimension(dm, &dim);
257:   DMPlexGetDepth(dm, &depth);
258:   PetscMalloc1((depth+1), &depthEnd);
259:   for (d = 0; d <= depth; ++d) {
260:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
261:   }
262:   /* Step 8: Convert coordinates */
263:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
264:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
265:   DMGetCoordinateSection(dm, &coordSection);
266:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
267:   PetscSectionSetNumFields(newCoordSection, 1);
268:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
269:   PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);
270:   for (v = vStartNew; v < vEndNew; ++v) {
271:     PetscSectionSetDof(newCoordSection, v, dim);
272:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
273:   }
274:   PetscSectionSetUp(newCoordSection);
275:   DMSetCoordinateSection(dmNew, newCoordSection);
276:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
277:   VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);
278:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
279:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
280:   VecSetType(newCoordinates,VECSTANDARD);
281:   DMSetCoordinatesLocal(dmNew, newCoordinates);
282:   DMGetCoordinatesLocal(dm, &coordinates);
283:   VecGetArray(coordinates, &coords);
284:   VecGetArray(newCoordinates, &newCoords);
285:   for (v = vStart; v < vEnd; ++v) {
286:     PetscInt dof, off, noff, d;

288:     PetscSectionGetDof(coordSection, v, &dof);
289:     PetscSectionGetOffset(coordSection, v, &off);
290:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthEnd, depthEnd, depthShift), &noff);
291:     for (d = 0; d < dof; ++d) {
292:       newCoords[noff+d] = coords[off+d];
293:     }
294:   }
295:   VecRestoreArray(coordinates, &coords);
296:   VecRestoreArray(newCoordinates, &newCoords);
297:   VecDestroy(&newCoordinates);
298:   PetscSectionDestroy(&newCoordSection);
299:   PetscFree(depthEnd);
300:   return(0);
301: }

305: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
306: {
307:   PetscInt          *depthMax, *depthEnd;
308:   PetscInt           depth = 0, d;
309:   PetscSF            sfPoint, sfPointNew;
310:   const PetscSFNode *remotePoints;
311:   PetscSFNode       *gremotePoints;
312:   const PetscInt    *localPoints;
313:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
314:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
315:   PetscErrorCode     ierr;

318:   DMPlexGetDepth(dm, &depth);
319:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
320:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
321:   for (d = 0; d <= depth; ++d) {
322:     totShift += depthShift[d];
323:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
324:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
325:   }
326:   /* Step 9: Convert pointSF */
327:   DMGetPointSF(dm, &sfPoint);
328:   DMGetPointSF(dmNew, &sfPointNew);
329:   DMPlexGetChart(dm, &pStart, &pEnd);
330:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
331:   if (numRoots >= 0) {
332:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
333:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthMax, depthEnd, depthShift);
334:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
335:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
336:     PetscMalloc1(numLeaves,    &glocalPoints);
337:     PetscMalloc1(numLeaves, &gremotePoints);
338:     for (l = 0; l < numLeaves; ++l) {
339:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthMax, depthEnd, depthShift);
340:       gremotePoints[l].rank  = remotePoints[l].rank;
341:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
342:     }
343:     PetscFree2(newLocation,newRemoteLocation);
344:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
345:   }
346:   PetscFree2(depthMax,depthEnd);
347:   return(0);
348: }

352: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
353: {
354:   PetscSF            sfPoint;
355:   DMLabel            vtkLabel, ghostLabel;
356:   PetscInt          *depthMax, *depthEnd;
357:   const PetscSFNode *leafRemote;
358:   const PetscInt    *leafLocal;
359:   PetscInt           depth = 0, d, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
360:   PetscMPIInt        rank;
361:   PetscErrorCode     ierr;

364:   DMPlexGetDepth(dm, &depth);
365:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
366:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
367:   for (d = 0; d <= depth; ++d) {
368:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
369:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
370:   }
371:   /* Step 10: Convert labels */
372:   DMPlexGetNumLabels(dm, &numLabels);
373:   for (l = 0; l < numLabels; ++l) {
374:     DMLabel         label, newlabel;
375:     const char     *lname;
376:     PetscBool       isDepth;
377:     IS              valueIS;
378:     const PetscInt *values;
379:     PetscInt        numValues, val;

381:     DMPlexGetLabelName(dm, l, &lname);
382:     PetscStrcmp(lname, "depth", &isDepth);
383:     if (isDepth) continue;
384:     DMPlexCreateLabel(dmNew, lname);
385:     DMPlexGetLabel(dm, lname, &label);
386:     DMPlexGetLabel(dmNew, lname, &newlabel);
387:     DMLabelGetValueIS(label, &valueIS);
388:     ISGetLocalSize(valueIS, &numValues);
389:     ISGetIndices(valueIS, &values);
390:     for (val = 0; val < numValues; ++val) {
391:       IS              pointIS;
392:       const PetscInt *points;
393:       PetscInt        numPoints, p;

395:       DMLabelGetStratumIS(label, values[val], &pointIS);
396:       ISGetLocalSize(pointIS, &numPoints);
397:       ISGetIndices(pointIS, &points);
398:       for (p = 0; p < numPoints; ++p) {
399:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthMax, depthEnd, depthShift);

401:         DMLabelSetValue(newlabel, newpoint, values[val]);
402:       }
403:       ISRestoreIndices(pointIS, &points);
404:       ISDestroy(&pointIS);
405:     }
406:     ISRestoreIndices(valueIS, &values);
407:     ISDestroy(&valueIS);
408:   }
409:   PetscFree2(depthMax,depthEnd);
410:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
411:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
412:   DMGetPointSF(dm, &sfPoint);
413:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
414:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
415:   DMPlexCreateLabel(dmNew, "vtk");
416:   DMPlexCreateLabel(dmNew, "ghost");
417:   DMPlexGetLabel(dmNew, "vtk", &vtkLabel);
418:   DMPlexGetLabel(dmNew, "ghost", &ghostLabel);
419:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
420:     for (; c < leafLocal[l] && c < cEnd; ++c) {
421:       DMLabelSetValue(vtkLabel, c, 1);
422:     }
423:     if (leafLocal[l] >= cEnd) break;
424:     if (leafRemote[l].rank == rank) {
425:       DMLabelSetValue(vtkLabel, c, 1);
426:     } else {
427:       DMLabelSetValue(ghostLabel, c, 2);
428:     }
429:   }
430:   for (; c < cEnd; ++c) {
431:     DMLabelSetValue(vtkLabel, c, 1);
432:   }
433:   if (0) {
434:     PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
435:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
436:     PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
437:   }
438:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
439:   for (f = fStart; f < fEnd; ++f) {
440:     PetscInt numCells;

442:     DMPlexGetSupportSize(dmNew, f, &numCells);
443:     if (numCells < 2) {
444:       DMLabelSetValue(ghostLabel, f, 1);
445:     } else {
446:       const PetscInt *cells = NULL;
447:       PetscInt        vA, vB;

449:       DMPlexGetSupport(dmNew, f, &cells);
450:       DMLabelGetValue(vtkLabel, cells[0], &vA);
451:       DMLabelGetValue(vtkLabel, cells[1], &vB);
452:       if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
453:     }
454:   }
455:   if (0) {
456:     PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
457:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
458:     PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
459:   }
460:   return(0);
461: }

465: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
466: {
467:   IS              valueIS;
468:   const PetscInt *values;
469:   PetscInt       *depthShift;
470:   PetscInt        depth = 0, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
471:   PetscErrorCode  ierr;

474:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
475:   /* Count ghost cells */
476:   DMLabelGetValueIS(label, &valueIS);
477:   ISGetLocalSize(valueIS, &numFS);
478:   ISGetIndices(valueIS, &values);
479:   *numGhostCells = 0;
480:   for (fs = 0; fs < numFS; ++fs) {
481:     IS              faceIS;
482:     const PetscInt *faces;
483:     PetscInt        numFaces, f, numBdFaces = 0;

485:     DMLabelGetStratumIS(label, values[fs], &faceIS);
486:     ISGetLocalSize(faceIS, &numFaces);
487:     ISGetIndices(faceIS, &faces);
488:     for (f = 0; f < numFaces; ++f) {
489:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
490:     }
491:     *numGhostCells += numBdFaces;
492:     ISDestroy(&faceIS);
493:   }
494:   DMPlexGetDepth(dm, &depth);
495:   PetscMalloc1((depth+1), &depthShift);
496:   PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
497:   if (depth >= 0) depthShift[depth] = *numGhostCells;
498:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
499:   /* Step 3: Set cone/support sizes for new points */
500:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
501:   for (c = cEnd; c < cEnd + *numGhostCells; ++c) {
502:     DMPlexSetConeSize(gdm, c, 1);
503:   }
504:   for (fs = 0; fs < numFS; ++fs) {
505:     IS              faceIS;
506:     const PetscInt *faces;
507:     PetscInt        numFaces, f;

509:     DMLabelGetStratumIS(label, values[fs], &faceIS);
510:     ISGetLocalSize(faceIS, &numFaces);
511:     ISGetIndices(faceIS, &faces);
512:     for (f = 0; f < numFaces; ++f) {
513:       PetscInt size;

515:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
516:       DMPlexGetSupportSize(dm, faces[f], &size);
517:       if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
518:       DMPlexSetSupportSize(gdm, faces[f] + *numGhostCells, 2);
519:     }
520:     ISRestoreIndices(faceIS, &faces);
521:     ISDestroy(&faceIS);
522:   }
523:   /* Step 4: Setup ghosted DM */
524:   DMSetUp(gdm);
525:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
526:   /* Step 6: Set cones and supports for new points */
527:   ghostCell = cEnd;
528:   for (fs = 0; fs < numFS; ++fs) {
529:     IS              faceIS;
530:     const PetscInt *faces;
531:     PetscInt        numFaces, f;

533:     DMLabelGetStratumIS(label, values[fs], &faceIS);
534:     ISGetLocalSize(faceIS, &numFaces);
535:     ISGetIndices(faceIS, &faces);
536:     for (f = 0; f < numFaces; ++f) {
537:       PetscInt newFace = faces[f] + *numGhostCells;

539:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
540:       DMPlexSetCone(gdm, ghostCell, &newFace);
541:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
542:       ++ghostCell;
543:     }
544:     ISRestoreIndices(faceIS, &faces);
545:     ISDestroy(&faceIS);
546:   }
547:   ISRestoreIndices(valueIS, &values);
548:   ISDestroy(&valueIS);
549:   /* Step 7: Stratify */
550:   DMPlexStratify(gdm);
551:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
552:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
553:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
554:   PetscFree(depthShift);
555:   return(0);
556: }

560: /*@C
561:   DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face

563:   Collective on dm

565:   Input Parameters:
566: + dm - The original DM
567: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL

569:   Output Parameters:
570: + numGhostCells - The number of ghost cells added to the DM
571: - dmGhosted - The new DM

573:   Note: If no label exists of that name, one will be created marking all boundary faces

575:   Level: developer

577: .seealso: DMCreate()
578: @*/
579: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
580: {
581:   DM             gdm;
582:   DMLabel        label;
583:   const char    *name = labelName ? labelName : "Face Sets";
584:   PetscInt       dim;
585:   PetscBool      flag;

592:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
593:   DMSetType(gdm, DMPLEX);
594:   DMPlexGetDimension(dm, &dim);
595:   DMPlexSetDimension(gdm, dim);
596:   DMPlexGetAdjacencyUseCone(dm, &flag);
597:   DMPlexSetAdjacencyUseCone(gdm, flag);
598:   DMPlexGetAdjacencyUseClosure(dm, &flag);
599:   DMPlexSetAdjacencyUseClosure(gdm, flag);
600:   DMPlexGetLabel(dm, name, &label);
601:   if (!label) {
602:     /* Get label for boundary faces */
603:     DMPlexCreateLabel(dm, name);
604:     DMPlexGetLabel(dm, name, &label);
605:     DMPlexMarkBoundaryFaces(dm, label);
606:   }
607:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
608:   *dmGhosted = gdm;
609:   return(0);
610: }

614: /*
615:   We are adding three kinds of points here:
616:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
617:     Non-replicated: Points which exist on the fault, but are not replicated
618:     Hybrid:         Entirely new points, such as cohesive cells

620:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
621:   each depth so that the new split/hybrid points can be inserted as a block.
622: */
623: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
624: {
625:   MPI_Comm         comm;
626:   IS               valueIS;
627:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
628:   const PetscInt  *values;          /* List of depths for which we have replicated points */
629:   IS              *splitIS;
630:   IS              *unsplitIS;
631:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
632:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
633:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
634:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
635:   const PetscInt **splitPoints;        /* Replicated points for each depth */
636:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
637:   PetscSection     coordSection;
638:   Vec              coordinates;
639:   PetscScalar     *coords;
640:   PetscInt         depths[4];          /* Depths in the order that plex points are numbered */
641:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
642:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
643:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
644:   PetscInt        *depthOffset;        /* Prefix sums of depthShift */
645:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
646:   PetscInt        *coneNew, *coneONew, *supportNew;
647:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
648:   PetscErrorCode   ierr;

651:   PetscObjectGetComm((PetscObject)dm,&comm);
652:   DMPlexGetDimension(dm, &dim);
653:   DMPlexGetDepth(dm, &depth);
654:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
655:   depths[0] = depth;
656:   depths[1] = 0;
657:   depths[2] = depth-1;
658:   depths[3] = 1;
659:   /* Count split points and add cohesive cells */
660:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
661:   PetscMalloc6(depth+1,&depthMax,depth+1,&depthEnd,depth+1,&depthShift,depth+1,&depthOffset,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
662:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
663:   PetscMemzero(depthShift,  (depth+1) * sizeof(PetscInt));
664:   PetscMemzero(depthOffset, (depth+1) * sizeof(PetscInt));
665:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
666:   for (d = 0; d <= depth; ++d) {
667:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
668:     depthEnd[d]           = pMaxNew[d];
669:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
670:     numSplitPoints[d]     = 0;
671:     numUnsplitPoints[d]   = 0;
672:     numHybridPoints[d]    = 0;
673:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
674:     splitPoints[d]        = NULL;
675:     unsplitPoints[d]      = NULL;
676:     splitIS[d]            = NULL;
677:     unsplitIS[d]          = NULL;
678:   }
679:   if (label) {
680:     DMLabelGetValueIS(label, &valueIS);
681:     ISGetLocalSize(valueIS, &numSP);
682:     ISGetIndices(valueIS, &values);
683:   }
684:   for (sp = 0; sp < numSP; ++sp) {
685:     const PetscInt dep = values[sp];

687:     if ((dep < 0) || (dep > depth)) continue;
688:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
689:     if (splitIS[dep]) {
690:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
691:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
692:     }
693:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
694:     if (unsplitIS[dep]) {
695:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
696:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
697:     }
698:   }
699:   /* Calculate number of hybrid points */
700:   for (d = 1; d <= depth; ++d) numHybridPoints[d]     = numSplitPoints[d-1] + numUnsplitPoints[d-1]; /* There is a hybrid cell/face/edge for every split face/edge/vertex   */
701:   for (d = 0; d <= depth; ++d) depthShift[d]          = numSplitPoints[d] + numHybridPoints[d];
702:   for (d = 1; d <= depth; ++d) depthOffset[depths[d]] = depthOffset[depths[d-1]] + depthShift[depths[d-1]];
703:   for (d = 0; d <= depth; ++d) pMaxNew[d]            += depthOffset[d] - numHybridPointsOld[d];
704:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
705:   /* Step 3: Set cone/support sizes for new points */
706:   for (dep = 0; dep <= depth; ++dep) {
707:     for (p = 0; p < numSplitPoints[dep]; ++p) {
708:       const PetscInt  oldp   = splitPoints[dep][p];
709:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
710:       const PetscInt  splitp = p    + pMaxNew[dep];
711:       const PetscInt *support;
712:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

714:       DMPlexGetConeSize(dm, oldp, &coneSize);
715:       DMPlexSetConeSize(sdm, splitp, coneSize);
716:       DMPlexGetSupportSize(dm, oldp, &supportSize);
717:       DMPlexSetSupportSize(sdm, splitp, supportSize);
718:       if (dep == depth-1) {
719:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

721:         /* Add cohesive cells, they are prisms */
722:         DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
723:       } else if (dep == 0) {
724:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

726:         DMPlexGetSupport(dm, oldp, &support);
727:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
728:           PetscInt val;

730:           DMLabelGetValue(label, support[e], &val);
731:           if (val == 1) ++qf;
732:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
733:           if ((val == 1) || (val == -(shift + 1))) ++qp;
734:         }
735:         /* Split old vertex: Edges into original vertex and new cohesive edge */
736:         DMPlexSetSupportSize(sdm, newp, qn+1);
737:         /* Split new vertex: Edges into split vertex and new cohesive edge */
738:         DMPlexSetSupportSize(sdm, splitp, qp+1);
739:         /* Add hybrid edge */
740:         DMPlexSetConeSize(sdm, hybedge, 2);
741:         DMPlexSetSupportSize(sdm, hybedge, qf);
742:       } else if (dep == dim-2) {
743:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

745:         DMPlexGetSupport(dm, oldp, &support);
746:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
747:           PetscInt val;

749:           DMLabelGetValue(label, support[e], &val);
750:           if (val == dim-1) ++qf;
751:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
752:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
753:         }
754:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
755:         DMPlexSetSupportSize(sdm, newp, qn+1);
756:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
757:         DMPlexSetSupportSize(sdm, splitp, qp+1);
758:         /* Add hybrid face */
759:         DMPlexSetConeSize(sdm, hybface, 4);
760:         DMPlexSetSupportSize(sdm, hybface, qf);
761:       }
762:     }
763:   }
764:   for (dep = 0; dep <= depth; ++dep) {
765:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
766:       const PetscInt  oldp   = unsplitPoints[dep][p];
767:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
768:       const PetscInt *support;
769:       PetscInt        coneSize, supportSize, qf, e, s;

771:       DMPlexGetConeSize(dm, oldp, &coneSize);
772:       DMPlexGetSupportSize(dm, oldp, &supportSize);
773:       DMPlexGetSupport(dm, oldp, &support);
774:       if (dep == 0) {
775:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

777:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
778:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
779:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
780:           if (e >= 0) ++qf;
781:         }
782:         DMPlexSetSupportSize(sdm, newp, qf+2);
783:         /* Add hybrid edge */
784:         DMPlexSetConeSize(sdm, hybedge, 2);
785:         for (e = 0, qf = 0; e < supportSize; ++e) {
786:           PetscInt val;

788:           DMLabelGetValue(label, support[e], &val);
789:           /* Split and unsplit edges produce hybrid faces */
790:           if (val == 1) ++qf;
791:           if (val == (shift2 + 1)) ++qf;
792:         }
793:         DMPlexSetSupportSize(sdm, hybedge, qf);
794:       } else if (dep == dim-2) {
795:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
796:         PetscInt       val;

798:         for (e = 0, qf = 0; e < supportSize; ++e) {
799:           DMLabelGetValue(label, support[e], &val);
800:           if (val == dim-1) qf += 2;
801:           else              ++qf;
802:         }
803:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
804:         DMPlexSetSupportSize(sdm, newp, qf+2);
805:         /* Add hybrid face */
806:         for (e = 0, qf = 0; e < supportSize; ++e) {
807:           DMLabelGetValue(label, support[e], &val);
808:           if (val == dim-1) ++qf;
809:         }
810:         DMPlexSetConeSize(sdm, hybface, 4);
811:         DMPlexSetSupportSize(sdm, hybface, qf);
812:       }
813:     }
814:   }
815:   /* Step 4: Setup split DM */
816:   DMSetUp(sdm);
817:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
818:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
819:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
820:   /* Step 6: Set cones and supports for new points */
821:   for (dep = 0; dep <= depth; ++dep) {
822:     for (p = 0; p < numSplitPoints[dep]; ++p) {
823:       const PetscInt  oldp   = splitPoints[dep][p];
824:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
825:       const PetscInt  splitp = p    + pMaxNew[dep];
826:       const PetscInt *cone, *support, *ornt;
827:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

829:       DMPlexGetConeSize(dm, oldp, &coneSize);
830:       DMPlexGetCone(dm, oldp, &cone);
831:       DMPlexGetConeOrientation(dm, oldp, &ornt);
832:       DMPlexGetSupportSize(dm, oldp, &supportSize);
833:       DMPlexGetSupport(dm, oldp, &support);
834:       if (dep == depth-1) {
835:         PetscBool       hasUnsplit = PETSC_FALSE;
836:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
837:         const PetscInt *supportF;

839:         /* Split face:       copy in old face to new face to start */
840:         DMPlexGetSupport(sdm, newp,  &supportF);
841:         DMPlexSetSupport(sdm, splitp, supportF);
842:         /* Split old face:   old vertices/edges in cone so no change */
843:         /* Split new face:   new vertices/edges in cone */
844:         for (q = 0; q < coneSize; ++q) {
845:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
846:           if (v < 0) {
847:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
848:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
849:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
850:             hasUnsplit   = PETSC_TRUE;
851:           } else {
852:             coneNew[2+q] = v + pMaxNew[dep-1];
853:             if (dep > 1) {
854:               const PetscInt *econe;
855:               PetscInt        econeSize, r, vs, vu;

857:               DMPlexGetConeSize(dm, cone[q], &econeSize);
858:               DMPlexGetCone(dm, cone[q], &econe);
859:               for (r = 0; r < econeSize; ++r) {
860:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
861:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
862:                 if (vs >= 0) continue;
863:                 if (vu < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", econe[r], dep-2);
864:                 hasUnsplit   = PETSC_TRUE;
865:               }
866:             }
867:           }
868:         }
869:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
870:         DMPlexSetConeOrientation(sdm, splitp, ornt);
871:         /* Face support */
872:         for (s = 0; s < supportSize; ++s) {
873:           PetscInt val;

875:           DMLabelGetValue(label, support[s], &val);
876:           if (val < 0) {
877:             /* Split old face:   Replace negative side cell with cohesive cell */
878:              DMPlexInsertSupport(sdm, newp, s, hybcell);
879:           } else {
880:             /* Split new face:   Replace positive side cell with cohesive cell */
881:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
882:             /* Get orientation for cohesive face */
883:             {
884:               const PetscInt *ncone, *nconeO;
885:               PetscInt        nconeSize, nc;

887:               DMPlexGetConeSize(dm, support[s], &nconeSize);
888:               DMPlexGetCone(dm, support[s], &ncone);
889:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
890:               for (nc = 0; nc < nconeSize; ++nc) {
891:                 if (ncone[nc] == oldp) {
892:                   coneONew[0] = nconeO[nc];
893:                   break;
894:                 }
895:               }
896:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
897:             }
898:           }
899:         }
900:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
901:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
902:         coneNew[1]  = splitp;
903:         coneONew[1] = coneONew[0];
904:         for (q = 0; q < coneSize; ++q) {
905:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
906:           if (v < 0) {
907:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
908:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
909:             coneONew[2+q] = 0;
910:           } else {
911:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
912:           }
913:           coneONew[2+q] = 0;
914:         }
915:         DMPlexSetCone(sdm, hybcell, coneNew);
916:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
917:         /* Label the hybrid cells on the boundary of the split */
918:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
919:       } else if (dep == 0) {
920:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

922:         /* Split old vertex: Edges in old split faces and new cohesive edge */
923:         for (e = 0, qn = 0; e < supportSize; ++e) {
924:           PetscInt val;

926:           DMLabelGetValue(label, support[e], &val);
927:           if ((val == 1) || (val == (shift + 1))) {
928:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
929:           }
930:         }
931:         supportNew[qn] = hybedge;
932:         DMPlexSetSupport(sdm, newp, supportNew);
933:         /* Split new vertex: Edges in new split faces and new cohesive edge */
934:         for (e = 0, qp = 0; e < supportSize; ++e) {
935:           PetscInt val, edge;

937:           DMLabelGetValue(label, support[e], &val);
938:           if (val == 1) {
939:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
940:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
941:             supportNew[qp++] = edge + pMaxNew[dep+1];
942:           } else if (val == -(shift + 1)) {
943:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
944:           }
945:         }
946:         supportNew[qp] = hybedge;
947:         DMPlexSetSupport(sdm, splitp, supportNew);
948:         /* Hybrid edge:    Old and new split vertex */
949:         coneNew[0] = newp;
950:         coneNew[1] = splitp;
951:         DMPlexSetCone(sdm, hybedge, coneNew);
952:         for (e = 0, qf = 0; e < supportSize; ++e) {
953:           PetscInt val, edge;

955:           DMLabelGetValue(label, support[e], &val);
956:           if (val == 1) {
957:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
958:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
959:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
960:           }
961:         }
962:         DMPlexSetSupport(sdm, hybedge, supportNew);
963:       } else if (dep == dim-2) {
964:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

966:         /* Split old edge:   old vertices in cone so no change */
967:         /* Split new edge:   new vertices in cone */
968:         for (q = 0; q < coneSize; ++q) {
969:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
970:           if (v < 0) {
971:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
972:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
973:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
974:           } else {
975:             coneNew[q] = v + pMaxNew[dep-1];
976:           }
977:         }
978:         DMPlexSetCone(sdm, splitp, coneNew);
979:         /* Split old edge: Faces in positive side cells and old split faces */
980:         for (e = 0, q = 0; e < supportSize; ++e) {
981:           PetscInt val;

983:           DMLabelGetValue(label, support[e], &val);
984:           if (val == dim-1) {
985:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
986:           } else if (val == (shift + dim-1)) {
987:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
988:           }
989:         }
990:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
991:         DMPlexSetSupport(sdm, newp, supportNew);
992:         /* Split new edge: Faces in negative side cells and new split faces */
993:         for (e = 0, q = 0; e < supportSize; ++e) {
994:           PetscInt val, face;

996:           DMLabelGetValue(label, support[e], &val);
997:           if (val == dim-1) {
998:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
999:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1000:             supportNew[q++] = face + pMaxNew[dep+1];
1001:           } else if (val == -(shift + dim-1)) {
1002:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1003:           }
1004:         }
1005:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1006:         DMPlexSetSupport(sdm, splitp, supportNew);
1007:         /* Hybrid face */
1008:         coneNew[0] = newp;
1009:         coneNew[1] = splitp;
1010:         for (v = 0; v < coneSize; ++v) {
1011:           PetscInt vertex;
1012:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1013:           if (vertex < 0) {
1014:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1015:             if (vertex < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[v], dep-1);
1016:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1017:           } else {
1018:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1019:           }
1020:         }
1021:         DMPlexSetCone(sdm, hybface, coneNew);
1022:         for (e = 0, qf = 0; e < supportSize; ++e) {
1023:           PetscInt val, face;

1025:           DMLabelGetValue(label, support[e], &val);
1026:           if (val == dim-1) {
1027:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1028:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1029:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1030:           }
1031:         }
1032:         DMPlexSetSupport(sdm, hybface, supportNew);
1033:       }
1034:     }
1035:   }
1036:   for (dep = 0; dep <= depth; ++dep) {
1037:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1038:       const PetscInt  oldp   = unsplitPoints[dep][p];
1039:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
1040:       const PetscInt *cone, *support, *ornt;
1041:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1043:       DMPlexGetConeSize(dm, oldp, &coneSize);
1044:       DMPlexGetCone(dm, oldp, &cone);
1045:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1046:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1047:       DMPlexGetSupport(dm, oldp, &support);
1048:       if (dep == 0) {
1049:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1051:         /* Unsplit vertex */
1052:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1053:         for (s = 0, q = 0; s < supportSize; ++s) {
1054:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthMax, depthEnd, depthShift) /*support[s] + depthOffset[dep+1]*/;
1055:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1056:           if (e >= 0) {
1057:             supportNew[q++] = e + pMaxNew[dep+1];
1058:           }
1059:         }
1060:         supportNew[q++] = hybedge;
1061:         supportNew[q++] = hybedge;
1062:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1063:         DMPlexSetSupport(sdm, newp, supportNew);
1064:         /* Hybrid edge */
1065:         coneNew[0] = newp;
1066:         coneNew[1] = newp;
1067:         DMPlexSetCone(sdm, hybedge, coneNew);
1068:         for (e = 0, qf = 0; e < supportSize; ++e) {
1069:           PetscInt val, edge;

1071:           DMLabelGetValue(label, support[e], &val);
1072:           if (val == 1) {
1073:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1074:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1075:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1076:           } else if  (val ==  (shift2 + 1)) {
1077:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1078:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1079:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1080:           }
1081:         }
1082:         DMPlexSetSupport(sdm, hybedge, supportNew);
1083:       } else if (dep == dim-2) {
1084:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1086:         /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1087:         for (f = 0, qf = 0; f < supportSize; ++f) {
1088:           PetscInt val, face;

1090:           DMLabelGetValue(label, support[f], &val);
1091:           if (val == dim-1) {
1092:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1093:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1094:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1095:             supportNew[qf++] = face + pMaxNew[dep+1];
1096:           } else {
1097:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1098:           }
1099:         }
1100:         supportNew[qf++] = hybface;
1101:         supportNew[qf++] = hybface;
1102:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1103:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1104:         DMPlexSetSupport(sdm, newp, supportNew);
1105:         /* Add hybrid face */
1106:         coneNew[0] = newp;
1107:         coneNew[1] = newp;
1108:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1109:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1110:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1111:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1112:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1113:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1114:         DMPlexSetCone(sdm, hybface, coneNew);
1115:         for (f = 0, qf = 0; f < supportSize; ++f) {
1116:           PetscInt val, face;

1118:           DMLabelGetValue(label, support[f], &val);
1119:           if (val == dim-1) {
1120:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1121:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1122:           }
1123:         }
1124:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1125:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1126:         DMPlexSetSupport(sdm, hybface, supportNew);
1127:       }
1128:     }
1129:   }
1130:   /* Step 6b: Replace split points in negative side cones */
1131:   for (sp = 0; sp < numSP; ++sp) {
1132:     PetscInt        dep = values[sp];
1133:     IS              pIS;
1134:     PetscInt        numPoints;
1135:     const PetscInt *points;

1137:     if (dep >= 0) continue;
1138:     DMLabelGetStratumIS(label, dep, &pIS);
1139:     if (!pIS) continue;
1140:     dep  = -dep - shift;
1141:     ISGetLocalSize(pIS, &numPoints);
1142:     ISGetIndices(pIS, &points);
1143:     for (p = 0; p < numPoints; ++p) {
1144:       const PetscInt  oldp = points[p];
1145:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + oldp*/;
1146:       const PetscInt *cone;
1147:       PetscInt        coneSize, c;
1148:       /* PetscBool       replaced = PETSC_FALSE; */

1150:       /* Negative edge: replace split vertex */
1151:       /* Negative cell: replace split face */
1152:       DMPlexGetConeSize(sdm, newp, &coneSize);
1153:       DMPlexGetCone(sdm, newp, &cone);
1154:       for (c = 0; c < coneSize; ++c) {
1155:         const PetscInt coldp = cone[c] - depthOffset[dep-1];
1156:         PetscInt       csplitp, cp, val;

1158:         DMLabelGetValue(label, coldp, &val);
1159:         if (val == dep-1) {
1160:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1161:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1162:           csplitp  = pMaxNew[dep-1] + cp;
1163:           DMPlexInsertCone(sdm, newp, c, csplitp);
1164:           /* replaced = PETSC_TRUE; */
1165:         }
1166:       }
1167:       /* Cells with only a vertex or edge on the submesh have no replacement */
1168:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1169:     }
1170:     ISRestoreIndices(pIS, &points);
1171:     ISDestroy(&pIS);
1172:   }
1173:   /* Step 7: Stratify */
1174:   DMPlexStratify(sdm);
1175:   /* Step 8: Coordinates */
1176:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1177:   DMGetCoordinateSection(sdm, &coordSection);
1178:   DMGetCoordinatesLocal(sdm, &coordinates);
1179:   VecGetArray(coordinates, &coords);
1180:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1181:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthMax, depthEnd, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1182:     const PetscInt splitp = pMaxNew[0] + v;
1183:     PetscInt       dof, off, soff, d;

1185:     PetscSectionGetDof(coordSection, newp, &dof);
1186:     PetscSectionGetOffset(coordSection, newp, &off);
1187:     PetscSectionGetOffset(coordSection, splitp, &soff);
1188:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1189:   }
1190:   VecRestoreArray(coordinates, &coords);
1191:   /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1192:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1193:   /* Step 10: Labels */
1194:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1195:   DMPlexGetNumLabels(sdm, &numLabels);
1196:   for (dep = 0; dep <= depth; ++dep) {
1197:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1198:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1199:       const PetscInt splitp = pMaxNew[dep] + p;
1200:       PetscInt       l;

1202:       for (l = 0; l < numLabels; ++l) {
1203:         DMLabel     mlabel;
1204:         const char *lname;
1205:         PetscInt    val;
1206:         PetscBool   isDepth;

1208:         DMPlexGetLabelName(sdm, l, &lname);
1209:         PetscStrcmp(lname, "depth", &isDepth);
1210:         if (isDepth) continue;
1211:         DMPlexGetLabel(sdm, lname, &mlabel);
1212:         DMLabelGetValue(mlabel, newp, &val);
1213:         if (val >= 0) {
1214:           DMLabelSetValue(mlabel, splitp, val);
1215: #if 0
1216:           /* Do not put cohesive edges into the label */
1217:           if (dep == 0) {
1218:             const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1219:             DMLabelSetValue(mlabel, cedge, val);
1220:           } else if (dep == dim-2) {
1221:             const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1222:             DMLabelSetValue(mlabel, cface, val);
1223:           }
1224:           /* Do not put cohesive faces into the label */
1225: #endif
1226:         }
1227:       }
1228:     }
1229:   }
1230:   for (sp = 0; sp < numSP; ++sp) {
1231:     const PetscInt dep = values[sp];

1233:     if ((dep < 0) || (dep > depth)) continue;
1234:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1235:     ISDestroy(&splitIS[dep]);
1236:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1237:     ISDestroy(&unsplitIS[dep]);
1238:   }
1239:   if (label) {
1240:     ISRestoreIndices(valueIS, &values);
1241:     ISDestroy(&valueIS);
1242:   }
1243:   for (d = 0; d <= depth; ++d) {
1244:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1245:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1246:   }
1247:   DMPlexSetHybridBounds(sdm, depth >= 0 ? pMaxNew[depth] : PETSC_DETERMINE, depth>1 ? pMaxNew[depth-1] : PETSC_DETERMINE, depth>2 ? pMaxNew[1] : PETSC_DETERMINE, depth >= 0 ? pMaxNew[0] : PETSC_DETERMINE);
1248:   PetscFree3(coneNew, coneONew, supportNew);
1249:   PetscFree6(depthMax, depthEnd, depthShift, depthOffset, pMaxNew, numHybridPointsOld);
1250:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1251:   return(0);
1252: }

1256: /*@C
1257:   DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface

1259:   Collective on dm

1261:   Input Parameters:
1262: + dm - The original DM
1263: - label - The label specifying the boundary faces (this could be auto-generated)

1265:   Output Parameters:
1266: - dmSplit - The new DM

1268:   Level: developer

1270: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1271: @*/
1272: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1273: {
1274:   DM             sdm;
1275:   PetscInt       dim;

1281:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1282:   DMSetType(sdm, DMPLEX);
1283:   DMPlexGetDimension(dm, &dim);
1284:   DMPlexSetDimension(sdm, dim);
1285:   switch (dim) {
1286:   case 2:
1287:   case 3:
1288:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1289:     break;
1290:   default:
1291:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1292:   }
1293:   *dmSplit = sdm;
1294:   return(0);
1295: }

1299: /*@
1300:   DMPlexLabelCohesiveComplete - Starting with a label marking vertices on an internal surface, we add all other mesh pieces
1301:   to complete the surface

1303:   Input Parameters:
1304: + dm - The DM
1305: . label - A DMLabel marking the surface vertices
1306: . flip  - Flag to flip the submesh normal and replace points on the other side
1307: - subdm - The subDM associated with the label, or NULL

1309:   Output Parameter:
1310: . label - A DMLabel marking all surface points

1312:   Level: developer

1314: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1315: @*/
1316: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, PetscBool flip, DM subdm)
1317: {
1318:   DMLabel         depthLabel;
1319:   IS              dimIS, subpointIS, facePosIS, faceNegIS;
1320:   const PetscInt *points, *subpoints;
1321:   const PetscInt  rev   = flip ? -1 : 1;
1322:   PetscInt        shift = 100, shift2 = 200, dim, dep, cStart, cEnd, numPoints, numSubpoints, p, val;
1323:   PetscErrorCode  ierr;

1326:   DMPlexGetDepthLabel(dm, &depthLabel);
1327:   DMPlexGetDimension(dm, &dim);
1328:   if (subdm) {
1329:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1330:     if (subpointIS) {
1331:       ISGetLocalSize(subpointIS, &numSubpoints);
1332:       ISGetIndices(subpointIS, &subpoints);
1333:     }
1334:   }
1335:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1336:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1337:   if (!dimIS) return(0);
1338:   ISGetLocalSize(dimIS, &numPoints);
1339:   ISGetIndices(dimIS, &points);
1340:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1341:     const PetscInt *support;
1342:     PetscInt        supportSize, s;

1344:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1345:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1346:     DMPlexGetSupport(dm, points[p], &support);
1347:     for (s = 0; s < supportSize; ++s) {
1348:       const PetscInt *cone, *ornt;
1349:       PetscInt        coneSize, c;
1350:       PetscBool       pos = PETSC_TRUE;

1352:       DMPlexGetConeSize(dm, support[s], &coneSize);
1353:       DMPlexGetCone(dm, support[s], &cone);
1354:       DMPlexGetConeOrientation(dm, support[s], &ornt);
1355:       for (c = 0; c < coneSize; ++c) {
1356:         if (cone[c] == points[p]) {
1357:           PetscInt o = ornt[c];

1359:           if (subdm) {
1360:             const PetscInt *subcone, *subornt;
1361:             PetscInt        subpoint, subface, subconeSize, sc;

1363:             PetscFindInt(support[s], numSubpoints, subpoints, &subpoint);
1364:             PetscFindInt(points[p],  numSubpoints, subpoints, &subface);
1365:             DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1366:             DMPlexGetCone(subdm, subpoint, &subcone);
1367:             DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1368:             for (sc = 0; sc < subconeSize; ++sc) {
1369:               if (subcone[sc] == subface) {
1370:                 o = subornt[0];
1371:                 break;
1372:               }
1373:             }
1374:             if (sc >= subconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find point %d in cone for subpoint %d", points[p], subpoint);
1375:           }
1376:           if (o >= 0) {
1377:             DMLabelSetValue(label, support[s],  rev*(shift+dim));
1378:             pos  = rev > 0 ? PETSC_TRUE : PETSC_FALSE;
1379:           } else {
1380:             DMLabelSetValue(label, support[s], -rev*(shift+dim));
1381:             pos  = rev > 0 ? PETSC_FALSE : PETSC_TRUE;
1382:           }
1383:           break;
1384:         }
1385:       }
1386:       if (c == coneSize) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell split face %d support does not have it in the cone", points[p]);
1387:       /* Put faces touching the fault in the label */
1388:       for (c = 0; c < coneSize; ++c) {
1389:         const PetscInt point = cone[c];

1391:         DMLabelGetValue(label, point, &val);
1392:         if (val == -1) {
1393:           PetscInt *closure = NULL;
1394:           PetscInt  closureSize, cl;

1396:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1397:           for (cl = 0; cl < closureSize*2; cl += 2) {
1398:             const PetscInt clp = closure[cl];

1400:             DMLabelGetValue(label, clp, &val);
1401:             if ((val >= 0) && (val < dim-1)) {
1402:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1403:               break;
1404:             }
1405:           }
1406:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1407:         }
1408:       }
1409:     }
1410:   }
1411:   if (subdm) {
1412:     if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1413:     ISDestroy(&subpointIS);
1414:   }
1415:   ISRestoreIndices(dimIS, &points);
1416:   ISDestroy(&dimIS);
1417:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1418:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1419:   DMLabelGetStratumIS(label, 0, &dimIS);
1420:   if (!dimIS) return(0);
1421:   ISGetLocalSize(dimIS, &numPoints);
1422:   ISGetIndices(dimIS, &points);
1423:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1424:     PetscInt *star = NULL;
1425:     PetscInt  starSize, s;
1426:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1428:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1429:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1430:     while (again) {
1431:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1432:       again = 0;
1433:       for (s = 0; s < starSize*2; s += 2) {
1434:         const PetscInt  point = star[s];
1435:         const PetscInt *cone;
1436:         PetscInt        coneSize, c;

1438:         if ((point < cStart) || (point >= cEnd)) continue;
1439:         DMLabelGetValue(label, point, &val);
1440:         if (val != -1) continue;
1441:         again = again == 1 ? 1 : 2;
1442:         DMPlexGetConeSize(dm, point, &coneSize);
1443:         DMPlexGetCone(dm, point, &cone);
1444:         for (c = 0; c < coneSize; ++c) {
1445:           DMLabelGetValue(label, cone[c], &val);
1446:           if (val != -1) {
1447:             const PetscInt *ccone;
1448:             PetscInt        cconeSize, cc, side;

1450:             if (abs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
1451:             if (val > 0) side =  1;
1452:             else         side = -1;
1453:             DMLabelSetValue(label, point, side*(shift+dim));
1454:             /* Mark cell faces which touch the fault */
1455:             DMPlexGetConeSize(dm, point, &cconeSize);
1456:             DMPlexGetCone(dm, point, &ccone);
1457:             for (cc = 0; cc < cconeSize; ++cc) {
1458:               PetscInt *closure = NULL;
1459:               PetscInt  closureSize, cl;

1461:               DMLabelGetValue(label, ccone[cc], &val);
1462:               if (val != -1) continue;
1463:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1464:               for (cl = 0; cl < closureSize*2; cl += 2) {
1465:                 const PetscInt clp = closure[cl];

1467:                 DMLabelGetValue(label, clp, &val);
1468:                 if (val == -1) continue;
1469:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1470:                 break;
1471:               }
1472:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1473:             }
1474:             again = 1;
1475:             break;
1476:           }
1477:         }
1478:       }
1479:     }
1480:     /* Classify the rest by cell membership */
1481:     for (s = 0; s < starSize*2; s += 2) {
1482:       const PetscInt point = star[s];

1484:       DMLabelGetValue(label, point, &val);
1485:       if (val == -1) {
1486:         PetscInt  *sstar = NULL;
1487:         PetscInt   sstarSize, ss;
1488:         PetscBool  marked = PETSC_FALSE;

1490:         DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1491:         for (ss = 0; ss < sstarSize*2; ss += 2) {
1492:           const PetscInt spoint = sstar[ss];

1494:           if ((spoint < cStart) || (spoint >= cEnd)) continue;
1495:           DMLabelGetValue(label, spoint, &val);
1496:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1497:           DMLabelGetValue(depthLabel, point, &dep);
1498:           if (val > 0) {
1499:             DMLabelSetValue(label, point,   shift+dep);
1500:           } else {
1501:             DMLabelSetValue(label, point, -(shift+dep));
1502:           }
1503:           marked = PETSC_TRUE;
1504:           break;
1505:         }
1506:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1507:         if (!marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1508:       }
1509:     }
1510:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1511:   }
1512:   ISRestoreIndices(dimIS, &points);
1513:   ISDestroy(&dimIS);
1514:   /* If any faces touching the fault divide cells on either side, split them */
1515:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1516:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1517:   ISExpand(facePosIS, faceNegIS, &dimIS);
1518:   ISDestroy(&facePosIS);
1519:   ISDestroy(&faceNegIS);
1520:   ISGetLocalSize(dimIS, &numPoints);
1521:   ISGetIndices(dimIS, &points);
1522:   for (p = 0; p < numPoints; ++p) {
1523:     const PetscInt  point = points[p];
1524:     const PetscInt *support;
1525:     PetscInt        supportSize, valA, valB;

1527:     DMPlexGetSupportSize(dm, point, &supportSize);
1528:     if (supportSize != 2) continue;
1529:     DMPlexGetSupport(dm, point, &support);
1530:     DMLabelGetValue(label, support[0], &valA);
1531:     DMLabelGetValue(label, support[1], &valB);
1532:     if ((valA == -1) || (valB == -1)) continue;
1533:     if (valA*valB > 0) continue;
1534:     /* Split the face */
1535:     DMLabelGetValue(label, point, &valA);
1536:     DMLabelClearValue(label, point, valA);
1537:     DMLabelSetValue(label, point, dim-1);
1538:     /* Label its closure:
1539:       unmarked: label as unsplit
1540:       incident: relabel as split
1541:       split:    do nothing
1542:     */
1543:     {
1544:       PetscInt *closure = NULL;
1545:       PetscInt  closureSize, cl;

1547:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1548:       for (cl = 0; cl < closureSize*2; cl += 2) {
1549:         DMLabelGetValue(label, closure[cl], &valA);
1550:         if (valA == -1) { /* Mark as unsplit */
1551:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1552:           DMLabelSetValue(label, closure[cl], shift2+dep);
1553:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1554:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1555:           DMLabelClearValue(label, closure[cl], valA);
1556:           DMLabelSetValue(label, closure[cl], dep);
1557:         }
1558:       }
1559:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1560:     }
1561:   }
1562:   ISRestoreIndices(dimIS, &points);
1563:   ISDestroy(&dimIS);
1564:   return(0);
1565: }

1569: /*@C
1570:   DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface

1572:   Collective on dm

1574:   Input Parameters:
1575: + dm - The original DM
1576: - labelName - The label specifying the interface vertices

1578:   Output Parameters:
1579: + hybridLabel - The label fully marking the interface
1580: - dmHybrid - The new DM

1582:   Level: developer

1584: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1585: @*/
1586: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1587: {
1588:   DM             idm;
1589:   DMLabel        subpointMap, hlabel;
1590:   PetscInt       dim;

1597:   DMPlexGetDimension(dm, &dim);
1598:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1599:   DMPlexOrient(idm);
1600:   DMPlexGetSubpointMap(idm, &subpointMap);
1601:   DMLabelDuplicate(subpointMap, &hlabel);
1602:   DMLabelClearStratum(hlabel, dim);
1603:   DMPlexLabelCohesiveComplete(dm, hlabel, PETSC_FALSE, idm);
1604:   DMDestroy(&idm);
1605:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1606:   if (hybridLabel) *hybridLabel = hlabel;
1607:   else             {DMLabelDestroy(&hlabel);}
1608:   return(0);
1609: }

1613: /* Here we need the explicit assumption that:

1615:      For any marked cell, the marked vertices constitute a single face
1616: */
1617: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1618: {
1619:   IS               subvertexIS = NULL;
1620:   const PetscInt  *subvertices;
1621:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1622:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1623:   PetscErrorCode   ierr;

1626:   *numFaces = 0;
1627:   *nFV      = 0;
1628:   DMPlexGetDepth(dm, &depth);
1629:   DMPlexGetDimension(dm, &dim);
1630:   pSize = PetscMax(depth, dim) + 1;
1631:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1632:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1633:   for (d = 0; d <= depth; ++d) {
1634:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1635:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1636:   }
1637:   /* Loop over initial vertices and mark all faces in the collective star() */
1638:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1639:   if (subvertexIS) {
1640:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1641:     ISGetIndices(subvertexIS, &subvertices);
1642:   }
1643:   for (v = 0; v < numSubVerticesInitial; ++v) {
1644:     const PetscInt vertex = subvertices[v];
1645:     PetscInt      *star   = NULL;
1646:     PetscInt       starSize, s, numCells = 0, c;

1648:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1649:     for (s = 0; s < starSize*2; s += 2) {
1650:       const PetscInt point = star[s];
1651:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1652:     }
1653:     for (c = 0; c < numCells; ++c) {
1654:       const PetscInt cell    = star[c];
1655:       PetscInt      *closure = NULL;
1656:       PetscInt       closureSize, cl;
1657:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1659:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1660:       if (cellLoc == 2) continue;
1661:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1662:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1663:       for (cl = 0; cl < closureSize*2; cl += 2) {
1664:         const PetscInt point = closure[cl];
1665:         PetscInt       vertexLoc;

1667:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1668:           ++numCorners;
1669:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1670:           if (vertexLoc == value) closure[faceSize++] = point;
1671:         }
1672:       }
1673:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1674:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1675:       if (faceSize == *nFV) {
1676:         const PetscInt *cells = NULL;
1677:         PetscInt        numCells, nc;

1679:         ++(*numFaces);
1680:         for (cl = 0; cl < faceSize; ++cl) {
1681:           DMLabelSetValue(subpointMap, closure[cl], 0);
1682:         }
1683:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1684:         for (nc = 0; nc < numCells; ++nc) {
1685:           DMLabelSetValue(subpointMap, cells[nc], 2);
1686:         }
1687:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1688:       }
1689:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1690:     }
1691:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1692:   }
1693:   if (subvertexIS) {
1694:     ISRestoreIndices(subvertexIS, &subvertices);
1695:   }
1696:   ISDestroy(&subvertexIS);
1697:   PetscFree3(pStart,pEnd,pMax);
1698:   return(0);
1699: }

1703: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1704: {
1705:   IS               subvertexIS;
1706:   const PetscInt  *subvertices;
1707:   PetscInt        *pStart, *pEnd, *pMax;
1708:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
1709:   PetscErrorCode   ierr;

1712:   DMPlexGetDimension(dm, &dim);
1713:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1714:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1715:   for (d = 0; d <= dim; ++d) {
1716:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1717:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1718:   }
1719:   /* Loop over initial vertices and mark all faces in the collective star() */
1720:   DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1721:   if (subvertexIS) {
1722:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1723:     ISGetIndices(subvertexIS, &subvertices);
1724:   }
1725:   for (v = 0; v < numSubVerticesInitial; ++v) {
1726:     const PetscInt vertex = subvertices[v];
1727:     PetscInt      *star   = NULL;
1728:     PetscInt       starSize, s, numFaces = 0, f;

1730:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1731:     for (s = 0; s < starSize*2; s += 2) {
1732:       const PetscInt point = star[s];
1733:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1734:     }
1735:     for (f = 0; f < numFaces; ++f) {
1736:       const PetscInt face    = star[f];
1737:       PetscInt      *closure = NULL;
1738:       PetscInt       closureSize, c;
1739:       PetscInt       faceLoc;

1741:       DMLabelGetValue(subpointMap, face, &faceLoc);
1742:       if (faceLoc == dim-1) continue;
1743:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1744:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1745:       for (c = 0; c < closureSize*2; c += 2) {
1746:         const PetscInt point = closure[c];
1747:         PetscInt       vertexLoc;

1749:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1750:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1751:           if (vertexLoc != value) break;
1752:         }
1753:       }
1754:       if (c == closureSize*2) {
1755:         const PetscInt *support;
1756:         PetscInt        supportSize, s;

1758:         for (c = 0; c < closureSize*2; c += 2) {
1759:           const PetscInt point = closure[c];

1761:           for (d = 0; d < dim; ++d) {
1762:             if ((point >= pStart[d]) && (point < pEnd[d])) {
1763:               DMLabelSetValue(subpointMap, point, d);
1764:               break;
1765:             }
1766:           }
1767:         }
1768:         DMPlexGetSupportSize(dm, face, &supportSize);
1769:         DMPlexGetSupport(dm, face, &support);
1770:         for (s = 0; s < supportSize; ++s) {
1771:           DMLabelSetValue(subpointMap, support[s], dim);
1772:         }
1773:       }
1774:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1775:     }
1776:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1777:   }
1778:   if (subvertexIS) {
1779:     ISRestoreIndices(subvertexIS, &subvertices);
1780:   }
1781:   ISDestroy(&subvertexIS);
1782:   PetscFree3(pStart,pEnd,pMax);
1783:   return(0);
1784: }

1788: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
1789: {
1790:   DMLabel         label = NULL;
1791:   const PetscInt *cone;
1792:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize;
1793:   PetscErrorCode  ierr;

1796:   *numFaces = 0;
1797:   *nFV = 0;
1798:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1799:   *subCells = NULL;
1800:   DMPlexGetDimension(dm, &dim);
1801:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1802:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1803:   if (cMax < 0) return(0);
1804:   if (label) {
1805:     for (c = cMax; c < cEnd; ++c) {
1806:       PetscInt val;

1808:       DMLabelGetValue(label, c, &val);
1809:       if (val == value) {
1810:         ++(*numFaces);
1811:         DMPlexGetConeSize(dm, c, &coneSize);
1812:       }
1813:     }
1814:   } else {
1815:     *numFaces = cEnd - cMax;
1816:     DMPlexGetConeSize(dm, cMax, &coneSize);
1817:   }
1818:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
1819:   PetscMalloc1(*numFaces *2, subCells);
1820:   for (c = cMax; c < cEnd; ++c) {
1821:     const PetscInt *cells;
1822:     PetscInt        numCells;

1824:     if (label) {
1825:       PetscInt val;

1827:       DMLabelGetValue(label, c, &val);
1828:       if (val != value) continue;
1829:     }
1830:     DMPlexGetCone(dm, c, &cone);
1831:     for (p = 0; p < *nFV; ++p) {
1832:       DMLabelSetValue(subpointMap, cone[p], 0);
1833:     }
1834:     /* Negative face */
1835:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
1836:     /* Not true in parallel
1837:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
1838:     for (p = 0; p < numCells; ++p) {
1839:       DMLabelSetValue(subpointMap, cells[p], 2);
1840:       (*subCells)[subc++] = cells[p];
1841:     }
1842:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
1843:     /* Positive face is not included */
1844:   }
1845:   return(0);
1846: }

1850: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DMLabel subpointMap, DM subdm)
1851: {
1852:   DMLabel        label = NULL;
1853:   PetscInt      *pStart, *pEnd;
1854:   PetscInt       dim, cMax, cEnd, c, d;

1858:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1859:   DMPlexGetDimension(dm, &dim);
1860:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1861:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1862:   if (cMax < 0) return(0);
1863:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
1864:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1865:   for (c = cMax; c < cEnd; ++c) {
1866:     const PetscInt *cone;
1867:     PetscInt       *closure = NULL;
1868:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

1870:     if (label) {
1871:       DMLabelGetValue(label, c, &val);
1872:       if (val != value) continue;
1873:     }
1874:     DMPlexGetConeSize(dm, c, &coneSize);
1875:     DMPlexGetCone(dm, c, &cone);
1876:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
1877:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
1878:     /* Negative face */
1879:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1880:     for (cl = 0; cl < closureSize*2; cl += 2) {
1881:       const PetscInt point = closure[cl];

1883:       for (d = 0; d <= dim; ++d) {
1884:         if ((point >= pStart[d]) && (point < pEnd[d])) {
1885:           DMLabelSetValue(subpointMap, point, d);
1886:           break;
1887:         }
1888:       }
1889:     }
1890:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1891:     /* Cells -- positive face is not included */
1892:     for (cl = 0; cl < 1; ++cl) {
1893:       const PetscInt *support;
1894:       PetscInt        supportSize, s;

1896:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
1897:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
1898:       DMPlexGetSupport(dm, cone[cl], &support);
1899:       for (s = 0; s < supportSize; ++s) {
1900:         DMLabelSetValue(subpointMap, support[s], dim);
1901:       }
1902:     }
1903:   }
1904:   PetscFree2(pStart, pEnd);
1905:   return(0);
1906: }

1910: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
1911: {
1912:   MPI_Comm       comm;
1913:   PetscBool      posOrient = PETSC_FALSE;
1914:   const PetscInt debug     = 0;
1915:   PetscInt       cellDim, faceSize, f;

1919:   PetscObjectGetComm((PetscObject)dm,&comm);
1920:   DMPlexGetDimension(dm, &cellDim);
1921:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

1923:   if (cellDim == 1 && numCorners == 2) {
1924:     /* Triangle */
1925:     faceSize  = numCorners-1;
1926:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
1927:   } else if (cellDim == 2 && numCorners == 3) {
1928:     /* Triangle */
1929:     faceSize  = numCorners-1;
1930:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
1931:   } else if (cellDim == 3 && numCorners == 4) {
1932:     /* Tetrahedron */
1933:     faceSize  = numCorners-1;
1934:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
1935:   } else if (cellDim == 1 && numCorners == 3) {
1936:     /* Quadratic line */
1937:     faceSize  = 1;
1938:     posOrient = PETSC_TRUE;
1939:   } else if (cellDim == 2 && numCorners == 4) {
1940:     /* Quads */
1941:     faceSize = 2;
1942:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
1943:       posOrient = PETSC_TRUE;
1944:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
1945:       posOrient = PETSC_TRUE;
1946:     } else {
1947:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
1948:         posOrient = PETSC_FALSE;
1949:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
1950:     }
1951:   } else if (cellDim == 2 && numCorners == 6) {
1952:     /* Quadratic triangle (I hate this) */
1953:     /* Edges are determined by the first 2 vertices (corners of edges) */
1954:     const PetscInt faceSizeTri = 3;
1955:     PetscInt       sortedIndices[3], i, iFace;
1956:     PetscBool      found                    = PETSC_FALSE;
1957:     PetscInt       faceVerticesTriSorted[9] = {
1958:       0, 3,  4, /* bottom */
1959:       1, 4,  5, /* right */
1960:       2, 3,  5, /* left */
1961:     };
1962:     PetscInt       faceVerticesTri[9] = {
1963:       0, 3,  4, /* bottom */
1964:       1, 4,  5, /* right */
1965:       2, 5,  3, /* left */
1966:     };

1968:     faceSize = faceSizeTri;
1969:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
1970:     PetscSortInt(faceSizeTri, sortedIndices);
1971:     for (iFace = 0; iFace < 3; ++iFace) {
1972:       const PetscInt ii = iFace*faceSizeTri;
1973:       PetscInt       fVertex, cVertex;

1975:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
1976:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
1977:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
1978:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
1979:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
1980:               faceVertices[fVertex] = origVertices[cVertex];
1981:               break;
1982:             }
1983:           }
1984:         }
1985:         found = PETSC_TRUE;
1986:         break;
1987:       }
1988:     }
1989:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
1990:     if (posOriented) *posOriented = PETSC_TRUE;
1991:     return(0);
1992:   } else if (cellDim == 2 && numCorners == 9) {
1993:     /* Quadratic quad (I hate this) */
1994:     /* Edges are determined by the first 2 vertices (corners of edges) */
1995:     const PetscInt faceSizeQuad = 3;
1996:     PetscInt       sortedIndices[3], i, iFace;
1997:     PetscBool      found                      = PETSC_FALSE;
1998:     PetscInt       faceVerticesQuadSorted[12] = {
1999:       0, 1,  4, /* bottom */
2000:       1, 2,  5, /* right */
2001:       2, 3,  6, /* top */
2002:       0, 3,  7, /* left */
2003:     };
2004:     PetscInt       faceVerticesQuad[12] = {
2005:       0, 1,  4, /* bottom */
2006:       1, 2,  5, /* right */
2007:       2, 3,  6, /* top */
2008:       3, 0,  7, /* left */
2009:     };

2011:     faceSize = faceSizeQuad;
2012:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2013:     PetscSortInt(faceSizeQuad, sortedIndices);
2014:     for (iFace = 0; iFace < 4; ++iFace) {
2015:       const PetscInt ii = iFace*faceSizeQuad;
2016:       PetscInt       fVertex, cVertex;

2018:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2019:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2020:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2021:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2022:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2023:               faceVertices[fVertex] = origVertices[cVertex];
2024:               break;
2025:             }
2026:           }
2027:         }
2028:         found = PETSC_TRUE;
2029:         break;
2030:       }
2031:     }
2032:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2033:     if (posOriented) *posOriented = PETSC_TRUE;
2034:     return(0);
2035:   } else if (cellDim == 3 && numCorners == 8) {
2036:     /* Hexes
2037:        A hex is two oriented quads with the normal of the first
2038:        pointing up at the second.

2040:           7---6
2041:          /|  /|
2042:         4---5 |
2043:         | 1-|-2
2044:         |/  |/
2045:         0---3

2047:         Faces are determined by the first 4 vertices (corners of faces) */
2048:     const PetscInt faceSizeHex = 4;
2049:     PetscInt       sortedIndices[4], i, iFace;
2050:     PetscBool      found                     = PETSC_FALSE;
2051:     PetscInt       faceVerticesHexSorted[24] = {
2052:       0, 1, 2, 3,  /* bottom */
2053:       4, 5, 6, 7,  /* top */
2054:       0, 3, 4, 5,  /* front */
2055:       2, 3, 5, 6,  /* right */
2056:       1, 2, 6, 7,  /* back */
2057:       0, 1, 4, 7,  /* left */
2058:     };
2059:     PetscInt       faceVerticesHex[24] = {
2060:       1, 2, 3, 0,  /* bottom */
2061:       4, 5, 6, 7,  /* top */
2062:       0, 3, 5, 4,  /* front */
2063:       3, 2, 6, 5,  /* right */
2064:       2, 1, 7, 6,  /* back */
2065:       1, 0, 4, 7,  /* left */
2066:     };

2068:     faceSize = faceSizeHex;
2069:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2070:     PetscSortInt(faceSizeHex, sortedIndices);
2071:     for (iFace = 0; iFace < 6; ++iFace) {
2072:       const PetscInt ii = iFace*faceSizeHex;
2073:       PetscInt       fVertex, cVertex;

2075:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2076:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2077:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2078:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2079:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2080:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2081:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2082:               faceVertices[fVertex] = origVertices[cVertex];
2083:               break;
2084:             }
2085:           }
2086:         }
2087:         found = PETSC_TRUE;
2088:         break;
2089:       }
2090:     }
2091:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2092:     if (posOriented) *posOriented = PETSC_TRUE;
2093:     return(0);
2094:   } else if (cellDim == 3 && numCorners == 10) {
2095:     /* Quadratic tet */
2096:     /* Faces are determined by the first 3 vertices (corners of faces) */
2097:     const PetscInt faceSizeTet = 6;
2098:     PetscInt       sortedIndices[6], i, iFace;
2099:     PetscBool      found                     = PETSC_FALSE;
2100:     PetscInt       faceVerticesTetSorted[24] = {
2101:       0, 1, 2,  6, 7, 8, /* bottom */
2102:       0, 3, 4,  6, 7, 9,  /* front */
2103:       1, 4, 5,  7, 8, 9,  /* right */
2104:       2, 3, 5,  6, 8, 9,  /* left */
2105:     };
2106:     PetscInt       faceVerticesTet[24] = {
2107:       0, 1, 2,  6, 7, 8, /* bottom */
2108:       0, 4, 3,  6, 7, 9,  /* front */
2109:       1, 5, 4,  7, 8, 9,  /* right */
2110:       2, 3, 5,  8, 6, 9,  /* left */
2111:     };

2113:     faceSize = faceSizeTet;
2114:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2115:     PetscSortInt(faceSizeTet, sortedIndices);
2116:     for (iFace=0; iFace < 4; ++iFace) {
2117:       const PetscInt ii = iFace*faceSizeTet;
2118:       PetscInt       fVertex, cVertex;

2120:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2121:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2122:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2123:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2124:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2125:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2126:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2127:               faceVertices[fVertex] = origVertices[cVertex];
2128:               break;
2129:             }
2130:           }
2131:         }
2132:         found = PETSC_TRUE;
2133:         break;
2134:       }
2135:     }
2136:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2137:     if (posOriented) *posOriented = PETSC_TRUE;
2138:     return(0);
2139:   } else if (cellDim == 3 && numCorners == 27) {
2140:     /* Quadratic hexes (I hate this)
2141:        A hex is two oriented quads with the normal of the first
2142:        pointing up at the second.

2144:          7---6
2145:         /|  /|
2146:        4---5 |
2147:        | 3-|-2
2148:        |/  |/
2149:        0---1

2151:        Faces are determined by the first 4 vertices (corners of faces) */
2152:     const PetscInt faceSizeQuadHex = 9;
2153:     PetscInt       sortedIndices[9], i, iFace;
2154:     PetscBool      found                         = PETSC_FALSE;
2155:     PetscInt       faceVerticesQuadHexSorted[54] = {
2156:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2157:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2158:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2159:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2160:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2161:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2162:     };
2163:     PetscInt       faceVerticesQuadHex[54] = {
2164:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2165:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2166:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2167:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2168:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2169:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2170:     };

2172:     faceSize = faceSizeQuadHex;
2173:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2174:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2175:     for (iFace = 0; iFace < 6; ++iFace) {
2176:       const PetscInt ii = iFace*faceSizeQuadHex;
2177:       PetscInt       fVertex, cVertex;

2179:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2180:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2181:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2182:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2183:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2184:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2185:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2186:               faceVertices[fVertex] = origVertices[cVertex];
2187:               break;
2188:             }
2189:           }
2190:         }
2191:         found = PETSC_TRUE;
2192:         break;
2193:       }
2194:     }
2195:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2196:     if (posOriented) *posOriented = PETSC_TRUE;
2197:     return(0);
2198:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2199:   if (!posOrient) {
2200:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2201:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2202:   } else {
2203:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2204:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2205:   }
2206:   if (posOriented) *posOriented = posOrient;
2207:   return(0);
2208: }

2212: /*
2213:     Given a cell and a face, as a set of vertices,
2214:       return the oriented face, as a set of vertices, in faceVertices
2215:     The orientation is such that the face normal points out of the cell
2216: */
2217: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2218: {
2219:   const PetscInt *cone = NULL;
2220:   PetscInt        coneSize, v, f, v2;
2221:   PetscInt        oppositeVertex = -1;
2222:   PetscErrorCode  ierr;

2225:   DMPlexGetConeSize(dm, cell, &coneSize);
2226:   DMPlexGetCone(dm, cell, &cone);
2227:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2228:     PetscBool found = PETSC_FALSE;

2230:     for (f = 0; f < faceSize; ++f) {
2231:       if (face[f] == cone[v]) {
2232:         found = PETSC_TRUE; break;
2233:       }
2234:     }
2235:     if (found) {
2236:       indices[v2]      = v;
2237:       origVertices[v2] = cone[v];
2238:       ++v2;
2239:     } else {
2240:       oppositeVertex = v;
2241:     }
2242:   }
2243:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2244:   return(0);
2245: }

2249: /*
2250:   DMPlexInsertFace_Internal - Puts a face into the mesh

2252:   Not collective

2254:   Input Parameters:
2255:   + dm              - The DMPlex
2256:   . numFaceVertex   - The number of vertices in the face
2257:   . faceVertices    - The vertices in the face for dm
2258:   . subfaceVertices - The vertices in the face for subdm
2259:   . numCorners      - The number of vertices in the cell
2260:   . cell            - A cell in dm containing the face
2261:   . subcell         - A cell in subdm containing the face
2262:   . firstFace       - First face in the mesh
2263:   - newFacePoint    - Next face in the mesh

2265:   Output Parameters:
2266:   . newFacePoint - Contains next face point number on input, updated on output

2268:   Level: developer
2269: */
2270: static PetscErrorCode DMPlexInsertFace_Internal(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint)
2271: {
2272:   MPI_Comm        comm;
2273:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2274:   const PetscInt *faces;
2275:   PetscInt        numFaces, coneSize;
2276:   PetscErrorCode  ierr;

2279:   PetscObjectGetComm((PetscObject)dm,&comm);
2280:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2281:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2282: #if 0
2283:   /* Cannot use this because support() has not been constructed yet */
2284:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2285: #else
2286:   {
2287:     PetscInt f;

2289:     numFaces = 0;
2290:     DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2291:     for (f = firstFace; f < *newFacePoint; ++f) {
2292:       PetscInt dof, off, d;

2294:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2295:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2296:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2297:       for (d = 0; d < dof; ++d) {
2298:         const PetscInt p = submesh->cones[off+d];
2299:         PetscInt       v;

2301:         for (v = 0; v < numFaceVertices; ++v) {
2302:           if (subfaceVertices[v] == p) break;
2303:         }
2304:         if (v == numFaceVertices) break;
2305:       }
2306:       if (d == dof) {
2307:         numFaces               = 1;
2308:         ((PetscInt*) faces)[0] = f;
2309:       }
2310:     }
2311:   }
2312: #endif
2313:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2314:   else if (numFaces == 1) {
2315:     /* Add the other cell neighbor for this face */
2316:     DMPlexSetCone(subdm, subcell, faces);
2317:   } else {
2318:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2319:     PetscBool posOriented;

2321:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2322:     origVertices        = &orientedVertices[numFaceVertices];
2323:     indices             = &orientedVertices[numFaceVertices*2];
2324:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2325:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2326:     /* TODO: I know that routine should return a permutation, not the indices */
2327:     for (v = 0; v < numFaceVertices; ++v) {
2328:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2329:       for (ov = 0; ov < numFaceVertices; ++ov) {
2330:         if (orientedVertices[ov] == vertex) {
2331:           orientedSubVertices[ov] = subvertex;
2332:           break;
2333:         }
2334:       }
2335:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2336:     }
2337:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2338:     DMPlexSetCone(subdm, subcell, newFacePoint);
2339:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2340:     ++(*newFacePoint);
2341:   }
2342: #if 0
2343:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2344: #else
2345:   DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2346: #endif
2347:   return(0);
2348: }

2352: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2353: {
2354:   MPI_Comm        comm;
2355:   DMLabel         subpointMap;
2356:   IS              subvertexIS,  subcellIS;
2357:   const PetscInt *subVertices, *subCells;
2358:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2359:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2360:   PetscInt        vStart, vEnd, c, f;
2361:   PetscErrorCode  ierr;

2364:   PetscObjectGetComm((PetscObject)dm,&comm);
2365:   /* Create subpointMap which marks the submesh */
2366:   DMLabelCreate("subpoint_map", &subpointMap);
2367:   DMPlexSetSubpointMap(subdm, subpointMap);
2368:   DMLabelDestroy(&subpointMap);
2369:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2370:   /* Setup chart */
2371:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2372:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2373:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2374:   DMPlexSetVTKCellHeight(subdm, 1);
2375:   /* Set cone sizes */
2376:   firstSubVertex = numSubCells;
2377:   firstSubFace   = numSubCells+numSubVertices;
2378:   newFacePoint   = firstSubFace;
2379:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2380:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2381:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2382:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2383:   for (c = 0; c < numSubCells; ++c) {
2384:     DMPlexSetConeSize(subdm, c, 1);
2385:   }
2386:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2387:     DMPlexSetConeSize(subdm, f, nFV);
2388:   }
2389:   DMSetUp(subdm);
2390:   /* Create face cones */
2391:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2392:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2393:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2394:   for (c = 0; c < numSubCells; ++c) {
2395:     const PetscInt cell    = subCells[c];
2396:     const PetscInt subcell = c;
2397:     PetscInt      *closure = NULL;
2398:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2400:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2401:     for (cl = 0; cl < closureSize*2; cl += 2) {
2402:       const PetscInt point = closure[cl];
2403:       PetscInt       subVertex;

2405:       if ((point >= vStart) && (point < vEnd)) {
2406:         ++numCorners;
2407:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2408:         if (subVertex >= 0) {
2409:           closure[faceSize] = point;
2410:           subface[faceSize] = firstSubVertex+subVertex;
2411:           ++faceSize;
2412:         }
2413:       }
2414:     }
2415:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2416:     if (faceSize == nFV) {
2417:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2418:     }
2419:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2420:   }
2421:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2422:   DMPlexSymmetrize(subdm);
2423:   DMPlexStratify(subdm);
2424:   /* Build coordinates */
2425:   {
2426:     PetscSection coordSection, subCoordSection;
2427:     Vec          coordinates, subCoordinates;
2428:     PetscScalar *coords, *subCoords;
2429:     PetscInt     numComp, coordSize, v;

2431:     DMGetCoordinateSection(dm, &coordSection);
2432:     DMGetCoordinatesLocal(dm, &coordinates);
2433:     DMGetCoordinateSection(subdm, &subCoordSection);
2434:     PetscSectionSetNumFields(subCoordSection, 1);
2435:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2436:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2437:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2438:     for (v = 0; v < numSubVertices; ++v) {
2439:       const PetscInt vertex    = subVertices[v];
2440:       const PetscInt subvertex = firstSubVertex+v;
2441:       PetscInt       dof;

2443:       PetscSectionGetDof(coordSection, vertex, &dof);
2444:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2445:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2446:     }
2447:     PetscSectionSetUp(subCoordSection);
2448:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2449:     VecCreate(comm, &subCoordinates);
2450:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2451:     VecSetType(subCoordinates,VECSTANDARD);
2452:     if (coordSize) {
2453:       VecGetArray(coordinates,    &coords);
2454:       VecGetArray(subCoordinates, &subCoords);
2455:       for (v = 0; v < numSubVertices; ++v) {
2456:         const PetscInt vertex    = subVertices[v];
2457:         const PetscInt subvertex = firstSubVertex+v;
2458:         PetscInt       dof, off, sdof, soff, d;

2460:         PetscSectionGetDof(coordSection, vertex, &dof);
2461:         PetscSectionGetOffset(coordSection, vertex, &off);
2462:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2463:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2464:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2465:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2466:       }
2467:       VecRestoreArray(coordinates,    &coords);
2468:       VecRestoreArray(subCoordinates, &subCoords);
2469:     }
2470:     DMSetCoordinatesLocal(subdm, subCoordinates);
2471:     VecDestroy(&subCoordinates);
2472:   }
2473:   /* Cleanup */
2474:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2475:   ISDestroy(&subvertexIS);
2476:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2477:   ISDestroy(&subcellIS);
2478:   return(0);
2479: }

2483: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2484: {
2485:   MPI_Comm         comm;
2486:   DMLabel          subpointMap;
2487:   IS              *subpointIS;
2488:   const PetscInt **subpoints;
2489:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *coneONew;
2490:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2491:   PetscErrorCode   ierr;

2494:   PetscObjectGetComm((PetscObject)dm,&comm);
2495:   /* Create subpointMap which marks the submesh */
2496:   DMLabelCreate("subpoint_map", &subpointMap);
2497:   DMPlexSetSubpointMap(subdm, subpointMap);
2498:   DMLabelDestroy(&subpointMap);
2499:   if (vertexLabel) {DMPlexMarkSubmesh_Interpolated(dm, vertexLabel, value, subpointMap, subdm);}
2500:   /* Setup chart */
2501:   DMPlexGetDimension(dm, &dim);
2502:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2503:   for (d = 0; d <= dim; ++d) {
2504:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2505:     totSubPoints += numSubPoints[d];
2506:   }
2507:   DMPlexSetChart(subdm, 0, totSubPoints);
2508:   DMPlexSetVTKCellHeight(subdm, 1);
2509:   /* Set cone sizes */
2510:   firstSubPoint[dim] = 0;
2511:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2512:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2513:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2514:   for (d = 0; d <= dim; ++d) {
2515:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2516:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2517:   }
2518:   for (d = 0; d <= dim; ++d) {
2519:     for (p = 0; p < numSubPoints[d]; ++p) {
2520:       const PetscInt  point    = subpoints[d][p];
2521:       const PetscInt  subpoint = firstSubPoint[d] + p;
2522:       const PetscInt *cone;
2523:       PetscInt        coneSize, coneSizeNew, c, val;

2525:       DMPlexGetConeSize(dm, point, &coneSize);
2526:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2527:       if (d == dim) {
2528:         DMPlexGetCone(dm, point, &cone);
2529:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2530:           DMLabelGetValue(subpointMap, cone[c], &val);
2531:           if (val >= 0) coneSizeNew++;
2532:         }
2533:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2534:       }
2535:     }
2536:   }
2537:   DMSetUp(subdm);
2538:   /* Set cones */
2539:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2540:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&coneONew);
2541:   for (d = 0; d <= dim; ++d) {
2542:     for (p = 0; p < numSubPoints[d]; ++p) {
2543:       const PetscInt  point    = subpoints[d][p];
2544:       const PetscInt  subpoint = firstSubPoint[d] + p;
2545:       const PetscInt *cone, *ornt;
2546:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

2548:       DMPlexGetConeSize(dm, point, &coneSize);
2549:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2550:       DMPlexGetCone(dm, point, &cone);
2551:       DMPlexGetConeOrientation(dm, point, &ornt);
2552:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2553:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2554:         if (subc >= 0) {
2555:           coneNew[coneSizeNew]  = firstSubPoint[d-1] + subc;
2556:           coneONew[coneSizeNew] = ornt[c];
2557:           ++coneSizeNew;
2558:         }
2559:       }
2560:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2561:       DMPlexSetCone(subdm, subpoint, coneNew);
2562:       DMPlexSetConeOrientation(subdm, subpoint, coneONew);
2563:     }
2564:   }
2565:   PetscFree2(coneNew,coneONew);
2566:   DMPlexSymmetrize(subdm);
2567:   DMPlexStratify(subdm);
2568:   /* Build coordinates */
2569:   {
2570:     PetscSection coordSection, subCoordSection;
2571:     Vec          coordinates, subCoordinates;
2572:     PetscScalar *coords, *subCoords;
2573:     PetscInt     numComp, coordSize;

2575:     DMGetCoordinateSection(dm, &coordSection);
2576:     DMGetCoordinatesLocal(dm, &coordinates);
2577:     DMGetCoordinateSection(subdm, &subCoordSection);
2578:     PetscSectionSetNumFields(subCoordSection, 1);
2579:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2580:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2581:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2582:     for (v = 0; v < numSubPoints[0]; ++v) {
2583:       const PetscInt vertex    = subpoints[0][v];
2584:       const PetscInt subvertex = firstSubPoint[0]+v;
2585:       PetscInt       dof;

2587:       PetscSectionGetDof(coordSection, vertex, &dof);
2588:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2589:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2590:     }
2591:     PetscSectionSetUp(subCoordSection);
2592:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2593:     VecCreate(comm, &subCoordinates);
2594:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2595:     VecSetType(subCoordinates,VECSTANDARD);
2596:     VecGetArray(coordinates,    &coords);
2597:     VecGetArray(subCoordinates, &subCoords);
2598:     for (v = 0; v < numSubPoints[0]; ++v) {
2599:       const PetscInt vertex    = subpoints[0][v];
2600:       const PetscInt subvertex = firstSubPoint[0]+v;
2601:       PetscInt dof, off, sdof, soff, d;

2603:       PetscSectionGetDof(coordSection, vertex, &dof);
2604:       PetscSectionGetOffset(coordSection, vertex, &off);
2605:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2606:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2607:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2608:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2609:     }
2610:     VecRestoreArray(coordinates,    &coords);
2611:     VecRestoreArray(subCoordinates, &subCoords);
2612:     DMSetCoordinatesLocal(subdm, subCoordinates);
2613:     VecDestroy(&subCoordinates);
2614:   }
2615:   /* Cleanup */
2616:   for (d = 0; d <= dim; ++d) {
2617:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
2618:     ISDestroy(&subpointIS[d]);
2619:   }
2620:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2621:   return(0);
2622: }

2626: /*@
2627:   DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label

2629:   Input Parameters:
2630: + dm           - The original mesh
2631: . vertexLabel  - The DMLabel marking vertices contained in the surface
2632: - value        - The label value to use

2634:   Output Parameter:
2635: . subdm - The surface mesh

2637:   Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().

2639:   Level: developer

2641: .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue()
2642: @*/
2643: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
2644: {
2645:   PetscInt       dim, depth;

2651:   DMPlexGetDimension(dm, &dim);
2652:   DMPlexGetDepth(dm, &depth);
2653:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2654:   DMSetType(*subdm, DMPLEX);
2655:   DMPlexSetDimension(*subdm, dim-1);
2656:   if (depth == dim) {
2657:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
2658:   } else {
2659:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
2660:   }
2661:   return(0);
2662: }

2666: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2667: {
2668:   PetscInt       subPoint;

2671:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2672:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2673: }

2677: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
2678: {
2679:   MPI_Comm        comm;
2680:   DMLabel         subpointMap;
2681:   IS              subvertexIS;
2682:   const PetscInt *subVertices;
2683:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
2684:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2685:   PetscInt        cMax, c, f;
2686:   PetscErrorCode  ierr;

2689:   PetscObjectGetComm((PetscObject)dm, &comm);
2690:   /* Create subpointMap which marks the submesh */
2691:   DMLabelCreate("subpoint_map", &subpointMap);
2692:   DMPlexSetSubpointMap(subdm, subpointMap);
2693:   DMLabelDestroy(&subpointMap);
2694:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
2695:   /* Setup chart */
2696:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2697:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2698:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2699:   DMPlexSetVTKCellHeight(subdm, 1);
2700:   /* Set cone sizes */
2701:   firstSubVertex = numSubCells;
2702:   firstSubFace   = numSubCells+numSubVertices;
2703:   newFacePoint   = firstSubFace;
2704:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2705:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2706:   for (c = 0; c < numSubCells; ++c) {
2707:     DMPlexSetConeSize(subdm, c, 1);
2708:   }
2709:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2710:     DMPlexSetConeSize(subdm, f, nFV);
2711:   }
2712:   DMSetUp(subdm);
2713:   /* Create face cones */
2714:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2715:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2716:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2717:   for (c = 0; c < numSubCells; ++c) {
2718:     const PetscInt  cell    = subCells[c];
2719:     const PetscInt  subcell = c;
2720:     const PetscInt *cone, *cells;
2721:     PetscInt        numCells, subVertex, p, v;

2723:     if (cell < cMax) continue;
2724:     DMPlexGetCone(dm, cell, &cone);
2725:     for (v = 0; v < nFV; ++v) {
2726:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
2727:       subface[v] = firstSubVertex+subVertex;
2728:     }
2729:     DMPlexSetCone(subdm, newFacePoint, subface);
2730:     DMPlexSetCone(subdm, subcell, &newFacePoint);
2731:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
2732:     /* Not true in parallel
2733:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2734:     for (p = 0; p < numCells; ++p) {
2735:       PetscInt negsubcell;

2737:       if (cells[p] >= cMax) continue;
2738:       /* I know this is a crap search */
2739:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2740:         if (subCells[negsubcell] == cells[p]) break;
2741:       }
2742:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2743:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
2744:     }
2745:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
2746:     ++newFacePoint;
2747:   }
2748:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2749:   DMPlexSymmetrize(subdm);
2750:   DMPlexStratify(subdm);
2751:   /* Build coordinates */
2752:   {
2753:     PetscSection coordSection, subCoordSection;
2754:     Vec          coordinates, subCoordinates;
2755:     PetscScalar *coords, *subCoords;
2756:     PetscInt     numComp, coordSize, v;

2758:     DMGetCoordinateSection(dm, &coordSection);
2759:     DMGetCoordinatesLocal(dm, &coordinates);
2760:     DMGetCoordinateSection(subdm, &subCoordSection);
2761:     PetscSectionSetNumFields(subCoordSection, 1);
2762:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2763:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2764:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2765:     for (v = 0; v < numSubVertices; ++v) {
2766:       const PetscInt vertex    = subVertices[v];
2767:       const PetscInt subvertex = firstSubVertex+v;
2768:       PetscInt       dof;

2770:       PetscSectionGetDof(coordSection, vertex, &dof);
2771:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2772:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2773:     }
2774:     PetscSectionSetUp(subCoordSection);
2775:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2776:     VecCreate(comm, &subCoordinates);
2777:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2778:     VecSetType(subCoordinates,VECSTANDARD);
2779:     VecGetArray(coordinates,    &coords);
2780:     VecGetArray(subCoordinates, &subCoords);
2781:     for (v = 0; v < numSubVertices; ++v) {
2782:       const PetscInt vertex    = subVertices[v];
2783:       const PetscInt subvertex = firstSubVertex+v;
2784:       PetscInt       dof, off, sdof, soff, d;

2786:       PetscSectionGetDof(coordSection, vertex, &dof);
2787:       PetscSectionGetOffset(coordSection, vertex, &off);
2788:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2789:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2790:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2791:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2792:     }
2793:     VecRestoreArray(coordinates,    &coords);
2794:     VecRestoreArray(subCoordinates, &subCoords);
2795:     DMSetCoordinatesLocal(subdm, subCoordinates);
2796:     VecDestroy(&subCoordinates);
2797:   }
2798:   /* Build SF */
2799:   CHKMEMQ;
2800:   {
2801:     PetscSF            sfPoint, sfPointSub;
2802:     const PetscSFNode *remotePoints;
2803:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
2804:     const PetscInt    *localPoints;
2805:     PetscInt          *slocalPoints;
2806:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
2807:     PetscMPIInt        rank;

2809:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2810:     DMGetPointSF(dm, &sfPoint);
2811:     DMGetPointSF(subdm, &sfPointSub);
2812:     DMPlexGetChart(dm, &pStart, &pEnd);
2813:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2814:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2815:     if (numRoots >= 0) {
2816:       /* Only vertices should be shared */
2817:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2818:       for (p = 0; p < pEnd-pStart; ++p) {
2819:         newLocalPoints[p].rank  = -2;
2820:         newLocalPoints[p].index = -2;
2821:       }
2822:       /* Set subleaves */
2823:       for (l = 0; l < numLeaves; ++l) {
2824:         const PetscInt point    = localPoints[l];
2825:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2827:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
2828:         if (subPoint < 0) continue;
2829:         newLocalPoints[point-pStart].rank  = rank;
2830:         newLocalPoints[point-pStart].index = subPoint;
2831:         ++numSubLeaves;
2832:       }
2833:       /* Must put in owned subpoints */
2834:       for (p = pStart; p < pEnd; ++p) {
2835:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

2837:         if (subPoint < 0) {
2838:           newOwners[p-pStart].rank  = -3;
2839:           newOwners[p-pStart].index = -3;
2840:         } else {
2841:           newOwners[p-pStart].rank  = rank;
2842:           newOwners[p-pStart].index = subPoint;
2843:         }
2844:       }
2845:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2846:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2847:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2848:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2849:       PetscMalloc1(numSubLeaves,    &slocalPoints);
2850:       PetscMalloc1(numSubLeaves, &sremotePoints);
2851:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
2852:         const PetscInt point    = localPoints[l];
2853:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2855:         if (subPoint < 0) continue;
2856:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
2857:         slocalPoints[sl]        = subPoint;
2858:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
2859:         sremotePoints[sl].index = newLocalPoints[point].index;
2860:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
2861:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
2862:         ++sl;
2863:       }
2864:       PetscFree2(newLocalPoints,newOwners);
2865:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
2866:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
2867:     }
2868:   }
2869:   CHKMEMQ;
2870:   /* Cleanup */
2871:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2872:   ISDestroy(&subvertexIS);
2873:   PetscFree(subCells);
2874:   return(0);
2875: }

2879: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char label[], PetscInt value, DM subdm)
2880: {
2881:   MPI_Comm         comm;
2882:   DMLabel          subpointMap;
2883:   IS              *subpointIS;
2884:   const PetscInt **subpoints;
2885:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2886:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2887:   PetscErrorCode   ierr;

2890:   PetscObjectGetComm((PetscObject)dm,&comm);
2891:   /* Create subpointMap which marks the submesh */
2892:   DMLabelCreate("subpoint_map", &subpointMap);
2893:   DMPlexSetSubpointMap(subdm, subpointMap);
2894:   DMLabelDestroy(&subpointMap);
2895:   DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
2896:   /* Setup chart */
2897:   DMPlexGetDimension(dm, &dim);
2898:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2899:   for (d = 0; d <= dim; ++d) {
2900:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2901:     totSubPoints += numSubPoints[d];
2902:   }
2903:   DMPlexSetChart(subdm, 0, totSubPoints);
2904:   DMPlexSetVTKCellHeight(subdm, 1);
2905:   /* Set cone sizes */
2906:   firstSubPoint[dim] = 0;
2907:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2908:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2909:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2910:   for (d = 0; d <= dim; ++d) {
2911:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2912:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2913:   }
2914:   for (d = 0; d <= dim; ++d) {
2915:     for (p = 0; p < numSubPoints[d]; ++p) {
2916:       const PetscInt  point    = subpoints[d][p];
2917:       const PetscInt  subpoint = firstSubPoint[d] + p;
2918:       const PetscInt *cone;
2919:       PetscInt        coneSize, coneSizeNew, c, val;

2921:       DMPlexGetConeSize(dm, point, &coneSize);
2922:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2923:       if (d == dim) {
2924:         DMPlexGetCone(dm, point, &cone);
2925:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2926:           DMLabelGetValue(subpointMap, cone[c], &val);
2927:           if (val >= 0) coneSizeNew++;
2928:         }
2929:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2930:       }
2931:     }
2932:   }
2933:   DMSetUp(subdm);
2934:   /* Set cones */
2935:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2936:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2937:   for (d = 0; d <= dim; ++d) {
2938:     for (p = 0; p < numSubPoints[d]; ++p) {
2939:       const PetscInt  point    = subpoints[d][p];
2940:       const PetscInt  subpoint = firstSubPoint[d] + p;
2941:       const PetscInt *cone, *ornt;
2942:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

2944:       DMPlexGetConeSize(dm, point, &coneSize);
2945:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2946:       DMPlexGetCone(dm, point, &cone);
2947:       DMPlexGetConeOrientation(dm, point, &ornt);
2948:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2949:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2950:         if (subc >= 0) {
2951:           coneNew[coneSizeNew]   = firstSubPoint[d-1] + subc;
2952:           orntNew[coneSizeNew++] = ornt[c];
2953:         }
2954:       }
2955:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2956:       DMPlexSetCone(subdm, subpoint, coneNew);
2957:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2958:     }
2959:   }
2960:   PetscFree2(coneNew,orntNew);
2961:   DMPlexSymmetrize(subdm);
2962:   DMPlexStratify(subdm);
2963:   /* Build coordinates */
2964:   {
2965:     PetscSection coordSection, subCoordSection;
2966:     Vec          coordinates, subCoordinates;
2967:     PetscScalar *coords, *subCoords;
2968:     PetscInt     numComp, coordSize;

2970:     DMGetCoordinateSection(dm, &coordSection);
2971:     DMGetCoordinatesLocal(dm, &coordinates);
2972:     DMGetCoordinateSection(subdm, &subCoordSection);
2973:     PetscSectionSetNumFields(subCoordSection, 1);
2974:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2975:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2976:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2977:     for (v = 0; v < numSubPoints[0]; ++v) {
2978:       const PetscInt vertex    = subpoints[0][v];
2979:       const PetscInt subvertex = firstSubPoint[0]+v;
2980:       PetscInt       dof;

2982:       PetscSectionGetDof(coordSection, vertex, &dof);
2983:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2984:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2985:     }
2986:     PetscSectionSetUp(subCoordSection);
2987:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2988:     VecCreate(comm, &subCoordinates);
2989:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2990:     VecSetType(subCoordinates,VECSTANDARD);
2991:     VecGetArray(coordinates,    &coords);
2992:     VecGetArray(subCoordinates, &subCoords);
2993:     for (v = 0; v < numSubPoints[0]; ++v) {
2994:       const PetscInt vertex    = subpoints[0][v];
2995:       const PetscInt subvertex = firstSubPoint[0]+v;
2996:       PetscInt dof, off, sdof, soff, d;

2998:       PetscSectionGetDof(coordSection, vertex, &dof);
2999:       PetscSectionGetOffset(coordSection, vertex, &off);
3000:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3001:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3002:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3003:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3004:     }
3005:     VecRestoreArray(coordinates,    &coords);
3006:     VecRestoreArray(subCoordinates, &subCoords);
3007:     DMSetCoordinatesLocal(subdm, subCoordinates);
3008:     VecDestroy(&subCoordinates);
3009:   }
3010:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3011:   {
3012:     PetscSF            sfPoint, sfPointSub;
3013:     IS                 subpIS;
3014:     const PetscSFNode *remotePoints;
3015:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3016:     const PetscInt    *localPoints, *subpoints;
3017:     PetscInt          *slocalPoints;
3018:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3019:     PetscMPIInt        rank;

3021:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3022:     DMGetPointSF(dm, &sfPoint);
3023:     DMGetPointSF(subdm, &sfPointSub);
3024:     DMPlexGetChart(dm, &pStart, &pEnd);
3025:     DMPlexGetChart(subdm, NULL, &numSubroots);
3026:     DMPlexCreateSubpointIS(subdm, &subpIS);
3027:     if (subpIS) {
3028:       ISGetIndices(subpIS, &subpoints);
3029:       ISGetLocalSize(subpIS, &numSubpoints);
3030:     }
3031:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3032:     if (numRoots >= 0) {
3033:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3034:       for (p = 0; p < pEnd-pStart; ++p) {
3035:         newLocalPoints[p].rank  = -2;
3036:         newLocalPoints[p].index = -2;
3037:       }
3038:       /* Set subleaves */
3039:       for (l = 0; l < numLeaves; ++l) {
3040:         const PetscInt point    = localPoints[l];
3041:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3043:         if (subpoint < 0) continue;
3044:         newLocalPoints[point-pStart].rank  = rank;
3045:         newLocalPoints[point-pStart].index = subpoint;
3046:         ++numSubleaves;
3047:       }
3048:       /* Must put in owned subpoints */
3049:       for (p = pStart; p < pEnd; ++p) {
3050:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3052:         if (subpoint < 0) {
3053:           newOwners[p-pStart].rank  = -3;
3054:           newOwners[p-pStart].index = -3;
3055:         } else {
3056:           newOwners[p-pStart].rank  = rank;
3057:           newOwners[p-pStart].index = subpoint;
3058:         }
3059:       }
3060:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3061:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3062:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3063:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3064:       PetscMalloc(numSubleaves * sizeof(PetscInt),    &slocalPoints);
3065:       PetscMalloc(numSubleaves * sizeof(PetscSFNode), &sremotePoints);
3066:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3067:         const PetscInt point    = localPoints[l];
3068:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3070:         if (subpoint < 0) continue;
3071:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3072:         slocalPoints[sl]        = subpoint;
3073:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3074:         sremotePoints[sl].index = newLocalPoints[point].index;
3075:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3076:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3077:         ++sl;
3078:       }
3079:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3080:       PetscFree2(newLocalPoints,newOwners);
3081:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3082:     }
3083:     if (subpIS) {
3084:       ISRestoreIndices(subpIS, &subpoints);
3085:       ISDestroy(&subpIS);
3086:     }
3087:   }
3088:   /* Cleanup */
3089:   for (d = 0; d <= dim; ++d) {
3090:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3091:     ISDestroy(&subpointIS[d]);
3092:   }
3093:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3094:   return(0);
3095: }

3099: /*
3100:   DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label an be given to restrict the cells.

3102:   Input Parameters:
3103: + dm          - The original mesh
3104: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3105: . label       - A label name, or NULL
3106: - value  - A label value

3108:   Output Parameter:
3109: . subdm - The surface mesh

3111:   Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().

3113:   Level: developer

3115: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3116: */
3117: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3118: {
3119:   PetscInt       dim, depth;

3125:   DMPlexGetDimension(dm, &dim);
3126:   DMPlexGetDepth(dm, &depth);
3127:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3128:   DMSetType(*subdm, DMPLEX);
3129:   DMPlexSetDimension(*subdm, dim-1);
3130:   if (depth == dim) {
3131:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3132:   } else {
3133:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3134:   }
3135:   return(0);
3136: }

3140: /*@
3141:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3143:   Input Parameter:
3144: . dm - The submesh DM

3146:   Output Parameter:
3147: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh

3149:   Level: developer

3151: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3152: @*/
3153: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3154: {
3155:   DM_Plex *mesh = (DM_Plex*) dm->data;

3160:   *subpointMap = mesh->subpointMap;
3161:   return(0);
3162: }

3166: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3167: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3168: {
3169:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3170:   DMLabel        tmp;

3175:   tmp  = mesh->subpointMap;
3176:   mesh->subpointMap = subpointMap;
3177:   ++mesh->subpointMap->refct;
3178:   DMLabelDestroy(&tmp);
3179:   return(0);
3180: }

3184: /*@
3185:   DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data

3187:   Input Parameter:
3188: . dm - The submesh DM

3190:   Output Parameter:
3191: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh

3193:   Note: This is IS is guaranteed to be sorted by the construction of the submesh

3195:   Level: developer

3197: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3198: @*/
3199: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3200: {
3201:   MPI_Comm        comm;
3202:   DMLabel         subpointMap;
3203:   IS              is;
3204:   const PetscInt *opoints;
3205:   PetscInt       *points, *depths;
3206:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3207:   PetscErrorCode  ierr;

3212:   PetscObjectGetComm((PetscObject)dm,&comm);
3213:   *subpointIS = NULL;
3214:   DMPlexGetSubpointMap(dm, &subpointMap);
3215:   DMPlexGetDepth(dm, &depth);
3216:   if (subpointMap && depth >= 0) {
3217:     DMPlexGetChart(dm, &pStart, &pEnd);
3218:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3219:     DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3220:     depths[0] = depth;
3221:     depths[1] = 0;
3222:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3223:     PetscMalloc1(pEnd, &points);
3224:     for(d = 0, off = 0; d <= depth; ++d) {
3225:       const PetscInt dep = depths[d];

3227:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3228:       DMLabelGetStratumSize(subpointMap, dep, &n);
3229:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3230:         if (n != depEnd-depStart) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
3231:       } else {
3232:         if (!n) {
3233:           if (d == 0) {
3234:             /* Missing cells */
3235:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3236:           } else {
3237:             /* Missing faces */
3238:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3239:           }
3240:         }
3241:       }
3242:       if (n) {
3243:         DMLabelGetStratumIS(subpointMap, dep, &is);
3244:         ISGetIndices(is, &opoints);
3245:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3246:         ISRestoreIndices(is, &opoints);
3247:         ISDestroy(&is);
3248:       }
3249:     }
3250:     DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3251:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3252:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3253:   }
3254:   return(0);
3255: }