Actual source code: plexsubmesh.c

petsc-master 2017-07-20
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petscsf.h>

  5: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label)
  6: {
  7:   PetscInt       fStart, fEnd, f;

 12:   DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
 13:   for (f = fStart; f < fEnd; ++f) {
 14:     PetscInt supportSize;

 16:     DMPlexGetSupportSize(dm, f, &supportSize);
 17:     if (supportSize == 1) {DMLabelSetValue(label, f, 1);}
 18:   }
 19:   return(0);
 20: }

 22: /*@
 23:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 25:   Not Collective

 27:   Input Parameter:
 28: . dm - The original DM

 30:   Output Parameter:
 31: . label - The DMLabel marking boundary faces with value 1

 33:   Level: developer

 35: .seealso: DMLabelCreate(), DMCreateLabel()
 36: @*/
 37: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
 38: {

 42:   DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
 43:   return(0);
 44: }

 46: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
 47: {
 48:   IS              valueIS;
 49:   const PetscInt *values;
 50:   PetscInt        numValues, v, cStart, cEnd;
 51:   PetscErrorCode  ierr;

 54:   DMLabelGetNumValues(label, &numValues);
 55:   DMLabelGetValueIS(label, &valueIS);
 56:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 57:   ISGetIndices(valueIS, &values);
 58:   for (v = 0; v < numValues; ++v) {
 59:     IS              pointIS;
 60:     const PetscInt *points;
 61:     PetscInt        numPoints, p;

 63:     DMLabelGetStratumSize(label, values[v], &numPoints);
 64:     DMLabelGetStratumIS(label, values[v], &pointIS);
 65:     ISGetIndices(pointIS, &points);
 66:     for (p = 0; p < numPoints; ++p) {
 67:       PetscInt  q = points[p];
 68:       PetscInt *closure = NULL;
 69:       PetscInt  closureSize, c;

 71:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
 72:         continue;
 73:       }
 74:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
 75:       for (c = 0; c < closureSize*2; c += 2) {
 76:         DMLabelSetValue(label, closure[c], values[v]);
 77:       }
 78:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
 79:     }
 80:     ISRestoreIndices(pointIS, &points);
 81:     ISDestroy(&pointIS);
 82:   }
 83:   ISRestoreIndices(valueIS, &values);
 84:   ISDestroy(&valueIS);
 85:   return(0);
 86: }

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

 91:   Input Parameters:
 92: + dm - The DM
 93: - label - A DMLabel marking the surface points

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

 98:   Level: developer

100: .seealso: DMPlexLabelCohesiveComplete()
101: @*/
102: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
103: {

107:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
108:   return(0);
109: }

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

114:   Input Parameters:
115: + dm - The DM
116: - label - A DMLabel marking the surface points

118:   Output Parameter:
119: . label - A DMLabel incorporating cells

121:   Level: developer

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

125: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
126: @*/
127: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
128: {
129:   IS              valueIS;
130:   const PetscInt *values;
131:   PetscInt        numValues, v, cStart, cEnd, cEndInterior;
132:   PetscErrorCode  ierr;

135:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
136:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
137:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
138:   DMLabelGetNumValues(label, &numValues);
139:   DMLabelGetValueIS(label, &valueIS);
140:   ISGetIndices(valueIS, &values);
141:   for (v = 0; v < numValues; ++v) {
142:     IS              pointIS;
143:     const PetscInt *points;
144:     PetscInt        numPoints, p;

146:     DMLabelGetStratumSize(label, values[v], &numPoints);
147:     DMLabelGetStratumIS(label, values[v], &pointIS);
148:     ISGetIndices(pointIS, &points);
149:     for (p = 0; p < numPoints; ++p) {
150:       PetscInt *closure = NULL;
151:       PetscInt  closureSize, point, cl;

153:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
154:       for (cl = closureSize-1; cl > 0; --cl) {
155:         point = closure[cl*2];
156:         if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
157:       }
158:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
159:     }
160:     ISRestoreIndices(pointIS, &points);
161:     ISDestroy(&pointIS);
162:   }
163:   ISRestoreIndices(valueIS, &values);
164:   ISDestroy(&valueIS);
165:   return(0);
166: }

168: /*@
169:   DMPlexLabelClearCells - Remove cells from a label

171:   Input Parameters:
172: + dm - The DM
173: - label - A DMLabel marking surface points and their adjacent cells

175:   Output Parameter:
176: . label - A DMLabel without cells

178:   Level: developer

180:   Note: This undoes DMPlexLabelAddCells()

182: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
183: @*/
184: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
185: {
186:   IS              valueIS;
187:   const PetscInt *values;
188:   PetscInt        numValues, v, cStart, cEnd, cEndInterior;
189:   PetscErrorCode  ierr;

192:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
193:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
194:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
195:   DMLabelGetNumValues(label, &numValues);
196:   DMLabelGetValueIS(label, &valueIS);
197:   ISGetIndices(valueIS, &values);
198:   for (v = 0; v < numValues; ++v) {
199:     IS              pointIS;
200:     const PetscInt *points;
201:     PetscInt        numPoints, p;

203:     DMLabelGetStratumSize(label, values[v], &numPoints);
204:     DMLabelGetStratumIS(label, values[v], &pointIS);
205:     ISGetIndices(pointIS, &points);
206:     for (p = 0; p < numPoints; ++p) {
207:       PetscInt point = points[p];

209:       if (point >= cStart && point < cEnd) {
210:         DMLabelClearValue(label,point,values[v]);
211:       }
212:     }
213:     ISRestoreIndices(pointIS, &points);
214:     ISDestroy(&pointIS);
215:   }
216:   ISRestoreIndices(valueIS, &values);
217:   ISDestroy(&valueIS);
218:   return(0);
219: }

221: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
222:  * index (skipping first, which is (0,0)) */
223: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
224: {
225:   PetscInt d, off = 0;

228:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
229:   for (d = 0; d < depth; d++) {
230:     PetscInt firstd = d;
231:     PetscInt firstStart = depthShift[2*d];
232:     PetscInt e;

234:     for (e = d+1; e <= depth; e++) {
235:       if (depthShift[2*e] < firstStart) {
236:         firstd = e;
237:         firstStart = depthShift[2*d];
238:       }
239:     }
240:     if (firstd != d) {
241:       PetscInt swap[2];

243:       e = firstd;
244:       swap[0] = depthShift[2*d];
245:       swap[1] = depthShift[2*d+1];
246:       depthShift[2*d]   = depthShift[2*e];
247:       depthShift[2*d+1] = depthShift[2*e+1];
248:       depthShift[2*e]   = swap[0];
249:       depthShift[2*e+1] = swap[1];
250:     }
251:   }
252:   /* convert (oldstart, added) to (oldstart, newstart) */
253:   for (d = 0; d <= depth; d++) {
254:     off += depthShift[2*d+1];
255:     depthShift[2*d+1] = depthShift[2*d] + off;
256:   }
257:   return(0);
258: }

260: /* depthShift is a list of (old, new) pairs */
261: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
262: {
263:   PetscInt d;
264:   PetscInt newOff = 0;

266:   for (d = 0; d <= depth; d++) {
267:     if (p < depthShift[2*d]) return p + newOff;
268:     else newOff = depthShift[2*d+1] - depthShift[2*d];
269:   }
270:   return p + newOff;
271: }

273: /* depthShift is a list of (old, new) pairs */
274: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
275: {
276:   PetscInt d;
277:   PetscInt newOff = 0;

279:   for (d = 0; d <= depth; d++) {
280:     if (p < depthShift[2*d+1]) return p + newOff;
281:     else newOff = depthShift[2*d] - depthShift[2*d+1];
282:   }
283:   return p + newOff;
284: }

286: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
287: {
288:   PetscInt       depth = 0, d, pStart, pEnd, p;

292:   DMPlexGetDepth(dm, &depth);
293:   if (depth < 0) return(0);
294:   /* Step 1: Expand chart */
295:   DMPlexGetChart(dm, &pStart, &pEnd);
296:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
297:   DMPlexSetChart(dmNew, pStart, pEnd);
298:   /* Step 2: Set cone and support sizes */
299:   for (d = 0; d <= depth; ++d) {
300:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
301:     for (p = pStart; p < pEnd; ++p) {
302:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
303:       PetscInt size;

305:       DMPlexGetConeSize(dm, p, &size);
306:       DMPlexSetConeSize(dmNew, newp, size);
307:       DMPlexGetSupportSize(dm, p, &size);
308:       DMPlexSetSupportSize(dmNew, newp, size);
309:     }
310:   }
311:   return(0);
312: }

314: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
315: {
316:   PetscInt      *newpoints;
317:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

321:   DMPlexGetDepth(dm, &depth);
322:   if (depth < 0) return(0);
323:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
324:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
325:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
326:   /* Step 5: Set cones and supports */
327:   DMPlexGetChart(dm, &pStart, &pEnd);
328:   for (p = pStart; p < pEnd; ++p) {
329:     const PetscInt *points = NULL, *orientations = NULL;
330:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

332:     DMPlexGetConeSize(dm, p, &size);
333:     DMPlexGetCone(dm, p, &points);
334:     DMPlexGetConeOrientation(dm, p, &orientations);
335:     for (i = 0; i < size; ++i) {
336:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
337:     }
338:     DMPlexSetCone(dmNew, newp, newpoints);
339:     DMPlexSetConeOrientation(dmNew, newp, orientations);
340:     DMPlexGetSupportSize(dm, p, &size);
341:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
342:     DMPlexGetSupport(dm, p, &points);
343:     for (i = 0; i < size; ++i) {
344:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
345:     }
346:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
347:     DMPlexSetSupport(dmNew, newp, newpoints);
348:   }
349:   PetscFree(newpoints);
350:   return(0);
351: }

353: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
354: {
355:   PetscSection   coordSection, newCoordSection;
356:   Vec            coordinates, newCoordinates;
357:   PetscScalar   *coords, *newCoords;
358:   PetscInt       coordSize, sStart, sEnd;
359:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
360:   PetscBool      hasCells;

364:   DMGetCoordinateDim(dm, &dim);
365:   DMSetCoordinateDim(dmNew, dim);
366:   DMPlexGetDepth(dm, &depth);
367:   /* Step 8: Convert coordinates */
368:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
369:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
370:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
371:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
372:   DMGetCoordinateSection(dm, &coordSection);
373:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
374:   PetscSectionSetNumFields(newCoordSection, 1);
375:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
376:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
377:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
378:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
379:   if (hasCells) {
380:     for (c = cStart; c < cEnd; ++c) {
381:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

383:       PetscSectionGetDof(coordSection, c, &dof);
384:       PetscSectionSetDof(newCoordSection, cNew, dof);
385:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
386:     }
387:   }
388:   for (v = vStartNew; v < vEndNew; ++v) {
389:     PetscSectionSetDof(newCoordSection, v, dim);
390:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
391:   }
392:   PetscSectionSetUp(newCoordSection);
393:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
394:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
395:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
396:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
397:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
398:   VecSetBlockSize(newCoordinates, dim);
399:   VecSetType(newCoordinates,VECSTANDARD);
400:   DMSetCoordinatesLocal(dmNew, newCoordinates);
401:   DMGetCoordinatesLocal(dm, &coordinates);
402:   VecGetArray(coordinates, &coords);
403:   VecGetArray(newCoordinates, &newCoords);
404:   if (hasCells) {
405:     for (c = cStart; c < cEnd; ++c) {
406:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

408:       PetscSectionGetDof(coordSection, c, &dof);
409:       PetscSectionGetOffset(coordSection, c, &off);
410:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
411:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
412:     }
413:   }
414:   for (v = vStart; v < vEnd; ++v) {
415:     PetscInt dof, off, noff, d;

417:     PetscSectionGetDof(coordSection, v, &dof);
418:     PetscSectionGetOffset(coordSection, v, &off);
419:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
420:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
421:   }
422:   VecRestoreArray(coordinates, &coords);
423:   VecRestoreArray(newCoordinates, &newCoords);
424:   VecDestroy(&newCoordinates);
425:   PetscSectionDestroy(&newCoordSection);
426:   return(0);
427: }

429: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
430: {
431:   PetscInt           depth = 0;
432:   PetscSF            sfPoint, sfPointNew;
433:   const PetscSFNode *remotePoints;
434:   PetscSFNode       *gremotePoints;
435:   const PetscInt    *localPoints;
436:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
437:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
438:   PetscErrorCode     ierr;

441:   DMPlexGetDepth(dm, &depth);
442:   /* Step 9: Convert pointSF */
443:   DMGetPointSF(dm, &sfPoint);
444:   DMGetPointSF(dmNew, &sfPointNew);
445:   DMPlexGetChart(dm, &pStart, &pEnd);
446:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
447:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
448:   if (numRoots >= 0) {
449:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
450:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
451:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
452:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
453:     PetscMalloc1(numLeaves,    &glocalPoints);
454:     PetscMalloc1(numLeaves, &gremotePoints);
455:     for (l = 0; l < numLeaves; ++l) {
456:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
457:       gremotePoints[l].rank  = remotePoints[l].rank;
458:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
459:     }
460:     PetscFree2(newLocation,newRemoteLocation);
461:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
462:   }
463:   return(0);
464: }

