Actual source code: plexsubmesh.c

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

  5: static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid)
  6: {
  7:   DMPolytopeType ct;

 11:   DMPlexGetCellType(dm, p, &ct);
 12:   switch (ct) {
 13:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
 14:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
 15:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
 16:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
 17:       *isHybrid = PETSC_TRUE;
 18:     default: *isHybrid = PETSC_FALSE;
 19:   }
 20:   return(0);
 21: }

 23: static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
 24: {
 25:   DMLabel        ctLabel;

 29:   if (cStart) *cStart = -1;
 30:   if (cEnd)   *cEnd   = -1;
 31:   DMPlexGetCellTypeLabel(dm, &ctLabel);
 32:   switch (dim) {
 33:     case 1: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd);break;
 34:     case 2: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd);break;
 35:     case 3:
 36:       DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd);
 37:       if (*cStart < 0) {DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd);}
 38:       break;
 39:     default: return(0);
 40:   }
 41:   return(0);
 42: }

 44: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
 45: {
 46:   PetscInt       fStart, fEnd, f;

 50:   DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
 51:   for (f = fStart; f < fEnd; ++f) {
 52:     PetscInt supportSize;

 54:     DMPlexGetSupportSize(dm, f, &supportSize);
 55:     if (supportSize == 1) {
 56:       if (val < 0) {
 57:         PetscInt *closure = NULL;
 58:         PetscInt  clSize, cl, cval;

 60:         DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 61:         for (cl = 0; cl < clSize*2; cl += 2) {
 62:           DMLabelGetValue(label, closure[cl], &cval);
 63:           if (cval < 0) continue;
 64:           DMLabelSetValue(label, f, cval);
 65:           break;
 66:         }
 67:         if (cl == clSize*2) {DMLabelSetValue(label, f, 1);}
 68:         DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 69:       } else {
 70:         DMLabelSetValue(label, f, val);
 71:       }
 72:     }
 73:   }
 74:   return(0);
 75: }

 77: /*@
 78:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 80:   Not Collective

 82:   Input Parameter:
 83: + dm - The original DM
 84: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)

 86:   Output Parameter:
 87: . label - The DMLabel marking boundary faces with the given value

 89:   Level: developer

 91: .seealso: DMLabelCreate(), DMCreateLabel()
 92: @*/
 93: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
 94: {
 95:   DMPlexInterpolatedFlag  flg;
 96:   PetscErrorCode          ierr;

100:   DMPlexIsInterpolated(dm, &flg);
101:   if (flg != DMPLEX_INTERPOLATED_FULL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not fully interpolated on this rank");
102:   DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
103:   return(0);
104: }

106: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
107: {
108:   IS              valueIS;
109:   PetscSF         sfPoint;
110:   const PetscInt *values;
111:   PetscInt        numValues, v, cStart, cEnd, nroots;
112:   PetscErrorCode  ierr;

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

124:     DMLabelGetStratumSize(label, values[v], &numPoints);
125:     DMLabelGetStratumIS(label, values[v], &pointIS);
126:     ISGetIndices(pointIS, &points);
127:     for (p = 0; p < numPoints; ++p) {
128:       PetscInt  q = points[p];
129:       PetscInt *closure = NULL;
130:       PetscInt  closureSize, c;

132:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
133:         continue;
134:       }
135:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
136:       for (c = 0; c < closureSize*2; c += 2) {
137:         DMLabelSetValue(label, closure[c], values[v]);
138:       }
139:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
140:     }
141:     ISRestoreIndices(pointIS, &points);
142:     ISDestroy(&pointIS);
143:   }
144:   ISRestoreIndices(valueIS, &values);
145:   ISDestroy(&valueIS);
146:   DMGetPointSF(dm, &sfPoint);
147:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
148:   if (nroots >= 0) {
149:     DMLabel         lblRoots, lblLeaves;
150:     IS              valueIS, pointIS;
151:     const PetscInt *values;
152:     PetscInt        numValues, v;
153:     PetscErrorCode  ierr;

155:     DMGetPointSF(dm, &sfPoint);
156:     /* Pull point contributions from remote leaves into local roots */
157:     DMLabelGather(label, sfPoint, &lblLeaves);
158:     DMLabelGetValueIS(lblLeaves, &valueIS);
159:     ISGetLocalSize(valueIS, &numValues);
160:     ISGetIndices(valueIS, &values);
161:     for (v = 0; v < numValues; ++v) {
162:       const PetscInt value = values[v];

164:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
165:       DMLabelInsertIS(label, pointIS, value);
166:       ISDestroy(&pointIS);
167:     }
168:     ISRestoreIndices(valueIS, &values);
169:     ISDestroy(&valueIS);
170:     DMLabelDestroy(&lblLeaves);
171:     /* Push point contributions from roots into remote leaves */
172:     DMLabelDistribute(label, sfPoint, &lblRoots);
173:     DMLabelGetValueIS(lblRoots, &valueIS);
174:     ISGetLocalSize(valueIS, &numValues);
175:     ISGetIndices(valueIS, &values);
176:     for (v = 0; v < numValues; ++v) {
177:       const PetscInt value = values[v];

179:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
180:       DMLabelInsertIS(label, pointIS, value);
181:       ISDestroy(&pointIS);
182:     }
183:     ISRestoreIndices(valueIS, &values);
184:     ISDestroy(&valueIS);
185:     DMLabelDestroy(&lblRoots);
186:   }
187:   return(0);
188: }

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

193:   Input Parameters:
194: + dm - The DM
195: - label - A DMLabel marking the surface points

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

200:   Level: developer

202: .seealso: DMPlexLabelCohesiveComplete()
203: @*/
204: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
205: {

209:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
210:   return(0);
211: }

213: /*@
214:   DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point

216:   Input Parameters:
217: + dm - The DM
218: - label - A DMLabel marking the surface points

220:   Output Parameter:
221: . label - A DMLabel incorporating cells

223:   Level: developer

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

227: .seealso: DMPlexLabelAddFaceCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
228: @*/
229: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
230: {
231:   IS              valueIS;
232:   const PetscInt *values;
233:   PetscInt        numValues, v, cStart, cEnd;
234:   PetscErrorCode  ierr;

237:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
238:   DMLabelGetNumValues(label, &numValues);
239:   DMLabelGetValueIS(label, &valueIS);
240:   ISGetIndices(valueIS, &values);
241:   for (v = 0; v < numValues; ++v) {
242:     IS              pointIS;
243:     const PetscInt *points;
244:     PetscInt        numPoints, p;

246:     DMLabelGetStratumSize(label, values[v], &numPoints);
247:     DMLabelGetStratumIS(label, values[v], &pointIS);
248:     ISGetIndices(pointIS, &points);
249:     for (p = 0; p < numPoints; ++p) {
250:       PetscInt *closure = NULL;
251:       PetscInt  closureSize, cl;

253:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
254:       for (cl = closureSize-1; cl > 0; --cl) {
255:         const PetscInt cell = closure[cl*2];
256:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
257:       }
258:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
259:     }
260:     ISRestoreIndices(pointIS, &points);
261:     ISDestroy(&pointIS);
262:   }
263:   ISRestoreIndices(valueIS, &values);
264:   ISDestroy(&valueIS);
265:   return(0);
266: }

268: /*@
269:   DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face

271:   Input Parameters:
272: + dm - The DM
273: - label - A DMLabel marking the surface points

275:   Output Parameter:
276: . label - A DMLabel incorporating cells

278:   Level: developer

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

282: .seealso: DMPlexLabelAddCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
283: @*/
284: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
285: {
286:   IS              valueIS;
287:   const PetscInt *values;
288:   PetscInt        numValues, v, cStart, cEnd, fStart, fEnd;
289:   PetscErrorCode  ierr;

292:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
293:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
294:   DMLabelGetNumValues(label, &numValues);
295:   DMLabelGetValueIS(label, &valueIS);
296:   ISGetIndices(valueIS, &values);
297:   for (v = 0; v < numValues; ++v) {
298:     IS              pointIS;
299:     const PetscInt *points;
300:     PetscInt        numPoints, p;

302:     DMLabelGetStratumSize(label, values[v], &numPoints);
303:     DMLabelGetStratumIS(label, values[v], &pointIS);
304:     ISGetIndices(pointIS, &points);
305:     for (p = 0; p < numPoints; ++p) {
306:       const PetscInt face = points[p];
307:       PetscInt      *closure = NULL;
308:       PetscInt       closureSize, cl;

310:       if ((face < fStart) || (face >= fEnd)) continue;
311:       DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
312:       for (cl = closureSize-1; cl > 0; --cl) {
313:         const PetscInt cell = closure[cl*2];
314:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
315:       }
316:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
317:     }
318:     ISRestoreIndices(pointIS, &points);
319:     ISDestroy(&pointIS);
320:   }
321:   ISRestoreIndices(valueIS, &values);
322:   ISDestroy(&valueIS);
323:   return(0);
324: }

326: /*@
327:   DMPlexLabelClearCells - Remove cells from a label

329:   Input Parameters:
330: + dm - The DM
331: - label - A DMLabel marking surface points and their adjacent cells

333:   Output Parameter:
334: . label - A DMLabel without cells

336:   Level: developer

338:   Note: This undoes DMPlexLabelAddCells() or DMPlexLabelAddFaceCells()

340: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
341: @*/
342: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
343: {
344:   IS              valueIS;
345:   const PetscInt *values;
346:   PetscInt        numValues, v, cStart, cEnd;
347:   PetscErrorCode  ierr;

350:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
351:   DMLabelGetNumValues(label, &numValues);
352:   DMLabelGetValueIS(label, &valueIS);
353:   ISGetIndices(valueIS, &values);
354:   for (v = 0; v < numValues; ++v) {
355:     IS              pointIS;
356:     const PetscInt *points;
357:     PetscInt        numPoints, p;

359:     DMLabelGetStratumSize(label, values[v], &numPoints);
360:     DMLabelGetStratumIS(label, values[v], &pointIS);
361:     ISGetIndices(pointIS, &points);
362:     for (p = 0; p < numPoints; ++p) {
363:       PetscInt point = points[p];

365:       if (point >= cStart && point < cEnd) {
366:         DMLabelClearValue(label,point,values[v]);
367:       }
368:     }
369:     ISRestoreIndices(pointIS, &points);
370:     ISDestroy(&pointIS);
371:   }
372:   ISRestoreIndices(valueIS, &values);
373:   ISDestroy(&valueIS);
374:   return(0);
375: }

377: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
378:  * index (skipping first, which is (0,0)) */
379: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
380: {
381:   PetscInt d, off = 0;

384:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
385:   for (d = 0; d < depth; d++) {
386:     PetscInt firstd = d;
387:     PetscInt firstStart = depthShift[2*d];
388:     PetscInt e;

390:     for (e = d+1; e <= depth; e++) {
391:       if (depthShift[2*e] < firstStart) {
392:         firstd = e;
393:         firstStart = depthShift[2*d];
394:       }
395:     }
396:     if (firstd != d) {
397:       PetscInt swap[2];

399:       e = firstd;
400:       swap[0] = depthShift[2*d];
401:       swap[1] = depthShift[2*d+1];
402:       depthShift[2*d]   = depthShift[2*e];
403:       depthShift[2*d+1] = depthShift[2*e+1];
404:       depthShift[2*e]   = swap[0];
405:       depthShift[2*e+1] = swap[1];
406:     }
407:   }
408:   /* convert (oldstart, added) to (oldstart, newstart) */
409:   for (d = 0; d <= depth; d++) {
410:     off += depthShift[2*d+1];
411:     depthShift[2*d+1] = depthShift[2*d] + off;
412:   }
413:   return(0);
414: }

416: /* depthShift is a list of (old, new) pairs */
417: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
418: {
419:   PetscInt d;
420:   PetscInt newOff = 0;

422:   for (d = 0; d <= depth; d++) {
423:     if (p < depthShift[2*d]) return p + newOff;
424:     else newOff = depthShift[2*d+1] - depthShift[2*d];
425:   }
426:   return p + newOff;
427: }

429: /* depthShift is a list of (old, new) pairs */
430: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
431: {
432:   PetscInt d;
433:   PetscInt newOff = 0;

435:   for (d = 0; d <= depth; d++) {
436:     if (p < depthShift[2*d+1]) return p + newOff;
437:     else newOff = depthShift[2*d] - depthShift[2*d+1];
438:   }
439:   return p + newOff;
440: }

442: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
443: {
444:   PetscInt       depth = 0, d, pStart, pEnd, p;
445:   DMLabel        depthLabel;

449:   DMPlexGetDepth(dm, &depth);
450:   if (depth < 0) return(0);
451:   /* Step 1: Expand chart */
452:   DMPlexGetChart(dm, &pStart, &pEnd);
453:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
454:   DMPlexSetChart(dmNew, pStart, pEnd);
455:   DMCreateLabel(dmNew,"depth");
456:   DMPlexGetDepthLabel(dmNew,&depthLabel);
457:   DMCreateLabel(dmNew, "celltype");
458:   /* Step 2: Set cone and support sizes */
459:   for (d = 0; d <= depth; ++d) {
460:     PetscInt pStartNew, pEndNew;
461:     IS pIS;

463:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
464:     pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
465:     pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
466:     ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
467:     DMLabelSetStratumIS(depthLabel, d, pIS);
468:     ISDestroy(&pIS);
469:     for (p = pStart; p < pEnd; ++p) {
470:       PetscInt       newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
471:       PetscInt       size;
472:       DMPolytopeType ct;

474:       DMPlexGetConeSize(dm, p, &size);
475:       DMPlexSetConeSize(dmNew, newp, size);
476:       DMPlexGetSupportSize(dm, p, &size);
477:       DMPlexSetSupportSize(dmNew, newp, size);
478:       DMPlexGetCellType(dm, p, &ct);
479:       DMPlexSetCellType(dmNew, newp, ct);
480:     }
481:   }
482:   return(0);
483: }

485: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
486: {
487:   PetscInt      *newpoints;
488:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

492:   DMPlexGetDepth(dm, &depth);
493:   if (depth < 0) return(0);
494:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
495:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
496:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
497:   /* Step 5: Set cones and supports */
498:   DMPlexGetChart(dm, &pStart, &pEnd);
499:   for (p = pStart; p < pEnd; ++p) {
500:     const PetscInt *points = NULL, *orientations = NULL;
501:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

503:     DMPlexGetConeSize(dm, p, &size);
504:     DMPlexGetCone(dm, p, &points);
505:     DMPlexGetConeOrientation(dm, p, &orientations);
506:     for (i = 0; i < size; ++i) {
507:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
508:     }
509:     DMPlexSetCone(dmNew, newp, newpoints);
510:     DMPlexSetConeOrientation(dmNew, newp, orientations);
511:     DMPlexGetSupportSize(dm, p, &size);
512:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
513:     DMPlexGetSupport(dm, p, &points);
514:     for (i = 0; i < size; ++i) {
515:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
516:     }
517:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
518:     DMPlexSetSupport(dmNew, newp, newpoints);
519:   }
520:   PetscFree(newpoints);
521:   return(0);
522: }

