Actual source code: plexsubmesh.c

petsc-main 2021-04-20
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:     /* Pull point contributions from remote leaves into local roots */
156:     DMLabelGather(label, sfPoint, &lblLeaves);
157:     DMLabelGetValueIS(lblLeaves, &valueIS);
158:     ISGetLocalSize(valueIS, &numValues);
159:     ISGetIndices(valueIS, &values);
160:     for (v = 0; v < numValues; ++v) {
161:       const PetscInt value = values[v];

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

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

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

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

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

199:   Level: developer

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

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

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

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

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

222:   Level: developer

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

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

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

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

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

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

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

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

277:   Level: developer

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

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

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

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

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

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

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

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

335:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

918:   Collective on dm

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

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

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

930:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1567:       if (splitLabel) {
1568:         const PetscInt val = 100 + dep;

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

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

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

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

1616:   Collective on dm

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

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

1626:   Level: developer

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

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

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

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

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

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

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

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

1709:   Output Parameter:
1710: . label - A DMLabel marking all surface points

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

1714:   Level: developer

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

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

1747:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1748: #if 0
1749:     if (supportSize != 2) {
1750:       const PetscInt *lp;
1751:       PetscInt        Nlp, pind;

1753:       /* Check that for a cell with a single support face, that face is in the SF */
1754:       /*   THis check only works for the remote side. We would need root side information */
1755:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1756:       PetscFindInt(points[p], Nlp, lp, &pind);
1757:       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);
1758:     }
1759: #endif
1760:     DMPlexGetSupport(dm, points[p], &support);
1761:     for (s = 0; s < supportSize; ++s) {
1762:       const PetscInt *cone;
1763:       PetscInt        coneSize, c;
1764:       PetscBool       pos;

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

1776:         DMLabelGetValue(label, point, &val);
1777:         if (val == -1) {
1778:           PetscInt *closure = NULL;
1779:           PetscInt  closureSize, cl;

1781:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1782:           for (cl = 0; cl < closureSize*2; cl += 2) {
1783:             const PetscInt clp  = closure[cl];
1784:             PetscInt       bval = -1;

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

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

1823:       DMLabelGetValue(blabel, point, &bval);
1824:       if (bval >= 0) {
1825:         const PetscInt *cone,    *support;
1826:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

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

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

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

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

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

1907:               DMLabelGetValue(label, ccone[cc], &val);
1908:               if (val != -1) continue;
1909:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1910:               for (cl = 0; cl < closureSize*2; cl += 2) {
1911:                 const PetscInt clp = closure[cl];

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

1930:       DMLabelGetValue(label, point, &val);
1931:       if (val == -1) {
1932:         PetscInt      *sstar = NULL;
1933:         PetscInt       sstarSize, ss;
1934:         PetscBool      marked = PETSC_FALSE, isHybrid;

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

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

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

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

2014: /* Check that no cell have all vertices on the fault */
2015: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2016: {
2017:   IS              subpointIS;
2018:   const PetscInt *dmpoints;
2019:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;
2020:   PetscErrorCode  ierr;

2023:   if (!label) return(0);
2024:   DMLabelGetDefaultValue(label, &defaultValue);
2025:   DMPlexGetSubpointIS(subdm, &subpointIS);
2026:   if (!subpointIS) return(0);
2027:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2028:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2029:   ISGetIndices(subpointIS, &dmpoints);
2030:   for (c = cStart; c < cEnd; ++c) {
2031:     PetscBool invalidCell = PETSC_TRUE;
2032:     PetscInt *closure     = NULL;
2033:     PetscInt  closureSize, cl;

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

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


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

2059:   Collective on dm

2061:   Input Parameters:
2062: + dm - The original DM
2063: . label - The label specifying the interface vertices
2064: - bdlabel - The optional label specifying the interface boundary vertices

2066:   Output Parameters:
2067: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2068: . splitLabel - The label containing the split points, or NULL if no output is desired
2069: . dmInterface - The new interface DM, or NULL
2070: - dmHybrid - The new DM with cohesive cells

2072:   Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
2073:   directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2074:   7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2075:   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
2076:   surface also hits vertex 3, it would be 9 (-101) in hybridLabel.

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

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

2085:   Level: developer

2087: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2088: @*/
2089: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2090: {
2091:   DM             idm;
2092:   DMLabel        subpointMap, hlabel, slabel = NULL;
2093:   PetscInt       dim;

2103:   DMGetDimension(dm, &dim);
2104:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2105:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
2106:   DMPlexOrient(idm);
2107:   DMPlexGetSubpointMap(idm, &subpointMap);
2108:   DMLabelDuplicate(subpointMap, &hlabel);
2109:   DMLabelClearStratum(hlabel, dim);
2110:   if (splitLabel) {
2111:     const char *name;
2112:     char        sname[PETSC_MAX_PATH_LEN];

2114:     PetscObjectGetName((PetscObject) hlabel, &name);
2115:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2116:     PetscStrcat(sname, " split");
2117:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2118:   }
2119:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2120:   if (dmInterface) {*dmInterface = idm;}
2121:   else             {DMDestroy(&idm);}
2122:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2123:   if (hybridLabel) *hybridLabel = hlabel;
2124:   else             {DMLabelDestroy(&hlabel);}
2125:   if (splitLabel)  *splitLabel  = slabel;
2126:   return(0);
2127: }

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

2131:      For any marked cell, the marked vertices constitute a single face
2132: */
2133: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2134: {
2135:   IS               subvertexIS = NULL;
2136:   const PetscInt  *subvertices;
2137:   PetscInt        *pStart, *pEnd, pSize;
2138:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
2139:   PetscErrorCode   ierr;

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

2162:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2163:     for (s = 0; s < starSize*2; s += 2) {
2164:       const PetscInt point = star[s];
2165:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2166:     }
2167:     for (c = 0; c < numCells; ++c) {
2168:       const PetscInt cell    = star[c];
2169:       PetscInt      *closure = NULL;
2170:       PetscInt       closureSize, cl;
2171:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

2173:       DMLabelGetValue(subpointMap, cell, &cellLoc);
2174:       if (cellLoc == 2) continue;
2175:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2176:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2177:       for (cl = 0; cl < closureSize*2; cl += 2) {
2178:         const PetscInt point = closure[cl];
2179:         PetscInt       vertexLoc;

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

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

2215: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2216: {
2217:   IS               subvertexIS = NULL;
2218:   const PetscInt  *subvertices;
2219:   PetscInt        *pStart, *pEnd;
2220:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2221:   PetscErrorCode   ierr;

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

2242:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2243:     for (s = 0; s < starSize*2; s += 2) {
2244:       const PetscInt point = star[s];
2245:       PetscInt       faceLoc;

2247:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2248:         if (markedFaces) {
2249:           DMLabelGetValue(vertexLabel, point, &faceLoc);
2250:           if (faceLoc < 0) continue;
2251:         }
2252:         star[numFaces++] = point;
2253:       }
2254:     }
2255:     for (f = 0; f < numFaces; ++f) {
2256:       const PetscInt face    = star[f];
2257:       PetscInt      *closure = NULL;
2258:       PetscInt       closureSize, c;
2259:       PetscInt       faceLoc;

2261:       DMLabelGetValue(subpointMap, face, &faceLoc);
2262:       if (faceLoc == dim-1) continue;
2263:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2264:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2265:       for (c = 0; c < closureSize*2; c += 2) {
2266:         const PetscInt point = closure[c];
2267:         PetscInt       vertexLoc;

2269:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2270:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2271:           if (vertexLoc != value) break;
2272:         }
2273:       }
2274:       if (c == closureSize*2) {
2275:         const PetscInt *support;
2276:         PetscInt        supportSize, s;

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

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

2304: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2305: {
2306:   DMLabel         label = NULL;
2307:   const PetscInt *cone;
2308:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2309:   PetscErrorCode  ierr;

2312:   *numFaces = 0;
2313:   *nFV = 0;
2314:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2315:   *subCells = NULL;
2316:   DMGetDimension(dm, &dim);
2317:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2318:   if (cMax < 0) return(0);
2319:   if (label) {
2320:     for (c = cMax; c < cEnd; ++c) {
2321:       PetscInt val;

2323:       DMLabelGetValue(label, c, &val);
2324:       if (val == value) {
2325:         ++(*numFaces);
2326:         DMPlexGetConeSize(dm, c, &coneSize);
2327:       }
2328:     }
2329:   } else {
2330:     *numFaces = cEnd - cMax;
2331:     DMPlexGetConeSize(dm, cMax, &coneSize);
2332:   }
2333:   PetscMalloc1(*numFaces *2, subCells);
2334:   if (!(*numFaces)) return(0);
2335:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2336:   for (c = cMax; c < cEnd; ++c) {
2337:     const PetscInt *cells;
2338:     PetscInt        numCells;

2340:     if (label) {
2341:       PetscInt val;

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

2364: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2365: {
2366:   PetscInt      *pStart, *pEnd;
2367:   PetscInt       dim, cMax, cEnd, c, d;

2371:   DMGetDimension(dm, &dim);
2372:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2373:   if (cMax < 0) return(0);
2374:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2375:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2376:   for (c = cMax; c < cEnd; ++c) {
2377:     const PetscInt *cone;
2378:     PetscInt       *closure = NULL;
2379:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

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

2394:       for (d = 0; d <= dim; ++d) {
2395:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2396:           DMLabelSetValue(subpointMap, point, d);
2397:           break;
2398:         }
2399:       }
2400:     }
2401:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2402:     /* Cells -- positive face is not included */
2403:     for (cl = 0; cl < 1; ++cl) {
2404:       const PetscInt *support;
2405:       PetscInt        supportSize, s;

2407:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2408:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2409:       DMPlexGetSupport(dm, cone[cl], &support);
2410:       for (s = 0; s < supportSize; ++s) {
2411:         DMLabelSetValue(subpointMap, support[s], dim);
2412:       }
2413:     }
2414:   }
2415:   PetscFree2(pStart, pEnd);
2416:   return(0);
2417: }

2419: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2420: {
2421:   MPI_Comm       comm;
2422:   PetscBool      posOrient = PETSC_FALSE;
2423:   const PetscInt debug     = 0;
2424:   PetscInt       cellDim, faceSize, f;

2428:   PetscObjectGetComm((PetscObject)dm,&comm);
2429:   DMGetDimension(dm, &cellDim);
2430:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

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

2477:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2478:     PetscSortInt(faceSizeTri, sortedIndices);
2479:     for (iFace = 0; iFace < 3; ++iFace) {
2480:       const PetscInt ii = iFace*faceSizeTri;
2481:       PetscInt       fVertex, cVertex;

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

2519:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2520:     PetscSortInt(faceSizeQuad, sortedIndices);
2521:     for (iFace = 0; iFace < 4; ++iFace) {
2522:       const PetscInt ii = iFace*faceSizeQuad;
2523:       PetscInt       fVertex, cVertex;

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

2547:           7---6
2548:          /|  /|
2549:         4---5 |
2550:         | 1-|-2
2551:         |/  |/
2552:         0---3

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

2575:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2576:     PetscSortInt(faceSizeHex, sortedIndices);
2577:     for (iFace = 0; iFace < 6; ++iFace) {
2578:       const PetscInt ii = iFace*faceSizeHex;
2579:       PetscInt       fVertex, cVertex;

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

2619:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2620:     PetscSortInt(faceSizeTet, sortedIndices);
2621:     for (iFace=0; iFace < 4; ++iFace) {
2622:       const PetscInt ii = iFace*faceSizeTet;
2623:       PetscInt       fVertex, cVertex;

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

2649:          7---6
2650:         /|  /|
2651:        4---5 |
2652:        | 3-|-2
2653:        |/  |/
2654:        0---1

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

2677:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2678:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2679:     for (iFace = 0; iFace < 6; ++iFace) {
2680:       const PetscInt ii = iFace*faceSizeQuadHex;
2681:       PetscInt       fVertex, cVertex;

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

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

2718:   Not collective

2720:   Input Parameters:
2721: + dm           - The original mesh
2722: . cell         - The cell mesh point
2723: . faceSize     - The number of vertices on the face
2724: . face         - The face vertices
2725: . numCorners   - The number of vertices on the cell
2726: . indices      - Local numbering of face vertices in cell cone
2727: - origVertices - Original face vertices

2729:   Output Parameter:
2730: + faceVertices - The face vertices properly oriented
2731: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2733:   Level: developer

2735: .seealso: DMPlexGetCone()
2736: @*/
2737: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2738: {
2739:   const PetscInt *cone = NULL;
2740:   PetscInt        coneSize, v, f, v2;
2741:   PetscInt        oppositeVertex = -1;
2742:   PetscErrorCode  ierr;

2745:   DMPlexGetConeSize(dm, cell, &coneSize);
2746:   DMPlexGetCone(dm, cell, &cone);
2747:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2748:     PetscBool found = PETSC_FALSE;

2750:     for (f = 0; f < faceSize; ++f) {
2751:       if (face[f] == cone[v]) {
2752:         found = PETSC_TRUE; break;
2753:       }
2754:     }
2755:     if (found) {
2756:       indices[v2]      = v;
2757:       origVertices[v2] = cone[v];
2758:       ++v2;
2759:     } else {
2760:       oppositeVertex = v;
2761:     }
2762:   }
2763:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2764:   return(0);
2765: }

2767: /*
2768:   DMPlexInsertFace_Internal - Puts a face into the mesh

2770:   Not collective

2772:   Input Parameters:
2773:   + dm              - The DMPlex
2774:   . numFaceVertex   - The number of vertices in the face
2775:   . faceVertices    - The vertices in the face for dm
2776:   . subfaceVertices - The vertices in the face for subdm
2777:   . numCorners      - The number of vertices in the cell
2778:   . cell            - A cell in dm containing the face
2779:   . subcell         - A cell in subdm containing the face
2780:   . firstFace       - First face in the mesh
2781:   - newFacePoint    - Next face in the mesh

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

2786:   Level: developer
2787: */
2788: 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)
2789: {
2790:   MPI_Comm        comm;
2791:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2792:   const PetscInt *faces;
2793:   PetscInt        numFaces, coneSize;
2794:   PetscErrorCode  ierr;

2797:   PetscObjectGetComm((PetscObject)dm,&comm);
2798:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2799:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2800: #if 0
2801:   /* Cannot use this because support() has not been constructed yet */
2802:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2803: #else
2804:   {
2805:     PetscInt f;

2807:     numFaces = 0;
2808:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2809:     for (f = firstFace; f < *newFacePoint; ++f) {
2810:       PetscInt dof, off, d;

2812:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2813:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2814:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2815:       for (d = 0; d < dof; ++d) {
2816:         const PetscInt p = submesh->cones[off+d];
2817:         PetscInt       v;

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

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

2868: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2869: {
2870:   MPI_Comm        comm;
2871:   DMLabel         subpointMap;
2872:   IS              subvertexIS,  subcellIS;
2873:   const PetscInt *subVertices, *subCells;
2874:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2875:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2876:   PetscInt        vStart, vEnd, c, f;
2877:   PetscErrorCode  ierr;

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

2916:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2917:     for (cl = 0; cl < closureSize*2; cl += 2) {
2918:       const PetscInt point = closure[cl];
2919:       PetscInt       subVertex;

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

2948:     DMGetCoordinateSection(dm, &coordSection);
2949:     DMGetCoordinatesLocal(dm, &coordinates);
2950:     DMGetCoordinateSection(subdm, &subCoordSection);
2951:     PetscSectionSetNumFields(subCoordSection, 1);
2952:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2953:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2954:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2955:     for (v = 0; v < numSubVertices; ++v) {
2956:       const PetscInt vertex    = subVertices[v];
2957:       const PetscInt subvertex = firstSubVertex+v;
2958:       PetscInt       dof;

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

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

3000: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3001: {
3002:   PetscInt       subPoint;

3005:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
3006:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
3007: }

3009: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3010: {
3011:   MPI_Comm         comm;
3012:   DMLabel          subpointMap;
3013:   IS              *subpointIS;
3014:   const PetscInt **subpoints;
3015:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3016:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
3017:   PetscMPIInt      rank;
3018:   PetscErrorCode   ierr;

3021:   PetscObjectGetComm((PetscObject)dm,&comm);
3022:   MPI_Comm_rank(comm, &rank);
3023:   /* Create subpointMap which marks the submesh */
3024:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3025:   DMPlexSetSubpointMap(subdm, subpointMap);
3026:   if (cellHeight) {
3027:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
3028:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
3029:   } else {
3030:     DMLabel         depth;
3031:     IS              pointIS;
3032:     const PetscInt *points;
3033:     PetscInt        numPoints=0;

3035:     DMPlexGetDepthLabel(dm, &depth);
3036:     DMLabelGetStratumIS(label, value, &pointIS);
3037:     if (pointIS) {
3038:       ISGetIndices(pointIS, &points);
3039:       ISGetLocalSize(pointIS, &numPoints);
3040:     }
3041:     for (p = 0; p < numPoints; ++p) {
3042:       PetscInt *closure = NULL;
3043:       PetscInt  closureSize, c, pdim;

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

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

3110:       if (d == dim-1) {
3111:         const PetscInt *support, *cone, *ornt;
3112:         PetscInt        supportSize, coneSize, s, subc;

3114:         DMPlexGetSupport(dm, point, &support);
3115:         DMPlexGetSupportSize(dm, point, &supportSize);
3116:         for (s = 0; s < supportSize; ++s) {
3117:           PetscBool isHybrid;

3119:           DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3120:           if (!isHybrid) continue;
3121:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3122:           if (subc >= 0) {
3123:             const PetscInt ccell = subpoints[d+1][subc];

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

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

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

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

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

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

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

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

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

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

3319:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3320:   return(0);
3321: }

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

3326:   Input Parameters:
3327: + dm           - The original mesh
3328: . vertexLabel  - The DMLabel marking points contained in the surface
3329: . value        - The label value to use
3330: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

3332:   Output Parameter:
3333: . subdm - The surface mesh

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

3337:   Level: developer

3339: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3340: @*/
3341: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3342: {
3343:   DMPlexInterpolatedFlag interpolated;
3344:   PetscInt       dim, cdim;

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

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

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

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

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

3451:     DMGetCoordinateDim(dm, &cdim);
3452:     DMGetCoordinateSection(dm, &coordSection);
3453:     DMGetCoordinatesLocal(dm, &coordinates);
3454:     DMGetCoordinateSection(subdm, &subCoordSection);
3455:     PetscSectionSetNumFields(subCoordSection, 1);
3456:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3457:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3458:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3459:     for (v = 0; v < numSubVertices; ++v) {
3460:       const PetscInt vertex    = subVertices[v];
3461:       const PetscInt subvertex = firstSubVertex+v;
3462:       PetscInt       dof;

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

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

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

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

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

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

3574: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3575: {
3576:   DMLabel        label = NULL;

3580:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3581:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3582:   return(0);
3583: }

3585: /*@C
3586:   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.

3588:   Input Parameters:
3589: + dm          - The original mesh
3590: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3591: . label       - A label name, or NULL
3592: - value  - A label value

3594:   Output Parameter:
3595: . subdm - The surface mesh

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

3599:   Level: developer

3601: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3602: @*/
3603: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3604: {
3605:   PetscInt       dim, cdim, depth;

3611:   DMGetDimension(dm, &dim);
3612:   DMPlexGetDepth(dm, &depth);
3613:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3614:   DMSetType(*subdm, DMPLEX);
3615:   DMSetDimension(*subdm, dim-1);
3616:   DMGetCoordinateDim(dm, &cdim);
3617:   DMSetCoordinateDim(*subdm, cdim);
3618:   if (depth == dim) {
3619:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3620:   } else {
3621:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3622:   }
3623:   return(0);
3624: }

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

3629:   Input Parameters:
3630: + dm        - The original mesh
3631: . cellLabel - The DMLabel marking cells contained in the new mesh
3632: - value     - The label value to use

3634:   Output Parameter:
3635: . subdm - The new mesh

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

3639:   Level: developer

3641: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3642: @*/
3643: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3644: {
3645:   PetscInt       dim;

3651:   DMGetDimension(dm, &dim);
3652:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3653:   DMSetType(*subdm, DMPLEX);
3654:   DMSetDimension(*subdm, dim);
3655:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3656:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3657:   return(0);
3658: }

3660: /*@
3661:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3663:   Input Parameter:
3664: . dm - The submesh DM

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

3669:   Level: developer

3671: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3672: @*/
3673: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3674: {
3678:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3679:   return(0);
3680: }

3682: /*@
3683:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3685:   Input Parameters:
3686: + dm - The submesh DM
3687: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3691:   Level: developer

3693: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS()
3694: @*/
3695: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3696: {
3697:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3698:   DMLabel        tmp;

3703:   tmp  = mesh->subpointMap;
3704:   mesh->subpointMap = subpointMap;
3705:   PetscObjectReference((PetscObject) mesh->subpointMap);
3706:   DMLabelDestroy(&tmp);
3707:   return(0);
3708: }

3710: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
3711: {
3712:   DMLabel        spmap;
3713:   PetscInt       depth, d;

3717:   DMPlexGetSubpointMap(dm, &spmap);
3718:   DMPlexGetDepth(dm, &depth);
3719:   if (spmap && depth >= 0) {
3720:     DM_Plex  *mesh = (DM_Plex *) dm->data;
3721:     PetscInt *points, *depths;
3722:     PetscInt  pStart, pEnd, p, off;

3724:     DMPlexGetChart(dm, &pStart, &pEnd);
3725:     if (pStart) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3726:     PetscMalloc1(pEnd, &points);
3727:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3728:     depths[0] = depth;
3729:     depths[1] = 0;
3730:     for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3731:     for (d = 0, off = 0; d <= depth; ++d) {
3732:       const PetscInt dep = depths[d];
3733:       PetscInt       depStart, depEnd, n;

3735:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3736:       DMLabelGetStratumSize(spmap, dep, &n);
3737:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3738:         if (n != depEnd-depStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
3739:       } else {
3740:         if (!n) {
3741:           if (d == 0) {
3742:             /* Missing cells */
3743:             for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3744:           } else {
3745:             /* Missing faces */
3746:             for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3747:           }
3748:         }
3749:       }
3750:       if (n) {
3751:         IS              is;
3752:         const PetscInt *opoints;

3754:         DMLabelGetStratumIS(spmap, dep, &is);
3755:         ISGetIndices(is, &opoints);
3756:         for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3757:         ISRestoreIndices(is, &opoints);
3758:         ISDestroy(&is);
3759:       }
3760:     }
3761:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3762:     if (off != pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3763:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3764:     PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState);
3765:   }
3766:   return(0);
3767: }

3769: /*@
3770:   DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data

3772:   Input Parameter:
3773: . dm - The submesh DM

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

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

3780:   Level: developer

3782: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3783: @*/
3784: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
3785: {
3786:   DM_Plex         *mesh = (DM_Plex *) dm->data;
3787:   DMLabel          spmap;
3788:   PetscObjectState state;
3789:   PetscErrorCode   ierr;

3794:   DMPlexGetSubpointMap(dm, &spmap);
3795:   PetscObjectStateGet((PetscObject) spmap, &state);
3796:   if (state != mesh->subpointState || !mesh->subpointIS) {DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);}
3797:   *subpointIS = mesh->subpointIS;
3798:   return(0);
3799: }

3801: /*@
3802:   DMGetEnclosureRelation - Get the relationship between dmA and dmB

3804:   Input Parameters:
3805: + dmA - The first DM
3806: - dmB - The second DM

3808:   Output Parameter:
3809: . rel - The relation of dmA to dmB

3811:   Level: intermediate

3813: .seealso: DMGetEnclosurePoint()
3814: @*/
3815: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3816: {
3817:   DM             plexA, plexB, sdm;
3818:   DMLabel        spmap;
3819:   PetscInt       pStartA, pEndA, pStartB, pEndB, NpA, NpB;

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

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

3860:   Input Parameters:
3861: + dmA   - The first DM
3862: . dmB   - The second DM
3863: . etype - The type of enclosure relation that dmA has to dmB
3864: - pB    - A point of dmB

3866:   Output Parameter:
3867: . pA    - The corresponding point of dmA

3869:   Level: intermediate

3871: .seealso: DMGetEnclosureRelation()
3872: @*/
3873: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3874: {
3875:   DM              sdm;
3876:   IS              subpointIS;
3877:   const PetscInt *subpoints;
3878:   PetscInt        numSubpoints;
3879:   PetscErrorCode  ierr;

3882:   /* TODO Cache the IS, making it look like an index */
3883:   switch (etype) {
3884:     case DM_ENC_SUPERMESH:
3885:     sdm  = dmB;
3886:     DMPlexGetSubpointIS(sdm, &subpointIS);
3887:     ISGetIndices(subpointIS, &subpoints);
3888:     *pA  = subpoints[pB];
3889:     ISRestoreIndices(subpointIS, &subpoints);
3890:     break;
3891:     case DM_ENC_SUBMESH:
3892:     sdm  = dmA;
3893:     DMPlexGetSubpointIS(sdm, &subpointIS);
3894:     ISGetLocalSize(subpointIS, &numSubpoints);
3895:     ISGetIndices(subpointIS, &subpoints);
3896:     PetscFindInt(pB, numSubpoints, subpoints, pA);
3897:     if (*pA < 0) {
3898:       DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3899:       DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3900:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3901:     }
3902:     ISRestoreIndices(subpointIS, &subpoints);
3903:     break;
3904:     case DM_ENC_EQUALITY:
3905:     case DM_ENC_NONE:
3906:     *pA = pB;break;
3907:     case DM_ENC_UNKNOWN:
3908:     {
3909:       DMEnclosureType enc;

3911:       DMGetEnclosureRelation(dmA, dmB, &enc);
3912:       DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3913:     }
3914:     break;
3915:     default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3916:   }
3917:   return(0);
3918: }