466: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
467: {
468:   PetscSF            sfPoint;
469:   DMLabel            vtkLabel, ghostLabel;
470:   const PetscSFNode *leafRemote;
471:   const PetscInt    *leafLocal;
472:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
473:   PetscMPIInt        rank;
474:   PetscErrorCode     ierr;

477:   DMPlexGetDepth(dm, &depth);
478:   /* Step 10: Convert labels */
479:   DMGetNumLabels(dm, &numLabels);
480:   for (l = 0; l < numLabels; ++l) {
481:     DMLabel         label, newlabel;
482:     const char     *lname;
483:     PetscBool       isDepth;
484:     IS              valueIS;
485:     const PetscInt *values;
486:     PetscInt        numValues, val;

488:     DMGetLabelName(dm, l, &lname);
489:     PetscStrcmp(lname, "depth", &isDepth);
490:     if (isDepth) continue;
491:     DMCreateLabel(dmNew, lname);
492:     DMGetLabel(dm, lname, &label);
493:     DMGetLabel(dmNew, lname, &newlabel);
494:     DMLabelGetDefaultValue(label,&val);
495:     DMLabelSetDefaultValue(newlabel,val);
496:     DMLabelGetValueIS(label, &valueIS);
497:     ISGetLocalSize(valueIS, &numValues);
498:     ISGetIndices(valueIS, &values);
499:     for (val = 0; val < numValues; ++val) {
500:       IS              pointIS;
501:       const PetscInt *points;
502:       PetscInt        numPoints, p;

504:       DMLabelGetStratumIS(label, values[val], &pointIS);
505:       ISGetLocalSize(pointIS, &numPoints);
506:       ISGetIndices(pointIS, &points);
507:       for (p = 0; p < numPoints; ++p) {
508:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

510:         DMLabelSetValue(newlabel, newpoint, values[val]);
511:       }
512:       ISRestoreIndices(pointIS, &points);
513:       ISDestroy(&pointIS);
514:     }
515:     ISRestoreIndices(valueIS, &values);
516:     ISDestroy(&valueIS);
517:   }
518:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
519:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
520:   DMGetPointSF(dm, &sfPoint);
521:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
522:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
523:   DMCreateLabel(dmNew, "vtk");
524:   DMCreateLabel(dmNew, "ghost");
525:   DMGetLabel(dmNew, "vtk", &vtkLabel);
526:   DMGetLabel(dmNew, "ghost", &ghostLabel);
527:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
528:     for (; c < leafLocal[l] && c < cEnd; ++c) {
529:       DMLabelSetValue(vtkLabel, c, 1);
530:     }
531:     if (leafLocal[l] >= cEnd) break;
532:     if (leafRemote[l].rank == rank) {
533:       DMLabelSetValue(vtkLabel, c, 1);
534:     } else {
535:       DMLabelSetValue(ghostLabel, c, 2);
536:     }
537:   }
538:   for (; c < cEnd; ++c) {
539:     DMLabelSetValue(vtkLabel, c, 1);
540:   }
541:   if (0) {
542:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
543:   }
544:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
545:   for (f = fStart; f < fEnd; ++f) {
546:     PetscInt numCells;

548:     DMPlexGetSupportSize(dmNew, f, &numCells);
549:     if (numCells < 2) {
550:       DMLabelSetValue(ghostLabel, f, 1);
551:     } else {
552:       const PetscInt *cells = NULL;
553:       PetscInt        vA, vB;

555:       DMPlexGetSupport(dmNew, f, &cells);
556:       DMLabelGetValue(vtkLabel, cells[0], &vA);
557:       DMLabelGetValue(vtkLabel, cells[1], &vB);
558:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
559:     }
560:   }
561:   if (0) {
562:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
563:   }
564:   return(0);
565: }

567: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
568: {
569:   DM             refTree;
570:   PetscSection   pSec;
571:   PetscInt       *parents, *childIDs;

575:   DMPlexGetReferenceTree(dm,&refTree);
576:   DMPlexSetReferenceTree(dmNew,refTree);
577:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
578:   if (pSec) {
579:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
580:     PetscInt *childIDsShifted;
581:     PetscSection pSecShifted;

583:     PetscSectionGetChart(pSec,&pStart,&pEnd);
584:     DMPlexGetDepth(dm,&depth);
585:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
586:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
587:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
588:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
589:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
590:     for (p = pStartShifted; p < pEndShifted; p++) {
591:       /* start off assuming no children */
592:       PetscSectionSetDof(pSecShifted,p,0);
593:     }
594:     for (p = pStart; p < pEnd; p++) {
595:       PetscInt dof;
596:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

598:       PetscSectionGetDof(pSec,p,&dof);
599:       PetscSectionSetDof(pSecShifted,pNew,dof);
600:     }
601:     PetscSectionSetUp(pSecShifted);
602:     for (p = pStart; p < pEnd; p++) {
603:       PetscInt dof;
604:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

606:       PetscSectionGetDof(pSec,p,&dof);
607:       if (dof) {
608:         PetscInt off, offNew;

610:         PetscSectionGetOffset(pSec,p,&off);
611:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
612:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
613:         childIDsShifted[offNew] = childIDs[off];
614:       }
615:     }
616:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
617:     PetscFree2(parentsShifted,childIDsShifted);
618:     PetscSectionDestroy(&pSecShifted);
619:   }
620:   return(0);
621: }

623: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
624: {
625:   PetscSF               sf;
626:   IS                    valueIS;
627:   const PetscInt       *values, *leaves;
628:   PetscInt             *depthShift;
629:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
630:   PetscBool             isper;
631:   const PetscReal      *maxCell, *L;
632:   const DMBoundaryType *bd;
633:   PetscErrorCode        ierr;

636:   DMGetPointSF(dm, &sf);
637:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
638:   nleaves = PetscMax(0, nleaves);
639:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
640:   /* Count ghost cells */
641:   DMLabelGetValueIS(label, &valueIS);
642:   ISGetLocalSize(valueIS, &numFS);
643:   ISGetIndices(valueIS, &values);
644:   Ng   = 0;
645:   for (fs = 0; fs < numFS; ++fs) {
646:     IS              faceIS;
647:     const PetscInt *faces;
648:     PetscInt        numFaces, f, numBdFaces = 0;

650:     DMLabelGetStratumIS(label, values[fs], &faceIS);
651:     ISGetLocalSize(faceIS, &numFaces);
652:     ISGetIndices(faceIS, &faces);
653:     for (f = 0; f < numFaces; ++f) {
654:       PetscInt numChildren;

656:       PetscFindInt(faces[f], nleaves, leaves, &loc);
657:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
658:       /* non-local and ancestors points don't get to register ghosts */
659:       if (loc >= 0 || numChildren) continue;
660:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
661:     }
662:     Ng += numBdFaces;
663:     ISDestroy(&faceIS);
664:   }
665:   DMPlexGetDepth(dm, &depth);
666:   PetscMalloc1(2*(depth+1), &depthShift);
667:   for (d = 0; d <= depth; d++) {
668:     PetscInt dEnd;

670:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
671:     depthShift[2*d]   = dEnd;
672:     depthShift[2*d+1] = 0;
673:   }
674:   if (depth >= 0) depthShift[2*depth+1] = Ng;
675:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
676:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
677:   /* Step 3: Set cone/support sizes for new points */
678:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
679:   DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
680:   for (c = cEnd; c < cEnd + Ng; ++c) {
681:     DMPlexSetConeSize(gdm, c, 1);
682:   }
683:   for (fs = 0; fs < numFS; ++fs) {
684:     IS              faceIS;
685:     const PetscInt *faces;
686:     PetscInt        numFaces, f;

688:     DMLabelGetStratumIS(label, values[fs], &faceIS);
689:     ISGetLocalSize(faceIS, &numFaces);
690:     ISGetIndices(faceIS, &faces);
691:     for (f = 0; f < numFaces; ++f) {
692:       PetscInt size, numChildren;

694:       PetscFindInt(faces[f], nleaves, leaves, &loc);
695:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
696:       if (loc >= 0 || numChildren) continue;
697:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
698:       DMPlexGetSupportSize(dm, faces[f], &size);
699:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
700:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
701:     }
702:     ISRestoreIndices(faceIS, &faces);
703:     ISDestroy(&faceIS);
704:   }
705:   /* Step 4: Setup ghosted DM */
706:   DMSetUp(gdm);
707:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
708:   /* Step 6: Set cones and supports for new points */
709:   ghostCell = cEnd;
710:   for (fs = 0; fs < numFS; ++fs) {
711:     IS              faceIS;
712:     const PetscInt *faces;
713:     PetscInt        numFaces, f;

715:     DMLabelGetStratumIS(label, values[fs], &faceIS);
716:     ISGetLocalSize(faceIS, &numFaces);
717:     ISGetIndices(faceIS, &faces);
718:     for (f = 0; f < numFaces; ++f) {
719:       PetscInt newFace = faces[f] + Ng, numChildren;

721:       PetscFindInt(faces[f], nleaves, leaves, &loc);
722:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
723:       if (loc >= 0 || numChildren) continue;
724:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
725:       DMPlexSetCone(gdm, ghostCell, &newFace);
726:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
727:       ++ghostCell;
728:     }
729:     ISRestoreIndices(faceIS, &faces);
730:     ISDestroy(&faceIS);
731:   }
732:   ISRestoreIndices(valueIS, &values);
733:   ISDestroy(&valueIS);
734:   /* Step 7: Stratify */
735:   DMPlexStratify(gdm);
736:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
737:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
738:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
739:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
740:   PetscFree(depthShift);
741:   /* Step 7: Periodicity */
742:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
743:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
744:   if (numGhostCells) *numGhostCells = Ng;
745:   return(0);
746: }

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

751:   Collective on dm

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

757:   Output Parameters:
758: + numGhostCells - The number of ghost cells added to the DM
759: - dmGhosted - The new DM

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

763:   Level: developer

765: .seealso: DMCreate()
766: @*/
767: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
768: {
769:   DM             gdm;
770:   DMLabel        label;
771:   const char    *name = labelName ? labelName : "Face Sets";
772:   PetscInt       dim;
773:   PetscBool      flag;

780:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
781:   DMSetType(gdm, DMPLEX);
782:   DMGetDimension(dm, &dim);
783:   DMSetDimension(gdm, dim);
784:   DMPlexGetAdjacencyUseCone(dm, &flag);
785:   DMPlexSetAdjacencyUseCone(gdm, flag);
786:   DMPlexGetAdjacencyUseClosure(dm, &flag);
787:   DMPlexSetAdjacencyUseClosure(gdm, flag);
788:   DMGetLabel(dm, name, &label);
789:   if (!label) {
790:     /* Get label for boundary faces */
791:     DMCreateLabel(dm, name);
792:     DMGetLabel(dm, name, &label);
793:     DMPlexMarkBoundaryFaces(dm, label);
794:   }
795:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
796:   DMCopyBoundary(dm, gdm);
797:   *dmGhosted = gdm;
798:   return(0);
799: }