524: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
525: {
526:   PetscSection   coordSection, newCoordSection;
527:   Vec            coordinates, newCoordinates;
528:   PetscScalar   *coords, *newCoords;
529:   PetscInt       coordSize, sStart, sEnd;
530:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
531:   PetscBool      hasCells;

535:   DMGetCoordinateDim(dm, &dim);
536:   DMSetCoordinateDim(dmNew, dim);
537:   DMPlexGetDepth(dm, &depth);
538:   /* Step 8: Convert coordinates */
539:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
540:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
541:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
542:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
543:   DMGetCoordinateSection(dm, &coordSection);
544:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
545:   PetscSectionSetNumFields(newCoordSection, 1);
546:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
547:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
548:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
549:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
550:   if (hasCells) {
551:     for (c = cStart; c < cEnd; ++c) {
552:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

554:       PetscSectionGetDof(coordSection, c, &dof);
555:       PetscSectionSetDof(newCoordSection, cNew, dof);
556:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
557:     }
558:   }
559:   for (v = vStartNew; v < vEndNew; ++v) {
560:     PetscSectionSetDof(newCoordSection, v, dim);
561:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
562:   }
563:   PetscSectionSetUp(newCoordSection);
564:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
565:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
566:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
567:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
568:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
569:   VecSetBlockSize(newCoordinates, dim);
570:   VecSetType(newCoordinates,VECSTANDARD);
571:   DMSetCoordinatesLocal(dmNew, newCoordinates);
572:   DMGetCoordinatesLocal(dm, &coordinates);
573:   VecGetArray(coordinates, &coords);
574:   VecGetArray(newCoordinates, &newCoords);
575:   if (hasCells) {
576:     for (c = cStart; c < cEnd; ++c) {
577:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

579:       PetscSectionGetDof(coordSection, c, &dof);
580:       PetscSectionGetOffset(coordSection, c, &off);
581:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
582:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
583:     }
584:   }
585:   for (v = vStart; v < vEnd; ++v) {
586:     PetscInt dof, off, noff, d;

588:     PetscSectionGetDof(coordSection, v, &dof);
589:     PetscSectionGetOffset(coordSection, v, &off);
590:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
591:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
592:   }
593:   VecRestoreArray(coordinates, &coords);
594:   VecRestoreArray(newCoordinates, &newCoords);
595:   VecDestroy(&newCoordinates);
596:   PetscSectionDestroy(&newCoordSection);
597:   return(0);
598: }

600: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
601: {
602:   PetscInt           depth = 0;
603:   PetscSF            sfPoint, sfPointNew;
604:   const PetscSFNode *remotePoints;
605:   PetscSFNode       *gremotePoints;
606:   const PetscInt    *localPoints;
607:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
608:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
609:   PetscErrorCode     ierr;

612:   DMPlexGetDepth(dm, &depth);
613:   /* Step 9: Convert pointSF */
614:   DMGetPointSF(dm, &sfPoint);
615:   DMGetPointSF(dmNew, &sfPointNew);
616:   DMPlexGetChart(dm, &pStart, &pEnd);
617:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
618:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
619:   if (numRoots >= 0) {
620:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
621:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
622:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
623:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
624:     PetscMalloc1(numLeaves,    &glocalPoints);
625:     PetscMalloc1(numLeaves, &gremotePoints);
626:     for (l = 0; l < numLeaves; ++l) {
627:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
628:       gremotePoints[l].rank  = remotePoints[l].rank;
629:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
630:     }
631:     PetscFree2(newLocation,newRemoteLocation);
632:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
633:   }
634:   return(0);
635: }

637: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
638: {
639:   PetscSF            sfPoint;
640:   DMLabel            vtkLabel, ghostLabel;
641:   const PetscSFNode *leafRemote;
642:   const PetscInt    *leafLocal;
643:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
644:   PetscMPIInt        rank;
645:   PetscErrorCode     ierr;

648:   DMPlexGetDepth(dm, &depth);
649:   /* Step 10: Convert labels */
650:   DMGetNumLabels(dm, &numLabels);
651:   for (l = 0; l < numLabels; ++l) {
652:     DMLabel         label, newlabel;
653:     const char     *lname;
654:     PetscBool       isDepth, isDim;
655:     IS              valueIS;
656:     const PetscInt *values;
657:     PetscInt        numValues, val;

659:     DMGetLabelName(dm, l, &lname);
660:     PetscStrcmp(lname, "depth", &isDepth);
661:     if (isDepth) continue;
662:     PetscStrcmp(lname, "dim", &isDim);
663:     if (isDim) continue;
664:     DMCreateLabel(dmNew, lname);
665:     DMGetLabel(dm, lname, &label);
666:     DMGetLabel(dmNew, lname, &newlabel);
667:     DMLabelGetDefaultValue(label,&val);
668:     DMLabelSetDefaultValue(newlabel,val);
669:     DMLabelGetValueIS(label, &valueIS);
670:     ISGetLocalSize(valueIS, &numValues);
671:     ISGetIndices(valueIS, &values);
672:     for (val = 0; val < numValues; ++val) {
673:       IS              pointIS;
674:       const PetscInt *points;
675:       PetscInt        numPoints, p;

677:       DMLabelGetStratumIS(label, values[val], &pointIS);
678:       ISGetLocalSize(pointIS, &numPoints);
679:       ISGetIndices(pointIS, &points);
680:       for (p = 0; p < numPoints; ++p) {
681:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

683:         DMLabelSetValue(newlabel, newpoint, values[val]);
684:       }
685:       ISRestoreIndices(pointIS, &points);
686:       ISDestroy(&pointIS);
687:     }
688:     ISRestoreIndices(valueIS, &values);
689:     ISDestroy(&valueIS);
690:   }
691:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
692:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
693:   DMGetPointSF(dm, &sfPoint);
694:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
695:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
696:   DMCreateLabel(dmNew, "vtk");
697:   DMCreateLabel(dmNew, "ghost");
698:   DMGetLabel(dmNew, "vtk", &vtkLabel);
699:   DMGetLabel(dmNew, "ghost", &ghostLabel);
700:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
701:     for (; c < leafLocal[l] && c < cEnd; ++c) {
702:       DMLabelSetValue(vtkLabel, c, 1);
703:     }
704:     if (leafLocal[l] >= cEnd) break;
705:     if (leafRemote[l].rank == rank) {
706:       DMLabelSetValue(vtkLabel, c, 1);
707:     } else {
708:       DMLabelSetValue(ghostLabel, c, 2);
709:     }
710:   }
711:   for (; c < cEnd; ++c) {
712:     DMLabelSetValue(vtkLabel, c, 1);
713:   }
714:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
715:   for (f = fStart; f < fEnd; ++f) {
716:     PetscInt numCells;

718:     DMPlexGetSupportSize(dmNew, f, &numCells);
719:     if (numCells < 2) {
720:       DMLabelSetValue(ghostLabel, f, 1);
721:     } else {
722:       const PetscInt *cells = NULL;
723:       PetscInt        vA, vB;

725:       DMPlexGetSupport(dmNew, f, &cells);
726:       DMLabelGetValue(vtkLabel, cells[0], &vA);
727:       DMLabelGetValue(vtkLabel, cells[1], &vB);
728:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
729:     }
730:   }
731:   return(0);
732: }

734: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
735: {
736:   DM             refTree;
737:   PetscSection   pSec;
738:   PetscInt       *parents, *childIDs;

742:   DMPlexGetReferenceTree(dm,&refTree);
743:   DMPlexSetReferenceTree(dmNew,refTree);
744:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
745:   if (pSec) {
746:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
747:     PetscInt *childIDsShifted;
748:     PetscSection pSecShifted;

750:     PetscSectionGetChart(pSec,&pStart,&pEnd);
751:     DMPlexGetDepth(dm,&depth);
752:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
753:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
754:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
755:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
756:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
757:     for (p = pStartShifted; p < pEndShifted; p++) {
758:       /* start off assuming no children */
759:       PetscSectionSetDof(pSecShifted,p,0);
760:     }
761:     for (p = pStart; p < pEnd; p++) {
762:       PetscInt dof;
763:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

765:       PetscSectionGetDof(pSec,p,&dof);
766:       PetscSectionSetDof(pSecShifted,pNew,dof);
767:     }
768:     PetscSectionSetUp(pSecShifted);
769:     for (p = pStart; p < pEnd; p++) {
770:       PetscInt dof;
771:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

773:       PetscSectionGetDof(pSec,p,&dof);
774:       if (dof) {
775:         PetscInt off, offNew;

777:         PetscSectionGetOffset(pSec,p,&off);
778:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
779:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
780:         childIDsShifted[offNew] = childIDs[off];
781:       }
782:     }
783:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
784:     PetscFree2(parentsShifted,childIDsShifted);
785:     PetscSectionDestroy(&pSecShifted);
786:   }
787:   return(0);
788: }

790: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
791: {
792:   PetscSF               sf;
793:   IS                    valueIS;
794:   const PetscInt       *values, *leaves;
795:   PetscInt             *depthShift;
796:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
797:   PetscBool             isper;
798:   const PetscReal      *maxCell, *L;
799:   const DMBoundaryType *bd;
800:   PetscErrorCode        ierr;

803:   DMGetPointSF(dm, &sf);
804:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
805:   nleaves = PetscMax(0, nleaves);
806:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
807:   /* Count ghost cells */
808:   DMLabelGetValueIS(label, &valueIS);
809:   ISGetLocalSize(valueIS, &numFS);
810:   ISGetIndices(valueIS, &values);
811:   Ng   = 0;
812:   for (fs = 0; fs < numFS; ++fs) {
813:     IS              faceIS;
814:     const PetscInt *faces;
815:     PetscInt        numFaces, f, numBdFaces = 0;

817:     DMLabelGetStratumIS(label, values[fs], &faceIS);
818:     ISGetLocalSize(faceIS, &numFaces);
819:     ISGetIndices(faceIS, &faces);
820:     for (f = 0; f < numFaces; ++f) {
821:       PetscInt numChildren;

823:       PetscFindInt(faces[f], nleaves, leaves, &loc);
824:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
825:       /* non-local and ancestors points don't get to register ghosts */
826:       if (loc >= 0 || numChildren) continue;
827:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
828:     }
829:     Ng += numBdFaces;
830:     ISRestoreIndices(faceIS, &faces);
831:     ISDestroy(&faceIS);
832:   }
833:   DMPlexGetDepth(dm, &depth);
834:   PetscMalloc1(2*(depth+1), &depthShift);
835:   for (d = 0; d <= depth; d++) {
836:     PetscInt dEnd;

838:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
839:     depthShift[2*d]   = dEnd;
840:     depthShift[2*d+1] = 0;
841:   }
842:   if (depth >= 0) depthShift[2*depth+1] = Ng;
843:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
844:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
845:   /* Step 3: Set cone/support sizes for new points */
846:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
847:   for (c = cEnd; c < cEnd + Ng; ++c) {
848:     DMPlexSetConeSize(gdm, c, 1);
849:   }
850:   for (fs = 0; fs < numFS; ++fs) {
851:     IS              faceIS;
852:     const PetscInt *faces;
853:     PetscInt        numFaces, f;

855:     DMLabelGetStratumIS(label, values[fs], &faceIS);
856:     ISGetLocalSize(faceIS, &numFaces);
857:     ISGetIndices(faceIS, &faces);
858:     for (f = 0; f < numFaces; ++f) {
859:       PetscInt size, numChildren;

861:       PetscFindInt(faces[f], nleaves, leaves, &loc);
862:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
863:       if (loc >= 0 || numChildren) continue;
864:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
865:       DMPlexGetSupportSize(dm, faces[f], &size);
866:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
867:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
868:     }
869:     ISRestoreIndices(faceIS, &faces);
870:     ISDestroy(&faceIS);
871:   }
872:   /* Step 4: Setup ghosted DM */
873:   DMSetUp(gdm);
874:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
875:   /* Step 6: Set cones and supports for new points */
876:   ghostCell = cEnd;
877:   for (fs = 0; fs < numFS; ++fs) {
878:     IS              faceIS;
879:     const PetscInt *faces;
880:     PetscInt        numFaces, f;

882:     DMLabelGetStratumIS(label, values[fs], &faceIS);
883:     ISGetLocalSize(faceIS, &numFaces);
884:     ISGetIndices(faceIS, &faces);
885:     for (f = 0; f < numFaces; ++f) {
886:       PetscInt newFace = faces[f] + Ng, numChildren;

888:       PetscFindInt(faces[f], nleaves, leaves, &loc);
889:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
890:       if (loc >= 0 || numChildren) continue;
891:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
892:       DMPlexSetCone(gdm, ghostCell, &newFace);
893:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
894:       ++ghostCell;
895:     }
896:     ISRestoreIndices(faceIS, &faces);
897:     ISDestroy(&faceIS);
898:   }
899:   ISRestoreIndices(valueIS, &values);
900:   ISDestroy(&valueIS);
901:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
902:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
903:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
904:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
905:   PetscFree(depthShift);
906:   for (c = cEnd; c < cEnd + Ng; ++c) {
907:     DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST);
908:   }
909:   /* Step 7: Periodicity */
910:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
911:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
912:   if (numGhostCells) *numGhostCells = Ng;
913:   return(0);
914: }

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

919:   Collective on dm

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

925:   Output Parameters:
926: + numGhostCells - The number of ghost cells added to the DM
927: - dmGhosted - The new DM

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

931:   Level: developer

933: .seealso: DMCreate()
934: @*/
935: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
936: {
937:   DM             gdm;
938:   DMLabel        label;
939:   const char    *name = labelName ? labelName : "Face Sets";
940:   PetscInt       dim, Ng = 0;
941:   PetscBool      useCone, useClosure;

948:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
949:   DMSetType(gdm, DMPLEX);
950:   DMGetDimension(dm, &dim);
951:   DMSetDimension(gdm, dim);
952:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
953:   DMSetBasicAdjacency(gdm, useCone,  useClosure);
954:   DMGetLabel(dm, name, &label);
955:   if (!label) {
956:     /* Get label for boundary faces */
957:     DMCreateLabel(dm, name);
958:     DMGetLabel(dm, name, &label);
959:     DMPlexMarkBoundaryFaces(dm, 1, label);
960:   }
961:   DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
962:   DMCopyBoundary(dm, gdm);
963:   DMCopyDisc(dm, gdm);
964:   gdm->setfromoptionscalled = dm->setfromoptionscalled;
965:   if (numGhostCells) *numGhostCells = Ng;
966:   *dmGhosted = gdm;
967:   return(0);
968: }

970: /*
971:   We are adding three kinds of points here:
972:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
973:     Non-replicated: Points which exist on the fault, but are not replicated
974:     Hybrid:         Entirely new points, such as cohesive cells

976:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
977:   each depth so that the new split/hybrid points can be inserted as a block.
978: */
979: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
980: {
981:   MPI_Comm         comm;
982:   IS               valueIS;
983:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
984:   const PetscInt  *values;          /* List of depths for which we have replicated points */
985:   IS              *splitIS;
986:   IS              *unsplitIS;
987:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
988:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
989:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
990:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
991:   const PetscInt **splitPoints;        /* Replicated points for each depth */
992:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
993:   PetscSection     coordSection;
994:   Vec              coordinates;
995:   PetscScalar     *coords;
996:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
997:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
998:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
999:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
1000:   PetscInt        *coneNew, *coneONew, *supportNew;
1001:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
1002:   PetscErrorCode   ierr;

1005:   PetscObjectGetComm((PetscObject)dm,&comm);
1006:   DMGetDimension(dm, &dim);
1007:   DMPlexGetDepth(dm, &depth);
1008:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1009:   /* We do not want this label automatically computed, instead we compute it here */
1010:   DMCreateLabel(sdm, "celltype");
1011:   /* Count split points and add cohesive cells */
1012:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
1013:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
1014:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
1015:   for (d = 0; d <= depth; ++d) {
1016:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
1017:     DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
1018:     depthEnd[d]           = pMaxNew[d];
1019:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1020:     numSplitPoints[d]     = 0;
1021:     numUnsplitPoints[d]   = 0;
1022:     numHybridPoints[d]    = 0;
1023:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1024:     splitPoints[d]        = NULL;
1025:     unsplitPoints[d]      = NULL;
1026:     splitIS[d]            = NULL;
1027:     unsplitIS[d]          = NULL;
1028:     /* we are shifting the existing hybrid points with the stratum behind them, so
1029:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
1030:     depthShift[2*d]       = depthMax[d];
1031:     depthShift[2*d+1]     = 0;
1032:   }
1033:   if (label) {
1034:     DMLabelGetValueIS(label, &valueIS);
1035:     ISGetLocalSize(valueIS, &numSP);
1036:     ISGetIndices(valueIS, &values);
1037:   }
1038:   for (sp = 0; sp < numSP; ++sp) {
1039:     const PetscInt dep = values[sp];

1041:     if ((dep < 0) || (dep > depth)) continue;
1042:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1043:     if (splitIS[dep]) {
1044:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1045:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
1046:     }
1047:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
1048:     if (unsplitIS[dep]) {
1049:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1050:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1051:     }
1052:   }
1053:   /* Calculate number of hybrid points */
1054:   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   */
1055:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
1056:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
1057:   /* the end of the points in this stratum that come before the new points:
1058:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1059:    * added points */
1060:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1061:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1062:   /* Step 3: Set cone/support sizes for new points */
1063:   for (dep = 0; dep <= depth; ++dep) {
1064:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1065:       const PetscInt  oldp   = splitPoints[dep][p];
1066:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1067:       const PetscInt  splitp = p    + pMaxNew[dep];
1068:       const PetscInt *support;
1069:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

1071:       DMPlexGetConeSize(dm, oldp, &coneSize);
1072:       DMPlexSetConeSize(sdm, splitp, coneSize);
1073:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1074:       DMPlexSetSupportSize(sdm, splitp, supportSize);
1075:       if (dep == depth-1) {
1076:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1078:         /* Add cohesive cells, they are prisms */
1079:         DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
1080:         switch (coneSize) {
1081:           case 2: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);break;
1082:           case 3: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);break;
1083:           case 4: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);break;
1084:         }
1085:       } else if (dep == 0) {
1086:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1088:         DMPlexGetSupport(dm, oldp, &support);
1089:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1090:           PetscInt val;

1092:           DMLabelGetValue(label, support[e], &val);
1093:           if (val == 1) ++qf;
1094:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
1095:           if ((val == 1) || (val == -(shift + 1))) ++qp;
1096:         }
1097:         /* Split old vertex: Edges into original vertex and new cohesive edge */
1098:         DMPlexSetSupportSize(sdm, newp, qn+1);
1099:         /* Split new vertex: Edges into split vertex and new cohesive edge */
1100:         DMPlexSetSupportSize(sdm, splitp, qp+1);
1101:         /* Add hybrid edge */
1102:         DMPlexSetConeSize(sdm, hybedge, 2);
1103:         DMPlexSetSupportSize(sdm, hybedge, qf);
1104:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1105:       } else if (dep == dim-2) {
1106:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1108:         DMPlexGetSupport(dm, oldp, &support);
1109:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1110:           PetscInt val;

1112:           DMLabelGetValue(label, support[e], &val);
1113:           if (val == dim-1) ++qf;
1114:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
1115:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1116:         }
1117:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1118:         DMPlexSetSupportSize(sdm, newp, qn+1);
1119:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1120:         DMPlexSetSupportSize(sdm, splitp, qp+1);
1121:         /* Add hybrid face */
1122:         DMPlexSetConeSize(sdm, hybface, 4);
1123:         DMPlexSetSupportSize(sdm, hybface, qf);
1124:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1125:       }
1126:     }
1127:   }
1128:   for (dep = 0; dep <= depth; ++dep) {
1129:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1130:       const PetscInt  oldp   = unsplitPoints[dep][p];
1131:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1132:       const PetscInt *support;
1133:       PetscInt        coneSize, supportSize, qf, e, s;

1135:       DMPlexGetConeSize(dm, oldp, &coneSize);
1136:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1137:       DMPlexGetSupport(dm, oldp, &support);
1138:       if (dep == 0) {
1139:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1141:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1142:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1143:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1144:           if (e >= 0) ++qf;
1145:         }
1146:         DMPlexSetSupportSize(sdm, newp, qf+2);
1147:         /* Add hybrid edge */
1148:         DMPlexSetConeSize(sdm, hybedge, 2);
1149:         for (e = 0, qf = 0; e < supportSize; ++e) {
1150:           PetscInt val;

1152:           DMLabelGetValue(label, support[e], &val);
1153:           /* Split and unsplit edges produce hybrid faces */
1154:           if (val == 1) ++qf;
1155:           if (val == (shift2 + 1)) ++qf;
1156:         }
1157:         DMPlexSetSupportSize(sdm, hybedge, qf);
1158:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1159:       } else if (dep == dim-2) {
1160:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1161:         PetscInt       val;

1163:         for (e = 0, qf = 0; e < supportSize; ++e) {
1164:           DMLabelGetValue(label, support[e], &val);
1165:           if (val == dim-1) qf += 2;
1166:           else              ++qf;
1167:         }
1168:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1169:         DMPlexSetSupportSize(sdm, newp, qf+2);
1170:         /* Add hybrid face */
1171:         for (e = 0, qf = 0; e < supportSize; ++e) {
1172:           DMLabelGetValue(label, support[e], &val);
1173:           if (val == dim-1) ++qf;
1174:         }
1175:         DMPlexSetConeSize(sdm, hybface, 4);
1176:         DMPlexSetSupportSize(sdm, hybface, qf);
1177:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1178:       }
1179:     }
1180:   }
1181:   /* Step 4: Setup split DM */
1182:   DMSetUp(sdm);
1183:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1184:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1185:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1186:   /* Step 6: Set cones and supports for new points */
1187:   for (dep = 0; dep <= depth; ++dep) {
1188:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1189:       const PetscInt  oldp   = splitPoints[dep][p];
1190:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1191:       const PetscInt  splitp = p    + pMaxNew[dep];
1192:       const PetscInt *cone, *support, *ornt;
1193:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1195:       DMPlexGetConeSize(dm, oldp, &coneSize);
1196:       DMPlexGetCone(dm, oldp, &cone);
1197:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1198:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1199:       DMPlexGetSupport(dm, oldp, &support);
1200:       if (dep == depth-1) {
1201:         PetscBool       hasUnsplit = PETSC_FALSE;
1202:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1203:         const PetscInt *supportF;

1205:         /* Split face:       copy in old face to new face to start */
1206:         DMPlexGetSupport(sdm, newp,  &supportF);
1207:         DMPlexSetSupport(sdm, splitp, supportF);
1208:         /* Split old face:   old vertices/edges in cone so no change */
1209:         /* Split new face:   new vertices/edges in cone */
1210:         for (q = 0; q < coneSize; ++q) {
1211:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1212:           if (v < 0) {
1213:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1214:             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);
1215:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1216:             hasUnsplit   = PETSC_TRUE;
1217:           } else {
1218:             coneNew[2+q] = v + pMaxNew[dep-1];
1219:             if (dep > 1) {
1220:               const PetscInt *econe;
1221:               PetscInt        econeSize, r, vs, vu;

1223:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1224:               DMPlexGetCone(dm, cone[q], &econe);
1225:               for (r = 0; r < econeSize; ++r) {
1226:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
1227:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1228:                 if (vs >= 0) continue;
1229:                 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);
1230:                 hasUnsplit   = PETSC_TRUE;
1231:               }
1232:             }
1233:           }
1234:         }
1235:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1236:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1237:         /* Face support */
1238:         for (s = 0; s < supportSize; ++s) {
1239:           PetscInt val;

1241:           DMLabelGetValue(label, support[s], &val);
1242:           if (val < 0) {
1243:             /* Split old face:   Replace negative side cell with cohesive cell */
1244:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1245:           } else {
1246:             /* Split new face:   Replace positive side cell with cohesive cell */
1247:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1248:             /* Get orientation for cohesive face */
1249:             {
1250:               const PetscInt *ncone, *nconeO;
1251:               PetscInt        nconeSize, nc;

1253:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1254:               DMPlexGetCone(dm, support[s], &ncone);
1255:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1256:               for (nc = 0; nc < nconeSize; ++nc) {
1257:                 if (ncone[nc] == oldp) {
1258:                   coneONew[0] = nconeO[nc];
1259:                   break;
1260:                 }
1261:               }
1262:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1263:             }
1264:           }
1265:         }
1266:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1267:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1268:         coneNew[1]  = splitp;
1269:         coneONew[1] = coneONew[0];
1270:         for (q = 0; q < coneSize; ++q) {
1271:           /* Hybrid faces must follow order from oriented end face */
1272:           const PetscInt o  = coneONew[0];
1273:           const PetscInt qo = o < 0 ? (-(o+1)+coneSize-q)%coneSize : (q+o)%coneSize;

1275:           PetscFindInt(cone[qo], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1276:           if (v < 0) {
1277:             PetscFindInt(cone[qo], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1278:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1279:           } else {
1280:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1281:           }
1282:           coneONew[2+q] = ((o < 0) + (ornt[qo] < 0))%2 ? -1 : 0;
1283:         }
1284:         DMPlexSetCone(sdm, hybcell, coneNew);
1285:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1286:         /* Label the hybrid cells on the boundary of the split */
1287:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1288:       } else if (dep == 0) {
1289:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1295:           DMLabelGetValue(label, support[e], &val);
1296:           if ((val == 1) || (val == (shift + 1))) {
1297:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1298:           }
1299:         }
1300:         supportNew[qn] = hybedge;
1301:         DMPlexSetSupport(sdm, newp, supportNew);
1302:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1303:         for (e = 0, qp = 0; e < supportSize; ++e) {
1304:           PetscInt val, edge;

1306:           DMLabelGetValue(label, support[e], &val);
1307:           if (val == 1) {
1308:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1309:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1310:             supportNew[qp++] = edge + pMaxNew[dep+1];
1311:           } else if (val == -(shift + 1)) {
1312:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1313:           }
1314:         }
1315:         supportNew[qp] = hybedge;
1316:         DMPlexSetSupport(sdm, splitp, supportNew);
1317:         /* Hybrid edge:    Old and new split vertex */
1318:         coneNew[0] = newp;
1319:         coneNew[1] = splitp;
1320:         DMPlexSetCone(sdm, hybedge, coneNew);
1321:         for (e = 0, qf = 0; e < supportSize; ++e) {
1322:           PetscInt val, edge;

1324:           DMLabelGetValue(label, support[e], &val);
1325:           if (val == 1) {
1326:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1327:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1328:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1329:           }
1330:         }
1331:         DMPlexSetSupport(sdm, hybedge, supportNew);
1332:       } else if (dep == dim-2) {
1333:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1335:         /* Split old edge:   old vertices in cone so no change */
1336:         /* Split new edge:   new vertices in cone */
1337:         for (q = 0; q < coneSize; ++q) {
1338:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1339:           if (v < 0) {
1340:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1341:             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);
1342:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1343:           } else {
1344:             coneNew[q] = v + pMaxNew[dep-1];
1345:           }
1346:         }
1347:         DMPlexSetCone(sdm, splitp, coneNew);
1348:         /* Split old edge: Faces in positive side cells and old split faces */
1349:         for (e = 0, q = 0; e < supportSize; ++e) {
1350:           PetscInt val;

1352:           DMLabelGetValue(label, support[e], &val);
1353:           if (val == dim-1) {
1354:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1355:           } else if (val == (shift + dim-1)) {
1356:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1357:           }
1358:         }
1359:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1360:         DMPlexSetSupport(sdm, newp, supportNew);
1361:         /* Split new edge: Faces in negative side cells and new split faces */
1362:         for (e = 0, q = 0; e < supportSize; ++e) {
1363:           PetscInt val, face;

1365:           DMLabelGetValue(label, support[e], &val);
1366:           if (val == dim-1) {
1367:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1368:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1369:             supportNew[q++] = face + pMaxNew[dep+1];
1370:           } else if (val == -(shift + dim-1)) {
1371:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1372:           }
1373:         }
1374:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1375:         DMPlexSetSupport(sdm, splitp, supportNew);
1376:         /* Hybrid face */
1377:         coneNew[0] = newp;
1378:         coneNew[1] = splitp;
1379:         for (v = 0; v < coneSize; ++v) {
1380:           PetscInt vertex;
1381:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1382:           if (vertex < 0) {
1383:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1384:             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);
1385:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1386:           } else {
1387:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1388:           }
1389:         }
1390:         DMPlexSetCone(sdm, hybface, coneNew);
1391:         for (e = 0, qf = 0; e < supportSize; ++e) {
1392:           PetscInt val, face;

1394:           DMLabelGetValue(label, support[e], &val);
1395:           if (val == dim-1) {
1396:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1397:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1398:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1399:           }
1400:         }
1401:         DMPlexSetSupport(sdm, hybface, supportNew);
1402:       }
1403:     }
1404:   }
1405:   for (dep = 0; dep <= depth; ++dep) {
1406:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1407:       const PetscInt  oldp   = unsplitPoints[dep][p];
1408:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1409:       const PetscInt *cone, *support, *ornt;
1410:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1412:       DMPlexGetConeSize(dm, oldp, &coneSize);
1413:       DMPlexGetCone(dm, oldp, &cone);
1414:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1415:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1416:       DMPlexGetSupport(dm, oldp, &support);
1417:       if (dep == 0) {
1418:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1420:         /* Unsplit vertex */
1421:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1422:         for (s = 0, q = 0; s < supportSize; ++s) {
1423:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1424:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1425:           if (e >= 0) {
1426:             supportNew[q++] = e + pMaxNew[dep+1];
1427:           }
1428:         }
1429:         supportNew[q++] = hybedge;
1430:         supportNew[q++] = hybedge;
1431:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1432:         DMPlexSetSupport(sdm, newp, supportNew);
1433:         /* Hybrid edge */
1434:         coneNew[0] = newp;
1435:         coneNew[1] = newp;
1436:         DMPlexSetCone(sdm, hybedge, coneNew);
1437:         for (e = 0, qf = 0; e < supportSize; ++e) {
1438:           PetscInt val, edge;

1440:           DMLabelGetValue(label, support[e], &val);
1441:           if (val == 1) {
1442:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1443:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1444:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1445:           } else if  (val ==  (shift2 + 1)) {
1446:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1447:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1448:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1449:           }
1450:         }
1451:         DMPlexSetSupport(sdm, hybedge, supportNew);
1452:       } else if (dep == dim-2) {
1453:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1459:           DMLabelGetValue(label, support[f], &val);
1460:           if (val == dim-1) {
1461:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1462:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1463:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1464:             supportNew[qf++] = face + pMaxNew[dep+1];
1465:           } else {
1466:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1467:           }
1468:         }
1469:         supportNew[qf++] = hybface;
1470:         supportNew[qf++] = hybface;
1471:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1472:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1473:         DMPlexSetSupport(sdm, newp, supportNew);
1474:         /* Add hybrid face */
1475:         coneNew[0] = newp;
1476:         coneNew[1] = newp;
1477:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1478:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1479:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1480:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1481:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1482:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1483:         DMPlexSetCone(sdm, hybface, coneNew);
1484:         for (f = 0, qf = 0; f < supportSize; ++f) {
1485:           PetscInt val, face;

1487:           DMLabelGetValue(label, support[f], &val);
1488:           if (val == dim-1) {
1489:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1490:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1491:           }
1492:         }
1493:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1494:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1495:         DMPlexSetSupport(sdm, hybface, supportNew);
1496:       }
1497:     }
1498:   }
1499:   /* Step 6b: Replace split points in negative side cones */
1500:   for (sp = 0; sp < numSP; ++sp) {
1501:     PetscInt        dep = values[sp];
1502:     IS              pIS;
1503:     PetscInt        numPoints;
1504:     const PetscInt *points;

1506:     if (dep >= 0) continue;
1507:     DMLabelGetStratumIS(label, dep, &pIS);
1508:     if (!pIS) continue;
1509:     dep  = -dep - shift;
1510:     ISGetLocalSize(pIS, &numPoints);
1511:     ISGetIndices(pIS, &points);
1512:     for (p = 0; p < numPoints; ++p) {
1513:       const PetscInt  oldp = points[p];
1514:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1515:       const PetscInt *cone;
1516:       PetscInt        coneSize, c;
1517:       /* PetscBool       replaced = PETSC_FALSE; */

1519:       /* Negative edge: replace split vertex */
1520:       /* Negative cell: replace split face */
1521:       DMPlexGetConeSize(sdm, newp, &coneSize);
1522:       DMPlexGetCone(sdm, newp, &cone);
1523:       for (c = 0; c < coneSize; ++c) {
1524:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1525:         PetscInt       csplitp, cp, val;

1527:         DMLabelGetValue(label, coldp, &val);
1528:         if (val == dep-1) {
1529:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1530:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1531:           csplitp  = pMaxNew[dep-1] + cp;
1532:           DMPlexInsertCone(sdm, newp, c, csplitp);
1533:           /* replaced = PETSC_TRUE; */
1534:         }
1535:       }
1536:       /* Cells with only a vertex or edge on the submesh have no replacement */
1537:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1538:     }
1539:     ISRestoreIndices(pIS, &points);
1540:     ISDestroy(&pIS);
1541:   }
1542:   /* Step 7: Coordinates */
1543:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1544:   DMGetCoordinateSection(sdm, &coordSection);
1545:   DMGetCoordinatesLocal(sdm, &coordinates);
1546:   VecGetArray(coordinates, &coords);
1547:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1548:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1549:     const PetscInt splitp = pMaxNew[0] + v;
1550:     PetscInt       dof, off, soff, d;