801: /*
802:   We are adding three kinds of points here:
803:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
804:     Non-replicated: Points which exist on the fault, but are not replicated
805:     Hybrid:         Entirely new points, such as cohesive cells

807:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
808:   each depth so that the new split/hybrid points can be inserted as a block.
809: */
810: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
811: {
812:   MPI_Comm         comm;
813:   IS               valueIS;
814:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
815:   const PetscInt  *values;          /* List of depths for which we have replicated points */
816:   IS              *splitIS;
817:   IS              *unsplitIS;
818:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
819:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
820:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
821:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
822:   const PetscInt **splitPoints;        /* Replicated points for each depth */
823:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
824:   PetscSection     coordSection;
825:   Vec              coordinates;
826:   PetscScalar     *coords;
827:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
828:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
829:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
830:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
831:   PetscInt        *coneNew, *coneONew, *supportNew;
832:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
833:   PetscErrorCode   ierr;

836:   PetscObjectGetComm((PetscObject)dm,&comm);
837:   DMGetDimension(dm, &dim);
838:   DMPlexGetDepth(dm, &depth);
839:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
840:   /* Count split points and add cohesive cells */
841:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
842:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
843:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
844:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
845:   for (d = 0; d <= depth; ++d) {
846:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
847:     depthEnd[d]           = pMaxNew[d];
848:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
849:     numSplitPoints[d]     = 0;
850:     numUnsplitPoints[d]   = 0;
851:     numHybridPoints[d]    = 0;
852:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
853:     splitPoints[d]        = NULL;
854:     unsplitPoints[d]      = NULL;
855:     splitIS[d]            = NULL;
856:     unsplitIS[d]          = NULL;
857:     /* we are shifting the existing hybrid points with the stratum behind them, so
858:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
859:     depthShift[2*d]       = depthMax[d];
860:     depthShift[2*d+1]     = 0;
861:   }
862:   if (label) {
863:     DMLabelGetValueIS(label, &valueIS);
864:     ISGetLocalSize(valueIS, &numSP);
865:     ISGetIndices(valueIS, &values);
866:   }
867:   for (sp = 0; sp < numSP; ++sp) {
868:     const PetscInt dep = values[sp];

870:     if ((dep < 0) || (dep > depth)) continue;
871:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
872:     if (splitIS[dep]) {
873:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
874:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
875:     }
876:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
877:     if (unsplitIS[dep]) {
878:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
879:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
880:     }
881:   }
882:   /* Calculate number of hybrid points */
883:   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   */
884:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
885:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
886:   /* the end of the points in this stratum that come before the new points:
887:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
888:    * added points */
889:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
890:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
891:   /* Step 3: Set cone/support sizes for new points */
892:   for (dep = 0; dep <= depth; ++dep) {
893:     for (p = 0; p < numSplitPoints[dep]; ++p) {
894:       const PetscInt  oldp   = splitPoints[dep][p];
895:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
896:       const PetscInt  splitp = p    + pMaxNew[dep];
897:       const PetscInt *support;
898:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

900:       DMPlexGetConeSize(dm, oldp, &coneSize);
901:       DMPlexSetConeSize(sdm, splitp, coneSize);
902:       DMPlexGetSupportSize(dm, oldp, &supportSize);
903:       DMPlexSetSupportSize(sdm, splitp, supportSize);
904:       if (dep == depth-1) {
905:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

912:         DMPlexGetSupport(dm, oldp, &support);
913:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
914:           PetscInt val;

916:           DMLabelGetValue(label, support[e], &val);
917:           if (val == 1) ++qf;
918:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
919:           if ((val == 1) || (val == -(shift + 1))) ++qp;
920:         }
921:         /* Split old vertex: Edges into original vertex and new cohesive edge */
922:         DMPlexSetSupportSize(sdm, newp, qn+1);
923:         /* Split new vertex: Edges into split vertex and new cohesive edge */
924:         DMPlexSetSupportSize(sdm, splitp, qp+1);
925:         /* Add hybrid edge */
926:         DMPlexSetConeSize(sdm, hybedge, 2);
927:         DMPlexSetSupportSize(sdm, hybedge, qf);
928:       } else if (dep == dim-2) {
929:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

931:         DMPlexGetSupport(dm, oldp, &support);
932:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
933:           PetscInt val;

935:           DMLabelGetValue(label, support[e], &val);
936:           if (val == dim-1) ++qf;
937:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
938:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
939:         }
940:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
941:         DMPlexSetSupportSize(sdm, newp, qn+1);
942:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
943:         DMPlexSetSupportSize(sdm, splitp, qp+1);
944:         /* Add hybrid face */
945:         DMPlexSetConeSize(sdm, hybface, 4);
946:         DMPlexSetSupportSize(sdm, hybface, qf);
947:       }
948:     }
949:   }
950:   for (dep = 0; dep <= depth; ++dep) {
951:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
952:       const PetscInt  oldp   = unsplitPoints[dep][p];
953:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
954:       const PetscInt *support;
955:       PetscInt        coneSize, supportSize, qf, e, s;

957:       DMPlexGetConeSize(dm, oldp, &coneSize);
958:       DMPlexGetSupportSize(dm, oldp, &supportSize);
959:       DMPlexGetSupport(dm, oldp, &support);
960:       if (dep == 0) {
961:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

963:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
964:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
965:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
966:           if (e >= 0) ++qf;
967:         }
968:         DMPlexSetSupportSize(sdm, newp, qf+2);
969:         /* Add hybrid edge */
970:         DMPlexSetConeSize(sdm, hybedge, 2);
971:         for (e = 0, qf = 0; e < supportSize; ++e) {
972:           PetscInt val;

974:           DMLabelGetValue(label, support[e], &val);
975:           /* Split and unsplit edges produce hybrid faces */
976:           if (val == 1) ++qf;
977:           if (val == (shift2 + 1)) ++qf;
978:         }
979:         DMPlexSetSupportSize(sdm, hybedge, qf);
980:       } else if (dep == dim-2) {
981:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
982:         PetscInt       val;

984:         for (e = 0, qf = 0; e < supportSize; ++e) {
985:           DMLabelGetValue(label, support[e], &val);
986:           if (val == dim-1) qf += 2;
987:           else              ++qf;
988:         }
989:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
990:         DMPlexSetSupportSize(sdm, newp, qf+2);
991:         /* Add hybrid face */
992:         for (e = 0, qf = 0; e < supportSize; ++e) {
993:           DMLabelGetValue(label, support[e], &val);
994:           if (val == dim-1) ++qf;
995:         }
996:         DMPlexSetConeSize(sdm, hybface, 4);
997:         DMPlexSetSupportSize(sdm, hybface, qf);
998:       }
999:     }
1000:   }
1001:   /* Step 4: Setup split DM */
1002:   DMSetUp(sdm);
1003:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1004:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1005:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1006:   /* Step 6: Set cones and supports for new points */
1007:   for (dep = 0; dep <= depth; ++dep) {
1008:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1009:       const PetscInt  oldp   = splitPoints[dep][p];
1010:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1011:       const PetscInt  splitp = p    + pMaxNew[dep];
1012:       const PetscInt *cone, *support, *ornt;
1013:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1015:       DMPlexGetConeSize(dm, oldp, &coneSize);
1016:       DMPlexGetCone(dm, oldp, &cone);
1017:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1018:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1019:       DMPlexGetSupport(dm, oldp, &support);
1020:       if (dep == depth-1) {
1021:         PetscBool       hasUnsplit = PETSC_FALSE;
1022:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1023:         const PetscInt *supportF;

1025:         /* Split face:       copy in old face to new face to start */
1026:         DMPlexGetSupport(sdm, newp,  &supportF);
1027:         DMPlexSetSupport(sdm, splitp, supportF);
1028:         /* Split old face:   old vertices/edges in cone so no change */
1029:         /* Split new face:   new vertices/edges in cone */
1030:         for (q = 0; q < coneSize; ++q) {
1031:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1032:           if (v < 0) {
1033:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1034:             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);
1035:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1036:             hasUnsplit   = PETSC_TRUE;
1037:           } else {
1038:             coneNew[2+q] = v + pMaxNew[dep-1];
1039:             if (dep > 1) {
1040:               const PetscInt *econe;
1041:               PetscInt        econeSize, r, vs, vu;

1043:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1044:               DMPlexGetCone(dm, cone[q], &econe);
1045:               for (r = 0; r < econeSize; ++r) {
1046:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
1047:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1048:                 if (vs >= 0) continue;
1049:                 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);
1050:                 hasUnsplit   = PETSC_TRUE;
1051:               }
1052:             }
1053:           }
1054:         }
1055:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1056:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1057:         /* Face support */
1058:         for (s = 0; s < supportSize; ++s) {
1059:           PetscInt val;

1061:           DMLabelGetValue(label, support[s], &val);
1062:           if (val < 0) {
1063:             /* Split old face:   Replace negative side cell with cohesive cell */
1064:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1065:           } else {
1066:             /* Split new face:   Replace positive side cell with cohesive cell */
1067:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1068:             /* Get orientation for cohesive face */
1069:             {
1070:               const PetscInt *ncone, *nconeO;
1071:               PetscInt        nconeSize, nc;

1073:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1074:               DMPlexGetCone(dm, support[s], &ncone);
1075:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1076:               for (nc = 0; nc < nconeSize; ++nc) {
1077:                 if (ncone[nc] == oldp) {
1078:                   coneONew[0] = nconeO[nc];
1079:                   break;
1080:                 }
1081:               }
1082:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1083:             }
1084:           }
1085:         }
1086:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1087:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1088:         coneNew[1]  = splitp;
1089:         coneONew[1] = coneONew[0];
1090:         for (q = 0; q < coneSize; ++q) {
1091:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1092:           if (v < 0) {
1093:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1094:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1095:             coneONew[2+q] = 0;
1096:           } else {
1097:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1098:           }
1099:           coneONew[2+q] = 0;
1100:         }
1101:         DMPlexSetCone(sdm, hybcell, coneNew);
1102:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1103:         /* Label the hybrid cells on the boundary of the split */
1104:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1105:       } else if (dep == 0) {
1106:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1112:           DMLabelGetValue(label, support[e], &val);
1113:           if ((val == 1) || (val == (shift + 1))) {
1114:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1115:           }
1116:         }
1117:         supportNew[qn] = hybedge;
1118:         DMPlexSetSupport(sdm, newp, supportNew);
1119:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1120:         for (e = 0, qp = 0; e < supportSize; ++e) {
1121:           PetscInt val, edge;

1123:           DMLabelGetValue(label, support[e], &val);
1124:           if (val == 1) {
1125:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1126:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1127:             supportNew[qp++] = edge + pMaxNew[dep+1];
1128:           } else if (val == -(shift + 1)) {
1129:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1130:           }
1131:         }
1132:         supportNew[qp] = hybedge;
1133:         DMPlexSetSupport(sdm, splitp, supportNew);
1134:         /* Hybrid edge:    Old and new split vertex */
1135:         coneNew[0] = newp;
1136:         coneNew[1] = splitp;
1137:         DMPlexSetCone(sdm, hybedge, coneNew);
1138:         for (e = 0, qf = 0; e < supportSize; ++e) {
1139:           PetscInt val, edge;

1141:           DMLabelGetValue(label, support[e], &val);
1142:           if (val == 1) {
1143:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1144:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1145:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1146:           }
1147:         }
1148:         DMPlexSetSupport(sdm, hybedge, supportNew);
1149:       } else if (dep == dim-2) {
1150:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1152:         /* Split old edge:   old vertices in cone so no change */
1153:         /* Split new edge:   new vertices in cone */
1154:         for (q = 0; q < coneSize; ++q) {
1155:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1156:           if (v < 0) {
1157:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1158:             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);
1159:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1160:           } else {
1161:             coneNew[q] = v + pMaxNew[dep-1];
1162:           }
1163:         }
1164:         DMPlexSetCone(sdm, splitp, coneNew);
1165:         /* Split old edge: Faces in positive side cells and old split faces */
1166:         for (e = 0, q = 0; e < supportSize; ++e) {
1167:           PetscInt val;

1169:           DMLabelGetValue(label, support[e], &val);
1170:           if (val == dim-1) {
1171:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1172:           } else if (val == (shift + dim-1)) {
1173:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1174:           }
1175:         }
1176:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1177:         DMPlexSetSupport(sdm, newp, supportNew);
1178:         /* Split new edge: Faces in negative side cells and new split faces */
1179:         for (e = 0, q = 0; e < supportSize; ++e) {
1180:           PetscInt val, face;

1182:           DMLabelGetValue(label, support[e], &val);
1183:           if (val == dim-1) {
1184:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1185:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1186:             supportNew[q++] = face + pMaxNew[dep+1];
1187:           } else if (val == -(shift + dim-1)) {
1188:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1189:           }
1190:         }
1191:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1192:         DMPlexSetSupport(sdm, splitp, supportNew);
1193:         /* Hybrid face */
1194:         coneNew[0] = newp;
1195:         coneNew[1] = splitp;
1196:         for (v = 0; v < coneSize; ++v) {
1197:           PetscInt vertex;
1198:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1199:           if (vertex < 0) {
1200:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1201:             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);
1202:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1203:           } else {
1204:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1205:           }
1206:         }
1207:         DMPlexSetCone(sdm, hybface, coneNew);
1208:         for (e = 0, qf = 0; e < supportSize; ++e) {
1209:           PetscInt val, face;

1211:           DMLabelGetValue(label, support[e], &val);
1212:           if (val == dim-1) {
1213:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1214:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1215:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1216:           }
1217:         }
1218:         DMPlexSetSupport(sdm, hybface, supportNew);
1219:       }
1220:     }
1221:   }
1222:   for (dep = 0; dep <= depth; ++dep) {
1223:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1224:       const PetscInt  oldp   = unsplitPoints[dep][p];
1225:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1226:       const PetscInt *cone, *support, *ornt;
1227:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1229:       DMPlexGetConeSize(dm, oldp, &coneSize);
1230:       DMPlexGetCone(dm, oldp, &cone);
1231:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1232:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1233:       DMPlexGetSupport(dm, oldp, &support);
1234:       if (dep == 0) {
1235:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1237:         /* Unsplit vertex */
1238:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1239:         for (s = 0, q = 0; s < supportSize; ++s) {
1240:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1241:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1242:           if (e >= 0) {
1243:             supportNew[q++] = e + pMaxNew[dep+1];
1244:           }
1245:         }
1246:         supportNew[q++] = hybedge;
1247:         supportNew[q++] = hybedge;
1248:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1249:         DMPlexSetSupport(sdm, newp, supportNew);
1250:         /* Hybrid edge */
1251:         coneNew[0] = newp;
1252:         coneNew[1] = newp;
1253:         DMPlexSetCone(sdm, hybedge, coneNew);
1254:         for (e = 0, qf = 0; e < supportSize; ++e) {
1255:           PetscInt val, edge;

1257:           DMLabelGetValue(label, support[e], &val);
1258:           if (val == 1) {
1259:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1260:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1261:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1262:           } else if  (val ==  (shift2 + 1)) {
1263:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1264:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1265:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1266:           }
1267:         }
1268:         DMPlexSetSupport(sdm, hybedge, supportNew);
1269:       } else if (dep == dim-2) {
1270:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1276:           DMLabelGetValue(label, support[f], &val);
1277:           if (val == dim-1) {
1278:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1279:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1280:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1281:             supportNew[qf++] = face + pMaxNew[dep+1];
1282:           } else {
1283:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1284:           }
1285:         }
1286:         supportNew[qf++] = hybface;
1287:         supportNew[qf++] = hybface;
1288:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1289:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1290:         DMPlexSetSupport(sdm, newp, supportNew);
1291:         /* Add hybrid face */
1292:         coneNew[0] = newp;
1293:         coneNew[1] = newp;
1294:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1295:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1296:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1297:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1298:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1299:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1300:         DMPlexSetCone(sdm, hybface, coneNew);
1301:         for (f = 0, qf = 0; f < supportSize; ++f) {
1302:           PetscInt val, face;

1304:           DMLabelGetValue(label, support[f], &val);
1305:           if (val == dim-1) {
1306:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1307:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1308:           }
1309:         }
1310:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1311:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1312:         DMPlexSetSupport(sdm, hybface, supportNew);
1313:       }
1314:     }
1315:   }
1316:   /* Step 6b: Replace split points in negative side cones */
1317:   for (sp = 0; sp < numSP; ++sp) {
1318:     PetscInt        dep = values[sp];
1319:     IS              pIS;
1320:     PetscInt        numPoints;
1321:     const PetscInt *points;

1323:     if (dep >= 0) continue;
1324:     DMLabelGetStratumIS(label, dep, &pIS);
1325:     if (!pIS) continue;
1326:     dep  = -dep - shift;
1327:     ISGetLocalSize(pIS, &numPoints);
1328:     ISGetIndices(pIS, &points);
1329:     for (p = 0; p < numPoints; ++p) {
1330:       const PetscInt  oldp = points[p];
1331:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1332:       const PetscInt *cone;
1333:       PetscInt        coneSize, c;
1334:       /* PetscBool       replaced = PETSC_FALSE; */

1336:       /* Negative edge: replace split vertex */
1337:       /* Negative cell: replace split face */
1338:       DMPlexGetConeSize(sdm, newp, &coneSize);
1339:       DMPlexGetCone(sdm, newp, &cone);
1340:       for (c = 0; c < coneSize; ++c) {
1341:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1342:         PetscInt       csplitp, cp, val;

1344:         DMLabelGetValue(label, coldp, &val);
1345:         if (val == dep-1) {
1346:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1347:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1348:           csplitp  = pMaxNew[dep-1] + cp;
1349:           DMPlexInsertCone(sdm, newp, c, csplitp);
1350:           /* replaced = PETSC_TRUE; */
1351:         }
1352:       }
1353:       /* Cells with only a vertex or edge on the submesh have no replacement */
1354:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1355:     }
1356:     ISRestoreIndices(pIS, &points);
1357:     ISDestroy(&pIS);
1358:   }
1359:   /* Step 7: Stratify */
1360:   DMPlexStratify(sdm);
1361:   /* Step 8: Coordinates */
1362:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1363:   DMGetCoordinateSection(sdm, &coordSection);
1364:   DMGetCoordinatesLocal(sdm, &coordinates);
1365:   VecGetArray(coordinates, &coords);
1366:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1367:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1368:     const PetscInt splitp = pMaxNew[0] + v;
1369:     PetscInt       dof, off, soff, d;