1552:     PetscSectionGetDof(coordSection, newp, &dof);
1553:     PetscSectionGetOffset(coordSection, newp, &off);
1554:     PetscSectionGetOffset(coordSection, splitp, &soff);
1555:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1556:   }
1557:   VecRestoreArray(coordinates, &coords);
1558:   /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1559:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1560:   /* Step 9: Labels */
1561:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1562:   DMGetNumLabels(sdm, &numLabels);
1563:   for (dep = 0; dep <= depth; ++dep) {
1564:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1565:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1566:       const PetscInt splitp = pMaxNew[dep] + p;
1567:       PetscInt       l;

1569:       if (splitLabel) {
1570:         const PetscInt val = 100 + dep;

1572:         DMLabelSetValue(splitLabel, newp,    val);
1573:         DMLabelSetValue(splitLabel, splitp, -val);
1574:       }
1575:       for (l = 0; l < numLabels; ++l) {
1576:         DMLabel     mlabel;
1577:         const char *lname;
1578:         PetscInt    val;
1579:         PetscBool   isDepth;

1581:         DMGetLabelName(sdm, l, &lname);
1582:         PetscStrcmp(lname, "depth", &isDepth);
1583:         if (isDepth) continue;
1584:         DMGetLabel(sdm, lname, &mlabel);
1585:         DMLabelGetValue(mlabel, newp, &val);
1586:         if (val >= 0) {
1587:           DMLabelSetValue(mlabel, splitp, val);
1588:         }
1589:       }
1590:     }
1591:   }
1592:   for (sp = 0; sp < numSP; ++sp) {
1593:     const PetscInt dep = values[sp];

1595:     if ((dep < 0) || (dep > depth)) continue;
1596:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1597:     ISDestroy(&splitIS[dep]);
1598:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1599:     ISDestroy(&unsplitIS[dep]);
1600:   }
1601:   if (label) {
1602:     ISRestoreIndices(valueIS, &values);
1603:     ISDestroy(&valueIS);
1604:   }
1605:   for (d = 0; d <= depth; ++d) {
1606:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1607:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1608:   }
1609:   PetscFree3(coneNew, coneONew, supportNew);
1610:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1611:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1612:   return(0);
1613: }

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

1618:   Collective on dm

1620:   Input Parameters:
1621: + dm - The original DM
1622: - label - The label specifying the boundary faces (this could be auto-generated)

1624:   Output Parameters:
1625: + splitLabel - The label containing the split points, or NULL if no output is desired
1626: - dmSplit - The new DM

1628:   Level: developer

1630: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1631: @*/
1632: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1633: {
1634:   DM             sdm;
1635:   PetscInt       dim;

1641:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1642:   DMSetType(sdm, DMPLEX);
1643:   DMGetDimension(dm, &dim);
1644:   DMSetDimension(sdm, dim);
1645:   switch (dim) {
1646:   case 2:
1647:   case 3:
1648:     DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1649:     break;
1650:   default:
1651:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1652:   }
1653:   *dmSplit = sdm;
1654:   return(0);
1655: }

1657: /* Returns the side of the surface for a given cell with a face on the surface */
1658: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1659: {
1660:   const PetscInt *cone, *ornt;
1661:   PetscInt        dim, coneSize, c;
1662:   PetscErrorCode  ierr;

1665:   *pos = PETSC_TRUE;
1666:   DMGetDimension(dm, &dim);
1667:   DMPlexGetConeSize(dm, cell, &coneSize);
1668:   DMPlexGetCone(dm, cell, &cone);
1669:   DMPlexGetConeOrientation(dm, cell, &ornt);
1670:   for (c = 0; c < coneSize; ++c) {
1671:     if (cone[c] == face) {
1672:       PetscInt o = ornt[c];

1674:       if (subdm) {
1675:         const PetscInt *subcone, *subornt;
1676:         PetscInt        subpoint, subface, subconeSize, sc;

1678:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1679:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1680:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1681:         DMPlexGetCone(subdm, subpoint, &subcone);
1682:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1683:         for (sc = 0; sc < subconeSize; ++sc) {
1684:           if (subcone[sc] == subface) {
1685:             o = subornt[0];
1686:             break;
1687:           }
1688:         }
1689:         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);
1690:       }
1691:       if (o >= 0) *pos = PETSC_TRUE;
1692:       else        *pos = PETSC_FALSE;
1693:       break;
1694:     }
1695:   }
1696:   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);
1697:   return(0);
1698: }

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

1704:   Input Parameters:
1705: + dm     - The DM
1706: . label  - A DMLabel marking the surface
1707: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1708: . flip   - Flag to flip the submesh normal and replace points on the other side
1709: - subdm  - The subDM associated with the label, or NULL

1711:   Output Parameter:
1712: . label - A DMLabel marking all surface points

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

1716:   Level: developer

1718: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1719: @*/
1720: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1721: {
1722:   DMLabel         depthLabel;
1723:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1724:   const PetscInt *points, *subpoints;
1725:   const PetscInt  rev   = flip ? -1 : 1;
1726:   PetscInt        shift = 100, shift2 = 200, dim, depth, dep, cStart, cEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1727:   PetscErrorCode  ierr;

1730:   DMPlexGetDepth(dm, &depth);
1731:   DMGetDimension(dm, &dim);
1732:   DMPlexGetDepthLabel(dm, &depthLabel);
1733:   if (subdm) {
1734:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1735:     if (subpointIS) {
1736:       ISGetLocalSize(subpointIS, &numSubpoints);
1737:       ISGetIndices(subpointIS, &subpoints);
1738:     }
1739:   }
1740:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1741:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1742:   if (!dimIS) {
1743:     ISDestroy(&subpointIS);
1744:     return(0);
1745:   }
1746:   ISGetLocalSize(dimIS, &numPoints);
1747:   ISGetIndices(dimIS, &points);
1748:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1749:     const PetscInt *support;
1750:     PetscInt        supportSize, s;

1752:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1753: #if 0
1754:     if (supportSize != 2) {
1755:       const PetscInt *lp;
1756:       PetscInt        Nlp, pind;

1758:       /* Check that for a cell with a single support face, that face is in the SF */
1759:       /*   THis check only works for the remote side. We would need root side information */
1760:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1761:       PetscFindInt(points[p], Nlp, lp, &pind);
1762:       if (pind < 0) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports, and the face is not shared with another process", points[p], supportSize);
1763:     }
1764: #endif
1765:     DMPlexGetSupport(dm, points[p], &support);
1766:     for (s = 0; s < supportSize; ++s) {
1767:       const PetscInt *cone;
1768:       PetscInt        coneSize, c;
1769:       PetscBool       pos;

1771:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1772:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1773:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1774:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1775:       /* Put faces touching the fault in the label */
1776:       DMPlexGetConeSize(dm, support[s], &coneSize);
1777:       DMPlexGetCone(dm, support[s], &cone);
1778:       for (c = 0; c < coneSize; ++c) {
1779:         const PetscInt point = cone[c];

1781:         DMLabelGetValue(label, point, &val);
1782:         if (val == -1) {
1783:           PetscInt *closure = NULL;
1784:           PetscInt  closureSize, cl;

1786:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1787:           for (cl = 0; cl < closureSize*2; cl += 2) {
1788:             const PetscInt clp  = closure[cl];
1789:             PetscInt       bval = -1;

1791:             DMLabelGetValue(label, clp, &val);
1792:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1793:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1794:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1795:               break;
1796:             }
1797:           }
1798:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1799:         }
1800:       }
1801:     }
1802:   }
1803:   ISRestoreIndices(dimIS, &points);
1804:   ISDestroy(&dimIS);
1805:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1806:   ISDestroy(&subpointIS);
1807:   /* Mark boundary points as unsplit */
1808:   if (blabel) {
1809:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1810:     ISGetLocalSize(dimIS, &numPoints);
1811:     ISGetIndices(dimIS, &points);
1812:     for (p = 0; p < numPoints; ++p) {
1813:       const PetscInt point = points[p];
1814:       PetscInt       val, bval;

1816:       DMLabelGetValue(blabel, point, &bval);
1817:       if (bval >= 0) {
1818:         DMLabelGetValue(label, point, &val);
1819:         if ((val < 0) || (val > dim)) {
1820:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1821:           DMLabelClearValue(blabel, point, bval);
1822:         }
1823:       }
1824:     }
1825:     for (p = 0; p < numPoints; ++p) {
1826:       const PetscInt point = points[p];
1827:       PetscInt       val, bval;

1829:       DMLabelGetValue(blabel, point, &bval);
1830:       if (bval >= 0) {
1831:         const PetscInt *cone,    *support;
1832:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1834:         /* Mark as unsplit */
1835:         DMLabelGetValue(label, point, &val);
1836:         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);
1837:         DMLabelClearValue(label, point, val);
1838:         DMLabelSetValue(label, point, shift2+val);
1839:         /* Check for cross-edge
1840:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1841:         if (val != 0) continue;
1842:         DMPlexGetSupport(dm, point, &support);
1843:         DMPlexGetSupportSize(dm, point, &supportSize);
1844:         for (s = 0; s < supportSize; ++s) {
1845:           DMPlexGetCone(dm, support[s], &cone);
1846:           DMPlexGetConeSize(dm, support[s], &coneSize);
1847:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1848:           DMLabelGetValue(blabel, cone[0], &valA);
1849:           DMLabelGetValue(blabel, cone[1], &valB);
1850:           DMLabelGetValue(blabel, support[s], &valE);
1851:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1852:         }
1853:       }
1854:     }
1855:     ISRestoreIndices(dimIS, &points);
1856:     ISDestroy(&dimIS);
1857:   }
1858:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1859:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1860:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1861:   DMLabelGetStratumIS(label, 0, &dimIS);
1862:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1863:   if (dimIS && crossEdgeIS) {
1864:     IS vertIS = dimIS;

1866:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1867:     ISDestroy(&crossEdgeIS);
1868:     ISDestroy(&vertIS);
1869:   }
1870:   if (!dimIS) {
1871:     return(0);
1872:   }
1873:   ISGetLocalSize(dimIS, &numPoints);
1874:   ISGetIndices(dimIS, &points);
1875:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1876:     PetscInt *star = NULL;
1877:     PetscInt  starSize, s;
1878:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1880:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1881:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1882:     while (again) {
1883:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1884:       again = 0;
1885:       for (s = 0; s < starSize*2; s += 2) {
1886:         const PetscInt  point = star[s];
1887:         const PetscInt *cone;
1888:         PetscInt        coneSize, c;

1890:         if ((point < cStart) || (point >= cEnd)) continue;
1891:         DMLabelGetValue(label, point, &val);
1892:         if (val != -1) continue;
1893:         again = again == 1 ? 1 : 2;
1894:         DMPlexGetConeSize(dm, point, &coneSize);
1895:         DMPlexGetCone(dm, point, &cone);
1896:         for (c = 0; c < coneSize; ++c) {
1897:           DMLabelGetValue(label, cone[c], &val);
1898:           if (val != -1) {
1899:             const PetscInt *ccone;
1900:             PetscInt        cconeSize, cc, side;

1902:             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);
1903:             if (val > 0) side =  1;
1904:             else         side = -1;
1905:             DMLabelSetValue(label, point, side*(shift+dim));
1906:             /* Mark cell faces which touch the fault */
1907:             DMPlexGetConeSize(dm, point, &cconeSize);
1908:             DMPlexGetCone(dm, point, &ccone);
1909:             for (cc = 0; cc < cconeSize; ++cc) {
1910:               PetscInt *closure = NULL;
1911:               PetscInt  closureSize, cl;

1913:               DMLabelGetValue(label, ccone[cc], &val);
1914:               if (val != -1) continue;
1915:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1916:               for (cl = 0; cl < closureSize*2; cl += 2) {
1917:                 const PetscInt clp = closure[cl];

1919:                 DMLabelGetValue(label, clp, &val);
1920:                 if (val == -1) continue;
1921:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1922:                 break;
1923:               }
1924:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1925:             }
1926:             again = 1;
1927:             break;
1928:           }
1929:         }
1930:       }
1931:     }
1932:     /* Classify the rest by cell membership */
1933:     for (s = 0; s < starSize*2; s += 2) {
1934:       const PetscInt point = star[s];

1936:       DMLabelGetValue(label, point, &val);
1937:       if (val == -1) {
1938:         PetscInt      *sstar = NULL;
1939:         PetscInt       sstarSize, ss;
1940:         PetscBool      marked = PETSC_FALSE, isHybrid;

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

1946:           if ((spoint < cStart) || (spoint >= cEnd)) continue;
1947:           DMLabelGetValue(label, spoint, &val);
1948:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1949:           DMLabelGetValue(depthLabel, point, &dep);
1950:           if (val > 0) {
1951:             DMLabelSetValue(label, point,   shift+dep);
1952:           } else {
1953:             DMLabelSetValue(label, point, -(shift+dep));
1954:           }
1955:           marked = PETSC_TRUE;
1956:           break;
1957:         }
1958:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1959:         DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);
1960:         if (!isHybrid && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1961:       }
1962:     }
1963:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1964:   }
1965:   ISRestoreIndices(dimIS, &points);
1966:   ISDestroy(&dimIS);
1967:   /* If any faces touching the fault divide cells on either side, split them */
1968:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1969:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1970:   ISExpand(facePosIS, faceNegIS, &dimIS);
1971:   ISDestroy(&facePosIS);
1972:   ISDestroy(&faceNegIS);
1973:   ISGetLocalSize(dimIS, &numPoints);
1974:   ISGetIndices(dimIS, &points);
1975:   for (p = 0; p < numPoints; ++p) {
1976:     const PetscInt  point = points[p];
1977:     const PetscInt *support;
1978:     PetscInt        supportSize, valA, valB;

1980:     DMPlexGetSupportSize(dm, point, &supportSize);
1981:     if (supportSize != 2) continue;
1982:     DMPlexGetSupport(dm, point, &support);
1983:     DMLabelGetValue(label, support[0], &valA);
1984:     DMLabelGetValue(label, support[1], &valB);
1985:     if ((valA == -1) || (valB == -1)) continue;
1986:     if (valA*valB > 0) continue;
1987:     /* Split the face */
1988:     DMLabelGetValue(label, point, &valA);
1989:     DMLabelClearValue(label, point, valA);
1990:     DMLabelSetValue(label, point, dim-1);
1991:     /* Label its closure:
1992:       unmarked: label as unsplit
1993:       incident: relabel as split
1994:       split:    do nothing
1995:     */
1996:     {
1997:       PetscInt *closure = NULL;
1998:       PetscInt  closureSize, cl;

2000:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2001:       for (cl = 0; cl < closureSize*2; cl += 2) {
2002:         DMLabelGetValue(label, closure[cl], &valA);
2003:         if (valA == -1) { /* Mark as unsplit */
2004:           DMLabelGetValue(depthLabel, closure[cl], &dep);
2005:           DMLabelSetValue(label, closure[cl], shift2+dep);
2006:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
2007:           DMLabelGetValue(depthLabel, closure[cl], &dep);
2008:           DMLabelClearValue(label, closure[cl], valA);
2009:           DMLabelSetValue(label, closure[cl], dep);
2010:         }
2011:       }
2012:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2013:     }
2014:   }
2015:   ISRestoreIndices(dimIS, &points);
2016:   ISDestroy(&dimIS);
2017:   return(0);
2018: }

2020: /* Check that no cell have all vertices on the fault */
2021: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2022: {
2023:   IS              subpointIS;
2024:   const PetscInt *dmpoints;
2025:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;
2026:   PetscErrorCode  ierr;

2029:   if (!label) return(0);
2030:   DMLabelGetDefaultValue(label, &defaultValue);
2031:   DMPlexCreateSubpointIS(subdm, &subpointIS);
2032:   if (!subpointIS) return(0);
2033:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2034:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2035:   ISGetIndices(subpointIS, &dmpoints);
2036:   for (c = cStart; c < cEnd; ++c) {
2037:     PetscBool invalidCell = PETSC_TRUE;
2038:     PetscInt *closure     = NULL;
2039:     PetscInt  closureSize, cl;

2041:     DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2042:     for (cl = 0; cl < closureSize*2; cl += 2) {
2043:       PetscInt value = 0;

2045:       if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2046:       DMLabelGetValue(label, closure[cl], &value);
2047:       if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
2048:     }
2049:     DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2050:     if (invalidCell) {
2051:       ISRestoreIndices(subpointIS, &dmpoints);
2052:       ISDestroy(&subpointIS);
2053:       DMDestroy(&subdm);
2054:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
2055:     }
2056:   }
2057:   ISRestoreIndices(subpointIS, &dmpoints);
2058:   ISDestroy(&subpointIS);
2059:   return(0);
2060: }


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

2066:   Collective on dm

2068:   Input Parameters:
2069: + dm - The original DM
2070: . label - The label specifying the interface vertices
2071: - bdlabel - The optional label specifying the interface boundary vertices

2073:   Output Parameters:
2074: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2075: . splitLabel - The label containing the split points, or NULL if no output is desired
2076: . dmInterface - The new interface DM, or NULL
2077: - dmHybrid - The new DM with cohesive cells

2079:   Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
2080:   directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2081:   7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2082:   one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
2083:   surface also hits vertex 3, it would be 9 (-101) in hybridLabel.

2085:   The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2086:   mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2087:   splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.

2089:   The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
2090:   DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.

2092:   Level: developer

2094: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2095: @*/
2096: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2097: {
2098:   DM             idm;
2099:   DMLabel        subpointMap, hlabel, slabel = NULL;
2100:   PetscInt       dim;

2110:   DMGetDimension(dm, &dim);
2111:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2112:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
2113:   DMPlexOrient(idm);
2114:   DMPlexGetSubpointMap(idm, &subpointMap);
2115:   DMLabelDuplicate(subpointMap, &hlabel);
2116:   DMLabelClearStratum(hlabel, dim);
2117:   if (splitLabel) {
2118:     const char *name;
2119:     char        sname[PETSC_MAX_PATH_LEN];

2121:     PetscObjectGetName((PetscObject) hlabel, &name);
2122:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2123:     PetscStrcat(sname, " split");
2124:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2125:   }
2126:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2127:   if (dmInterface) {*dmInterface = idm;}
2128:   else             {DMDestroy(&idm);}
2129:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2130:   if (hybridLabel) *hybridLabel = hlabel;
2131:   else             {DMLabelDestroy(&hlabel);}
2132:   if (splitLabel)  *splitLabel  = slabel;
2133:   return(0);
2134: }

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

2138:      For any marked cell, the marked vertices constitute a single face
2139: */
2140: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2141: {
2142:   IS               subvertexIS = NULL;
2143:   const PetscInt  *subvertices;
2144:   PetscInt        *pStart, *pEnd, pSize;
2145:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
2146:   PetscErrorCode   ierr;

2149:   *numFaces = 0;
2150:   *nFV      = 0;
2151:   DMPlexGetDepth(dm, &depth);
2152:   DMGetDimension(dm, &dim);
2153:   pSize = PetscMax(depth, dim) + 1;
2154:   PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2155:   for (d = 0; d <= depth; ++d) {
2156:     DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);
2157:   }
2158:   /* Loop over initial vertices and mark all faces in the collective star() */
2159:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2160:   if (subvertexIS) {
2161:     ISGetSize(subvertexIS, &numSubVerticesInitial);
2162:     ISGetIndices(subvertexIS, &subvertices);
2163:   }
2164:   for (v = 0; v < numSubVerticesInitial; ++v) {
2165:     const PetscInt vertex = subvertices[v];
2166:     PetscInt      *star   = NULL;
2167:     PetscInt       starSize, s, numCells = 0, c;

2169:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2170:     for (s = 0; s < starSize*2; s += 2) {
2171:       const PetscInt point = star[s];
2172:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2173:     }
2174:     for (c = 0; c < numCells; ++c) {
2175:       const PetscInt cell    = star[c];
2176:       PetscInt      *closure = NULL;
2177:       PetscInt       closureSize, cl;
2178:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

2180:       DMLabelGetValue(subpointMap, cell, &cellLoc);
2181:       if (cellLoc == 2) continue;
2182:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2183:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2184:       for (cl = 0; cl < closureSize*2; cl += 2) {
2185:         const PetscInt point = closure[cl];
2186:         PetscInt       vertexLoc;

2188:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2189:           ++numCorners;
2190:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2191:           if (vertexLoc == value) closure[faceSize++] = point;
2192:         }
2193:       }
2194:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2195:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2196:       if (faceSize == *nFV) {
2197:         const PetscInt *cells = NULL;
2198:         PetscInt        numCells, nc;

2200:         ++(*numFaces);
2201:         for (cl = 0; cl < faceSize; ++cl) {
2202:           DMLabelSetValue(subpointMap, closure[cl], 0);
2203:         }
2204:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2205:         for (nc = 0; nc < numCells; ++nc) {
2206:           DMLabelSetValue(subpointMap, cells[nc], 2);
2207:         }
2208:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2209:       }
2210:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2211:     }
2212:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2213:   }
2214:   if (subvertexIS) {
2215:     ISRestoreIndices(subvertexIS, &subvertices);
2216:   }
2217:   ISDestroy(&subvertexIS);
2218:   PetscFree2(pStart, pEnd);
2219:   return(0);
2220: }

2222: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2223: {
2224:   IS               subvertexIS = NULL;
2225:   const PetscInt  *subvertices;
2226:   PetscInt        *pStart, *pEnd;
2227:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2228:   PetscErrorCode   ierr;

2231:   DMGetDimension(dm, &dim);
2232:   PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);
2233:   for (d = 0; d <= dim; ++d) {
2234:     DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);
2235:   }
2236:   /* Loop over initial vertices and mark all faces in the collective star() */
2237:   if (vertexLabel) {
2238:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2239:     if (subvertexIS) {
2240:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2241:       ISGetIndices(subvertexIS, &subvertices);
2242:     }
2243:   }
2244:   for (v = 0; v < numSubVerticesInitial; ++v) {
2245:     const PetscInt vertex = subvertices[v];
2246:     PetscInt      *star   = NULL;
2247:     PetscInt       starSize, s, numFaces = 0, f;

2249:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2250:     for (s = 0; s < starSize*2; s += 2) {
2251:       const PetscInt point = star[s];
2252:       PetscInt       faceLoc;

2254:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2255:         if (markedFaces) {
2256:           DMLabelGetValue(vertexLabel, point, &faceLoc);
2257:           if (faceLoc < 0) continue;
2258:         }
2259:         star[numFaces++] = point;
2260:       }
2261:     }
2262:     for (f = 0; f < numFaces; ++f) {
2263:       const PetscInt face    = star[f];
2264:       PetscInt      *closure = NULL;
2265:       PetscInt       closureSize, c;
2266:       PetscInt       faceLoc;

2268:       DMLabelGetValue(subpointMap, face, &faceLoc);
2269:       if (faceLoc == dim-1) continue;
2270:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2271:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2272:       for (c = 0; c < closureSize*2; c += 2) {
2273:         const PetscInt point = closure[c];
2274:         PetscInt       vertexLoc;

2276:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2277:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2278:           if (vertexLoc != value) break;
2279:         }
2280:       }
2281:       if (c == closureSize*2) {
2282:         const PetscInt *support;
2283:         PetscInt        supportSize, s;

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

2288:           for (d = 0; d < dim; ++d) {
2289:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2290:               DMLabelSetValue(subpointMap, point, d);
2291:               break;
2292:             }
2293:           }
2294:         }
2295:         DMPlexGetSupportSize(dm, face, &supportSize);
2296:         DMPlexGetSupport(dm, face, &support);
2297:         for (s = 0; s < supportSize; ++s) {
2298:           DMLabelSetValue(subpointMap, support[s], dim);
2299:         }
2300:       }
2301:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2302:     }
2303:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2304:   }
2305:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2306:   ISDestroy(&subvertexIS);
2307:   PetscFree2(pStart, pEnd);
2308:   return(0);
2309: }

2311: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2312: {
2313:   DMLabel         label = NULL;
2314:   const PetscInt *cone;
2315:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2316:   PetscErrorCode  ierr;

2319:   *numFaces = 0;
2320:   *nFV = 0;
2321:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2322:   *subCells = NULL;
2323:   DMGetDimension(dm, &dim);
2324:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2325:   if (cMax < 0) return(0);
2326:   if (label) {
2327:     for (c = cMax; c < cEnd; ++c) {
2328:       PetscInt val;

2330:       DMLabelGetValue(label, c, &val);
2331:       if (val == value) {
2332:         ++(*numFaces);
2333:         DMPlexGetConeSize(dm, c, &coneSize);
2334:       }
2335:     }
2336:   } else {
2337:     *numFaces = cEnd - cMax;
2338:     DMPlexGetConeSize(dm, cMax, &coneSize);
2339:   }
2340:   PetscMalloc1(*numFaces *2, subCells);
2341:   if (!(*numFaces)) return(0);
2342:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2343:   for (c = cMax; c < cEnd; ++c) {
2344:     const PetscInt *cells;
2345:     PetscInt        numCells;

2347:     if (label) {
2348:       PetscInt val;

2350:       DMLabelGetValue(label, c, &val);
2351:       if (val != value) continue;
2352:     }
2353:     DMPlexGetCone(dm, c, &cone);
2354:     for (p = 0; p < *nFV; ++p) {
2355:       DMLabelSetValue(subpointMap, cone[p], 0);
2356:     }
2357:     /* Negative face */
2358:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2359:     /* Not true in parallel
2360:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2361:     for (p = 0; p < numCells; ++p) {
2362:       DMLabelSetValue(subpointMap, cells[p], 2);
2363:       (*subCells)[subc++] = cells[p];
2364:     }
2365:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2366:     /* Positive face is not included */
2367:   }
2368:   return(0);
2369: }

2371: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2372: {
2373:   PetscInt      *pStart, *pEnd;
2374:   PetscInt       dim, cMax, cEnd, c, d;

2378:   DMGetDimension(dm, &dim);
2379:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2380:   if (cMax < 0) return(0);
2381:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2382:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2383:   for (c = cMax; c < cEnd; ++c) {
2384:     const PetscInt *cone;
2385:     PetscInt       *closure = NULL;
2386:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2388:     if (label) {
2389:       DMLabelGetValue(label, c, &val);
2390:       if (val != value) continue;
2391:     }
2392:     DMPlexGetConeSize(dm, c, &coneSize);
2393:     DMPlexGetCone(dm, c, &cone);
2394:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2395:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2396:     /* Negative face */
2397:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2398:     for (cl = 0; cl < closureSize*2; cl += 2) {
2399:       const PetscInt point = closure[cl];

2401:       for (d = 0; d <= dim; ++d) {
2402:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2403:           DMLabelSetValue(subpointMap, point, d);
2404:           break;
2405:         }
2406:       }
2407:     }
2408:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2409:     /* Cells -- positive face is not included */
2410:     for (cl = 0; cl < 1; ++cl) {
2411:       const PetscInt *support;
2412:       PetscInt        supportSize, s;

2414:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2415:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2416:       DMPlexGetSupport(dm, cone[cl], &support);
2417:       for (s = 0; s < supportSize; ++s) {
2418:         DMLabelSetValue(subpointMap, support[s], dim);
2419:       }
2420:     }
2421:   }
2422:   PetscFree2(pStart, pEnd);
2423:   return(0);
2424: }