1371:     PetscSectionGetDof(coordSection, newp, &dof);
1372:     PetscSectionGetOffset(coordSection, newp, &off);
1373:     PetscSectionGetOffset(coordSection, splitp, &soff);
1374:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1375:   }
1376:   VecRestoreArray(coordinates, &coords);
1377:   /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1378:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1379:   /* Step 10: Labels */
1380:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1381:   DMGetNumLabels(sdm, &numLabels);
1382:   for (dep = 0; dep <= depth; ++dep) {
1383:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1384:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1385:       const PetscInt splitp = pMaxNew[dep] + p;
1386:       PetscInt       l;

1388:       for (l = 0; l < numLabels; ++l) {
1389:         DMLabel     mlabel;
1390:         const char *lname;
1391:         PetscInt    val;
1392:         PetscBool   isDepth;

1394:         DMGetLabelName(sdm, l, &lname);
1395:         PetscStrcmp(lname, "depth", &isDepth);
1396:         if (isDepth) continue;
1397:         DMGetLabel(sdm, lname, &mlabel);
1398:         DMLabelGetValue(mlabel, newp, &val);
1399:         if (val >= 0) {
1400:           DMLabelSetValue(mlabel, splitp, val);
1401:         }
1402:       }
1403:     }
1404:   }
1405:   for (sp = 0; sp < numSP; ++sp) {
1406:     const PetscInt dep = values[sp];

1408:     if ((dep < 0) || (dep > depth)) continue;
1409:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1410:     ISDestroy(&splitIS[dep]);
1411:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1412:     ISDestroy(&unsplitIS[dep]);
1413:   }
1414:   if (label) {
1415:     ISRestoreIndices(valueIS, &values);
1416:     ISDestroy(&valueIS);
1417:   }
1418:   for (d = 0; d <= depth; ++d) {
1419:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1420:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1421:   }
1422:   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);
1423:   PetscFree3(coneNew, coneONew, supportNew);
1424:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1425:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1426:   return(0);
1427: }

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

1432:   Collective on dm

1434:   Input Parameters:
1435: + dm - The original DM
1436: - label - The label specifying the boundary faces (this could be auto-generated)

1438:   Output Parameters:
1439: - dmSplit - The new DM

1441:   Level: developer

1443: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1444: @*/
1445: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1446: {
1447:   DM             sdm;
1448:   PetscInt       dim;

1454:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1455:   DMSetType(sdm, DMPLEX);
1456:   DMGetDimension(dm, &dim);
1457:   DMSetDimension(sdm, dim);
1458:   switch (dim) {
1459:   case 2:
1460:   case 3:
1461:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1462:     break;
1463:   default:
1464:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1465:   }
1466:   *dmSplit = sdm;
1467:   return(0);
1468: }

1470: /* Returns the side of the surface for a given cell with a face on the surface */
1471: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1472: {
1473:   const PetscInt *cone, *ornt;
1474:   PetscInt        dim, coneSize, c;
1475:   PetscErrorCode  ierr;

1478:   *pos = PETSC_TRUE;
1479:   DMGetDimension(dm, &dim);
1480:   DMPlexGetConeSize(dm, cell, &coneSize);
1481:   DMPlexGetCone(dm, cell, &cone);
1482:   DMPlexGetConeOrientation(dm, cell, &ornt);
1483:   for (c = 0; c < coneSize; ++c) {
1484:     if (cone[c] == face) {
1485:       PetscInt o = ornt[c];

1487:       if (subdm) {
1488:         const PetscInt *subcone, *subornt;
1489:         PetscInt        subpoint, subface, subconeSize, sc;

1491:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1492:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1493:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1494:         DMPlexGetCone(subdm, subpoint, &subcone);
1495:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1496:         for (sc = 0; sc < subconeSize; ++sc) {
1497:           if (subcone[sc] == subface) {
1498:             o = subornt[0];
1499:             break;
1500:           }
1501:         }
1502:         if (sc >= subconeSize) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find subpoint %d (%d) in cone for subpoint %d (%d)", subface, face, subpoint, cell);
1503:       }
1504:       if (o >= 0) *pos = PETSC_TRUE;
1505:       else        *pos = PETSC_FALSE;
1506:       break;
1507:     }
1508:   }
1509:   if (c == coneSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell %d in split face %d support does not have it in the cone", cell, face);
1510:   return(0);
1511: }

1513: /*@
1514:   DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1515:   to complete the surface

1517:   Input Parameters:
1518: + dm     - The DM
1519: . label  - A DMLabel marking the surface
1520: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1521: . flip   - Flag to flip the submesh normal and replace points on the other side
1522: - subdm  - The subDM associated with the label, or NULL

1524:   Output Parameter:
1525: . label - A DMLabel marking all surface points

1527:   Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.

1529:   Level: developer

1531: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1532: @*/
1533: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1534: {
1535:   DMLabel         depthLabel;
1536:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1537:   const PetscInt *points, *subpoints;
1538:   const PetscInt  rev   = flip ? -1 : 1;
1539:   PetscInt       *pMax;
1540:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1541:   PetscErrorCode  ierr;

1544:   DMPlexGetDepth(dm, &depth);
1545:   DMGetDimension(dm, &dim);
1546:   pSize = PetscMax(depth, dim) + 1;
1547:   PetscMalloc1(pSize,&pMax);
1548:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1549:   DMPlexGetDepthLabel(dm, &depthLabel);
1550:   DMGetDimension(dm, &dim);
1551:   if (subdm) {
1552:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1553:     if (subpointIS) {
1554:       ISGetLocalSize(subpointIS, &numSubpoints);
1555:       ISGetIndices(subpointIS, &subpoints);
1556:     }
1557:   }
1558:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1559:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1560:   if (!dimIS) {
1561:     PetscFree(pMax);
1562:     ISDestroy(&subpointIS);
1563:     return(0);
1564:   }
1565:   ISGetLocalSize(dimIS, &numPoints);
1566:   ISGetIndices(dimIS, &points);
1567:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1568:     const PetscInt *support;
1569:     PetscInt        supportSize, s;

1571:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1572:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1573:     DMPlexGetSupport(dm, points[p], &support);
1574:     for (s = 0; s < supportSize; ++s) {
1575:       const PetscInt *cone;
1576:       PetscInt        coneSize, c;
1577:       PetscBool       pos;

1579:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1580:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1581:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1582:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1583:       /* Put faces touching the fault in the label */
1584:       DMPlexGetConeSize(dm, support[s], &coneSize);
1585:       DMPlexGetCone(dm, support[s], &cone);
1586:       for (c = 0; c < coneSize; ++c) {
1587:         const PetscInt point = cone[c];

1589:         DMLabelGetValue(label, point, &val);
1590:         if (val == -1) {
1591:           PetscInt *closure = NULL;
1592:           PetscInt  closureSize, cl;

1594:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1595:           for (cl = 0; cl < closureSize*2; cl += 2) {
1596:             const PetscInt clp  = closure[cl];
1597:             PetscInt       bval = -1;

1599:             DMLabelGetValue(label, clp, &val);
1600:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1601:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1602:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1603:               break;
1604:             }
1605:           }
1606:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1607:         }
1608:       }
1609:     }
1610:   }
1611:   ISRestoreIndices(dimIS, &points);
1612:   ISDestroy(&dimIS);
1613:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1614:   ISDestroy(&subpointIS);
1615:   /* Mark boundary points as unsplit */
1616:   if (blabel) {
1617:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1618:     ISGetLocalSize(dimIS, &numPoints);
1619:     ISGetIndices(dimIS, &points);
1620:     for (p = 0; p < numPoints; ++p) {
1621:       const PetscInt point = points[p];
1622:       PetscInt       val, bval;

1624:       DMLabelGetValue(blabel, point, &bval);
1625:       if (bval >= 0) {
1626:         DMLabelGetValue(label, point, &val);
1627:         if ((val < 0) || (val > dim)) {
1628:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1629:           DMLabelClearValue(blabel, point, bval);
1630:         }
1631:       }
1632:     }
1633:     for (p = 0; p < numPoints; ++p) {
1634:       const PetscInt point = points[p];
1635:       PetscInt       val, bval;

1637:       DMLabelGetValue(blabel, point, &bval);
1638:       if (bval >= 0) {
1639:         const PetscInt *cone,    *support;
1640:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1642:         /* Mark as unsplit */
1643:         DMLabelGetValue(label, point, &val);
1644:         if ((val < 0) || (val > dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d has label value %d, should be part of the fault", point, val);
1645:         DMLabelClearValue(label, point, val);
1646:         DMLabelSetValue(label, point, shift2+val);
1647:         /* Check for cross-edge
1648:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1649:         if (val != 0) continue;
1650:         DMPlexGetSupport(dm, point, &support);
1651:         DMPlexGetSupportSize(dm, point, &supportSize);
1652:         for (s = 0; s < supportSize; ++s) {
1653:           DMPlexGetCone(dm, support[s], &cone);
1654:           DMPlexGetConeSize(dm, support[s], &coneSize);
1655:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1656:           DMLabelGetValue(blabel, cone[0], &valA);
1657:           DMLabelGetValue(blabel, cone[1], &valB);
1658:           DMLabelGetValue(blabel, support[s], &valE);
1659:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1660:         }
1661:       }
1662:     }
1663:     ISRestoreIndices(dimIS, &points);
1664:     ISDestroy(&dimIS);
1665:   }
1666:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1667:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1668:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1669:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1670:   cMax = cMax < 0 ? cEnd : cMax;
1671:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1672:   DMLabelGetStratumIS(label, 0, &dimIS);
1673:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1674:   if (dimIS && crossEdgeIS) {
1675:     IS vertIS = dimIS;

1677:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1678:     ISDestroy(&crossEdgeIS);
1679:     ISDestroy(&vertIS);
1680:   }
1681:   if (!dimIS) {
1682:     PetscFree(pMax);
1683:     return(0);
1684:   }
1685:   ISGetLocalSize(dimIS, &numPoints);
1686:   ISGetIndices(dimIS, &points);
1687:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1688:     PetscInt *star = NULL;
1689:     PetscInt  starSize, s;
1690:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1692:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1693:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1694:     while (again) {
1695:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1696:       again = 0;
1697:       for (s = 0; s < starSize*2; s += 2) {
1698:         const PetscInt  point = star[s];
1699:         const PetscInt *cone;
1700:         PetscInt        coneSize, c;

1702:         if ((point < cStart) || (point >= cMax)) continue;
1703:         DMLabelGetValue(label, point, &val);
1704:         if (val != -1) continue;
1705:         again = again == 1 ? 1 : 2;
1706:         DMPlexGetConeSize(dm, point, &coneSize);
1707:         DMPlexGetCone(dm, point, &cone);
1708:         for (c = 0; c < coneSize; ++c) {
1709:           DMLabelGetValue(label, cone[c], &val);
1710:           if (val != -1) {
1711:             const PetscInt *ccone;
1712:             PetscInt        cconeSize, cc, side;

1714:             if (PetscAbs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
1715:             if (val > 0) side =  1;
1716:             else         side = -1;
1717:             DMLabelSetValue(label, point, side*(shift+dim));
1718:             /* Mark cell faces which touch the fault */
1719:             DMPlexGetConeSize(dm, point, &cconeSize);
1720:             DMPlexGetCone(dm, point, &ccone);
1721:             for (cc = 0; cc < cconeSize; ++cc) {
1722:               PetscInt *closure = NULL;
1723:               PetscInt  closureSize, cl;

1725:               DMLabelGetValue(label, ccone[cc], &val);
1726:               if (val != -1) continue;
1727:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1728:               for (cl = 0; cl < closureSize*2; cl += 2) {
1729:                 const PetscInt clp = closure[cl];

1731:                 DMLabelGetValue(label, clp, &val);
1732:                 if (val == -1) continue;
1733:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1734:                 break;
1735:               }
1736:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1737:             }
1738:             again = 1;
1739:             break;
1740:           }
1741:         }
1742:       }
1743:     }
1744:     /* Classify the rest by cell membership */
1745:     for (s = 0; s < starSize*2; s += 2) {
1746:       const PetscInt point = star[s];

1748:       DMLabelGetValue(label, point, &val);
1749:       if (val == -1) {
1750:         PetscInt  *sstar = NULL;
1751:         PetscInt   sstarSize, ss;
1752:         PetscBool  marked = PETSC_FALSE;

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

1758:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1759:           DMLabelGetValue(label, spoint, &val);
1760:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1761:           DMLabelGetValue(depthLabel, point, &dep);
1762:           if (val > 0) {
1763:             DMLabelSetValue(label, point,   shift+dep);
1764:           } else {
1765:             DMLabelSetValue(label, point, -(shift+dep));
1766:           }
1767:           marked = PETSC_TRUE;
1768:           break;
1769:         }
1770:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1771:         DMLabelGetValue(depthLabel, point, &dep);
1772:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1773:       }
1774:     }
1775:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1776:   }
1777:   ISRestoreIndices(dimIS, &points);
1778:   ISDestroy(&dimIS);
1779:   /* If any faces touching the fault divide cells on either side, split them */
1780:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1781:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1782:   ISExpand(facePosIS, faceNegIS, &dimIS);
1783:   ISDestroy(&facePosIS);
1784:   ISDestroy(&faceNegIS);
1785:   ISGetLocalSize(dimIS, &numPoints);
1786:   ISGetIndices(dimIS, &points);
1787:   for (p = 0; p < numPoints; ++p) {
1788:     const PetscInt  point = points[p];
1789:     const PetscInt *support;
1790:     PetscInt        supportSize, valA, valB;

1792:     DMPlexGetSupportSize(dm, point, &supportSize);
1793:     if (supportSize != 2) continue;
1794:     DMPlexGetSupport(dm, point, &support);
1795:     DMLabelGetValue(label, support[0], &valA);
1796:     DMLabelGetValue(label, support[1], &valB);
1797:     if ((valA == -1) || (valB == -1)) continue;
1798:     if (valA*valB > 0) continue;
1799:     /* Split the face */
1800:     DMLabelGetValue(label, point, &valA);
1801:     DMLabelClearValue(label, point, valA);
1802:     DMLabelSetValue(label, point, dim-1);
1803:     /* Label its closure:
1804:       unmarked: label as unsplit
1805:       incident: relabel as split
1806:       split:    do nothing
1807:     */
1808:     {
1809:       PetscInt *closure = NULL;
1810:       PetscInt  closureSize, cl;

1812:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1813:       for (cl = 0; cl < closureSize*2; cl += 2) {
1814:         DMLabelGetValue(label, closure[cl], &valA);
1815:         if (valA == -1) { /* Mark as unsplit */
1816:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1817:           DMLabelSetValue(label, closure[cl], shift2+dep);
1818:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1819:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1820:           DMLabelClearValue(label, closure[cl], valA);
1821:           DMLabelSetValue(label, closure[cl], dep);
1822:         }
1823:       }
1824:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1825:     }
1826:   }
1827:   ISRestoreIndices(dimIS, &points);
1828:   ISDestroy(&dimIS);
1829:   PetscFree(pMax);
1830:   return(0);
1831: }

1833: /*@
1834:   DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface

1836:   Collective on dm

1838:   Input Parameters:
1839: + dm - The original DM
1840: - labelName - The label specifying the interface vertices

1842:   Output Parameters:
1843: + hybridLabel - The label fully marking the interface
1844: - dmHybrid - The new DM

1846:   Level: developer

1848: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1849: @*/
1850: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1851: {
1852:   DM             idm;
1853:   DMLabel        subpointMap, hlabel;
1854:   PetscInt       dim;

1861:   DMGetDimension(dm, &dim);
1862:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1863:   DMPlexOrient(idm);
1864:   DMPlexGetSubpointMap(idm, &subpointMap);
1865:   DMLabelDuplicate(subpointMap, &hlabel);
1866:   DMLabelClearStratum(hlabel, dim);
1867:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1868:   DMDestroy(&idm);
1869:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1870:   if (hybridLabel) *hybridLabel = hlabel;
1871:   else             {DMLabelDestroy(&hlabel);}
1872:   return(0);
1873: }

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

1877:      For any marked cell, the marked vertices constitute a single face
1878: */
1879: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1880: {
1881:   IS               subvertexIS = NULL;
1882:   const PetscInt  *subvertices;
1883:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1884:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1885:   PetscErrorCode   ierr;

1888:   *numFaces = 0;
1889:   *nFV      = 0;
1890:   DMPlexGetDepth(dm, &depth);
1891:   DMGetDimension(dm, &dim);
1892:   pSize = PetscMax(depth, dim) + 1;
1893:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1894:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1895:   for (d = 0; d <= depth; ++d) {
1896:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1897:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1898:   }
1899:   /* Loop over initial vertices and mark all faces in the collective star() */
1900:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1901:   if (subvertexIS) {
1902:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1903:     ISGetIndices(subvertexIS, &subvertices);
1904:   }
1905:   for (v = 0; v < numSubVerticesInitial; ++v) {
1906:     const PetscInt vertex = subvertices[v];
1907:     PetscInt      *star   = NULL;
1908:     PetscInt       starSize, s, numCells = 0, c;

1910:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1911:     for (s = 0; s < starSize*2; s += 2) {
1912:       const PetscInt point = star[s];
1913:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1914:     }
1915:     for (c = 0; c < numCells; ++c) {
1916:       const PetscInt cell    = star[c];
1917:       PetscInt      *closure = NULL;
1918:       PetscInt       closureSize, cl;
1919:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1921:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1922:       if (cellLoc == 2) continue;
1923:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1924:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1925:       for (cl = 0; cl < closureSize*2; cl += 2) {
1926:         const PetscInt point = closure[cl];
1927:         PetscInt       vertexLoc;

1929:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1930:           ++numCorners;
1931:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1932:           if (vertexLoc == value) closure[faceSize++] = point;
1933:         }
1934:       }
1935:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1936:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1937:       if (faceSize == *nFV) {
1938:         const PetscInt *cells = NULL;
1939:         PetscInt        numCells, nc;

1941:         ++(*numFaces);
1942:         for (cl = 0; cl < faceSize; ++cl) {
1943:           DMLabelSetValue(subpointMap, closure[cl], 0);
1944:         }
1945:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1946:         for (nc = 0; nc < numCells; ++nc) {
1947:           DMLabelSetValue(subpointMap, cells[nc], 2);
1948:         }
1949:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1950:       }
1951:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1952:     }
1953:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1954:   }
1955:   if (subvertexIS) {
1956:     ISRestoreIndices(subvertexIS, &subvertices);
1957:   }
1958:   ISDestroy(&subvertexIS);
1959:   PetscFree3(pStart,pEnd,pMax);
1960:   return(0);
1961: }

1963: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1964: {
1965:   IS               subvertexIS = NULL;
1966:   const PetscInt  *subvertices;
1967:   PetscInt        *pStart, *pEnd, *pMax;
1968:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
1969:   PetscErrorCode   ierr;

1972:   DMGetDimension(dm, &dim);
1973:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1974:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1975:   for (d = 0; d <= dim; ++d) {
1976:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1977:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1978:   }
1979:   /* Loop over initial vertices and mark all faces in the collective star() */
1980:   if (vertexLabel) {
1981:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1982:     if (subvertexIS) {
1983:       ISGetSize(subvertexIS, &numSubVerticesInitial);
1984:       ISGetIndices(subvertexIS, &subvertices);
1985:     }
1986:   }
1987:   for (v = 0; v < numSubVerticesInitial; ++v) {
1988:     const PetscInt vertex = subvertices[v];
1989:     PetscInt      *star   = NULL;
1990:     PetscInt       starSize, s, numFaces = 0, f;

1992:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1993:     for (s = 0; s < starSize*2; s += 2) {
1994:       const PetscInt point = star[s];
1995:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1996:     }
1997:     for (f = 0; f < numFaces; ++f) {
1998:       const PetscInt face    = star[f];
1999:       PetscInt      *closure = NULL;
2000:       PetscInt       closureSize, c;
2001:       PetscInt       faceLoc;

2003:       DMLabelGetValue(subpointMap, face, &faceLoc);
2004:       if (faceLoc == dim-1) continue;
2005:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2006:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2007:       for (c = 0; c < closureSize*2; c += 2) {
2008:         const PetscInt point = closure[c];
2009:         PetscInt       vertexLoc;

2011:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2012:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2013:           if (vertexLoc != value) break;
2014:         }
2015:       }
2016:       if (c == closureSize*2) {
2017:         const PetscInt *support;
2018:         PetscInt        supportSize, s;

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

2023:           for (d = 0; d < dim; ++d) {
2024:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2025:               DMLabelSetValue(subpointMap, point, d);
2026:               break;
2027:             }
2028:           }
2029:         }
2030:         DMPlexGetSupportSize(dm, face, &supportSize);
2031:         DMPlexGetSupport(dm, face, &support);
2032:         for (s = 0; s < supportSize; ++s) {
2033:           DMLabelSetValue(subpointMap, support[s], dim);
2034:         }
2035:       }
2036:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2037:     }
2038:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2039:   }
2040:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2041:   ISDestroy(&subvertexIS);
2042:   PetscFree3(pStart,pEnd,pMax);
2043:   return(0);
2044: }

2046: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2047: {
2048:   DMLabel         label = NULL;
2049:   const PetscInt *cone;
2050:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2051:   PetscErrorCode  ierr;

2054:   *numFaces = 0;
2055:   *nFV = 0;
2056:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2057:   *subCells = NULL;
2058:   DMGetDimension(dm, &dim);
2059:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2060:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2061:   if (cMax < 0) return(0);
2062:   if (label) {
2063:     for (c = cMax; c < cEnd; ++c) {
2064:       PetscInt val;

2066:       DMLabelGetValue(label, c, &val);
2067:       if (val == value) {
2068:         ++(*numFaces);
2069:         DMPlexGetConeSize(dm, c, &coneSize);
2070:       }
2071:     }
2072:   } else {
2073:     *numFaces = cEnd - cMax;
2074:     DMPlexGetConeSize(dm, cMax, &coneSize);
2075:   }
2076:   PetscMalloc1(*numFaces *2, subCells);
2077:   if (!(*numFaces)) return(0);
2078:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2079:   for (c = cMax; c < cEnd; ++c) {
2080:     const PetscInt *cells;
2081:     PetscInt        numCells;

2083:     if (label) {
2084:       PetscInt val;

2086:       DMLabelGetValue(label, c, &val);
2087:       if (val != value) continue;
2088:     }
2089:     DMPlexGetCone(dm, c, &cone);
2090:     for (p = 0; p < *nFV; ++p) {
2091:       DMLabelSetValue(subpointMap, cone[p], 0);
2092:     }
2093:     /* Negative face */
2094:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2095:     /* Not true in parallel
2096:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2097:     for (p = 0; p < numCells; ++p) {
2098:       DMLabelSetValue(subpointMap, cells[p], 2);
2099:       (*subCells)[subc++] = cells[p];
2100:     }
2101:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2102:     /* Positive face is not included */
2103:   }
2104:   return(0);
2105: }

2107: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2108: {
2109:   PetscInt      *pStart, *pEnd;
2110:   PetscInt       dim, cMax, cEnd, c, d;

2114:   DMGetDimension(dm, &dim);
2115:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2116:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2117:   if (cMax < 0) return(0);
2118:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2119:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2120:   for (c = cMax; c < cEnd; ++c) {
2121:     const PetscInt *cone;
2122:     PetscInt       *closure = NULL;
2123:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2125:     if (label) {
2126:       DMLabelGetValue(label, c, &val);
2127:       if (val != value) continue;
2128:     }
2129:     DMPlexGetConeSize(dm, c, &coneSize);
2130:     DMPlexGetCone(dm, c, &cone);
2131:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2132:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2133:     /* Negative face */
2134:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2135:     for (cl = 0; cl < closureSize*2; cl += 2) {
2136:       const PetscInt point = closure[cl];

2138:       for (d = 0; d <= dim; ++d) {
2139:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2140:           DMLabelSetValue(subpointMap, point, d);
2141:           break;
2142:         }
2143:       }
2144:     }
2145:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2146:     /* Cells -- positive face is not included */
2147:     for (cl = 0; cl < 1; ++cl) {
2148:       const PetscInt *support;
2149:       PetscInt        supportSize, s;

2151:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2152:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2153:       DMPlexGetSupport(dm, cone[cl], &support);
2154:       for (s = 0; s < supportSize; ++s) {
2155:         DMLabelSetValue(subpointMap, support[s], dim);
2156:       }
2157:     }
2158:   }
2159:   PetscFree2(pStart, pEnd);
2160:   return(0);
2161: }