2426: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2427: {
2428:   MPI_Comm       comm;
2429:   PetscBool      posOrient = PETSC_FALSE;
2430:   const PetscInt debug     = 0;
2431:   PetscInt       cellDim, faceSize, f;

2435:   PetscObjectGetComm((PetscObject)dm,&comm);
2436:   DMGetDimension(dm, &cellDim);
2437:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2439:   if (cellDim == 1 && numCorners == 2) {
2440:     /* Triangle */
2441:     faceSize  = numCorners-1;
2442:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2443:   } else if (cellDim == 2 && numCorners == 3) {
2444:     /* Triangle */
2445:     faceSize  = numCorners-1;
2446:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2447:   } else if (cellDim == 3 && numCorners == 4) {
2448:     /* Tetrahedron */
2449:     faceSize  = numCorners-1;
2450:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2451:   } else if (cellDim == 1 && numCorners == 3) {
2452:     /* Quadratic line */
2453:     faceSize  = 1;
2454:     posOrient = PETSC_TRUE;
2455:   } else if (cellDim == 2 && numCorners == 4) {
2456:     /* Quads */
2457:     faceSize = 2;
2458:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2459:       posOrient = PETSC_TRUE;
2460:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2461:       posOrient = PETSC_TRUE;
2462:     } else {
2463:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2464:         posOrient = PETSC_FALSE;
2465:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2466:     }
2467:   } else if (cellDim == 2 && numCorners == 6) {
2468:     /* Quadratic triangle (I hate this) */
2469:     /* Edges are determined by the first 2 vertices (corners of edges) */
2470:     const PetscInt faceSizeTri = 3;
2471:     PetscInt       sortedIndices[3], i, iFace;
2472:     PetscBool      found                    = PETSC_FALSE;
2473:     PetscInt       faceVerticesTriSorted[9] = {
2474:       0, 3,  4, /* bottom */
2475:       1, 4,  5, /* right */
2476:       2, 3,  5, /* left */
2477:     };
2478:     PetscInt       faceVerticesTri[9] = {
2479:       0, 3,  4, /* bottom */
2480:       1, 4,  5, /* right */
2481:       2, 5,  3, /* left */
2482:     };

2484:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2485:     PetscSortInt(faceSizeTri, sortedIndices);
2486:     for (iFace = 0; iFace < 3; ++iFace) {
2487:       const PetscInt ii = iFace*faceSizeTri;
2488:       PetscInt       fVertex, cVertex;

2490:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2491:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2492:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2493:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2494:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2495:               faceVertices[fVertex] = origVertices[cVertex];
2496:               break;
2497:             }
2498:           }
2499:         }
2500:         found = PETSC_TRUE;
2501:         break;
2502:       }
2503:     }
2504:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2505:     if (posOriented) *posOriented = PETSC_TRUE;
2506:     return(0);
2507:   } else if (cellDim == 2 && numCorners == 9) {
2508:     /* Quadratic quad (I hate this) */
2509:     /* Edges are determined by the first 2 vertices (corners of edges) */
2510:     const PetscInt faceSizeQuad = 3;
2511:     PetscInt       sortedIndices[3], i, iFace;
2512:     PetscBool      found                      = PETSC_FALSE;
2513:     PetscInt       faceVerticesQuadSorted[12] = {
2514:       0, 1,  4, /* bottom */
2515:       1, 2,  5, /* right */
2516:       2, 3,  6, /* top */
2517:       0, 3,  7, /* left */
2518:     };
2519:     PetscInt       faceVerticesQuad[12] = {
2520:       0, 1,  4, /* bottom */
2521:       1, 2,  5, /* right */
2522:       2, 3,  6, /* top */
2523:       3, 0,  7, /* left */
2524:     };

2526:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2527:     PetscSortInt(faceSizeQuad, sortedIndices);
2528:     for (iFace = 0; iFace < 4; ++iFace) {
2529:       const PetscInt ii = iFace*faceSizeQuad;
2530:       PetscInt       fVertex, cVertex;

2532:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2533:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2534:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2535:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2536:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2537:               faceVertices[fVertex] = origVertices[cVertex];
2538:               break;
2539:             }
2540:           }
2541:         }
2542:         found = PETSC_TRUE;
2543:         break;
2544:       }
2545:     }
2546:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2547:     if (posOriented) *posOriented = PETSC_TRUE;
2548:     return(0);
2549:   } else if (cellDim == 3 && numCorners == 8) {
2550:     /* Hexes
2551:        A hex is two oriented quads with the normal of the first
2552:        pointing up at the second.

2554:           7---6
2555:          /|  /|
2556:         4---5 |
2557:         | 1-|-2
2558:         |/  |/
2559:         0---3

2561:         Faces are determined by the first 4 vertices (corners of faces) */
2562:     const PetscInt faceSizeHex = 4;
2563:     PetscInt       sortedIndices[4], i, iFace;
2564:     PetscBool      found                     = PETSC_FALSE;
2565:     PetscInt       faceVerticesHexSorted[24] = {
2566:       0, 1, 2, 3,  /* bottom */
2567:       4, 5, 6, 7,  /* top */
2568:       0, 3, 4, 5,  /* front */
2569:       2, 3, 5, 6,  /* right */
2570:       1, 2, 6, 7,  /* back */
2571:       0, 1, 4, 7,  /* left */
2572:     };
2573:     PetscInt       faceVerticesHex[24] = {
2574:       1, 2, 3, 0,  /* bottom */
2575:       4, 5, 6, 7,  /* top */
2576:       0, 3, 5, 4,  /* front */
2577:       3, 2, 6, 5,  /* right */
2578:       2, 1, 7, 6,  /* back */
2579:       1, 0, 4, 7,  /* left */
2580:     };

2582:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2583:     PetscSortInt(faceSizeHex, sortedIndices);
2584:     for (iFace = 0; iFace < 6; ++iFace) {
2585:       const PetscInt ii = iFace*faceSizeHex;
2586:       PetscInt       fVertex, cVertex;

2588:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2589:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2590:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2591:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2592:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2593:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2594:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2595:               faceVertices[fVertex] = origVertices[cVertex];
2596:               break;
2597:             }
2598:           }
2599:         }
2600:         found = PETSC_TRUE;
2601:         break;
2602:       }
2603:     }
2604:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2605:     if (posOriented) *posOriented = PETSC_TRUE;
2606:     return(0);
2607:   } else if (cellDim == 3 && numCorners == 10) {
2608:     /* Quadratic tet */
2609:     /* Faces are determined by the first 3 vertices (corners of faces) */
2610:     const PetscInt faceSizeTet = 6;
2611:     PetscInt       sortedIndices[6], i, iFace;
2612:     PetscBool      found                     = PETSC_FALSE;
2613:     PetscInt       faceVerticesTetSorted[24] = {
2614:       0, 1, 2,  6, 7, 8, /* bottom */
2615:       0, 3, 4,  6, 7, 9,  /* front */
2616:       1, 4, 5,  7, 8, 9,  /* right */
2617:       2, 3, 5,  6, 8, 9,  /* left */
2618:     };
2619:     PetscInt       faceVerticesTet[24] = {
2620:       0, 1, 2,  6, 7, 8, /* bottom */
2621:       0, 4, 3,  6, 7, 9,  /* front */
2622:       1, 5, 4,  7, 8, 9,  /* right */
2623:       2, 3, 5,  8, 6, 9,  /* left */
2624:     };

2626:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2627:     PetscSortInt(faceSizeTet, sortedIndices);
2628:     for (iFace=0; iFace < 4; ++iFace) {
2629:       const PetscInt ii = iFace*faceSizeTet;
2630:       PetscInt       fVertex, cVertex;

2632:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2633:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2634:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2635:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2636:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2637:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2638:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2639:               faceVertices[fVertex] = origVertices[cVertex];
2640:               break;
2641:             }
2642:           }
2643:         }
2644:         found = PETSC_TRUE;
2645:         break;
2646:       }
2647:     }
2648:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2649:     if (posOriented) *posOriented = PETSC_TRUE;
2650:     return(0);
2651:   } else if (cellDim == 3 && numCorners == 27) {
2652:     /* Quadratic hexes (I hate this)
2653:        A hex is two oriented quads with the normal of the first
2654:        pointing up at the second.

2656:          7---6
2657:         /|  /|
2658:        4---5 |
2659:        | 3-|-2
2660:        |/  |/
2661:        0---1

2663:        Faces are determined by the first 4 vertices (corners of faces) */
2664:     const PetscInt faceSizeQuadHex = 9;
2665:     PetscInt       sortedIndices[9], i, iFace;
2666:     PetscBool      found                         = PETSC_FALSE;
2667:     PetscInt       faceVerticesQuadHexSorted[54] = {
2668:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2669:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2670:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2671:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2672:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2673:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2674:     };
2675:     PetscInt       faceVerticesQuadHex[54] = {
2676:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2677:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2678:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2679:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2680:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2681:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2682:     };

2684:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2685:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2686:     for (iFace = 0; iFace < 6; ++iFace) {
2687:       const PetscInt ii = iFace*faceSizeQuadHex;
2688:       PetscInt       fVertex, cVertex;

2690:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2691:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2692:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2693:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2694:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2695:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2696:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2697:               faceVertices[fVertex] = origVertices[cVertex];
2698:               break;
2699:             }
2700:           }
2701:         }
2702:         found = PETSC_TRUE;
2703:         break;
2704:       }
2705:     }
2706:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2707:     if (posOriented) *posOriented = PETSC_TRUE;
2708:     return(0);
2709:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2710:   if (!posOrient) {
2711:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2712:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2713:   } else {
2714:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2715:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2716:   }
2717:   if (posOriented) *posOriented = posOrient;
2718:   return(0);
2719: }

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

2725:   Not collective

2727:   Input Parameters:
2728: + dm           - The original mesh
2729: . cell         - The cell mesh point
2730: . faceSize     - The number of vertices on the face
2731: . face         - The face vertices
2732: . numCorners   - The number of vertices on the cell
2733: . indices      - Local numbering of face vertices in cell cone
2734: - origVertices - Original face vertices

2736:   Output Parameter:
2737: + faceVertices - The face vertices properly oriented
2738: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2740:   Level: developer

2742: .seealso: DMPlexGetCone()
2743: @*/
2744: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2745: {
2746:   const PetscInt *cone = NULL;
2747:   PetscInt        coneSize, v, f, v2;
2748:   PetscInt        oppositeVertex = -1;
2749:   PetscErrorCode  ierr;

2752:   DMPlexGetConeSize(dm, cell, &coneSize);
2753:   DMPlexGetCone(dm, cell, &cone);
2754:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2755:     PetscBool found = PETSC_FALSE;

2757:     for (f = 0; f < faceSize; ++f) {
2758:       if (face[f] == cone[v]) {
2759:         found = PETSC_TRUE; break;
2760:       }
2761:     }
2762:     if (found) {
2763:       indices[v2]      = v;
2764:       origVertices[v2] = cone[v];
2765:       ++v2;
2766:     } else {
2767:       oppositeVertex = v;
2768:     }
2769:   }
2770:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2771:   return(0);
2772: }

2774: /*
2775:   DMPlexInsertFace_Internal - Puts a face into the mesh

2777:   Not collective

2779:   Input Parameters:
2780:   + dm              - The DMPlex
2781:   . numFaceVertex   - The number of vertices in the face
2782:   . faceVertices    - The vertices in the face for dm
2783:   . subfaceVertices - The vertices in the face for subdm
2784:   . numCorners      - The number of vertices in the cell
2785:   . cell            - A cell in dm containing the face
2786:   . subcell         - A cell in subdm containing the face
2787:   . firstFace       - First face in the mesh
2788:   - newFacePoint    - Next face in the mesh

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

2793:   Level: developer
2794: */
2795: 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)
2796: {
2797:   MPI_Comm        comm;
2798:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2799:   const PetscInt *faces;
2800:   PetscInt        numFaces, coneSize;
2801:   PetscErrorCode  ierr;

2804:   PetscObjectGetComm((PetscObject)dm,&comm);
2805:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2806:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2807: #if 0
2808:   /* Cannot use this because support() has not been constructed yet */
2809:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2810: #else
2811:   {
2812:     PetscInt f;

2814:     numFaces = 0;
2815:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2816:     for (f = firstFace; f < *newFacePoint; ++f) {
2817:       PetscInt dof, off, d;

2819:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2820:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2821:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2822:       for (d = 0; d < dof; ++d) {
2823:         const PetscInt p = submesh->cones[off+d];
2824:         PetscInt       v;

2826:         for (v = 0; v < numFaceVertices; ++v) {
2827:           if (subfaceVertices[v] == p) break;
2828:         }
2829:         if (v == numFaceVertices) break;
2830:       }
2831:       if (d == dof) {
2832:         numFaces               = 1;
2833:         ((PetscInt*) faces)[0] = f;
2834:       }
2835:     }
2836:   }
2837: #endif
2838:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2839:   else if (numFaces == 1) {
2840:     /* Add the other cell neighbor for this face */
2841:     DMPlexSetCone(subdm, subcell, faces);
2842:   } else {
2843:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2844:     PetscBool posOriented;

2846:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2847:     origVertices        = &orientedVertices[numFaceVertices];
2848:     indices             = &orientedVertices[numFaceVertices*2];
2849:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2850:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2851:     /* TODO: I know that routine should return a permutation, not the indices */
2852:     for (v = 0; v < numFaceVertices; ++v) {
2853:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2854:       for (ov = 0; ov < numFaceVertices; ++ov) {
2855:         if (orientedVertices[ov] == vertex) {
2856:           orientedSubVertices[ov] = subvertex;
2857:           break;
2858:         }
2859:       }
2860:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2861:     }
2862:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2863:     DMPlexSetCone(subdm, subcell, newFacePoint);
2864:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2865:     ++(*newFacePoint);
2866:   }
2867: #if 0
2868:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2869: #else
2870:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2871: #endif
2872:   return(0);
2873: }