2163: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2164: {
2165:   MPI_Comm       comm;
2166:   PetscBool      posOrient = PETSC_FALSE;
2167:   const PetscInt debug     = 0;
2168:   PetscInt       cellDim, faceSize, f;

2172:   PetscObjectGetComm((PetscObject)dm,&comm);
2173:   DMGetDimension(dm, &cellDim);
2174:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2176:   if (cellDim == 1 && numCorners == 2) {
2177:     /* Triangle */
2178:     faceSize  = numCorners-1;
2179:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2180:   } else if (cellDim == 2 && numCorners == 3) {
2181:     /* Triangle */
2182:     faceSize  = numCorners-1;
2183:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2184:   } else if (cellDim == 3 && numCorners == 4) {
2185:     /* Tetrahedron */
2186:     faceSize  = numCorners-1;
2187:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2188:   } else if (cellDim == 1 && numCorners == 3) {
2189:     /* Quadratic line */
2190:     faceSize  = 1;
2191:     posOrient = PETSC_TRUE;
2192:   } else if (cellDim == 2 && numCorners == 4) {
2193:     /* Quads */
2194:     faceSize = 2;
2195:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2196:       posOrient = PETSC_TRUE;
2197:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2198:       posOrient = PETSC_TRUE;
2199:     } else {
2200:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2201:         posOrient = PETSC_FALSE;
2202:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2203:     }
2204:   } else if (cellDim == 2 && numCorners == 6) {
2205:     /* Quadratic triangle (I hate this) */
2206:     /* Edges are determined by the first 2 vertices (corners of edges) */
2207:     const PetscInt faceSizeTri = 3;
2208:     PetscInt       sortedIndices[3], i, iFace;
2209:     PetscBool      found                    = PETSC_FALSE;
2210:     PetscInt       faceVerticesTriSorted[9] = {
2211:       0, 3,  4, /* bottom */
2212:       1, 4,  5, /* right */
2213:       2, 3,  5, /* left */
2214:     };
2215:     PetscInt       faceVerticesTri[9] = {
2216:       0, 3,  4, /* bottom */
2217:       1, 4,  5, /* right */
2218:       2, 5,  3, /* left */
2219:     };

2221:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2222:     PetscSortInt(faceSizeTri, sortedIndices);
2223:     for (iFace = 0; iFace < 3; ++iFace) {
2224:       const PetscInt ii = iFace*faceSizeTri;
2225:       PetscInt       fVertex, cVertex;

2227:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2228:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2229:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2230:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2231:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2232:               faceVertices[fVertex] = origVertices[cVertex];
2233:               break;
2234:             }
2235:           }
2236:         }
2237:         found = PETSC_TRUE;
2238:         break;
2239:       }
2240:     }
2241:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2242:     if (posOriented) *posOriented = PETSC_TRUE;
2243:     return(0);
2244:   } else if (cellDim == 2 && numCorners == 9) {
2245:     /* Quadratic quad (I hate this) */
2246:     /* Edges are determined by the first 2 vertices (corners of edges) */
2247:     const PetscInt faceSizeQuad = 3;
2248:     PetscInt       sortedIndices[3], i, iFace;
2249:     PetscBool      found                      = PETSC_FALSE;
2250:     PetscInt       faceVerticesQuadSorted[12] = {
2251:       0, 1,  4, /* bottom */
2252:       1, 2,  5, /* right */
2253:       2, 3,  6, /* top */
2254:       0, 3,  7, /* left */
2255:     };
2256:     PetscInt       faceVerticesQuad[12] = {
2257:       0, 1,  4, /* bottom */
2258:       1, 2,  5, /* right */
2259:       2, 3,  6, /* top */
2260:       3, 0,  7, /* left */
2261:     };

2263:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2264:     PetscSortInt(faceSizeQuad, sortedIndices);
2265:     for (iFace = 0; iFace < 4; ++iFace) {
2266:       const PetscInt ii = iFace*faceSizeQuad;
2267:       PetscInt       fVertex, cVertex;

2269:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2270:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2271:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2272:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2273:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2274:               faceVertices[fVertex] = origVertices[cVertex];
2275:               break;
2276:             }
2277:           }
2278:         }
2279:         found = PETSC_TRUE;
2280:         break;
2281:       }
2282:     }
2283:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2284:     if (posOriented) *posOriented = PETSC_TRUE;
2285:     return(0);
2286:   } else if (cellDim == 3 && numCorners == 8) {
2287:     /* Hexes
2288:        A hex is two oriented quads with the normal of the first
2289:        pointing up at the second.

2291:           7---6
2292:          /|  /|
2293:         4---5 |
2294:         | 1-|-2
2295:         |/  |/
2296:         0---3

2298:         Faces are determined by the first 4 vertices (corners of faces) */
2299:     const PetscInt faceSizeHex = 4;
2300:     PetscInt       sortedIndices[4], i, iFace;
2301:     PetscBool      found                     = PETSC_FALSE;
2302:     PetscInt       faceVerticesHexSorted[24] = {
2303:       0, 1, 2, 3,  /* bottom */
2304:       4, 5, 6, 7,  /* top */
2305:       0, 3, 4, 5,  /* front */
2306:       2, 3, 5, 6,  /* right */
2307:       1, 2, 6, 7,  /* back */
2308:       0, 1, 4, 7,  /* left */
2309:     };
2310:     PetscInt       faceVerticesHex[24] = {
2311:       1, 2, 3, 0,  /* bottom */
2312:       4, 5, 6, 7,  /* top */
2313:       0, 3, 5, 4,  /* front */
2314:       3, 2, 6, 5,  /* right */
2315:       2, 1, 7, 6,  /* back */
2316:       1, 0, 4, 7,  /* left */
2317:     };

2319:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2320:     PetscSortInt(faceSizeHex, sortedIndices);
2321:     for (iFace = 0; iFace < 6; ++iFace) {
2322:       const PetscInt ii = iFace*faceSizeHex;
2323:       PetscInt       fVertex, cVertex;

2325:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2326:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2327:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2328:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2329:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2330:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2331:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2332:               faceVertices[fVertex] = origVertices[cVertex];
2333:               break;
2334:             }
2335:           }
2336:         }
2337:         found = PETSC_TRUE;
2338:         break;
2339:       }
2340:     }
2341:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2342:     if (posOriented) *posOriented = PETSC_TRUE;
2343:     return(0);
2344:   } else if (cellDim == 3 && numCorners == 10) {
2345:     /* Quadratic tet */
2346:     /* Faces are determined by the first 3 vertices (corners of faces) */
2347:     const PetscInt faceSizeTet = 6;
2348:     PetscInt       sortedIndices[6], i, iFace;
2349:     PetscBool      found                     = PETSC_FALSE;
2350:     PetscInt       faceVerticesTetSorted[24] = {
2351:       0, 1, 2,  6, 7, 8, /* bottom */
2352:       0, 3, 4,  6, 7, 9,  /* front */
2353:       1, 4, 5,  7, 8, 9,  /* right */
2354:       2, 3, 5,  6, 8, 9,  /* left */
2355:     };
2356:     PetscInt       faceVerticesTet[24] = {
2357:       0, 1, 2,  6, 7, 8, /* bottom */
2358:       0, 4, 3,  6, 7, 9,  /* front */
2359:       1, 5, 4,  7, 8, 9,  /* right */
2360:       2, 3, 5,  8, 6, 9,  /* left */
2361:     };

2363:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2364:     PetscSortInt(faceSizeTet, sortedIndices);
2365:     for (iFace=0; iFace < 4; ++iFace) {
2366:       const PetscInt ii = iFace*faceSizeTet;
2367:       PetscInt       fVertex, cVertex;

2369:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2370:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2371:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2372:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2373:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2374:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2375:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2376:               faceVertices[fVertex] = origVertices[cVertex];
2377:               break;
2378:             }
2379:           }
2380:         }
2381:         found = PETSC_TRUE;
2382:         break;
2383:       }
2384:     }
2385:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2386:     if (posOriented) *posOriented = PETSC_TRUE;
2387:     return(0);
2388:   } else if (cellDim == 3 && numCorners == 27) {
2389:     /* Quadratic hexes (I hate this)
2390:        A hex is two oriented quads with the normal of the first
2391:        pointing up at the second.

2393:          7---6
2394:         /|  /|
2395:        4---5 |
2396:        | 3-|-2
2397:        |/  |/
2398:        0---1

2400:        Faces are determined by the first 4 vertices (corners of faces) */
2401:     const PetscInt faceSizeQuadHex = 9;
2402:     PetscInt       sortedIndices[9], i, iFace;
2403:     PetscBool      found                         = PETSC_FALSE;
2404:     PetscInt       faceVerticesQuadHexSorted[54] = {
2405:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2406:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2407:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2408:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2409:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2410:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2411:     };
2412:     PetscInt       faceVerticesQuadHex[54] = {
2413:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2414:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2415:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2416:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2417:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2418:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2419:     };

2421:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2422:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2423:     for (iFace = 0; iFace < 6; ++iFace) {
2424:       const PetscInt ii = iFace*faceSizeQuadHex;
2425:       PetscInt       fVertex, cVertex;

2427:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2428:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2429:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2430:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2431:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2432:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2433:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2434:               faceVertices[fVertex] = origVertices[cVertex];
2435:               break;
2436:             }
2437:           }
2438:         }
2439:         found = PETSC_TRUE;
2440:         break;
2441:       }
2442:     }
2443:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2444:     if (posOriented) *posOriented = PETSC_TRUE;
2445:     return(0);
2446:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2447:   if (!posOrient) {
2448:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2449:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2450:   } else {
2451:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2452:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2453:   }
2454:   if (posOriented) *posOriented = posOrient;
2455:   return(0);
2456: }

2458: /*@
2459:   DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2460:   in faceVertices. The orientation is such that the face normal points out of the cell

2462:   Not collective

2464:   Input Parameters:
2465: + dm           - The original mesh
2466: . cell         - The cell mesh point
2467: . faceSize     - The number of vertices on the face
2468: . face         - The face vertices
2469: . numCorners   - The number of vertices on the cell
2470: . indices      - Local numbering of face vertices in cell cone
2471: - origVertices - Original face vertices

2473:   Output Parameter:
2474: + faceVertices - The face vertices properly oriented
2475: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2477:   Level: developer

2479: .seealso: DMPlexGetCone()
2480: @*/
2481: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2482: {
2483:   const PetscInt *cone = NULL;
2484:   PetscInt        coneSize, v, f, v2;
2485:   PetscInt        oppositeVertex = -1;
2486:   PetscErrorCode  ierr;

2489:   DMPlexGetConeSize(dm, cell, &coneSize);
2490:   DMPlexGetCone(dm, cell, &cone);
2491:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2492:     PetscBool found = PETSC_FALSE;

2494:     for (f = 0; f < faceSize; ++f) {
2495:       if (face[f] == cone[v]) {
2496:         found = PETSC_TRUE; break;
2497:       }
2498:     }
2499:     if (found) {
2500:       indices[v2]      = v;
2501:       origVertices[v2] = cone[v];
2502:       ++v2;
2503:     } else {
2504:       oppositeVertex = v;
2505:     }
2506:   }
2507:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2508:   return(0);
2509: }

2511: /*
2512:   DMPlexInsertFace_Internal - Puts a face into the mesh

2514:   Not collective

2516:   Input Parameters:
2517:   + dm              - The DMPlex
2518:   . numFaceVertex   - The number of vertices in the face
2519:   . faceVertices    - The vertices in the face for dm
2520:   . subfaceVertices - The vertices in the face for subdm
2521:   . numCorners      - The number of vertices in the cell
2522:   . cell            - A cell in dm containing the face
2523:   . subcell         - A cell in subdm containing the face
2524:   . firstFace       - First face in the mesh
2525:   - newFacePoint    - Next face in the mesh

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

2530:   Level: developer
2531: */
2532: 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)
2533: {
2534:   MPI_Comm        comm;
2535:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2536:   const PetscInt *faces;
2537:   PetscInt        numFaces, coneSize;
2538:   PetscErrorCode  ierr;

2541:   PetscObjectGetComm((PetscObject)dm,&comm);
2542:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2543:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2544: #if 0
2545:   /* Cannot use this because support() has not been constructed yet */
2546:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2547: #else
2548:   {
2549:     PetscInt f;

2551:     numFaces = 0;
2552:     DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2553:     for (f = firstFace; f < *newFacePoint; ++f) {
2554:       PetscInt dof, off, d;

2556:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2557:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2558:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2559:       for (d = 0; d < dof; ++d) {
2560:         const PetscInt p = submesh->cones[off+d];
2561:         PetscInt       v;

2563:         for (v = 0; v < numFaceVertices; ++v) {
2564:           if (subfaceVertices[v] == p) break;
2565:         }
2566:         if (v == numFaceVertices) break;
2567:       }
2568:       if (d == dof) {
2569:         numFaces               = 1;
2570:         ((PetscInt*) faces)[0] = f;
2571:       }
2572:     }
2573:   }
2574: #endif
2575:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2576:   else if (numFaces == 1) {
2577:     /* Add the other cell neighbor for this face */
2578:     DMPlexSetCone(subdm, subcell, faces);
2579:   } else {
2580:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2581:     PetscBool posOriented;

2583:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2584:     origVertices        = &orientedVertices[numFaceVertices];
2585:     indices             = &orientedVertices[numFaceVertices*2];
2586:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2587:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2588:     /* TODO: I know that routine should return a permutation, not the indices */
2589:     for (v = 0; v < numFaceVertices; ++v) {
2590:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2591:       for (ov = 0; ov < numFaceVertices; ++ov) {
2592:         if (orientedVertices[ov] == vertex) {
2593:           orientedSubVertices[ov] = subvertex;
2594:           break;
2595:         }
2596:       }
2597:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2598:     }
2599:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2600:     DMPlexSetCone(subdm, subcell, newFacePoint);
2601:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2602:     ++(*newFacePoint);
2603:   }
2604: #if 0
2605:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2606: #else
2607:   DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2608: #endif
2609:   return(0);
2610: }