2875: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2876: {
2877:   MPI_Comm        comm;
2878:   DMLabel         subpointMap;
2879:   IS              subvertexIS,  subcellIS;
2880:   const PetscInt *subVertices, *subCells;
2881:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2882:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2883:   PetscInt        vStart, vEnd, c, f;
2884:   PetscErrorCode  ierr;

2887:   PetscObjectGetComm((PetscObject)dm,&comm);
2888:   /* Create subpointMap which marks the submesh */
2889:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2890:   DMPlexSetSubpointMap(subdm, subpointMap);
2891:   DMLabelDestroy(&subpointMap);
2892:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2893:   /* Setup chart */
2894:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2895:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2896:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2897:   DMPlexSetVTKCellHeight(subdm, 1);
2898:   /* Set cone sizes */
2899:   firstSubVertex = numSubCells;
2900:   firstSubFace   = numSubCells+numSubVertices;
2901:   newFacePoint   = firstSubFace;
2902:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2903:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2904:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2905:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2906:   for (c = 0; c < numSubCells; ++c) {
2907:     DMPlexSetConeSize(subdm, c, 1);
2908:   }
2909:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2910:     DMPlexSetConeSize(subdm, f, nFV);
2911:   }
2912:   DMSetUp(subdm);
2913:   /* Create face cones */
2914:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2915:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2916:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2917:   for (c = 0; c < numSubCells; ++c) {
2918:     const PetscInt cell    = subCells[c];
2919:     const PetscInt subcell = c;
2920:     PetscInt      *closure = NULL;
2921:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2923:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2924:     for (cl = 0; cl < closureSize*2; cl += 2) {
2925:       const PetscInt point = closure[cl];
2926:       PetscInt       subVertex;

2928:       if ((point >= vStart) && (point < vEnd)) {
2929:         ++numCorners;
2930:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2931:         if (subVertex >= 0) {
2932:           closure[faceSize] = point;
2933:           subface[faceSize] = firstSubVertex+subVertex;
2934:           ++faceSize;
2935:         }
2936:       }
2937:     }
2938:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2939:     if (faceSize == nFV) {
2940:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2941:     }
2942:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2943:   }
2944:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2945:   DMPlexSymmetrize(subdm);
2946:   DMPlexStratify(subdm);
2947:   /* Build coordinates */
2948:   {
2949:     PetscSection coordSection, subCoordSection;
2950:     Vec          coordinates, subCoordinates;
2951:     PetscScalar *coords, *subCoords;
2952:     PetscInt     numComp, coordSize, v;
2953:     const char  *name;

2955:     DMGetCoordinateSection(dm, &coordSection);
2956:     DMGetCoordinatesLocal(dm, &coordinates);
2957:     DMGetCoordinateSection(subdm, &subCoordSection);
2958:     PetscSectionSetNumFields(subCoordSection, 1);
2959:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2960:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2961:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2962:     for (v = 0; v < numSubVertices; ++v) {
2963:       const PetscInt vertex    = subVertices[v];
2964:       const PetscInt subvertex = firstSubVertex+v;
2965:       PetscInt       dof;

2967:       PetscSectionGetDof(coordSection, vertex, &dof);
2968:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2969:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2970:     }
2971:     PetscSectionSetUp(subCoordSection);
2972:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2973:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2974:     PetscObjectGetName((PetscObject)coordinates,&name);
2975:     PetscObjectSetName((PetscObject)subCoordinates,name);
2976:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2977:     VecSetType(subCoordinates,VECSTANDARD);
2978:     if (coordSize) {
2979:       VecGetArray(coordinates,    &coords);
2980:       VecGetArray(subCoordinates, &subCoords);
2981:       for (v = 0; v < numSubVertices; ++v) {
2982:         const PetscInt vertex    = subVertices[v];
2983:         const PetscInt subvertex = firstSubVertex+v;
2984:         PetscInt       dof, off, sdof, soff, d;

2986:         PetscSectionGetDof(coordSection, vertex, &dof);
2987:         PetscSectionGetOffset(coordSection, vertex, &off);
2988:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2989:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2990:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2991:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2992:       }
2993:       VecRestoreArray(coordinates,    &coords);
2994:       VecRestoreArray(subCoordinates, &subCoords);
2995:     }
2996:     DMSetCoordinatesLocal(subdm, subCoordinates);
2997:     VecDestroy(&subCoordinates);
2998:   }
2999:   /* Cleanup */
3000:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3001:   ISDestroy(&subvertexIS);
3002:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
3003:   ISDestroy(&subcellIS);
3004:   return(0);
3005: }

3007: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3008: {
3009:   PetscInt       subPoint;

3012:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
3013:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
3014: }

3016: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3017: {
3018:   MPI_Comm         comm;
3019:   DMLabel          subpointMap;
3020:   IS              *subpointIS;
3021:   const PetscInt **subpoints;
3022:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3023:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
3024:   PetscMPIInt      rank;
3025:   PetscErrorCode   ierr;

3028:   PetscObjectGetComm((PetscObject)dm,&comm);
3029:   MPI_Comm_rank(comm, &rank);
3030:   /* Create subpointMap which marks the submesh */
3031:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3032:   DMPlexSetSubpointMap(subdm, subpointMap);
3033:   if (cellHeight) {
3034:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
3035:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
3036:   } else {
3037:     DMLabel         depth;
3038:     IS              pointIS;
3039:     const PetscInt *points;
3040:     PetscInt        numPoints;

3042:     DMPlexGetDepthLabel(dm, &depth);
3043:     DMLabelGetStratumSize(label, value, &numPoints);
3044:     DMLabelGetStratumIS(label, value, &pointIS);
3045:     ISGetIndices(pointIS, &points);
3046:     for (p = 0; p < numPoints; ++p) {
3047:       PetscInt *closure = NULL;
3048:       PetscInt  closureSize, c, pdim;

3050:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3051:       for (c = 0; c < closureSize*2; c += 2) {
3052:         DMLabelGetValue(depth, closure[c], &pdim);
3053:         DMLabelSetValue(subpointMap, closure[c], pdim);
3054:       }
3055:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3056:     }
3057:     ISRestoreIndices(pointIS, &points);
3058:     ISDestroy(&pointIS);
3059:   }
3060:   /* Setup chart */
3061:   DMGetDimension(dm, &dim);
3062:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3063:   for (d = 0; d <= dim; ++d) {
3064:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3065:     totSubPoints += numSubPoints[d];
3066:   }
3067:   DMPlexSetChart(subdm, 0, totSubPoints);
3068:   DMPlexSetVTKCellHeight(subdm, cellHeight);
3069:   /* Set cone sizes */
3070:   firstSubPoint[dim] = 0;
3071:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
3072:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
3073:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3074:   for (d = 0; d <= dim; ++d) {
3075:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3076:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3077:   }
3078:   /* We do not want this label automatically computed, instead we compute it here */
3079:   DMCreateLabel(subdm, "celltype");
3080:   for (d = 0; d <= dim; ++d) {
3081:     for (p = 0; p < numSubPoints[d]; ++p) {
3082:       const PetscInt  point    = subpoints[d][p];
3083:       const PetscInt  subpoint = firstSubPoint[d] + p;
3084:       const PetscInt *cone;
3085:       PetscInt        coneSize, coneSizeNew, c, val;
3086:       DMPolytopeType  ct;

3088:       DMPlexGetConeSize(dm, point, &coneSize);
3089:       DMPlexSetConeSize(subdm, subpoint, coneSize);
3090:       DMPlexGetCellType(dm, point, &ct);
3091:       DMPlexSetCellType(subdm, subpoint, ct);
3092:       if (cellHeight && (d == dim)) {
3093:         DMPlexGetCone(dm, point, &cone);
3094:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3095:           DMLabelGetValue(subpointMap, cone[c], &val);
3096:           if (val >= 0) coneSizeNew++;
3097:         }
3098:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3099:         DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3100:       }
3101:     }
3102:   }
3103:   DMLabelDestroy(&subpointMap);
3104:   DMSetUp(subdm);
3105:   /* Set cones */
3106:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3107:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3108:   for (d = 0; d <= dim; ++d) {
3109:     for (p = 0; p < numSubPoints[d]; ++p) {
3110:       const PetscInt  point    = subpoints[d][p];
3111:       const PetscInt  subpoint = firstSubPoint[d] + p;
3112:       const PetscInt *cone, *ornt;
3113:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

3115:       if (d == dim-1) {
3116:         const PetscInt *support, *cone, *ornt;
3117:         PetscInt        supportSize, coneSize, s, subc;

3119:         DMPlexGetSupport(dm, point, &support);
3120:         DMPlexGetSupportSize(dm, point, &supportSize);
3121:         for (s = 0; s < supportSize; ++s) {
3122:           PetscBool isHybrid;

3124:           DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3125:           if (!isHybrid) continue;
3126:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3127:           if (subc >= 0) {
3128:             const PetscInt ccell = subpoints[d+1][subc];

3130:             DMPlexGetCone(dm, ccell, &cone);
3131:             DMPlexGetConeSize(dm, ccell, &coneSize);
3132:             DMPlexGetConeOrientation(dm, ccell, &ornt);
3133:             for (c = 0; c < coneSize; ++c) {
3134:               if (cone[c] == point) {
3135:                 fornt = ornt[c];
3136:                 break;
3137:               }
3138:             }
3139:             break;
3140:           }
3141:         }
3142:       }
3143:       DMPlexGetConeSize(dm, point, &coneSize);
3144:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3145:       DMPlexGetCone(dm, point, &cone);
3146:       DMPlexGetConeOrientation(dm, point, &ornt);
3147:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3148:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3149:         if (subc >= 0) {
3150:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3151:           orntNew[coneSizeNew] = ornt[c];
3152:           ++coneSizeNew;
3153:         }
3154:       }
3155:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3156:       if (fornt < 0) {
3157:         /* This should be replaced by a call to DMPlexReverseCell() */
3158: #if 0
3159:         DMPlexReverseCell(subdm, subpoint);
3160: #else
3161:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3162:           PetscInt faceSize, tmp;

3164:           tmp        = coneNew[c];
3165:           coneNew[c] = coneNew[coneSizeNew-1-c];
3166:           coneNew[coneSizeNew-1-c] = tmp;
3167:           DMPlexGetConeSize(dm, cone[c], &faceSize);
3168:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3169:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3170:           orntNew[coneSizeNew-1-c] = tmp;
3171:         }
3172:       }
3173:       DMPlexSetCone(subdm, subpoint, coneNew);
3174:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3175: #endif
3176:     }
3177:   }
3178:   PetscFree2(coneNew,orntNew);
3179:   DMPlexSymmetrize(subdm);
3180:   DMPlexStratify(subdm);
3181:   /* Build coordinates */
3182:   {
3183:     PetscSection coordSection, subCoordSection;
3184:     Vec          coordinates, subCoordinates;
3185:     PetscScalar *coords, *subCoords;
3186:     PetscInt     cdim, numComp, coordSize;
3187:     const char  *name;

3189:     DMGetCoordinateDim(dm, &cdim);
3190:     DMGetCoordinateSection(dm, &coordSection);
3191:     DMGetCoordinatesLocal(dm, &coordinates);
3192:     DMGetCoordinateSection(subdm, &subCoordSection);
3193:     PetscSectionSetNumFields(subCoordSection, 1);
3194:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3195:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3196:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3197:     for (v = 0; v < numSubPoints[0]; ++v) {
3198:       const PetscInt vertex    = subpoints[0][v];
3199:       const PetscInt subvertex = firstSubPoint[0]+v;
3200:       PetscInt       dof;

3202:       PetscSectionGetDof(coordSection, vertex, &dof);
3203:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3204:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3205:     }
3206:     PetscSectionSetUp(subCoordSection);
3207:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3208:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3209:     PetscObjectGetName((PetscObject)coordinates,&name);
3210:     PetscObjectSetName((PetscObject)subCoordinates,name);
3211:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3212:     VecSetBlockSize(subCoordinates, cdim);
3213:     VecSetType(subCoordinates,VECSTANDARD);
3214:     VecGetArray(coordinates,    &coords);
3215:     VecGetArray(subCoordinates, &subCoords);
3216:     for (v = 0; v < numSubPoints[0]; ++v) {
3217:       const PetscInt vertex    = subpoints[0][v];
3218:       const PetscInt subvertex = firstSubPoint[0]+v;
3219:       PetscInt dof, off, sdof, soff, d;

3221:       PetscSectionGetDof(coordSection, vertex, &dof);
3222:       PetscSectionGetOffset(coordSection, vertex, &off);
3223:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3224:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3225:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3226:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3227:     }
3228:     VecRestoreArray(coordinates,    &coords);
3229:     VecRestoreArray(subCoordinates, &subCoords);
3230:     DMSetCoordinatesLocal(subdm, subCoordinates);
3231:     VecDestroy(&subCoordinates);
3232:   }
3233:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3234:   {
3235:     PetscSF            sfPoint, sfPointSub;
3236:     IS                 subpIS;
3237:     const PetscSFNode *remotePoints;
3238:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3239:     const PetscInt    *localPoints, *subpoints;
3240:     PetscInt          *slocalPoints;
3241:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3242:     PetscMPIInt        rank;

3244:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3245:     DMGetPointSF(dm, &sfPoint);
3246:     DMGetPointSF(subdm, &sfPointSub);
3247:     DMPlexGetChart(dm, &pStart, &pEnd);
3248:     DMPlexGetChart(subdm, NULL, &numSubroots);
3249:     DMPlexCreateSubpointIS(subdm, &subpIS);
3250:     if (subpIS) {
3251:       ISGetIndices(subpIS, &subpoints);
3252:       ISGetLocalSize(subpIS, &numSubpoints);
3253:     }
3254:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3255:     if (numRoots >= 0) {
3256:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3257:       for (p = 0; p < pEnd-pStart; ++p) {
3258:         newLocalPoints[p].rank  = -2;
3259:         newLocalPoints[p].index = -2;
3260:       }
3261:       /* Set subleaves */
3262:       for (l = 0; l < numLeaves; ++l) {
3263:         const PetscInt point    = localPoints[l];
3264:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3266:         if (subpoint < 0) continue;
3267:         newLocalPoints[point-pStart].rank  = rank;
3268:         newLocalPoints[point-pStart].index = subpoint;
3269:         ++numSubleaves;
3270:       }
3271:       /* Must put in owned subpoints */
3272:       for (p = pStart; p < pEnd; ++p) {
3273:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3275:         if (subpoint < 0) {
3276:           newOwners[p-pStart].rank  = -3;
3277:           newOwners[p-pStart].index = -3;
3278:         } else {
3279:           newOwners[p-pStart].rank  = rank;
3280:           newOwners[p-pStart].index = subpoint;
3281:         }
3282:       }
3283:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3284:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3285:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3286:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3287:       PetscMalloc1(numSubleaves, &slocalPoints);
3288:       PetscMalloc1(numSubleaves, &sremotePoints);
3289:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3290:         const PetscInt point    = localPoints[l];
3291:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3293:         if (subpoint < 0) continue;
3294:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3295:         slocalPoints[sl]        = subpoint;
3296:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3297:         sremotePoints[sl].index = newLocalPoints[point].index;
3298:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3299:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3300:         ++sl;
3301:       }
3302:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3303:       PetscFree2(newLocalPoints,newOwners);
3304:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3305:     }
3306:     if (subpIS) {
3307:       ISRestoreIndices(subpIS, &subpoints);
3308:       ISDestroy(&subpIS);
3309:     }
3310:   }
3311:   /* Cleanup */
3312:   for (d = 0; d <= dim; ++d) {
3313:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3314:     ISDestroy(&subpointIS[d]);
3315:   }
3316:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3317:   return(0);
3318: }

3320: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3321: {

3325:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3326:   return(0);
3327: }

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

3332:   Input Parameters:
3333: + dm           - The original mesh
3334: . vertexLabel  - The DMLabel marking points contained in the surface
3335: . value        - The label value to use
3336: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

3338:   Output Parameter:
3339: . subdm - The surface mesh

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

3343:   Level: developer

3345: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3346: @*/
3347: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3348: {
3349:   DMPlexInterpolatedFlag interpolated;
3350:   PetscInt       dim, cdim;

3356:   DMGetDimension(dm, &dim);
3357:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3358:   DMSetType(*subdm, DMPLEX);
3359:   DMSetDimension(*subdm, dim-1);
3360:   DMGetCoordinateDim(dm, &cdim);
3361:   DMSetCoordinateDim(*subdm, cdim);
3362:   DMPlexIsInterpolated(dm, &interpolated);
3363:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
3364:   if (interpolated) {
3365:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3366:   } else {
3367:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3368:   }
3369:   return(0);
3370: }

3372: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3373: {
3374:   MPI_Comm        comm;
3375:   DMLabel         subpointMap;
3376:   IS              subvertexIS;
3377:   const PetscInt *subVertices;
3378:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3379:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3380:   PetscInt        c, f;
3381:   PetscErrorCode  ierr;

3384:   PetscObjectGetComm((PetscObject)dm, &comm);
3385:   /* Create subpointMap which marks the submesh */
3386:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3387:   DMPlexSetSubpointMap(subdm, subpointMap);
3388:   DMLabelDestroy(&subpointMap);
3389:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3390:   /* Setup chart */
3391:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3392:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3393:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3394:   DMPlexSetVTKCellHeight(subdm, 1);
3395:   /* Set cone sizes */
3396:   firstSubVertex = numSubCells;
3397:   firstSubFace   = numSubCells+numSubVertices;
3398:   newFacePoint   = firstSubFace;
3399:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3400:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3401:   for (c = 0; c < numSubCells; ++c) {
3402:     DMPlexSetConeSize(subdm, c, 1);
3403:   }
3404:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3405:     DMPlexSetConeSize(subdm, f, nFV);
3406:   }
3407:   DMSetUp(subdm);
3408:   /* Create face cones */
3409:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3410:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3411:   for (c = 0; c < numSubCells; ++c) {
3412:     const PetscInt  cell    = subCells[c];
3413:     const PetscInt  subcell = c;
3414:     const PetscInt *cone, *cells;
3415:     PetscBool       isHybrid;
3416:     PetscInt        numCells, subVertex, p, v;

3418:     DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3419:     if (!isHybrid) continue;
3420:     DMPlexGetCone(dm, cell, &cone);
3421:     for (v = 0; v < nFV; ++v) {
3422:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3423:       subface[v] = firstSubVertex+subVertex;
3424:     }
3425:     DMPlexSetCone(subdm, newFacePoint, subface);
3426:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3427:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3428:     /* Not true in parallel
3429:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3430:     for (p = 0; p < numCells; ++p) {
3431:       PetscInt  negsubcell;
3432:       PetscBool isHybrid;

3434:       DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3435:       if (isHybrid) continue;
3436:       /* I know this is a crap search */
3437:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3438:         if (subCells[negsubcell] == cells[p]) break;
3439:       }
3440:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3441:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3442:     }
3443:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3444:     ++newFacePoint;
3445:   }
3446:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3447:   DMPlexSymmetrize(subdm);
3448:   DMPlexStratify(subdm);
3449:   /* Build coordinates */
3450:   {
3451:     PetscSection coordSection, subCoordSection;
3452:     Vec          coordinates, subCoordinates;
3453:     PetscScalar *coords, *subCoords;
3454:     PetscInt     cdim, numComp, coordSize, v;
3455:     const char  *name;

3457:     DMGetCoordinateDim(dm, &cdim);
3458:     DMGetCoordinateSection(dm, &coordSection);
3459:     DMGetCoordinatesLocal(dm, &coordinates);
3460:     DMGetCoordinateSection(subdm, &subCoordSection);
3461:     PetscSectionSetNumFields(subCoordSection, 1);
3462:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3463:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3464:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3465:     for (v = 0; v < numSubVertices; ++v) {
3466:       const PetscInt vertex    = subVertices[v];
3467:       const PetscInt subvertex = firstSubVertex+v;
3468:       PetscInt       dof;

3470:       PetscSectionGetDof(coordSection, vertex, &dof);
3471:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3472:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3473:     }
3474:     PetscSectionSetUp(subCoordSection);
3475:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3476:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3477:     PetscObjectGetName((PetscObject)coordinates,&name);
3478:     PetscObjectSetName((PetscObject)subCoordinates,name);
3479:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3480:     VecSetBlockSize(subCoordinates, cdim);
3481:     VecSetType(subCoordinates,VECSTANDARD);
3482:     VecGetArray(coordinates,    &coords);
3483:     VecGetArray(subCoordinates, &subCoords);
3484:     for (v = 0; v < numSubVertices; ++v) {
3485:       const PetscInt vertex    = subVertices[v];
3486:       const PetscInt subvertex = firstSubVertex+v;
3487:       PetscInt       dof, off, sdof, soff, d;

3489:       PetscSectionGetDof(coordSection, vertex, &dof);
3490:       PetscSectionGetOffset(coordSection, vertex, &off);
3491:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3492:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3493:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3494:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3495:     }
3496:     VecRestoreArray(coordinates,    &coords);
3497:     VecRestoreArray(subCoordinates, &subCoords);
3498:     DMSetCoordinatesLocal(subdm, subCoordinates);
3499:     VecDestroy(&subCoordinates);
3500:   }
3501:   /* Build SF */
3502:   CHKMEMQ;
3503:   {
3504:     PetscSF            sfPoint, sfPointSub;
3505:     const PetscSFNode *remotePoints;
3506:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3507:     const PetscInt    *localPoints;
3508:     PetscInt          *slocalPoints;
3509:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3510:     PetscMPIInt        rank;

3512:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3513:     DMGetPointSF(dm, &sfPoint);
3514:     DMGetPointSF(subdm, &sfPointSub);
3515:     DMPlexGetChart(dm, &pStart, &pEnd);
3516:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3517:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3518:     if (numRoots >= 0) {
3519:       /* Only vertices should be shared */
3520:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3521:       for (p = 0; p < pEnd-pStart; ++p) {
3522:         newLocalPoints[p].rank  = -2;
3523:         newLocalPoints[p].index = -2;
3524:       }
3525:       /* Set subleaves */
3526:       for (l = 0; l < numLeaves; ++l) {
3527:         const PetscInt point    = localPoints[l];
3528:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3530:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3531:         if (subPoint < 0) continue;
3532:         newLocalPoints[point-pStart].rank  = rank;
3533:         newLocalPoints[point-pStart].index = subPoint;
3534:         ++numSubLeaves;
3535:       }
3536:       /* Must put in owned subpoints */
3537:       for (p = pStart; p < pEnd; ++p) {
3538:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3540:         if (subPoint < 0) {
3541:           newOwners[p-pStart].rank  = -3;
3542:           newOwners[p-pStart].index = -3;
3543:         } else {
3544:           newOwners[p-pStart].rank  = rank;
3545:           newOwners[p-pStart].index = subPoint;
3546:         }
3547:       }
3548:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3549:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3550:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3551:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3552:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3553:       PetscMalloc1(numSubLeaves, &sremotePoints);
3554:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3555:         const PetscInt point    = localPoints[l];
3556:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3558:         if (subPoint < 0) continue;
3559:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3560:         slocalPoints[sl]        = subPoint;
3561:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3562:         sremotePoints[sl].index = newLocalPoints[point].index;
3563:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3564:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3565:         ++sl;
3566:       }
3567:       PetscFree2(newLocalPoints,newOwners);
3568:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3569:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3570:     }
3571:   }
3572:   CHKMEMQ;
3573:   /* Cleanup */
3574:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3575:   ISDestroy(&subvertexIS);
3576:   PetscFree(subCells);
3577:   return(0);
3578: }

3580: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3581: {
3582:   DMLabel        label = NULL;

3586:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3587:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3588:   return(0);
3589: }

3591: /*@C
3592:   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.

3594:   Input Parameters:
3595: + dm          - The original mesh
3596: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3597: . label       - A label name, or NULL
3598: - value  - A label value

3600:   Output Parameter:
3601: . subdm - The surface mesh

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

3605:   Level: developer

3607: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3608: @*/
3609: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3610: {
3611:   PetscInt       dim, cdim, depth;

3617:   DMGetDimension(dm, &dim);
3618:   DMPlexGetDepth(dm, &depth);
3619:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3620:   DMSetType(*subdm, DMPLEX);
3621:   DMSetDimension(*subdm, dim-1);
3622:   DMGetCoordinateDim(dm, &cdim);
3623:   DMSetCoordinateDim(*subdm, cdim);
3624:   if (depth == dim) {
3625:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3626:   } else {
3627:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3628:   }
3629:   return(0);
3630: }

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

3635:   Input Parameters:
3636: + dm        - The original mesh
3637: . cellLabel - The DMLabel marking cells contained in the new mesh
3638: - value     - The label value to use

3640:   Output Parameter:
3641: . subdm - The new mesh

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

3645:   Level: developer

3647: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3648: @*/
3649: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3650: {
3651:   PetscInt       dim;

3657:   DMGetDimension(dm, &dim);
3658:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3659:   DMSetType(*subdm, DMPLEX);
3660:   DMSetDimension(*subdm, dim);
3661:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3662:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3663:   return(0);
3664: }

3666: /*@
3667:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3669:   Input Parameter:
3670: . dm - The submesh DM

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

3675:   Level: developer

3677: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3678: @*/
3679: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3680: {
3684:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3685:   return(0);
3686: }

3688: /*@
3689:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3691:   Input Parameters:
3692: + dm - The submesh DM
3693: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3697:   Level: developer

3699: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3700: @*/
3701: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3702: {
3703:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3704:   DMLabel        tmp;

3709:   tmp  = mesh->subpointMap;
3710:   mesh->subpointMap = subpointMap;
3711:   PetscObjectReference((PetscObject) mesh->subpointMap);
3712:   DMLabelDestroy(&tmp);
3713:   return(0);
3714: }

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

3719:   Input Parameter:
3720: . dm - The submesh DM

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

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

3727:   Level: developer

3729: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3730: @*/
3731: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3732: {
3733:   MPI_Comm        comm;
3734:   DMLabel         subpointMap;
3735:   IS              is;
3736:   const PetscInt *opoints;
3737:   PetscInt       *points, *depths;
3738:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3739:   PetscErrorCode  ierr;

3744:   PetscObjectGetComm((PetscObject)dm,&comm);
3745:   *subpointIS = NULL;
3746:   DMPlexGetSubpointMap(dm, &subpointMap);
3747:   DMPlexGetDepth(dm, &depth);
3748:   if (subpointMap && depth >= 0) {
3749:     DMPlexGetChart(dm, &pStart, &pEnd);
3750:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3751:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3752:     depths[0] = depth;
3753:     depths[1] = 0;
3754:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3755:     PetscMalloc1(pEnd, &points);
3756:     for(d = 0, off = 0; d <= depth; ++d) {
3757:       const PetscInt dep = depths[d];

3759:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3760:       DMLabelGetStratumSize(subpointMap, dep, &n);
3761:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3762:         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);
3763:       } else {
3764:         if (!n) {
3765:           if (d == 0) {
3766:             /* Missing cells */
3767:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3768:           } else {
3769:             /* Missing faces */
3770:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3771:           }
3772:         }
3773:       }
3774:       if (n) {
3775:         DMLabelGetStratumIS(subpointMap, dep, &is);
3776:         ISGetIndices(is, &opoints);
3777:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3778:         ISRestoreIndices(is, &opoints);
3779:         ISDestroy(&is);
3780:       }
3781:     }
3782:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3783:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3784:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3785:   }
3786:   return(0);
3787: }

3789: /*@
3790:   DMGetEnclosureRelation - Get the relationship between dmA and dmB

3792:   Input Parameters:
3793: + dmA - The first DM
3794: - dmB - The second DM

3796:   Output Parameter:
3797: . rel - The relation of dmA to dmB

3799:   Level: intermediate

3801: .seealso: DMPlexGetEnclosurePoint()
3802: @*/
3803: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3804: {
3805:   DM             plexA, plexB, sdm;
3806:   DMLabel        spmap;
3807:   PetscInt       pStartA, pEndA, pStartB, pEndB, NpA, NpB;

3812:   *rel = DM_ENC_NONE;
3813:   if (!dmA || !dmB) return(0);
3816:   if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return(0);}
3817:   DMConvert(dmA, DMPLEX, &plexA);
3818:   DMConvert(dmB, DMPLEX, &plexB);
3819:   DMPlexGetChart(plexA, &pStartA, &pEndA);
3820:   DMPlexGetChart(plexB, &pStartB, &pEndB);
3821:   /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
3822:     - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
3823:   if ((pStartA == pStartB) && (pEndA == pEndB)) {
3824:     *rel = DM_ENC_EQUALITY;
3825:     goto end;
3826:   }
3827:   NpA = pEndA - pStartA;
3828:   NpB = pEndB - pStartB;
3829:   if (NpA == NpB) goto end;
3830:   sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
3831:   DMPlexGetSubpointMap(sdm, &spmap);
3832:   if (!spmap) goto end;
3833:   /* TODO Check the space mapped to by subpointMap is same size as dm */
3834:   if (NpA > NpB) {
3835:     *rel = DM_ENC_SUPERMESH;
3836:   } else {
3837:     *rel = DM_ENC_SUBMESH;
3838:   }
3839:   end:
3840:   DMDestroy(&plexA);
3841:   DMDestroy(&plexB);
3842:   return(0);
3843: }

3845: /*@
3846:   DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB

3848:   Input Parameters:
3849: + dmA   - The first DM
3850: . dmB   - The second DM
3851: . etype - The type of enclosure relation that dmA has to dmB
3852: - pB    - A point of dmB

3854:   Output Parameter:
3855: . pA    - The corresponding point of dmA

3857:   Level: intermediate

3859: .seealso: DMGetEnclosureRelation()
3860: @*/
3861: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3862: {
3863:   DM              sdm;
3864:   DMLabel         spmap;
3865:   IS              subpointIS;
3866:   const PetscInt *subpoints;
3867:   PetscInt        numSubpoints;
3868:   PetscErrorCode  ierr;

3871:   /* TODO Cache the IS, making it look like an index */
3872:   switch (etype) {
3873:     case DM_ENC_SUPERMESH:
3874:     sdm  = dmB;
3875:     DMPlexGetSubpointMap(sdm, &spmap);
3876:     DMPlexCreateSubpointIS(sdm, &subpointIS);
3877:     ISGetIndices(subpointIS, &subpoints);
3878:     *pA  = subpoints[pB];
3879:     ISRestoreIndices(subpointIS, &subpoints);
3880:     ISDestroy(&subpointIS);
3881:     break;
3882:     case DM_ENC_SUBMESH:
3883:     sdm  = dmA;
3884:     DMPlexGetSubpointMap(sdm, &spmap);
3885:     DMPlexCreateSubpointIS(sdm, &subpointIS);
3886:     ISGetLocalSize(subpointIS, &numSubpoints);
3887:     ISGetIndices(subpointIS, &subpoints);
3888:     PetscFindInt(pB, numSubpoints, subpoints, pA);
3889:     if (*pA < 0) {
3890:       DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3891:       DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3892:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3893:     }
3894:     ISRestoreIndices(subpointIS, &subpoints);
3895:     ISDestroy(&subpointIS);
3896:     break;
3897:     case DM_ENC_EQUALITY:
3898:     case DM_ENC_NONE:
3899:     *pA = pB;break;
3900:     case DM_ENC_UNKNOWN:
3901:     {
3902:       DMEnclosureType enc;

3904:       DMGetEnclosureRelation(dmA, dmB, &enc);
3905:       DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3906:     }
3907:     break;
3908:     default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3909:   }
3910:   return(0);
3911: }