2612: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2613: {
2614:   MPI_Comm        comm;
2615:   DMLabel         subpointMap;
2616:   IS              subvertexIS,  subcellIS;
2617:   const PetscInt *subVertices, *subCells;
2618:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2619:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2620:   PetscInt        vStart, vEnd, c, f;
2621:   PetscErrorCode  ierr;

2624:   PetscObjectGetComm((PetscObject)dm,&comm);
2625:   /* Create subpointMap which marks the submesh */
2626:   DMLabelCreate("subpoint_map", &subpointMap);
2627:   DMPlexSetSubpointMap(subdm, subpointMap);
2628:   DMLabelDestroy(&subpointMap);
2629:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2630:   /* Setup chart */
2631:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2632:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2633:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2634:   DMPlexSetVTKCellHeight(subdm, 1);
2635:   /* Set cone sizes */
2636:   firstSubVertex = numSubCells;
2637:   firstSubFace   = numSubCells+numSubVertices;
2638:   newFacePoint   = firstSubFace;
2639:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2640:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2641:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2642:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2643:   for (c = 0; c < numSubCells; ++c) {
2644:     DMPlexSetConeSize(subdm, c, 1);
2645:   }
2646:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2647:     DMPlexSetConeSize(subdm, f, nFV);
2648:   }
2649:   DMSetUp(subdm);
2650:   /* Create face cones */
2651:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2652:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2653:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2654:   for (c = 0; c < numSubCells; ++c) {
2655:     const PetscInt cell    = subCells[c];
2656:     const PetscInt subcell = c;
2657:     PetscInt      *closure = NULL;
2658:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2660:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2661:     for (cl = 0; cl < closureSize*2; cl += 2) {
2662:       const PetscInt point = closure[cl];
2663:       PetscInt       subVertex;

2665:       if ((point >= vStart) && (point < vEnd)) {
2666:         ++numCorners;
2667:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2668:         if (subVertex >= 0) {
2669:           closure[faceSize] = point;
2670:           subface[faceSize] = firstSubVertex+subVertex;
2671:           ++faceSize;
2672:         }
2673:       }
2674:     }
2675:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2676:     if (faceSize == nFV) {
2677:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2678:     }
2679:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2680:   }
2681:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2682:   DMPlexSymmetrize(subdm);
2683:   DMPlexStratify(subdm);
2684:   /* Build coordinates */
2685:   {
2686:     PetscSection coordSection, subCoordSection;
2687:     Vec          coordinates, subCoordinates;
2688:     PetscScalar *coords, *subCoords;
2689:     PetscInt     numComp, coordSize, v;
2690:     const char  *name;

2692:     DMGetCoordinateSection(dm, &coordSection);
2693:     DMGetCoordinatesLocal(dm, &coordinates);
2694:     DMGetCoordinateSection(subdm, &subCoordSection);
2695:     PetscSectionSetNumFields(subCoordSection, 1);
2696:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2697:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2698:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2699:     for (v = 0; v < numSubVertices; ++v) {
2700:       const PetscInt vertex    = subVertices[v];
2701:       const PetscInt subvertex = firstSubVertex+v;
2702:       PetscInt       dof;

2704:       PetscSectionGetDof(coordSection, vertex, &dof);
2705:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2706:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2707:     }
2708:     PetscSectionSetUp(subCoordSection);
2709:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2710:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2711:     PetscObjectGetName((PetscObject)coordinates,&name);
2712:     PetscObjectSetName((PetscObject)subCoordinates,name);
2713:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2714:     VecSetType(subCoordinates,VECSTANDARD);
2715:     if (coordSize) {
2716:       VecGetArray(coordinates,    &coords);
2717:       VecGetArray(subCoordinates, &subCoords);
2718:       for (v = 0; v < numSubVertices; ++v) {
2719:         const PetscInt vertex    = subVertices[v];
2720:         const PetscInt subvertex = firstSubVertex+v;
2721:         PetscInt       dof, off, sdof, soff, d;

2723:         PetscSectionGetDof(coordSection, vertex, &dof);
2724:         PetscSectionGetOffset(coordSection, vertex, &off);
2725:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2726:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2727:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2728:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2729:       }
2730:       VecRestoreArray(coordinates,    &coords);
2731:       VecRestoreArray(subCoordinates, &subCoords);
2732:     }
2733:     DMSetCoordinatesLocal(subdm, subCoordinates);
2734:     VecDestroy(&subCoordinates);
2735:   }
2736:   /* Cleanup */
2737:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2738:   ISDestroy(&subvertexIS);
2739:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2740:   ISDestroy(&subcellIS);
2741:   return(0);
2742: }

2744: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2745: {
2746:   PetscInt       subPoint;

2749:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2750:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2751: }

2753: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2754: {
2755:   MPI_Comm         comm;
2756:   DMLabel          subpointMap;
2757:   IS              *subpointIS;
2758:   const PetscInt **subpoints;
2759:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2760:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2761:   PetscMPIInt      rank;
2762:   PetscErrorCode   ierr;

2765:   PetscObjectGetComm((PetscObject)dm,&comm);
2766:   MPI_Comm_rank(comm, &rank);
2767:   /* Create subpointMap which marks the submesh */
2768:   DMLabelCreate("subpoint_map", &subpointMap);
2769:   DMPlexSetSubpointMap(subdm, subpointMap);
2770:   if (cellHeight) {
2771:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2772:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2773:   } else {
2774:     DMLabel         depth;
2775:     IS              pointIS;
2776:     const PetscInt *points;
2777:     PetscInt        numPoints;

2779:     DMPlexGetDepthLabel(dm, &depth);
2780:     DMLabelGetStratumSize(label, value, &numPoints);
2781:     DMLabelGetStratumIS(label, value, &pointIS);
2782:     ISGetIndices(pointIS, &points);
2783:     for (p = 0; p < numPoints; ++p) {
2784:       PetscInt *closure = NULL;
2785:       PetscInt  closureSize, c, pdim;

2787:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2788:       for (c = 0; c < closureSize*2; c += 2) {
2789:         DMLabelGetValue(depth, closure[c], &pdim);
2790:         DMLabelSetValue(subpointMap, closure[c], pdim);
2791:       }
2792:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2793:     }
2794:     ISRestoreIndices(pointIS, &points);
2795:     ISDestroy(&pointIS);
2796:   }
2797:   DMLabelDestroy(&subpointMap);
2798:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2799:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2800:   cMax = (cMax < 0) ? cEnd : cMax;
2801:   /* Setup chart */
2802:   DMGetDimension(dm, &dim);
2803:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2804:   for (d = 0; d <= dim; ++d) {
2805:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2806:     totSubPoints += numSubPoints[d];
2807:   }
2808:   DMPlexSetChart(subdm, 0, totSubPoints);
2809:   DMPlexSetVTKCellHeight(subdm, cellHeight);
2810:   /* Set cone sizes */
2811:   firstSubPoint[dim] = 0;
2812:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2813:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2814:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2815:   for (d = 0; d <= dim; ++d) {
2816:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2817:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2818:   }
2819:   for (d = 0; d <= dim; ++d) {
2820:     for (p = 0; p < numSubPoints[d]; ++p) {
2821:       const PetscInt  point    = subpoints[d][p];
2822:       const PetscInt  subpoint = firstSubPoint[d] + p;
2823:       const PetscInt *cone;
2824:       PetscInt        coneSize, coneSizeNew, c, val;

2826:       DMPlexGetConeSize(dm, point, &coneSize);
2827:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2828:       if (cellHeight && (d == dim)) {
2829:         DMPlexGetCone(dm, point, &cone);
2830:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2831:           DMLabelGetValue(subpointMap, cone[c], &val);
2832:           if (val >= 0) coneSizeNew++;
2833:         }
2834:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2835:       }
2836:     }
2837:   }
2838:   DMSetUp(subdm);
2839:   /* Set cones */
2840:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2841:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2842:   for (d = 0; d <= dim; ++d) {
2843:     for (p = 0; p < numSubPoints[d]; ++p) {
2844:       const PetscInt  point    = subpoints[d][p];
2845:       const PetscInt  subpoint = firstSubPoint[d] + p;
2846:       const PetscInt *cone, *ornt;
2847:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

2849:       if (d == dim-1) {
2850:         const PetscInt *support, *cone, *ornt;
2851:         PetscInt        supportSize, coneSize, s, subc;

2853:         DMPlexGetSupport(dm, point, &support);
2854:         DMPlexGetSupportSize(dm, point, &supportSize);
2855:         for (s = 0; s < supportSize; ++s) {
2856:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2857:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2858:           if (subc >= 0) {
2859:             const PetscInt ccell = subpoints[d+1][subc];

2861:             DMPlexGetCone(dm, ccell, &cone);
2862:             DMPlexGetConeSize(dm, ccell, &coneSize);
2863:             DMPlexGetConeOrientation(dm, ccell, &ornt);
2864:             for (c = 0; c < coneSize; ++c) {
2865:               if (cone[c] == point) {
2866:                 fornt = ornt[c];
2867:                 break;
2868:               }
2869:             }
2870:             break;
2871:           }
2872:         }
2873:       }
2874:       DMPlexGetConeSize(dm, point, &coneSize);
2875:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2876:       DMPlexGetCone(dm, point, &cone);
2877:       DMPlexGetConeOrientation(dm, point, &ornt);
2878:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2879:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2880:         if (subc >= 0) {
2881:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2882:           orntNew[coneSizeNew] = ornt[c];
2883:           ++coneSizeNew;
2884:         }
2885:       }
2886:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2887:       if (fornt < 0) {
2888:         /* This should be replaced by a call to DMPlexReverseCell() */
2889: #if 0
2890:         DMPlexReverseCell(subdm, subpoint);
2891: #else
2892:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
2893:           PetscInt faceSize, tmp;

2895:           tmp        = coneNew[c];
2896:           coneNew[c] = coneNew[coneSizeNew-1-c];
2897:           coneNew[coneSizeNew-1-c] = tmp;
2898:           DMPlexGetConeSize(dm, cone[c], &faceSize);
2899:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2900:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2901:           orntNew[coneSizeNew-1-c] = tmp;
2902:         }
2903:       }
2904:       DMPlexSetCone(subdm, subpoint, coneNew);
2905:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2906: #endif
2907:     }
2908:   }
2909:   PetscFree2(coneNew,orntNew);
2910:   DMPlexSymmetrize(subdm);
2911:   DMPlexStratify(subdm);
2912:   /* Build coordinates */
2913:   {
2914:     PetscSection coordSection, subCoordSection;
2915:     Vec          coordinates, subCoordinates;
2916:     PetscScalar *coords, *subCoords;
2917:     PetscInt     cdim, numComp, coordSize;
2918:     const char  *name;

2920:     DMGetCoordinateDim(dm, &cdim);
2921:     DMGetCoordinateSection(dm, &coordSection);
2922:     DMGetCoordinatesLocal(dm, &coordinates);
2923:     DMGetCoordinateSection(subdm, &subCoordSection);
2924:     PetscSectionSetNumFields(subCoordSection, 1);
2925:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2926:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2927:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2928:     for (v = 0; v < numSubPoints[0]; ++v) {
2929:       const PetscInt vertex    = subpoints[0][v];
2930:       const PetscInt subvertex = firstSubPoint[0]+v;
2931:       PetscInt       dof;

2933:       PetscSectionGetDof(coordSection, vertex, &dof);
2934:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2935:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2936:     }
2937:     PetscSectionSetUp(subCoordSection);
2938:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2939:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2940:     PetscObjectGetName((PetscObject)coordinates,&name);
2941:     PetscObjectSetName((PetscObject)subCoordinates,name);
2942:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2943:     VecSetBlockSize(subCoordinates, cdim);
2944:     VecSetType(subCoordinates,VECSTANDARD);
2945:     VecGetArray(coordinates,    &coords);
2946:     VecGetArray(subCoordinates, &subCoords);
2947:     for (v = 0; v < numSubPoints[0]; ++v) {
2948:       const PetscInt vertex    = subpoints[0][v];
2949:       const PetscInt subvertex = firstSubPoint[0]+v;
2950:       PetscInt dof, off, sdof, soff, d;

2952:       PetscSectionGetDof(coordSection, vertex, &dof);
2953:       PetscSectionGetOffset(coordSection, vertex, &off);
2954:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2955:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2956:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2957:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2958:     }
2959:     VecRestoreArray(coordinates,    &coords);
2960:     VecRestoreArray(subCoordinates, &subCoords);
2961:     DMSetCoordinatesLocal(subdm, subCoordinates);
2962:     VecDestroy(&subCoordinates);
2963:   }
2964:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
2965:   {
2966:     PetscSF            sfPoint, sfPointSub;
2967:     IS                 subpIS;
2968:     const PetscSFNode *remotePoints;
2969:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
2970:     const PetscInt    *localPoints, *subpoints;
2971:     PetscInt          *slocalPoints;
2972:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
2973:     PetscMPIInt        rank;

2975:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2976:     DMGetPointSF(dm, &sfPoint);
2977:     DMGetPointSF(subdm, &sfPointSub);
2978:     DMPlexGetChart(dm, &pStart, &pEnd);
2979:     DMPlexGetChart(subdm, NULL, &numSubroots);
2980:     DMPlexCreateSubpointIS(subdm, &subpIS);
2981:     if (subpIS) {
2982:       ISGetIndices(subpIS, &subpoints);
2983:       ISGetLocalSize(subpIS, &numSubpoints);
2984:     }
2985:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2986:     if (numRoots >= 0) {
2987:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2988:       for (p = 0; p < pEnd-pStart; ++p) {
2989:         newLocalPoints[p].rank  = -2;
2990:         newLocalPoints[p].index = -2;
2991:       }
2992:       /* Set subleaves */
2993:       for (l = 0; l < numLeaves; ++l) {
2994:         const PetscInt point    = localPoints[l];
2995:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

2997:         if (subpoint < 0) continue;
2998:         newLocalPoints[point-pStart].rank  = rank;
2999:         newLocalPoints[point-pStart].index = subpoint;
3000:         ++numSubleaves;
3001:       }
3002:       /* Must put in owned subpoints */
3003:       for (p = pStart; p < pEnd; ++p) {
3004:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3006:         if (subpoint < 0) {
3007:           newOwners[p-pStart].rank  = -3;
3008:           newOwners[p-pStart].index = -3;
3009:         } else {
3010:           newOwners[p-pStart].rank  = rank;
3011:           newOwners[p-pStart].index = subpoint;
3012:         }
3013:       }
3014:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3015:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3016:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3017:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3018:       PetscMalloc1(numSubleaves, &slocalPoints);
3019:       PetscMalloc1(numSubleaves, &sremotePoints);
3020:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3021:         const PetscInt point    = localPoints[l];
3022:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3024:         if (subpoint < 0) continue;
3025:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3026:         slocalPoints[sl]        = subpoint;
3027:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3028:         sremotePoints[sl].index = newLocalPoints[point].index;
3029:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3030:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3031:         ++sl;
3032:       }
3033:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3034:       PetscFree2(newLocalPoints,newOwners);
3035:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3036:     }
3037:     if (subpIS) {
3038:       ISRestoreIndices(subpIS, &subpoints);
3039:       ISDestroy(&subpIS);
3040:     }
3041:   }
3042:   /* Cleanup */
3043:   for (d = 0; d <= dim; ++d) {
3044:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3045:     ISDestroy(&subpointIS[d]);
3046:   }
3047:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3048:   return(0);
3049: }

3051: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3052: {

3056:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3057:   return(0);
3058: }

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

3063:   Input Parameters:
3064: + dm           - The original mesh
3065: . vertexLabel  - The DMLabel marking vertices contained in the surface
3066: - value        - The label value to use

3068:   Output Parameter:
3069: . subdm - The surface mesh

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

3073:   Level: developer

3075: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3076: @*/
3077: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3078: {
3079:   PetscInt       dim, cdim, depth;

3085:   DMGetDimension(dm, &dim);
3086:   DMPlexGetDepth(dm, &depth);
3087:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3088:   DMSetType(*subdm, DMPLEX);
3089:   DMSetDimension(*subdm, dim-1);
3090:   DMGetCoordinateDim(dm, &cdim);
3091:   DMSetCoordinateDim(*subdm, cdim);
3092:   if (depth == dim) {
3093:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3094:   } else {
3095:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3096:   }
3097:   return(0);
3098: }

3100: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3101: {
3102:   MPI_Comm        comm;
3103:   DMLabel         subpointMap;
3104:   IS              subvertexIS;
3105:   const PetscInt *subVertices;
3106:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3107:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3108:   PetscInt        cMax, c, f;
3109:   PetscErrorCode  ierr;

3112:   PetscObjectGetComm((PetscObject)dm, &comm);
3113:   /* Create subpointMap which marks the submesh */
3114:   DMLabelCreate("subpoint_map", &subpointMap);
3115:   DMPlexSetSubpointMap(subdm, subpointMap);
3116:   DMLabelDestroy(&subpointMap);
3117:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3118:   /* Setup chart */
3119:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3120:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3121:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3122:   DMPlexSetVTKCellHeight(subdm, 1);
3123:   /* Set cone sizes */
3124:   firstSubVertex = numSubCells;
3125:   firstSubFace   = numSubCells+numSubVertices;
3126:   newFacePoint   = firstSubFace;
3127:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3128:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3129:   for (c = 0; c < numSubCells; ++c) {
3130:     DMPlexSetConeSize(subdm, c, 1);
3131:   }
3132:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3133:     DMPlexSetConeSize(subdm, f, nFV);
3134:   }
3135:   DMSetUp(subdm);
3136:   /* Create face cones */
3137:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3138:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3139:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3140:   for (c = 0; c < numSubCells; ++c) {
3141:     const PetscInt  cell    = subCells[c];
3142:     const PetscInt  subcell = c;
3143:     const PetscInt *cone, *cells;
3144:     PetscInt        numCells, subVertex, p, v;

3146:     if (cell < cMax) continue;
3147:     DMPlexGetCone(dm, cell, &cone);
3148:     for (v = 0; v < nFV; ++v) {
3149:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3150:       subface[v] = firstSubVertex+subVertex;
3151:     }
3152:     DMPlexSetCone(subdm, newFacePoint, subface);
3153:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3154:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3155:     /* Not true in parallel
3156:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3157:     for (p = 0; p < numCells; ++p) {
3158:       PetscInt negsubcell;

3160:       if (cells[p] >= cMax) continue;
3161:       /* I know this is a crap search */
3162:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3163:         if (subCells[negsubcell] == cells[p]) break;
3164:       }
3165:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3166:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3167:     }
3168:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3169:     ++newFacePoint;
3170:   }
3171:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3172:   DMPlexSymmetrize(subdm);
3173:   DMPlexStratify(subdm);
3174:   /* Build coordinates */
3175:   {
3176:     PetscSection coordSection, subCoordSection;
3177:     Vec          coordinates, subCoordinates;
3178:     PetscScalar *coords, *subCoords;
3179:     PetscInt     cdim, numComp, coordSize, v;
3180:     const char  *name;

3182:     DMGetCoordinateDim(dm, &cdim);
3183:     DMGetCoordinateSection(dm, &coordSection);
3184:     DMGetCoordinatesLocal(dm, &coordinates);
3185:     DMGetCoordinateSection(subdm, &subCoordSection);
3186:     PetscSectionSetNumFields(subCoordSection, 1);
3187:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3188:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3189:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3190:     for (v = 0; v < numSubVertices; ++v) {
3191:       const PetscInt vertex    = subVertices[v];
3192:       const PetscInt subvertex = firstSubVertex+v;
3193:       PetscInt       dof;

3195:       PetscSectionGetDof(coordSection, vertex, &dof);
3196:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3197:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3198:     }
3199:     PetscSectionSetUp(subCoordSection);
3200:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3201:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3202:     PetscObjectGetName((PetscObject)coordinates,&name);
3203:     PetscObjectSetName((PetscObject)subCoordinates,name);
3204:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3205:     VecSetBlockSize(subCoordinates, cdim);
3206:     VecSetType(subCoordinates,VECSTANDARD);
3207:     VecGetArray(coordinates,    &coords);
3208:     VecGetArray(subCoordinates, &subCoords);
3209:     for (v = 0; v < numSubVertices; ++v) {
3210:       const PetscInt vertex    = subVertices[v];
3211:       const PetscInt subvertex = firstSubVertex+v;
3212:       PetscInt       dof, off, sdof, soff, d;

3214:       PetscSectionGetDof(coordSection, vertex, &dof);
3215:       PetscSectionGetOffset(coordSection, vertex, &off);
3216:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3217:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3218:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3219:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3220:     }
3221:     VecRestoreArray(coordinates,    &coords);
3222:     VecRestoreArray(subCoordinates, &subCoords);
3223:     DMSetCoordinatesLocal(subdm, subCoordinates);
3224:     VecDestroy(&subCoordinates);
3225:   }
3226:   /* Build SF */
3227:   CHKMEMQ;
3228:   {
3229:     PetscSF            sfPoint, sfPointSub;
3230:     const PetscSFNode *remotePoints;
3231:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3232:     const PetscInt    *localPoints;
3233:     PetscInt          *slocalPoints;
3234:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3235:     PetscMPIInt        rank;

3237:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3238:     DMGetPointSF(dm, &sfPoint);
3239:     DMGetPointSF(subdm, &sfPointSub);
3240:     DMPlexGetChart(dm, &pStart, &pEnd);
3241:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3242:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3243:     if (numRoots >= 0) {
3244:       /* Only vertices should be shared */
3245:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3246:       for (p = 0; p < pEnd-pStart; ++p) {
3247:         newLocalPoints[p].rank  = -2;
3248:         newLocalPoints[p].index = -2;
3249:       }
3250:       /* Set subleaves */
3251:       for (l = 0; l < numLeaves; ++l) {
3252:         const PetscInt point    = localPoints[l];
3253:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3255:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3256:         if (subPoint < 0) continue;
3257:         newLocalPoints[point-pStart].rank  = rank;
3258:         newLocalPoints[point-pStart].index = subPoint;
3259:         ++numSubLeaves;
3260:       }
3261:       /* Must put in owned subpoints */
3262:       for (p = pStart; p < pEnd; ++p) {
3263:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3265:         if (subPoint < 0) {
3266:           newOwners[p-pStart].rank  = -3;
3267:           newOwners[p-pStart].index = -3;
3268:         } else {
3269:           newOwners[p-pStart].rank  = rank;
3270:           newOwners[p-pStart].index = subPoint;
3271:         }
3272:       }
3273:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3274:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3275:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3276:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3277:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3278:       PetscMalloc1(numSubLeaves, &sremotePoints);
3279:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3280:         const PetscInt point    = localPoints[l];
3281:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3283:         if (subPoint < 0) continue;
3284:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3285:         slocalPoints[sl]        = subPoint;
3286:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3287:         sremotePoints[sl].index = newLocalPoints[point].index;
3288:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3289:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3290:         ++sl;
3291:       }
3292:       PetscFree2(newLocalPoints,newOwners);
3293:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3294:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3295:     }
3296:   }
3297:   CHKMEMQ;
3298:   /* Cleanup */
3299:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3300:   ISDestroy(&subvertexIS);
3301:   PetscFree(subCells);
3302:   return(0);
3303: }

3305: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3306: {
3307:   DMLabel        label = NULL;

3311:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3312:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3313:   return(0);
3314: }

3316: /*@C
3317:   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.

3319:   Input Parameters:
3320: + dm          - The original mesh
3321: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3322: . label       - A label name, or NULL
3323: - value  - A label value

3325:   Output Parameter:
3326: . subdm - The surface mesh

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

3330:   Level: developer

3332: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3333: @*/
3334: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3335: {
3336:   PetscInt       dim, cdim, depth;

3342:   DMGetDimension(dm, &dim);
3343:   DMPlexGetDepth(dm, &depth);
3344:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3345:   DMSetType(*subdm, DMPLEX);
3346:   DMSetDimension(*subdm, dim-1);
3347:   DMGetCoordinateDim(dm, &cdim);
3348:   DMSetCoordinateDim(*subdm, cdim);
3349:   if (depth == dim) {
3350:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3351:   } else {
3352:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3353:   }
3354:   return(0);
3355: }

3357: /*@
3358:   DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh

3360:   Input Parameters:
3361: + dm        - The original mesh
3362: . cellLabel - The DMLabel marking cells contained in the new mesh
3363: - value     - The label value to use

3365:   Output Parameter:
3366: . subdm - The new mesh

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

3370:   Level: developer

3372: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3373: @*/
3374: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3375: {
3376:   PetscInt       dim;

3382:   DMGetDimension(dm, &dim);
3383:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3384:   DMSetType(*subdm, DMPLEX);
3385:   DMSetDimension(*subdm, dim);
3386:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3387:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3388:   return(0);
3389: }

3391: /*@
3392:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3394:   Input Parameter:
3395: . dm - The submesh DM

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

3400:   Level: developer

3402: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3403: @*/
3404: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3405: {
3409:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3410:   return(0);
3411: }

3413: /*@
3414:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3416:   Input Parameters:
3417: + dm - The submesh DM
3418: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

3420:   Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()

3422:   Level: developer

3424: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3425: @*/
3426: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3427: {
3428:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3429:   DMLabel        tmp;

3434:   tmp  = mesh->subpointMap;
3435:   mesh->subpointMap = subpointMap;
3436:   ++mesh->subpointMap->refct;
3437:   DMLabelDestroy(&tmp);
3438:   return(0);
3439: }

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

3444:   Input Parameter:
3445: . dm - The submesh DM

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

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

3452:   Level: developer

3454: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3455: @*/
3456: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3457: {
3458:   MPI_Comm        comm;
3459:   DMLabel         subpointMap;
3460:   IS              is;
3461:   const PetscInt *opoints;
3462:   PetscInt       *points, *depths;
3463:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3464:   PetscErrorCode  ierr;

3469:   PetscObjectGetComm((PetscObject)dm,&comm);
3470:   *subpointIS = NULL;
3471:   DMPlexGetSubpointMap(dm, &subpointMap);
3472:   DMPlexGetDepth(dm, &depth);
3473:   if (subpointMap && depth >= 0) {
3474:     DMPlexGetChart(dm, &pStart, &pEnd);
3475:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3476:     DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3477:     depths[0] = depth;
3478:     depths[1] = 0;
3479:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3480:     PetscMalloc1(pEnd, &points);
3481:     for(d = 0, off = 0; d <= depth; ++d) {
3482:       const PetscInt dep = depths[d];

3484:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3485:       DMLabelGetStratumSize(subpointMap, dep, &n);
3486:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3487:         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);
3488:       } else {
3489:         if (!n) {
3490:           if (d == 0) {
3491:             /* Missing cells */
3492:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3493:           } else {
3494:             /* Missing faces */
3495:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3496:           }
3497:         }
3498:       }
3499:       if (n) {
3500:         DMLabelGetStratumIS(subpointMap, dep, &is);
3501:         ISGetIndices(is, &opoints);
3502:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3503:         ISRestoreIndices(is, &opoints);
3504:         ISDestroy(&is);
3505:       }
3506:     }
3507:     DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3508:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3509:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3510:   }
3511:   return(0);
3512: }