Actual source code: plexsubmesh.c

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

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

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

 15:     DMPlexGetSupportSize(dm, f, &supportSize);
 16:     if (supportSize == 1) {
 17:       if (val < 0) {
 18:         PetscInt *closure = NULL;
 19:         PetscInt  clSize, cl, cval;

 21:         DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 22:         for (cl = 0; cl < clSize*2; cl += 2) {
 23:           DMLabelGetValue(label, closure[cl], &cval);
 24:           if (cval < 0) continue;
 25:           DMLabelSetValue(label, f, cval);
 26:           break;
 27:         }
 28:         if (cl == clSize*2) {DMLabelSetValue(label, f, 1);}
 29:         DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 30:       } else {
 31:         DMLabelSetValue(label, f, val);
 32:       }
 33:     }
 34:   }
 35:   return(0);
 36: }

 38: /*@
 39:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 41:   Not Collective

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

 47:   Output Parameter:
 48: . label - The DMLabel marking boundary faces with the given value

 50:   Level: developer

 52: .seealso: DMLabelCreate(), DMCreateLabel()
 53: @*/
 54: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
 55: {
 56:   DMPlexInterpolatedFlag  flg;
 57:   PetscErrorCode          ierr;

 61:   DMPlexIsInterpolated(dm, &flg);
 62:   if (flg != DMPLEX_INTERPOLATED_FULL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not fully interpolated on this rank");
 63:   DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
 64:   return(0);
 65: }

 67: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
 68: {
 69:   IS              valueIS;
 70:   PetscSF         sfPoint;
 71:   const PetscInt *values;
 72:   PetscInt        numValues, v, cStart, cEnd, nroots;
 73:   PetscErrorCode  ierr;

 76:   DMLabelGetNumValues(label, &numValues);
 77:   DMLabelGetValueIS(label, &valueIS);
 78:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 79:   ISGetIndices(valueIS, &values);
 80:   for (v = 0; v < numValues; ++v) {
 81:     IS              pointIS;
 82:     const PetscInt *points;
 83:     PetscInt        numPoints, p;

 85:     DMLabelGetStratumSize(label, values[v], &numPoints);
 86:     DMLabelGetStratumIS(label, values[v], &pointIS);
 87:     ISGetIndices(pointIS, &points);
 88:     for (p = 0; p < numPoints; ++p) {
 89:       PetscInt  q = points[p];
 90:       PetscInt *closure = NULL;
 91:       PetscInt  closureSize, c;

 93:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
 94:         continue;
 95:       }
 96:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
 97:       for (c = 0; c < closureSize*2; c += 2) {
 98:         DMLabelSetValue(label, closure[c], values[v]);
 99:       }
100:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
101:     }
102:     ISRestoreIndices(pointIS, &points);
103:     ISDestroy(&pointIS);
104:   }
105:   ISRestoreIndices(valueIS, &values);
106:   ISDestroy(&valueIS);
107:   DMGetPointSF(dm, &sfPoint);
108:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
109:   if (nroots >= 0) {
110:     DMLabel         lblRoots, lblLeaves;
111:     IS              valueIS, pointIS;
112:     const PetscInt *values;
113:     PetscInt        numValues, v;
114:     PetscErrorCode  ierr;

116:     DMGetPointSF(dm, &sfPoint);
117:     /* Pull point contributions from remote leaves into local roots */
118:     DMLabelGather(label, sfPoint, &lblLeaves);
119:     DMLabelGetValueIS(lblLeaves, &valueIS);
120:     ISGetLocalSize(valueIS, &numValues);
121:     ISGetIndices(valueIS, &values);
122:     for (v = 0; v < numValues; ++v) {
123:       const PetscInt value = values[v];

125:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
126:       DMLabelInsertIS(label, pointIS, value);
127:       ISDestroy(&pointIS);
128:     }
129:     ISRestoreIndices(valueIS, &values);
130:     ISDestroy(&valueIS);
131:     DMLabelDestroy(&lblLeaves);
132:     /* Push point contributions from roots into remote leaves */
133:     DMLabelDistribute(label, sfPoint, &lblRoots);
134:     DMLabelGetValueIS(lblRoots, &valueIS);
135:     ISGetLocalSize(valueIS, &numValues);
136:     ISGetIndices(valueIS, &values);
137:     for (v = 0; v < numValues; ++v) {
138:       const PetscInt value = values[v];

140:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
141:       DMLabelInsertIS(label, pointIS, value);
142:       ISDestroy(&pointIS);
143:     }
144:     ISRestoreIndices(valueIS, &values);
145:     ISDestroy(&valueIS);
146:     DMLabelDestroy(&lblRoots);
147:   }
148:   return(0);
149: }

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

154:   Input Parameters:
155: + dm - The DM
156: - label - A DMLabel marking the surface points

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

161:   Level: developer

163: .seealso: DMPlexLabelCohesiveComplete()
164: @*/
165: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
166: {

170:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
171:   return(0);
172: }

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

177:   Input Parameters:
178: + dm - The DM
179: - label - A DMLabel marking the surface points

181:   Output Parameter:
182: . label - A DMLabel incorporating cells

184:   Level: developer

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

188: .seealso: DMPlexLabelAddFaceCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
189: @*/
190: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
191: {
192:   IS              valueIS;
193:   const PetscInt *values;
194:   PetscInt        numValues, v, cStart, cEnd;
195:   PetscErrorCode  ierr;

198:   DMPlexGetInteriorCellStratum(dm, &cStart, &cEnd);
199:   DMLabelGetNumValues(label, &numValues);
200:   DMLabelGetValueIS(label, &valueIS);
201:   ISGetIndices(valueIS, &values);
202:   for (v = 0; v < numValues; ++v) {
203:     IS              pointIS;
204:     const PetscInt *points;
205:     PetscInt        numPoints, p;

207:     DMLabelGetStratumSize(label, values[v], &numPoints);
208:     DMLabelGetStratumIS(label, values[v], &pointIS);
209:     ISGetIndices(pointIS, &points);
210:     for (p = 0; p < numPoints; ++p) {
211:       PetscInt *closure = NULL;
212:       PetscInt  closureSize, cl;

214:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
215:       for (cl = closureSize-1; cl > 0; --cl) {
216:         const PetscInt cell = closure[cl*2];
217:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
218:       }
219:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
220:     }
221:     ISRestoreIndices(pointIS, &points);
222:     ISDestroy(&pointIS);
223:   }
224:   ISRestoreIndices(valueIS, &values);
225:   ISDestroy(&valueIS);
226:   return(0);
227: }

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

232:   Input Parameters:
233: + dm - The DM
234: - label - A DMLabel marking the surface points

236:   Output Parameter:
237: . label - A DMLabel incorporating cells

239:   Level: developer

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

243: .seealso: DMPlexLabelAddCells(), DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
244: @*/
245: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
246: {
247:   IS              valueIS;
248:   const PetscInt *values;
249:   PetscInt        numValues, v, cStart, cEnd, fStart, fEnd;
250:   PetscErrorCode  ierr;

253:   DMPlexGetInteriorCellStratum(dm, &cStart, &cEnd);
254:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
255:   DMLabelGetNumValues(label, &numValues);
256:   DMLabelGetValueIS(label, &valueIS);
257:   ISGetIndices(valueIS, &values);
258:   for (v = 0; v < numValues; ++v) {
259:     IS              pointIS;
260:     const PetscInt *points;
261:     PetscInt        numPoints, p;

263:     DMLabelGetStratumSize(label, values[v], &numPoints);
264:     DMLabelGetStratumIS(label, values[v], &pointIS);
265:     ISGetIndices(pointIS, &points);
266:     for (p = 0; p < numPoints; ++p) {
267:       const PetscInt face = points[p];
268:       PetscInt      *closure = NULL;
269:       PetscInt       closureSize, cl;

271:       if ((face < fStart) || (face >= fEnd)) continue;
272:       DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
273:       for (cl = closureSize-1; cl > 0; --cl) {
274:         const PetscInt cell = closure[cl*2];
275:         if ((cell >= cStart) && (cell < cEnd)) {DMLabelSetValue(label, cell, values[v]); break;}
276:       }
277:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
278:     }
279:     ISRestoreIndices(pointIS, &points);
280:     ISDestroy(&pointIS);
281:   }
282:   ISRestoreIndices(valueIS, &values);
283:   ISDestroy(&valueIS);
284:   return(0);
285: }

287: /*@
288:   DMPlexLabelClearCells - Remove cells from a label

290:   Input Parameters:
291: + dm - The DM
292: - label - A DMLabel marking surface points and their adjacent cells

294:   Output Parameter:
295: . label - A DMLabel without cells

297:   Level: developer

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

301: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
302: @*/
303: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
304: {
305:   IS              valueIS;
306:   const PetscInt *values;
307:   PetscInt        numValues, v, cStart, cEnd;
308:   PetscErrorCode  ierr;

311:   DMPlexGetInteriorCellStratum(dm, &cStart, &cEnd);
312:   DMLabelGetNumValues(label, &numValues);
313:   DMLabelGetValueIS(label, &valueIS);
314:   ISGetIndices(valueIS, &values);
315:   for (v = 0; v < numValues; ++v) {
316:     IS              pointIS;
317:     const PetscInt *points;
318:     PetscInt        numPoints, p;

320:     DMLabelGetStratumSize(label, values[v], &numPoints);
321:     DMLabelGetStratumIS(label, values[v], &pointIS);
322:     ISGetIndices(pointIS, &points);
323:     for (p = 0; p < numPoints; ++p) {
324:       PetscInt point = points[p];

326:       if (point >= cStart && point < cEnd) {
327:         DMLabelClearValue(label,point,values[v]);
328:       }
329:     }
330:     ISRestoreIndices(pointIS, &points);
331:     ISDestroy(&pointIS);
332:   }
333:   ISRestoreIndices(valueIS, &values);
334:   ISDestroy(&valueIS);
335:   return(0);
336: }

338: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
339:  * index (skipping first, which is (0,0)) */
340: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
341: {
342:   PetscInt d, off = 0;

345:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
346:   for (d = 0; d < depth; d++) {
347:     PetscInt firstd = d;
348:     PetscInt firstStart = depthShift[2*d];
349:     PetscInt e;

351:     for (e = d+1; e <= depth; e++) {
352:       if (depthShift[2*e] < firstStart) {
353:         firstd = e;
354:         firstStart = depthShift[2*d];
355:       }
356:     }
357:     if (firstd != d) {
358:       PetscInt swap[2];

360:       e = firstd;
361:       swap[0] = depthShift[2*d];
362:       swap[1] = depthShift[2*d+1];
363:       depthShift[2*d]   = depthShift[2*e];
364:       depthShift[2*d+1] = depthShift[2*e+1];
365:       depthShift[2*e]   = swap[0];
366:       depthShift[2*e+1] = swap[1];
367:     }
368:   }
369:   /* convert (oldstart, added) to (oldstart, newstart) */
370:   for (d = 0; d <= depth; d++) {
371:     off += depthShift[2*d+1];
372:     depthShift[2*d+1] = depthShift[2*d] + off;
373:   }
374:   return(0);
375: }

377: /* depthShift is a list of (old, new) pairs */
378: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
379: {
380:   PetscInt d;
381:   PetscInt newOff = 0;

383:   for (d = 0; d <= depth; d++) {
384:     if (p < depthShift[2*d]) return p + newOff;
385:     else newOff = depthShift[2*d+1] - depthShift[2*d];
386:   }
387:   return p + newOff;
388: }

390: /* depthShift is a list of (old, new) pairs */
391: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
392: {
393:   PetscInt d;
394:   PetscInt newOff = 0;

396:   for (d = 0; d <= depth; d++) {
397:     if (p < depthShift[2*d+1]) return p + newOff;
398:     else newOff = depthShift[2*d] - depthShift[2*d+1];
399:   }
400:   return p + newOff;
401: }

403: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
404: {
405:   PetscInt       depth = 0, d, pStart, pEnd, p;
406:   DMLabel        depthLabel;

410:   DMPlexGetDepth(dm, &depth);
411:   if (depth < 0) return(0);
412:   /* Step 1: Expand chart */
413:   DMPlexGetChart(dm, &pStart, &pEnd);
414:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
415:   DMPlexSetChart(dmNew, pStart, pEnd);
416:   DMCreateLabel(dmNew,"depth");
417:   DMPlexGetDepthLabel(dmNew,&depthLabel);
418:   /* Step 2: Set cone and support sizes */
419:   for (d = 0; d <= depth; ++d) {
420:     PetscInt pStartNew, pEndNew;
421:     IS pIS;

423:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
424:     pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
425:     pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
426:     ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
427:     DMLabelSetStratumIS(depthLabel, d, pIS);
428:     ISDestroy(&pIS);
429:     for (p = pStart; p < pEnd; ++p) {
430:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
431:       PetscInt size;

433:       DMPlexGetConeSize(dm, p, &size);
434:       DMPlexSetConeSize(dmNew, newp, size);
435:       DMPlexGetSupportSize(dm, p, &size);
436:       DMPlexSetSupportSize(dmNew, newp, size);
437:     }
438:   }
439:   return(0);
440: }

442: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
443: {
444:   PetscInt      *newpoints;
445:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

449:   DMPlexGetDepth(dm, &depth);
450:   if (depth < 0) return(0);
451:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
452:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
453:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
454:   /* Step 5: Set cones and supports */
455:   DMPlexGetChart(dm, &pStart, &pEnd);
456:   for (p = pStart; p < pEnd; ++p) {
457:     const PetscInt *points = NULL, *orientations = NULL;
458:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

460:     DMPlexGetConeSize(dm, p, &size);
461:     DMPlexGetCone(dm, p, &points);
462:     DMPlexGetConeOrientation(dm, p, &orientations);
463:     for (i = 0; i < size; ++i) {
464:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
465:     }
466:     DMPlexSetCone(dmNew, newp, newpoints);
467:     DMPlexSetConeOrientation(dmNew, newp, orientations);
468:     DMPlexGetSupportSize(dm, p, &size);
469:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
470:     DMPlexGetSupport(dm, p, &points);
471:     for (i = 0; i < size; ++i) {
472:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
473:     }
474:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
475:     DMPlexSetSupport(dmNew, newp, newpoints);
476:   }
477:   PetscFree(newpoints);
478:   return(0);
479: }

481: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
482: {
483:   PetscSection   coordSection, newCoordSection;
484:   Vec            coordinates, newCoordinates;
485:   PetscScalar   *coords, *newCoords;
486:   PetscInt       coordSize, sStart, sEnd;
487:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
488:   PetscBool      hasCells;

492:   DMGetCoordinateDim(dm, &dim);
493:   DMSetCoordinateDim(dmNew, dim);
494:   DMPlexGetDepth(dm, &depth);
495:   /* Step 8: Convert coordinates */
496:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
497:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
498:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
499:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
500:   DMGetCoordinateSection(dm, &coordSection);
501:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
502:   PetscSectionSetNumFields(newCoordSection, 1);
503:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
504:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
505:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
506:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
507:   if (hasCells) {
508:     for (c = cStart; c < cEnd; ++c) {
509:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

511:       PetscSectionGetDof(coordSection, c, &dof);
512:       PetscSectionSetDof(newCoordSection, cNew, dof);
513:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
514:     }
515:   }
516:   for (v = vStartNew; v < vEndNew; ++v) {
517:     PetscSectionSetDof(newCoordSection, v, dim);
518:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
519:   }
520:   PetscSectionSetUp(newCoordSection);
521:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
522:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
523:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
524:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
525:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
526:   VecSetBlockSize(newCoordinates, dim);
527:   VecSetType(newCoordinates,VECSTANDARD);
528:   DMSetCoordinatesLocal(dmNew, newCoordinates);
529:   DMGetCoordinatesLocal(dm, &coordinates);
530:   VecGetArray(coordinates, &coords);
531:   VecGetArray(newCoordinates, &newCoords);
532:   if (hasCells) {
533:     for (c = cStart; c < cEnd; ++c) {
534:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

536:       PetscSectionGetDof(coordSection, c, &dof);
537:       PetscSectionGetOffset(coordSection, c, &off);
538:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
539:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
540:     }
541:   }
542:   for (v = vStart; v < vEnd; ++v) {
543:     PetscInt dof, off, noff, d;

545:     PetscSectionGetDof(coordSection, v, &dof);
546:     PetscSectionGetOffset(coordSection, v, &off);
547:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
548:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
549:   }
550:   VecRestoreArray(coordinates, &coords);
551:   VecRestoreArray(newCoordinates, &newCoords);
552:   VecDestroy(&newCoordinates);
553:   PetscSectionDestroy(&newCoordSection);
554:   return(0);
555: }

557: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
558: {
559:   PetscInt           depth = 0;
560:   PetscSF            sfPoint, sfPointNew;
561:   const PetscSFNode *remotePoints;
562:   PetscSFNode       *gremotePoints;
563:   const PetscInt    *localPoints;
564:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
565:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
566:   PetscErrorCode     ierr;

569:   DMPlexGetDepth(dm, &depth);
570:   /* Step 9: Convert pointSF */
571:   DMGetPointSF(dm, &sfPoint);
572:   DMGetPointSF(dmNew, &sfPointNew);
573:   DMPlexGetChart(dm, &pStart, &pEnd);
574:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
575:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
576:   if (numRoots >= 0) {
577:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
578:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
579:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
580:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
581:     PetscMalloc1(numLeaves,    &glocalPoints);
582:     PetscMalloc1(numLeaves, &gremotePoints);
583:     for (l = 0; l < numLeaves; ++l) {
584:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
585:       gremotePoints[l].rank  = remotePoints[l].rank;
586:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
587:     }
588:     PetscFree2(newLocation,newRemoteLocation);
589:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
590:   }
591:   return(0);
592: }

594: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
595: {
596:   PetscSF            sfPoint;
597:   DMLabel            vtkLabel, ghostLabel;
598:   const PetscSFNode *leafRemote;
599:   const PetscInt    *leafLocal;
600:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
601:   PetscMPIInt        rank;
602:   PetscErrorCode     ierr;

605:   DMPlexGetDepth(dm, &depth);
606:   /* Step 10: Convert labels */
607:   DMGetNumLabels(dm, &numLabels);
608:   for (l = 0; l < numLabels; ++l) {
609:     DMLabel         label, newlabel;
610:     const char     *lname;
611:     PetscBool       isDepth, isDim;
612:     IS              valueIS;
613:     const PetscInt *values;
614:     PetscInt        numValues, val;

616:     DMGetLabelName(dm, l, &lname);
617:     PetscStrcmp(lname, "depth", &isDepth);
618:     if (isDepth) continue;
619:     PetscStrcmp(lname, "dim", &isDim);
620:     if (isDim) continue;
621:     DMCreateLabel(dmNew, lname);
622:     DMGetLabel(dm, lname, &label);
623:     DMGetLabel(dmNew, lname, &newlabel);
624:     DMLabelGetDefaultValue(label,&val);
625:     DMLabelSetDefaultValue(newlabel,val);
626:     DMLabelGetValueIS(label, &valueIS);
627:     ISGetLocalSize(valueIS, &numValues);
628:     ISGetIndices(valueIS, &values);
629:     for (val = 0; val < numValues; ++val) {
630:       IS              pointIS;
631:       const PetscInt *points;
632:       PetscInt        numPoints, p;

634:       DMLabelGetStratumIS(label, values[val], &pointIS);
635:       ISGetLocalSize(pointIS, &numPoints);
636:       ISGetIndices(pointIS, &points);
637:       for (p = 0; p < numPoints; ++p) {
638:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

640:         DMLabelSetValue(newlabel, newpoint, values[val]);
641:       }
642:       ISRestoreIndices(pointIS, &points);
643:       ISDestroy(&pointIS);
644:     }
645:     ISRestoreIndices(valueIS, &values);
646:     ISDestroy(&valueIS);
647:   }
648:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
649:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
650:   DMGetPointSF(dm, &sfPoint);
651:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
652:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
653:   DMCreateLabel(dmNew, "vtk");
654:   DMCreateLabel(dmNew, "ghost");
655:   DMGetLabel(dmNew, "vtk", &vtkLabel);
656:   DMGetLabel(dmNew, "ghost", &ghostLabel);
657:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
658:     for (; c < leafLocal[l] && c < cEnd; ++c) {
659:       DMLabelSetValue(vtkLabel, c, 1);
660:     }
661:     if (leafLocal[l] >= cEnd) break;
662:     if (leafRemote[l].rank == rank) {
663:       DMLabelSetValue(vtkLabel, c, 1);
664:     } else {
665:       DMLabelSetValue(ghostLabel, c, 2);
666:     }
667:   }
668:   for (; c < cEnd; ++c) {
669:     DMLabelSetValue(vtkLabel, c, 1);
670:   }
671:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
672:   for (f = fStart; f < fEnd; ++f) {
673:     PetscInt numCells;

675:     DMPlexGetSupportSize(dmNew, f, &numCells);
676:     if (numCells < 2) {
677:       DMLabelSetValue(ghostLabel, f, 1);
678:     } else {
679:       const PetscInt *cells = NULL;
680:       PetscInt        vA, vB;

682:       DMPlexGetSupport(dmNew, f, &cells);
683:       DMLabelGetValue(vtkLabel, cells[0], &vA);
684:       DMLabelGetValue(vtkLabel, cells[1], &vB);
685:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
686:     }
687:   }
688:   return(0);
689: }

691: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
692: {
693:   DM             refTree;
694:   PetscSection   pSec;
695:   PetscInt       *parents, *childIDs;

699:   DMPlexGetReferenceTree(dm,&refTree);
700:   DMPlexSetReferenceTree(dmNew,refTree);
701:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
702:   if (pSec) {
703:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
704:     PetscInt *childIDsShifted;
705:     PetscSection pSecShifted;

707:     PetscSectionGetChart(pSec,&pStart,&pEnd);
708:     DMPlexGetDepth(dm,&depth);
709:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
710:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
711:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
712:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
713:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
714:     for (p = pStartShifted; p < pEndShifted; p++) {
715:       /* start off assuming no children */
716:       PetscSectionSetDof(pSecShifted,p,0);
717:     }
718:     for (p = pStart; p < pEnd; p++) {
719:       PetscInt dof;
720:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

722:       PetscSectionGetDof(pSec,p,&dof);
723:       PetscSectionSetDof(pSecShifted,pNew,dof);
724:     }
725:     PetscSectionSetUp(pSecShifted);
726:     for (p = pStart; p < pEnd; p++) {
727:       PetscInt dof;
728:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

730:       PetscSectionGetDof(pSec,p,&dof);
731:       if (dof) {
732:         PetscInt off, offNew;

734:         PetscSectionGetOffset(pSec,p,&off);
735:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
736:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
737:         childIDsShifted[offNew] = childIDs[off];
738:       }
739:     }
740:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
741:     PetscFree2(parentsShifted,childIDsShifted);
742:     PetscSectionDestroy(&pSecShifted);
743:   }
744:   return(0);
745: }

747: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
748: {
749:   PetscSF               sf;
750:   IS                    valueIS;
751:   const PetscInt       *values, *leaves;
752:   PetscInt             *depthShift;
753:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
754:   PetscBool             isper;
755:   const PetscReal      *maxCell, *L;
756:   const DMBoundaryType *bd;
757:   PetscErrorCode        ierr;

760:   DMGetPointSF(dm, &sf);
761:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
762:   nleaves = PetscMax(0, nleaves);
763:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
764:   /* Count ghost cells */
765:   DMLabelGetValueIS(label, &valueIS);
766:   ISGetLocalSize(valueIS, &numFS);
767:   ISGetIndices(valueIS, &values);
768:   Ng   = 0;
769:   for (fs = 0; fs < numFS; ++fs) {
770:     IS              faceIS;
771:     const PetscInt *faces;
772:     PetscInt        numFaces, f, numBdFaces = 0;

774:     DMLabelGetStratumIS(label, values[fs], &faceIS);
775:     ISGetLocalSize(faceIS, &numFaces);
776:     ISGetIndices(faceIS, &faces);
777:     for (f = 0; f < numFaces; ++f) {
778:       PetscInt numChildren;

780:       PetscFindInt(faces[f], nleaves, leaves, &loc);
781:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
782:       /* non-local and ancestors points don't get to register ghosts */
783:       if (loc >= 0 || numChildren) continue;
784:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
785:     }
786:     Ng += numBdFaces;
787:     ISRestoreIndices(faceIS, &faces);
788:     ISDestroy(&faceIS);
789:   }
790:   DMPlexGetDepth(dm, &depth);
791:   PetscMalloc1(2*(depth+1), &depthShift);
792:   for (d = 0; d <= depth; d++) {
793:     PetscInt dEnd;

795:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
796:     depthShift[2*d]   = dEnd;
797:     depthShift[2*d+1] = 0;
798:   }
799:   if (depth >= 0) depthShift[2*depth+1] = Ng;
800:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
801:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
802:   /* Step 3: Set cone/support sizes for new points */
803:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
804:   DMPlexSetGhostCellStratum(gdm, cEnd, PETSC_DETERMINE);
805:   for (c = cEnd; c < cEnd + Ng; ++c) {
806:     DMPlexSetConeSize(gdm, c, 1);
807:   }
808:   for (fs = 0; fs < numFS; ++fs) {
809:     IS              faceIS;
810:     const PetscInt *faces;
811:     PetscInt        numFaces, f;

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

819:       PetscFindInt(faces[f], nleaves, leaves, &loc);
820:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
821:       if (loc >= 0 || numChildren) continue;
822:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
823:       DMPlexGetSupportSize(dm, faces[f], &size);
824:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
825:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
826:     }
827:     ISRestoreIndices(faceIS, &faces);
828:     ISDestroy(&faceIS);
829:   }
830:   /* Step 4: Setup ghosted DM */
831:   DMSetUp(gdm);
832:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
833:   /* Step 6: Set cones and supports for new points */
834:   ghostCell = cEnd;
835:   for (fs = 0; fs < numFS; ++fs) {
836:     IS              faceIS;
837:     const PetscInt *faces;
838:     PetscInt        numFaces, f;

840:     DMLabelGetStratumIS(label, values[fs], &faceIS);
841:     ISGetLocalSize(faceIS, &numFaces);
842:     ISGetIndices(faceIS, &faces);
843:     for (f = 0; f < numFaces; ++f) {
844:       PetscInt newFace = faces[f] + Ng, numChildren;

846:       PetscFindInt(faces[f], nleaves, leaves, &loc);
847:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
848:       if (loc >= 0 || numChildren) continue;
849:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
850:       DMPlexSetCone(gdm, ghostCell, &newFace);
851:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
852:       ++ghostCell;
853:     }
854:     ISRestoreIndices(faceIS, &faces);
855:     ISDestroy(&faceIS);
856:   }
857:   ISRestoreIndices(valueIS, &values);
858:   ISDestroy(&valueIS);
859:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
860:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
861:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
862:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
863:   PetscFree(depthShift);
864:   /* Step 7: Periodicity */
865:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
866:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
867:   if (numGhostCells) *numGhostCells = Ng;
868:   return(0);
869: }

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

874:   Collective on dm

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

880:   Output Parameters:
881: + numGhostCells - The number of ghost cells added to the DM
882: - dmGhosted - The new DM

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

886:   Level: developer

888: .seealso: DMCreate()
889: @*/
890: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
891: {
892:   DM             gdm;
893:   DMLabel        label;
894:   const char    *name = labelName ? labelName : "Face Sets";
895:   PetscInt       dim, Ng = 0, cMax, fMax, eMax, vMax;
896:   PetscBool      useCone, useClosure;

903:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
904:   DMSetType(gdm, DMPLEX);
905:   DMGetDimension(dm, &dim);
906:   DMSetDimension(gdm, dim);
907:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
908:   DMSetBasicAdjacency(gdm, useCone,  useClosure);
909:   DMGetLabel(dm, name, &label);
910:   if (!label) {
911:     /* Get label for boundary faces */
912:     DMCreateLabel(dm, name);
913:     DMGetLabel(dm, name, &label);
914:     DMPlexMarkBoundaryFaces(dm, 1, label);
915:   }
916:   DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
917:   DMCopyBoundary(dm, gdm);
918:   DMCopyDisc(dm, gdm);
919:   /* Copy the hybrid bounds */
920:   DMPlexGetHybridBounds(dm, &cMax, &fMax, &eMax, &vMax);
921:   if (fMax >= 0) fMax += Ng;
922:   if (eMax >= 0) eMax += Ng;
923:   if (vMax >= 0) vMax += Ng;
924:   DMPlexSetHybridBounds(gdm, cMax, fMax, eMax, vMax);
925:   gdm->setfromoptionscalled = dm->setfromoptionscalled;
926:   if (numGhostCells) *numGhostCells = Ng;
927:   *dmGhosted = gdm;
928:   return(0);
929: }

931: /*
932:   We are adding three kinds of points here:
933:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
934:     Non-replicated: Points which exist on the fault, but are not replicated
935:     Hybrid:         Entirely new points, such as cohesive cells

937:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
938:   each depth so that the new split/hybrid points can be inserted as a block.
939: */
940: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
941: {
942:   MPI_Comm         comm;
943:   IS               valueIS;
944:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
945:   const PetscInt  *values;          /* List of depths for which we have replicated points */
946:   IS              *splitIS;
947:   IS              *unsplitIS;
948:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
949:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
950:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
951:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
952:   const PetscInt **splitPoints;        /* Replicated points for each depth */
953:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
954:   PetscSection     coordSection;
955:   Vec              coordinates;
956:   PetscScalar     *coords;
957:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
958:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
959:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
960:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
961:   PetscInt        *coneNew, *coneONew, *supportNew;
962:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
963:   PetscErrorCode   ierr;

966:   PetscObjectGetComm((PetscObject)dm,&comm);
967:   DMGetDimension(dm, &dim);
968:   DMPlexGetDepth(dm, &depth);
969:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
970:   /* Count split points and add cohesive cells */
971:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
972:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
973:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
974:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
975:   for (d = 0; d <= depth; ++d) {
976:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
977:     depthEnd[d]           = pMaxNew[d];
978:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
979:     numSplitPoints[d]     = 0;
980:     numUnsplitPoints[d]   = 0;
981:     numHybridPoints[d]    = 0;
982:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
983:     splitPoints[d]        = NULL;
984:     unsplitPoints[d]      = NULL;
985:     splitIS[d]            = NULL;
986:     unsplitIS[d]          = NULL;
987:     /* we are shifting the existing hybrid points with the stratum behind them, so
988:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
989:     depthShift[2*d]       = depthMax[d];
990:     depthShift[2*d+1]     = 0;
991:   }
992:   if (label) {
993:     DMLabelGetValueIS(label, &valueIS);
994:     ISGetLocalSize(valueIS, &numSP);
995:     ISGetIndices(valueIS, &values);
996:   }
997:   for (sp = 0; sp < numSP; ++sp) {
998:     const PetscInt dep = values[sp];

1000:     if ((dep < 0) || (dep > depth)) continue;
1001:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1002:     if (splitIS[dep]) {
1003:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1004:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
1005:     }
1006:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
1007:     if (unsplitIS[dep]) {
1008:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1009:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1010:     }
1011:   }
1012:   /* Calculate number of hybrid points */
1013:   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   */
1014:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
1015:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
1016:   /* the end of the points in this stratum that come before the new points:
1017:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1018:    * added points */
1019:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1020:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1021:   /* Step 3: Set cone/support sizes for new points */
1022:   for (dep = 0; dep <= depth; ++dep) {
1023:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1024:       const PetscInt  oldp   = splitPoints[dep][p];
1025:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1026:       const PetscInt  splitp = p    + pMaxNew[dep];
1027:       const PetscInt *support;
1028:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

1030:       DMPlexGetConeSize(dm, oldp, &coneSize);
1031:       DMPlexSetConeSize(sdm, splitp, coneSize);
1032:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1033:       DMPlexSetSupportSize(sdm, splitp, supportSize);
1034:       if (dep == depth-1) {
1035:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1042:         DMPlexGetSupport(dm, oldp, &support);
1043:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1044:           PetscInt val;

1046:           DMLabelGetValue(label, support[e], &val);
1047:           if (val == 1) ++qf;
1048:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
1049:           if ((val == 1) || (val == -(shift + 1))) ++qp;
1050:         }
1051:         /* Split old vertex: Edges into original vertex and new cohesive edge */
1052:         DMPlexSetSupportSize(sdm, newp, qn+1);
1053:         /* Split new vertex: Edges into split vertex and new cohesive edge */
1054:         DMPlexSetSupportSize(sdm, splitp, qp+1);
1055:         /* Add hybrid edge */
1056:         DMPlexSetConeSize(sdm, hybedge, 2);
1057:         DMPlexSetSupportSize(sdm, hybedge, qf);
1058:       } else if (dep == dim-2) {
1059:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1061:         DMPlexGetSupport(dm, oldp, &support);
1062:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1063:           PetscInt val;

1065:           DMLabelGetValue(label, support[e], &val);
1066:           if (val == dim-1) ++qf;
1067:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
1068:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1069:         }
1070:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1071:         DMPlexSetSupportSize(sdm, newp, qn+1);
1072:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1073:         DMPlexSetSupportSize(sdm, splitp, qp+1);
1074:         /* Add hybrid face */
1075:         DMPlexSetConeSize(sdm, hybface, 4);
1076:         DMPlexSetSupportSize(sdm, hybface, qf);
1077:       }
1078:     }
1079:   }
1080:   for (dep = 0; dep <= depth; ++dep) {
1081:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1082:       const PetscInt  oldp   = unsplitPoints[dep][p];
1083:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1084:       const PetscInt *support;
1085:       PetscInt        coneSize, supportSize, qf, e, s;

1087:       DMPlexGetConeSize(dm, oldp, &coneSize);
1088:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1089:       DMPlexGetSupport(dm, oldp, &support);
1090:       if (dep == 0) {
1091:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1093:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1094:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1095:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1096:           if (e >= 0) ++qf;
1097:         }
1098:         DMPlexSetSupportSize(sdm, newp, qf+2);
1099:         /* Add hybrid edge */
1100:         DMPlexSetConeSize(sdm, hybedge, 2);
1101:         for (e = 0, qf = 0; e < supportSize; ++e) {
1102:           PetscInt val;

1104:           DMLabelGetValue(label, support[e], &val);
1105:           /* Split and unsplit edges produce hybrid faces */
1106:           if (val == 1) ++qf;
1107:           if (val == (shift2 + 1)) ++qf;
1108:         }
1109:         DMPlexSetSupportSize(sdm, hybedge, qf);
1110:       } else if (dep == dim-2) {
1111:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1112:         PetscInt       val;

1114:         for (e = 0, qf = 0; e < supportSize; ++e) {
1115:           DMLabelGetValue(label, support[e], &val);
1116:           if (val == dim-1) qf += 2;
1117:           else              ++qf;
1118:         }
1119:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1120:         DMPlexSetSupportSize(sdm, newp, qf+2);
1121:         /* Add hybrid face */
1122:         for (e = 0, qf = 0; e < supportSize; ++e) {
1123:           DMLabelGetValue(label, support[e], &val);
1124:           if (val == dim-1) ++qf;
1125:         }
1126:         DMPlexSetConeSize(sdm, hybface, 4);
1127:         DMPlexSetSupportSize(sdm, hybface, qf);
1128:       }
1129:     }
1130:   }
1131:   /* Step 4: Setup split DM */
1132:   DMSetUp(sdm);
1133:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1134:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1135:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1136:   /* Step 6: Set cones and supports for new points */
1137:   for (dep = 0; dep <= depth; ++dep) {
1138:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1139:       const PetscInt  oldp   = splitPoints[dep][p];
1140:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1141:       const PetscInt  splitp = p    + pMaxNew[dep];
1142:       const PetscInt *cone, *support, *ornt;
1143:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1145:       DMPlexGetConeSize(dm, oldp, &coneSize);
1146:       DMPlexGetCone(dm, oldp, &cone);
1147:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1148:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1149:       DMPlexGetSupport(dm, oldp, &support);
1150:       if (dep == depth-1) {
1151:         PetscBool       hasUnsplit = PETSC_FALSE;
1152:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1153:         const PetscInt *supportF;

1155:         /* Split face:       copy in old face to new face to start */
1156:         DMPlexGetSupport(sdm, newp,  &supportF);
1157:         DMPlexSetSupport(sdm, splitp, supportF);
1158:         /* Split old face:   old vertices/edges in cone so no change */
1159:         /* Split new face:   new vertices/edges in cone */
1160:         for (q = 0; q < coneSize; ++q) {
1161:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1162:           if (v < 0) {
1163:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1164:             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);
1165:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1166:             hasUnsplit   = PETSC_TRUE;
1167:           } else {
1168:             coneNew[2+q] = v + pMaxNew[dep-1];
1169:             if (dep > 1) {
1170:               const PetscInt *econe;
1171:               PetscInt        econeSize, r, vs, vu;

1173:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1174:               DMPlexGetCone(dm, cone[q], &econe);
1175:               for (r = 0; r < econeSize; ++r) {
1176:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
1177:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1178:                 if (vs >= 0) continue;
1179:                 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);
1180:                 hasUnsplit   = PETSC_TRUE;
1181:               }
1182:             }
1183:           }
1184:         }
1185:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1186:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1187:         /* Face support */
1188:         for (s = 0; s < supportSize; ++s) {
1189:           PetscInt val;

1191:           DMLabelGetValue(label, support[s], &val);
1192:           if (val < 0) {
1193:             /* Split old face:   Replace negative side cell with cohesive cell */
1194:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1195:           } else {
1196:             /* Split new face:   Replace positive side cell with cohesive cell */
1197:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1198:             /* Get orientation for cohesive face */
1199:             {
1200:               const PetscInt *ncone, *nconeO;
1201:               PetscInt        nconeSize, nc;

1203:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1204:               DMPlexGetCone(dm, support[s], &ncone);
1205:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1206:               for (nc = 0; nc < nconeSize; ++nc) {
1207:                 if (ncone[nc] == oldp) {
1208:                   coneONew[0] = nconeO[nc];
1209:                   break;
1210:                 }
1211:               }
1212:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1213:             }
1214:           }
1215:         }
1216:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1217:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1218:         coneNew[1]  = splitp;
1219:         coneONew[1] = coneONew[0];
1220:         for (q = 0; q < coneSize; ++q) {
1221:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1222:           if (v < 0) {
1223:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1224:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1225:             coneONew[2+q] = 0;
1226:           } else {
1227:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1228:           }
1229:           coneONew[2+q] = 0;
1230:         }
1231:         DMPlexSetCone(sdm, hybcell, coneNew);
1232:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1233:         /* Label the hybrid cells on the boundary of the split */
1234:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1235:       } else if (dep == 0) {
1236:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1242:           DMLabelGetValue(label, support[e], &val);
1243:           if ((val == 1) || (val == (shift + 1))) {
1244:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1245:           }
1246:         }
1247:         supportNew[qn] = hybedge;
1248:         DMPlexSetSupport(sdm, newp, supportNew);
1249:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1250:         for (e = 0, qp = 0; e < supportSize; ++e) {
1251:           PetscInt val, edge;

1253:           DMLabelGetValue(label, support[e], &val);
1254:           if (val == 1) {
1255:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1256:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1257:             supportNew[qp++] = edge + pMaxNew[dep+1];
1258:           } else if (val == -(shift + 1)) {
1259:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1260:           }
1261:         }
1262:         supportNew[qp] = hybedge;
1263:         DMPlexSetSupport(sdm, splitp, supportNew);
1264:         /* Hybrid edge:    Old and new split vertex */
1265:         coneNew[0] = newp;
1266:         coneNew[1] = splitp;
1267:         DMPlexSetCone(sdm, hybedge, coneNew);
1268:         for (e = 0, qf = 0; e < supportSize; ++e) {
1269:           PetscInt val, edge;

1271:           DMLabelGetValue(label, support[e], &val);
1272:           if (val == 1) {
1273:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1274:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1275:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1276:           }
1277:         }
1278:         DMPlexSetSupport(sdm, hybedge, supportNew);
1279:       } else if (dep == dim-2) {
1280:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1282:         /* Split old edge:   old vertices in cone so no change */
1283:         /* Split new edge:   new vertices in cone */
1284:         for (q = 0; q < coneSize; ++q) {
1285:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1286:           if (v < 0) {
1287:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1288:             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);
1289:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1290:           } else {
1291:             coneNew[q] = v + pMaxNew[dep-1];
1292:           }
1293:         }
1294:         DMPlexSetCone(sdm, splitp, coneNew);
1295:         /* Split old edge: Faces in positive side cells and old split faces */
1296:         for (e = 0, q = 0; e < supportSize; ++e) {
1297:           PetscInt val;

1299:           DMLabelGetValue(label, support[e], &val);
1300:           if (val == dim-1) {
1301:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1302:           } else if (val == (shift + dim-1)) {
1303:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1304:           }
1305:         }
1306:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1307:         DMPlexSetSupport(sdm, newp, supportNew);
1308:         /* Split new edge: Faces in negative side cells and new split faces */
1309:         for (e = 0, q = 0; e < supportSize; ++e) {
1310:           PetscInt val, face;

1312:           DMLabelGetValue(label, support[e], &val);
1313:           if (val == dim-1) {
1314:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1315:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1316:             supportNew[q++] = face + pMaxNew[dep+1];
1317:           } else if (val == -(shift + dim-1)) {
1318:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1319:           }
1320:         }
1321:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1322:         DMPlexSetSupport(sdm, splitp, supportNew);
1323:         /* Hybrid face */
1324:         coneNew[0] = newp;
1325:         coneNew[1] = splitp;
1326:         for (v = 0; v < coneSize; ++v) {
1327:           PetscInt vertex;
1328:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1329:           if (vertex < 0) {
1330:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1331:             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);
1332:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1333:           } else {
1334:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1335:           }
1336:         }
1337:         DMPlexSetCone(sdm, hybface, coneNew);
1338:         for (e = 0, qf = 0; e < supportSize; ++e) {
1339:           PetscInt val, face;

1341:           DMLabelGetValue(label, support[e], &val);
1342:           if (val == dim-1) {
1343:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1344:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1345:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1346:           }
1347:         }
1348:         DMPlexSetSupport(sdm, hybface, supportNew);
1349:       }
1350:     }
1351:   }
1352:   for (dep = 0; dep <= depth; ++dep) {
1353:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1354:       const PetscInt  oldp   = unsplitPoints[dep][p];
1355:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1356:       const PetscInt *cone, *support, *ornt;
1357:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1359:       DMPlexGetConeSize(dm, oldp, &coneSize);
1360:       DMPlexGetCone(dm, oldp, &cone);
1361:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1362:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1363:       DMPlexGetSupport(dm, oldp, &support);
1364:       if (dep == 0) {
1365:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1367:         /* Unsplit vertex */
1368:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1369:         for (s = 0, q = 0; s < supportSize; ++s) {
1370:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1371:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1372:           if (e >= 0) {
1373:             supportNew[q++] = e + pMaxNew[dep+1];
1374:           }
1375:         }
1376:         supportNew[q++] = hybedge;
1377:         supportNew[q++] = hybedge;
1378:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1379:         DMPlexSetSupport(sdm, newp, supportNew);
1380:         /* Hybrid edge */
1381:         coneNew[0] = newp;
1382:         coneNew[1] = newp;
1383:         DMPlexSetCone(sdm, hybedge, coneNew);
1384:         for (e = 0, qf = 0; e < supportSize; ++e) {
1385:           PetscInt val, edge;

1387:           DMLabelGetValue(label, support[e], &val);
1388:           if (val == 1) {
1389:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1390:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1391:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1392:           } else if  (val ==  (shift2 + 1)) {
1393:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1394:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1395:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1396:           }
1397:         }
1398:         DMPlexSetSupport(sdm, hybedge, supportNew);
1399:       } else if (dep == dim-2) {
1400:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1406:           DMLabelGetValue(label, support[f], &val);
1407:           if (val == dim-1) {
1408:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1409:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1410:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1411:             supportNew[qf++] = face + pMaxNew[dep+1];
1412:           } else {
1413:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1414:           }
1415:         }
1416:         supportNew[qf++] = hybface;
1417:         supportNew[qf++] = hybface;
1418:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1419:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1420:         DMPlexSetSupport(sdm, newp, supportNew);
1421:         /* Add hybrid face */
1422:         coneNew[0] = newp;
1423:         coneNew[1] = newp;
1424:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1425:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1426:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1427:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1428:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1429:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1430:         DMPlexSetCone(sdm, hybface, coneNew);
1431:         for (f = 0, qf = 0; f < supportSize; ++f) {
1432:           PetscInt val, face;

1434:           DMLabelGetValue(label, support[f], &val);
1435:           if (val == dim-1) {
1436:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1437:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1438:           }
1439:         }
1440:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1441:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1442:         DMPlexSetSupport(sdm, hybface, supportNew);
1443:       }
1444:     }
1445:   }
1446:   /* Step 6b: Replace split points in negative side cones */
1447:   for (sp = 0; sp < numSP; ++sp) {
1448:     PetscInt        dep = values[sp];
1449:     IS              pIS;
1450:     PetscInt        numPoints;
1451:     const PetscInt *points;

1453:     if (dep >= 0) continue;
1454:     DMLabelGetStratumIS(label, dep, &pIS);
1455:     if (!pIS) continue;
1456:     dep  = -dep - shift;
1457:     ISGetLocalSize(pIS, &numPoints);
1458:     ISGetIndices(pIS, &points);
1459:     for (p = 0; p < numPoints; ++p) {
1460:       const PetscInt  oldp = points[p];
1461:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1462:       const PetscInt *cone;
1463:       PetscInt        coneSize, c;
1464:       /* PetscBool       replaced = PETSC_FALSE; */

1466:       /* Negative edge: replace split vertex */
1467:       /* Negative cell: replace split face */
1468:       DMPlexGetConeSize(sdm, newp, &coneSize);
1469:       DMPlexGetCone(sdm, newp, &cone);
1470:       for (c = 0; c < coneSize; ++c) {
1471:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1472:         PetscInt       csplitp, cp, val;

1474:         DMLabelGetValue(label, coldp, &val);
1475:         if (val == dep-1) {
1476:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1477:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1478:           csplitp  = pMaxNew[dep-1] + cp;
1479:           DMPlexInsertCone(sdm, newp, c, csplitp);
1480:           /* replaced = PETSC_TRUE; */
1481:         }
1482:       }
1483:       /* Cells with only a vertex or edge on the submesh have no replacement */
1484:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1485:     }
1486:     ISRestoreIndices(pIS, &points);
1487:     ISDestroy(&pIS);
1488:   }
1489:   /* Step 7: Coordinates */
1490:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1491:   DMGetCoordinateSection(sdm, &coordSection);
1492:   DMGetCoordinatesLocal(sdm, &coordinates);
1493:   VecGetArray(coordinates, &coords);
1494:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1495:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1496:     const PetscInt splitp = pMaxNew[0] + v;
1497:     PetscInt       dof, off, soff, d;

1499:     PetscSectionGetDof(coordSection, newp, &dof);
1500:     PetscSectionGetOffset(coordSection, newp, &off);
1501:     PetscSectionGetOffset(coordSection, splitp, &soff);
1502:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1503:   }
1504:   VecRestoreArray(coordinates, &coords);
1505:   /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1506:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1507:   /* Step 9: Labels */
1508:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1509:   DMGetNumLabels(sdm, &numLabels);
1510:   for (dep = 0; dep <= depth; ++dep) {
1511:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1512:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1513:       const PetscInt splitp = pMaxNew[dep] + p;
1514:       PetscInt       l;

1516:       if (splitLabel) {
1517:         const PetscInt val = 100 + dep;

1519:         DMLabelSetValue(splitLabel, newp,    val);
1520:         DMLabelSetValue(splitLabel, splitp, -val);
1521:       }
1522:       for (l = 0; l < numLabels; ++l) {
1523:         DMLabel     mlabel;
1524:         const char *lname;
1525:         PetscInt    val;
1526:         PetscBool   isDepth;

1528:         DMGetLabelName(sdm, l, &lname);
1529:         PetscStrcmp(lname, "depth", &isDepth);
1530:         if (isDepth) continue;
1531:         DMGetLabel(sdm, lname, &mlabel);
1532:         DMLabelGetValue(mlabel, newp, &val);
1533:         if (val >= 0) {
1534:           DMLabelSetValue(mlabel, splitp, val);
1535:         }
1536:       }
1537:     }
1538:   }
1539:   for (sp = 0; sp < numSP; ++sp) {
1540:     const PetscInt dep = values[sp];

1542:     if ((dep < 0) || (dep > depth)) continue;
1543:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1544:     ISDestroy(&splitIS[dep]);
1545:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1546:     ISDestroy(&unsplitIS[dep]);
1547:   }
1548:   if (label) {
1549:     ISRestoreIndices(valueIS, &values);
1550:     ISDestroy(&valueIS);
1551:   }
1552:   for (d = 0; d <= depth; ++d) {
1553:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1554:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1555:   }
1556:   DMPlexSetHybridBounds(sdm, depth >= 0 ? pMaxNew[depth] : PETSC_DETERMINE, depth>1 ? pMaxNew[depth-1] : PETSC_DETERMINE, depth>2 ? pMaxNew[1] : PETSC_DETERMINE, depth >= 0 ? pMaxNew[0] : PETSC_DETERMINE);
1557:   PetscFree3(coneNew, coneONew, supportNew);
1558:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1559:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1560:   return(0);
1561: }

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

1566:   Collective on dm

1568:   Input Parameters:
1569: + dm - The original DM
1570: - label - The label specifying the boundary faces (this could be auto-generated)

1572:   Output Parameters:
1573: + splitLabel - The label containing the split points, or NULL if no output is desired
1574: - dmSplit - The new DM

1576:   Level: developer

1578: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1579: @*/
1580: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1581: {
1582:   DM             sdm;
1583:   PetscInt       dim;

1589:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1590:   DMSetType(sdm, DMPLEX);
1591:   DMGetDimension(dm, &dim);
1592:   DMSetDimension(sdm, dim);
1593:   switch (dim) {
1594:   case 2:
1595:   case 3:
1596:     DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1597:     break;
1598:   default:
1599:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1600:   }
1601:   *dmSplit = sdm;
1602:   return(0);
1603: }

1605: /* Returns the side of the surface for a given cell with a face on the surface */
1606: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1607: {
1608:   const PetscInt *cone, *ornt;
1609:   PetscInt        dim, coneSize, c;
1610:   PetscErrorCode  ierr;

1613:   *pos = PETSC_TRUE;
1614:   DMGetDimension(dm, &dim);
1615:   DMPlexGetConeSize(dm, cell, &coneSize);
1616:   DMPlexGetCone(dm, cell, &cone);
1617:   DMPlexGetConeOrientation(dm, cell, &ornt);
1618:   for (c = 0; c < coneSize; ++c) {
1619:     if (cone[c] == face) {
1620:       PetscInt o = ornt[c];

1622:       if (subdm) {
1623:         const PetscInt *subcone, *subornt;
1624:         PetscInt        subpoint, subface, subconeSize, sc;

1626:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1627:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1628:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1629:         DMPlexGetCone(subdm, subpoint, &subcone);
1630:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1631:         for (sc = 0; sc < subconeSize; ++sc) {
1632:           if (subcone[sc] == subface) {
1633:             o = subornt[0];
1634:             break;
1635:           }
1636:         }
1637:         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);
1638:       }
1639:       if (o >= 0) *pos = PETSC_TRUE;
1640:       else        *pos = PETSC_FALSE;
1641:       break;
1642:     }
1643:   }
1644:   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);
1645:   return(0);
1646: }

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

1652:   Input Parameters:
1653: + dm     - The DM
1654: . label  - A DMLabel marking the surface
1655: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1656: . flip   - Flag to flip the submesh normal and replace points on the other side
1657: - subdm  - The subDM associated with the label, or NULL

1659:   Output Parameter:
1660: . label - A DMLabel marking all surface points

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

1664:   Level: developer

1666: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1667: @*/
1668: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1669: {
1670:   DMLabel         depthLabel;
1671:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1672:   const PetscInt *points, *subpoints;
1673:   const PetscInt  rev   = flip ? -1 : 1;
1674:   PetscInt       *pMax;
1675:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1676:   PetscErrorCode  ierr;

1679:   DMPlexGetDepth(dm, &depth);
1680:   DMGetDimension(dm, &dim);
1681:   pSize = PetscMax(depth, dim) + 1;
1682:   PetscMalloc1(pSize,&pMax);
1683:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1684:   DMPlexGetDepthLabel(dm, &depthLabel);
1685:   DMGetDimension(dm, &dim);
1686:   if (subdm) {
1687:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1688:     if (subpointIS) {
1689:       ISGetLocalSize(subpointIS, &numSubpoints);
1690:       ISGetIndices(subpointIS, &subpoints);
1691:     }
1692:   }
1693:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1694:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1695:   if (!dimIS) {
1696:     PetscFree(pMax);
1697:     ISDestroy(&subpointIS);
1698:     return(0);
1699:   }
1700:   ISGetLocalSize(dimIS, &numPoints);
1701:   ISGetIndices(dimIS, &points);
1702:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1703:     const PetscInt *support;
1704:     PetscInt        supportSize, s;

1706:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1707: #if 0
1708:     if (supportSize != 2) {
1709:       const PetscInt *lp;
1710:       PetscInt        Nlp, pind;

1712:       /* Check that for a cell with a single support face, that face is in the SF */
1713:       /*   THis check only works for the remote side. We would need root side information */
1714:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1715:       PetscFindInt(points[p], Nlp, lp, &pind);
1716:       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);
1717:     }
1718: #endif
1719:     DMPlexGetSupport(dm, points[p], &support);
1720:     for (s = 0; s < supportSize; ++s) {
1721:       const PetscInt *cone;
1722:       PetscInt        coneSize, c;
1723:       PetscBool       pos;

1725:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1726:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1727:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1728:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1729:       /* Put faces touching the fault in the label */
1730:       DMPlexGetConeSize(dm, support[s], &coneSize);
1731:       DMPlexGetCone(dm, support[s], &cone);
1732:       for (c = 0; c < coneSize; ++c) {
1733:         const PetscInt point = cone[c];

1735:         DMLabelGetValue(label, point, &val);
1736:         if (val == -1) {
1737:           PetscInt *closure = NULL;
1738:           PetscInt  closureSize, cl;

1740:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1741:           for (cl = 0; cl < closureSize*2; cl += 2) {
1742:             const PetscInt clp  = closure[cl];
1743:             PetscInt       bval = -1;

1745:             DMLabelGetValue(label, clp, &val);
1746:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1747:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1748:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1749:               break;
1750:             }
1751:           }
1752:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1753:         }
1754:       }
1755:     }
1756:   }
1757:   ISRestoreIndices(dimIS, &points);
1758:   ISDestroy(&dimIS);
1759:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1760:   ISDestroy(&subpointIS);
1761:   /* Mark boundary points as unsplit */
1762:   if (blabel) {
1763:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1764:     ISGetLocalSize(dimIS, &numPoints);
1765:     ISGetIndices(dimIS, &points);
1766:     for (p = 0; p < numPoints; ++p) {
1767:       const PetscInt point = points[p];
1768:       PetscInt       val, bval;

1770:       DMLabelGetValue(blabel, point, &bval);
1771:       if (bval >= 0) {
1772:         DMLabelGetValue(label, point, &val);
1773:         if ((val < 0) || (val > dim)) {
1774:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1775:           DMLabelClearValue(blabel, point, bval);
1776:         }
1777:       }
1778:     }
1779:     for (p = 0; p < numPoints; ++p) {
1780:       const PetscInt point = points[p];
1781:       PetscInt       val, bval;

1783:       DMLabelGetValue(blabel, point, &bval);
1784:       if (bval >= 0) {
1785:         const PetscInt *cone,    *support;
1786:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1788:         /* Mark as unsplit */
1789:         DMLabelGetValue(label, point, &val);
1790:         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);
1791:         DMLabelClearValue(label, point, val);
1792:         DMLabelSetValue(label, point, shift2+val);
1793:         /* Check for cross-edge
1794:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1795:         if (val != 0) continue;
1796:         DMPlexGetSupport(dm, point, &support);
1797:         DMPlexGetSupportSize(dm, point, &supportSize);
1798:         for (s = 0; s < supportSize; ++s) {
1799:           DMPlexGetCone(dm, support[s], &cone);
1800:           DMPlexGetConeSize(dm, support[s], &coneSize);
1801:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1802:           DMLabelGetValue(blabel, cone[0], &valA);
1803:           DMLabelGetValue(blabel, cone[1], &valB);
1804:           DMLabelGetValue(blabel, support[s], &valE);
1805:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1806:         }
1807:       }
1808:     }
1809:     ISRestoreIndices(dimIS, &points);
1810:     ISDestroy(&dimIS);
1811:   }
1812:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1813:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1814:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1815:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1816:   cMax = cMax < 0 ? cEnd : cMax;
1817:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1818:   DMLabelGetStratumIS(label, 0, &dimIS);
1819:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1820:   if (dimIS && crossEdgeIS) {
1821:     IS vertIS = dimIS;

1823:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1824:     ISDestroy(&crossEdgeIS);
1825:     ISDestroy(&vertIS);
1826:   }
1827:   if (!dimIS) {
1828:     PetscFree(pMax);
1829:     return(0);
1830:   }
1831:   ISGetLocalSize(dimIS, &numPoints);
1832:   ISGetIndices(dimIS, &points);
1833:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1834:     PetscInt *star = NULL;
1835:     PetscInt  starSize, s;
1836:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1838:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1839:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1840:     while (again) {
1841:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1842:       again = 0;
1843:       for (s = 0; s < starSize*2; s += 2) {
1844:         const PetscInt  point = star[s];
1845:         const PetscInt *cone;
1846:         PetscInt        coneSize, c;

1848:         if ((point < cStart) || (point >= cMax)) continue;
1849:         DMLabelGetValue(label, point, &val);
1850:         if (val != -1) continue;
1851:         again = again == 1 ? 1 : 2;
1852:         DMPlexGetConeSize(dm, point, &coneSize);
1853:         DMPlexGetCone(dm, point, &cone);
1854:         for (c = 0; c < coneSize; ++c) {
1855:           DMLabelGetValue(label, cone[c], &val);
1856:           if (val != -1) {
1857:             const PetscInt *ccone;
1858:             PetscInt        cconeSize, cc, side;

1860:             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);
1861:             if (val > 0) side =  1;
1862:             else         side = -1;
1863:             DMLabelSetValue(label, point, side*(shift+dim));
1864:             /* Mark cell faces which touch the fault */
1865:             DMPlexGetConeSize(dm, point, &cconeSize);
1866:             DMPlexGetCone(dm, point, &ccone);
1867:             for (cc = 0; cc < cconeSize; ++cc) {
1868:               PetscInt *closure = NULL;
1869:               PetscInt  closureSize, cl;

1871:               DMLabelGetValue(label, ccone[cc], &val);
1872:               if (val != -1) continue;
1873:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1874:               for (cl = 0; cl < closureSize*2; cl += 2) {
1875:                 const PetscInt clp = closure[cl];

1877:                 DMLabelGetValue(label, clp, &val);
1878:                 if (val == -1) continue;
1879:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1880:                 break;
1881:               }
1882:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1883:             }
1884:             again = 1;
1885:             break;
1886:           }
1887:         }
1888:       }
1889:     }
1890:     /* Classify the rest by cell membership */
1891:     for (s = 0; s < starSize*2; s += 2) {
1892:       const PetscInt point = star[s];

1894:       DMLabelGetValue(label, point, &val);
1895:       if (val == -1) {
1896:         PetscInt  *sstar = NULL;
1897:         PetscInt   sstarSize, ss;
1898:         PetscBool  marked = PETSC_FALSE;

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

1904:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1905:           DMLabelGetValue(label, spoint, &val);
1906:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1907:           DMLabelGetValue(depthLabel, point, &dep);
1908:           if (val > 0) {
1909:             DMLabelSetValue(label, point,   shift+dep);
1910:           } else {
1911:             DMLabelSetValue(label, point, -(shift+dep));
1912:           }
1913:           marked = PETSC_TRUE;
1914:           break;
1915:         }
1916:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1917:         DMLabelGetValue(depthLabel, point, &dep);
1918:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1919:       }
1920:     }
1921:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1922:   }
1923:   ISRestoreIndices(dimIS, &points);
1924:   ISDestroy(&dimIS);
1925:   /* If any faces touching the fault divide cells on either side, split them */
1926:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1927:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1928:   ISExpand(facePosIS, faceNegIS, &dimIS);
1929:   ISDestroy(&facePosIS);
1930:   ISDestroy(&faceNegIS);
1931:   ISGetLocalSize(dimIS, &numPoints);
1932:   ISGetIndices(dimIS, &points);
1933:   for (p = 0; p < numPoints; ++p) {
1934:     const PetscInt  point = points[p];
1935:     const PetscInt *support;
1936:     PetscInt        supportSize, valA, valB;

1938:     DMPlexGetSupportSize(dm, point, &supportSize);
1939:     if (supportSize != 2) continue;
1940:     DMPlexGetSupport(dm, point, &support);
1941:     DMLabelGetValue(label, support[0], &valA);
1942:     DMLabelGetValue(label, support[1], &valB);
1943:     if ((valA == -1) || (valB == -1)) continue;
1944:     if (valA*valB > 0) continue;
1945:     /* Split the face */
1946:     DMLabelGetValue(label, point, &valA);
1947:     DMLabelClearValue(label, point, valA);
1948:     DMLabelSetValue(label, point, dim-1);
1949:     /* Label its closure:
1950:       unmarked: label as unsplit
1951:       incident: relabel as split
1952:       split:    do nothing
1953:     */
1954:     {
1955:       PetscInt *closure = NULL;
1956:       PetscInt  closureSize, cl;

1958:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1959:       for (cl = 0; cl < closureSize*2; cl += 2) {
1960:         DMLabelGetValue(label, closure[cl], &valA);
1961:         if (valA == -1) { /* Mark as unsplit */
1962:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1963:           DMLabelSetValue(label, closure[cl], shift2+dep);
1964:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1965:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1966:           DMLabelClearValue(label, closure[cl], valA);
1967:           DMLabelSetValue(label, closure[cl], dep);
1968:         }
1969:       }
1970:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1971:     }
1972:   }
1973:   ISRestoreIndices(dimIS, &points);
1974:   ISDestroy(&dimIS);
1975:   PetscFree(pMax);
1976:   return(0);
1977: }

1979: /* Check that no cell have all vertices on the fault */
1980: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
1981: {
1982:   IS              subpointIS;
1983:   const PetscInt *dmpoints;
1984:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;
1985:   PetscErrorCode  ierr;

1988:   if (!label) return(0);
1989:   DMLabelGetDefaultValue(label, &defaultValue);
1990:   DMPlexCreateSubpointIS(subdm, &subpointIS);
1991:   if (!subpointIS) return(0);
1992:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
1993:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1994:   ISGetIndices(subpointIS, &dmpoints);
1995:   for (c = cStart; c < cEnd; ++c) {
1996:     PetscBool invalidCell = PETSC_TRUE;
1997:     PetscInt *closure     = NULL;
1998:     PetscInt  closureSize, cl;

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

2004:       if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2005:       DMLabelGetValue(label, closure[cl], &value);
2006:       if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
2007:     }
2008:     DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2009:     if (invalidCell) {
2010:       ISRestoreIndices(subpointIS, &dmpoints);
2011:       ISDestroy(&subpointIS);
2012:       DMDestroy(&subdm);
2013:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
2014:     }
2015:   }
2016:   ISRestoreIndices(subpointIS, &dmpoints);
2017:   ISDestroy(&subpointIS);
2018:   return(0);
2019: }


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

2025:   Collective on dm

2027:   Input Parameters:
2028: + dm - The original DM
2029: . label - The label specifying the interface vertices
2030: - bdlabel - The optional label specifying the interface boundary vertices

2032:   Output Parameters:
2033: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2034: . splitLabel - The label containing the split points, or NULL if no output is desired
2035: . dmInterface - The new interface DM, or NULL
2036: - dmHybrid - The new DM with cohesive cells

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

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

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

2051:   Level: developer

2053: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
2054: @*/
2055: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2056: {
2057:   DM             idm;
2058:   DMLabel        subpointMap, hlabel, slabel = NULL;
2059:   PetscInt       dim;

2069:   DMGetDimension(dm, &dim);
2070:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2071:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
2072:   DMPlexOrient(idm);
2073:   DMPlexGetSubpointMap(idm, &subpointMap);
2074:   DMLabelDuplicate(subpointMap, &hlabel);
2075:   DMLabelClearStratum(hlabel, dim);
2076:   if (splitLabel) {
2077:     const char *name;
2078:     char        sname[PETSC_MAX_PATH_LEN];

2080:     PetscObjectGetName((PetscObject) hlabel, &name);
2081:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2082:     PetscStrcat(sname, " split");
2083:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2084:   }
2085:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2086:   if (dmInterface) {*dmInterface = idm;}
2087:   else             {DMDestroy(&idm);}
2088:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2089:   if (hybridLabel) *hybridLabel = hlabel;
2090:   else             {DMLabelDestroy(&hlabel);}
2091:   if (splitLabel)  *splitLabel  = slabel;
2092:   return(0);
2093: }

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

2097:      For any marked cell, the marked vertices constitute a single face
2098: */
2099: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2100: {
2101:   IS               subvertexIS = NULL;
2102:   const PetscInt  *subvertices;
2103:   PetscInt        *pStart, *pEnd, *pMax, pSize;
2104:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
2105:   PetscErrorCode   ierr;

2108:   *numFaces = 0;
2109:   *nFV      = 0;
2110:   DMPlexGetDepth(dm, &depth);
2111:   DMGetDimension(dm, &dim);
2112:   pSize = PetscMax(depth, dim) + 1;
2113:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
2114:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2115:   for (d = 0; d <= depth; ++d) {
2116:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2117:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2118:   }
2119:   /* Loop over initial vertices and mark all faces in the collective star() */
2120:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2121:   if (subvertexIS) {
2122:     ISGetSize(subvertexIS, &numSubVerticesInitial);
2123:     ISGetIndices(subvertexIS, &subvertices);
2124:   }
2125:   for (v = 0; v < numSubVerticesInitial; ++v) {
2126:     const PetscInt vertex = subvertices[v];
2127:     PetscInt      *star   = NULL;
2128:     PetscInt       starSize, s, numCells = 0, c;

2130:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2131:     for (s = 0; s < starSize*2; s += 2) {
2132:       const PetscInt point = star[s];
2133:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2134:     }
2135:     for (c = 0; c < numCells; ++c) {
2136:       const PetscInt cell    = star[c];
2137:       PetscInt      *closure = NULL;
2138:       PetscInt       closureSize, cl;
2139:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

2141:       DMLabelGetValue(subpointMap, cell, &cellLoc);
2142:       if (cellLoc == 2) continue;
2143:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2144:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2145:       for (cl = 0; cl < closureSize*2; cl += 2) {
2146:         const PetscInt point = closure[cl];
2147:         PetscInt       vertexLoc;

2149:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2150:           ++numCorners;
2151:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2152:           if (vertexLoc == value) closure[faceSize++] = point;
2153:         }
2154:       }
2155:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2156:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2157:       if (faceSize == *nFV) {
2158:         const PetscInt *cells = NULL;
2159:         PetscInt        numCells, nc;

2161:         ++(*numFaces);
2162:         for (cl = 0; cl < faceSize; ++cl) {
2163:           DMLabelSetValue(subpointMap, closure[cl], 0);
2164:         }
2165:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2166:         for (nc = 0; nc < numCells; ++nc) {
2167:           DMLabelSetValue(subpointMap, cells[nc], 2);
2168:         }
2169:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2170:       }
2171:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2172:     }
2173:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2174:   }
2175:   if (subvertexIS) {
2176:     ISRestoreIndices(subvertexIS, &subvertices);
2177:   }
2178:   ISDestroy(&subvertexIS);
2179:   PetscFree3(pStart,pEnd,pMax);
2180:   return(0);
2181: }

2183: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2184: {
2185:   IS               subvertexIS = NULL;
2186:   const PetscInt  *subvertices;
2187:   PetscInt        *pStart, *pEnd, *pMax;
2188:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2189:   PetscErrorCode   ierr;

2192:   DMGetDimension(dm, &dim);
2193:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2194:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2195:   for (d = 0; d <= dim; ++d) {
2196:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2197:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2198:   }
2199:   /* Loop over initial vertices and mark all faces in the collective star() */
2200:   if (vertexLabel) {
2201:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2202:     if (subvertexIS) {
2203:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2204:       ISGetIndices(subvertexIS, &subvertices);
2205:     }
2206:   }
2207:   for (v = 0; v < numSubVerticesInitial; ++v) {
2208:     const PetscInt vertex = subvertices[v];
2209:     PetscInt      *star   = NULL;
2210:     PetscInt       starSize, s, numFaces = 0, f;

2212:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2213:     for (s = 0; s < starSize*2; s += 2) {
2214:       const PetscInt point = star[s];
2215:       PetscInt       faceLoc;

2217:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2218:         if (markedFaces) {
2219:           DMLabelGetValue(vertexLabel, point, &faceLoc);
2220:           if (faceLoc < 0) continue;
2221:         }
2222:         star[numFaces++] = point;
2223:       }
2224:     }
2225:     for (f = 0; f < numFaces; ++f) {
2226:       const PetscInt face    = star[f];
2227:       PetscInt      *closure = NULL;
2228:       PetscInt       closureSize, c;
2229:       PetscInt       faceLoc;

2231:       DMLabelGetValue(subpointMap, face, &faceLoc);
2232:       if (faceLoc == dim-1) continue;
2233:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2234:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2235:       for (c = 0; c < closureSize*2; c += 2) {
2236:         const PetscInt point = closure[c];
2237:         PetscInt       vertexLoc;

2239:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2240:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2241:           if (vertexLoc != value) break;
2242:         }
2243:       }
2244:       if (c == closureSize*2) {
2245:         const PetscInt *support;
2246:         PetscInt        supportSize, s;

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

2251:           for (d = 0; d < dim; ++d) {
2252:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2253:               DMLabelSetValue(subpointMap, point, d);
2254:               break;
2255:             }
2256:           }
2257:         }
2258:         DMPlexGetSupportSize(dm, face, &supportSize);
2259:         DMPlexGetSupport(dm, face, &support);
2260:         for (s = 0; s < supportSize; ++s) {
2261:           DMLabelSetValue(subpointMap, support[s], dim);
2262:         }
2263:       }
2264:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2265:     }
2266:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2267:   }
2268:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2269:   ISDestroy(&subvertexIS);
2270:   PetscFree3(pStart,pEnd,pMax);
2271:   return(0);
2272: }

2274: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2275: {
2276:   DMLabel         label = NULL;
2277:   const PetscInt *cone;
2278:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2279:   PetscErrorCode  ierr;

2282:   *numFaces = 0;
2283:   *nFV = 0;
2284:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2285:   *subCells = NULL;
2286:   DMGetDimension(dm, &dim);
2287:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2288:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2289:   if (cMax < 0) return(0);
2290:   if (label) {
2291:     for (c = cMax; c < cEnd; ++c) {
2292:       PetscInt val;

2294:       DMLabelGetValue(label, c, &val);
2295:       if (val == value) {
2296:         ++(*numFaces);
2297:         DMPlexGetConeSize(dm, c, &coneSize);
2298:       }
2299:     }
2300:   } else {
2301:     *numFaces = cEnd - cMax;
2302:     DMPlexGetConeSize(dm, cMax, &coneSize);
2303:   }
2304:   PetscMalloc1(*numFaces *2, subCells);
2305:   if (!(*numFaces)) return(0);
2306:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2307:   for (c = cMax; c < cEnd; ++c) {
2308:     const PetscInt *cells;
2309:     PetscInt        numCells;

2311:     if (label) {
2312:       PetscInt val;

2314:       DMLabelGetValue(label, c, &val);
2315:       if (val != value) continue;
2316:     }
2317:     DMPlexGetCone(dm, c, &cone);
2318:     for (p = 0; p < *nFV; ++p) {
2319:       DMLabelSetValue(subpointMap, cone[p], 0);
2320:     }
2321:     /* Negative face */
2322:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2323:     /* Not true in parallel
2324:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2325:     for (p = 0; p < numCells; ++p) {
2326:       DMLabelSetValue(subpointMap, cells[p], 2);
2327:       (*subCells)[subc++] = cells[p];
2328:     }
2329:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2330:     /* Positive face is not included */
2331:   }
2332:   return(0);
2333: }

2335: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2336: {
2337:   PetscInt      *pStart, *pEnd;
2338:   PetscInt       dim, cMax, cEnd, c, d;

2342:   DMGetDimension(dm, &dim);
2343:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2344:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2345:   if (cMax < 0) return(0);
2346:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2347:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2348:   for (c = cMax; c < cEnd; ++c) {
2349:     const PetscInt *cone;
2350:     PetscInt       *closure = NULL;
2351:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2353:     if (label) {
2354:       DMLabelGetValue(label, c, &val);
2355:       if (val != value) continue;
2356:     }
2357:     DMPlexGetConeSize(dm, c, &coneSize);
2358:     DMPlexGetCone(dm, c, &cone);
2359:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2360:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2361:     /* Negative face */
2362:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2363:     for (cl = 0; cl < closureSize*2; cl += 2) {
2364:       const PetscInt point = closure[cl];

2366:       for (d = 0; d <= dim; ++d) {
2367:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2368:           DMLabelSetValue(subpointMap, point, d);
2369:           break;
2370:         }
2371:       }
2372:     }
2373:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2374:     /* Cells -- positive face is not included */
2375:     for (cl = 0; cl < 1; ++cl) {
2376:       const PetscInt *support;
2377:       PetscInt        supportSize, s;

2379:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2380:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2381:       DMPlexGetSupport(dm, cone[cl], &support);
2382:       for (s = 0; s < supportSize; ++s) {
2383:         DMLabelSetValue(subpointMap, support[s], dim);
2384:       }
2385:     }
2386:   }
2387:   PetscFree2(pStart, pEnd);
2388:   return(0);
2389: }

2391: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2392: {
2393:   MPI_Comm       comm;
2394:   PetscBool      posOrient = PETSC_FALSE;
2395:   const PetscInt debug     = 0;
2396:   PetscInt       cellDim, faceSize, f;

2400:   PetscObjectGetComm((PetscObject)dm,&comm);
2401:   DMGetDimension(dm, &cellDim);
2402:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2404:   if (cellDim == 1 && numCorners == 2) {
2405:     /* Triangle */
2406:     faceSize  = numCorners-1;
2407:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2408:   } else if (cellDim == 2 && numCorners == 3) {
2409:     /* Triangle */
2410:     faceSize  = numCorners-1;
2411:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2412:   } else if (cellDim == 3 && numCorners == 4) {
2413:     /* Tetrahedron */
2414:     faceSize  = numCorners-1;
2415:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2416:   } else if (cellDim == 1 && numCorners == 3) {
2417:     /* Quadratic line */
2418:     faceSize  = 1;
2419:     posOrient = PETSC_TRUE;
2420:   } else if (cellDim == 2 && numCorners == 4) {
2421:     /* Quads */
2422:     faceSize = 2;
2423:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2424:       posOrient = PETSC_TRUE;
2425:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2426:       posOrient = PETSC_TRUE;
2427:     } else {
2428:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2429:         posOrient = PETSC_FALSE;
2430:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2431:     }
2432:   } else if (cellDim == 2 && numCorners == 6) {
2433:     /* Quadratic triangle (I hate this) */
2434:     /* Edges are determined by the first 2 vertices (corners of edges) */
2435:     const PetscInt faceSizeTri = 3;
2436:     PetscInt       sortedIndices[3], i, iFace;
2437:     PetscBool      found                    = PETSC_FALSE;
2438:     PetscInt       faceVerticesTriSorted[9] = {
2439:       0, 3,  4, /* bottom */
2440:       1, 4,  5, /* right */
2441:       2, 3,  5, /* left */
2442:     };
2443:     PetscInt       faceVerticesTri[9] = {
2444:       0, 3,  4, /* bottom */
2445:       1, 4,  5, /* right */
2446:       2, 5,  3, /* left */
2447:     };

2449:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2450:     PetscSortInt(faceSizeTri, sortedIndices);
2451:     for (iFace = 0; iFace < 3; ++iFace) {
2452:       const PetscInt ii = iFace*faceSizeTri;
2453:       PetscInt       fVertex, cVertex;

2455:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2456:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2457:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2458:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2459:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2460:               faceVertices[fVertex] = origVertices[cVertex];
2461:               break;
2462:             }
2463:           }
2464:         }
2465:         found = PETSC_TRUE;
2466:         break;
2467:       }
2468:     }
2469:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2470:     if (posOriented) *posOriented = PETSC_TRUE;
2471:     return(0);
2472:   } else if (cellDim == 2 && numCorners == 9) {
2473:     /* Quadratic quad (I hate this) */
2474:     /* Edges are determined by the first 2 vertices (corners of edges) */
2475:     const PetscInt faceSizeQuad = 3;
2476:     PetscInt       sortedIndices[3], i, iFace;
2477:     PetscBool      found                      = PETSC_FALSE;
2478:     PetscInt       faceVerticesQuadSorted[12] = {
2479:       0, 1,  4, /* bottom */
2480:       1, 2,  5, /* right */
2481:       2, 3,  6, /* top */
2482:       0, 3,  7, /* left */
2483:     };
2484:     PetscInt       faceVerticesQuad[12] = {
2485:       0, 1,  4, /* bottom */
2486:       1, 2,  5, /* right */
2487:       2, 3,  6, /* top */
2488:       3, 0,  7, /* left */
2489:     };

2491:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2492:     PetscSortInt(faceSizeQuad, sortedIndices);
2493:     for (iFace = 0; iFace < 4; ++iFace) {
2494:       const PetscInt ii = iFace*faceSizeQuad;
2495:       PetscInt       fVertex, cVertex;

2497:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2498:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2499:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2500:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2501:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2502:               faceVertices[fVertex] = origVertices[cVertex];
2503:               break;
2504:             }
2505:           }
2506:         }
2507:         found = PETSC_TRUE;
2508:         break;
2509:       }
2510:     }
2511:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2512:     if (posOriented) *posOriented = PETSC_TRUE;
2513:     return(0);
2514:   } else if (cellDim == 3 && numCorners == 8) {
2515:     /* Hexes
2516:        A hex is two oriented quads with the normal of the first
2517:        pointing up at the second.

2519:           7---6
2520:          /|  /|
2521:         4---5 |
2522:         | 1-|-2
2523:         |/  |/
2524:         0---3

2526:         Faces are determined by the first 4 vertices (corners of faces) */
2527:     const PetscInt faceSizeHex = 4;
2528:     PetscInt       sortedIndices[4], i, iFace;
2529:     PetscBool      found                     = PETSC_FALSE;
2530:     PetscInt       faceVerticesHexSorted[24] = {
2531:       0, 1, 2, 3,  /* bottom */
2532:       4, 5, 6, 7,  /* top */
2533:       0, 3, 4, 5,  /* front */
2534:       2, 3, 5, 6,  /* right */
2535:       1, 2, 6, 7,  /* back */
2536:       0, 1, 4, 7,  /* left */
2537:     };
2538:     PetscInt       faceVerticesHex[24] = {
2539:       1, 2, 3, 0,  /* bottom */
2540:       4, 5, 6, 7,  /* top */
2541:       0, 3, 5, 4,  /* front */
2542:       3, 2, 6, 5,  /* right */
2543:       2, 1, 7, 6,  /* back */
2544:       1, 0, 4, 7,  /* left */
2545:     };

2547:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2548:     PetscSortInt(faceSizeHex, sortedIndices);
2549:     for (iFace = 0; iFace < 6; ++iFace) {
2550:       const PetscInt ii = iFace*faceSizeHex;
2551:       PetscInt       fVertex, cVertex;

2553:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2554:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2555:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2556:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2557:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2558:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2559:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2560:               faceVertices[fVertex] = origVertices[cVertex];
2561:               break;
2562:             }
2563:           }
2564:         }
2565:         found = PETSC_TRUE;
2566:         break;
2567:       }
2568:     }
2569:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2570:     if (posOriented) *posOriented = PETSC_TRUE;
2571:     return(0);
2572:   } else if (cellDim == 3 && numCorners == 10) {
2573:     /* Quadratic tet */
2574:     /* Faces are determined by the first 3 vertices (corners of faces) */
2575:     const PetscInt faceSizeTet = 6;
2576:     PetscInt       sortedIndices[6], i, iFace;
2577:     PetscBool      found                     = PETSC_FALSE;
2578:     PetscInt       faceVerticesTetSorted[24] = {
2579:       0, 1, 2,  6, 7, 8, /* bottom */
2580:       0, 3, 4,  6, 7, 9,  /* front */
2581:       1, 4, 5,  7, 8, 9,  /* right */
2582:       2, 3, 5,  6, 8, 9,  /* left */
2583:     };
2584:     PetscInt       faceVerticesTet[24] = {
2585:       0, 1, 2,  6, 7, 8, /* bottom */
2586:       0, 4, 3,  6, 7, 9,  /* front */
2587:       1, 5, 4,  7, 8, 9,  /* right */
2588:       2, 3, 5,  8, 6, 9,  /* left */
2589:     };

2591:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2592:     PetscSortInt(faceSizeTet, sortedIndices);
2593:     for (iFace=0; iFace < 4; ++iFace) {
2594:       const PetscInt ii = iFace*faceSizeTet;
2595:       PetscInt       fVertex, cVertex;

2597:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2598:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2599:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2600:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2601:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2602:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2603:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2604:               faceVertices[fVertex] = origVertices[cVertex];
2605:               break;
2606:             }
2607:           }
2608:         }
2609:         found = PETSC_TRUE;
2610:         break;
2611:       }
2612:     }
2613:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2614:     if (posOriented) *posOriented = PETSC_TRUE;
2615:     return(0);
2616:   } else if (cellDim == 3 && numCorners == 27) {
2617:     /* Quadratic hexes (I hate this)
2618:        A hex is two oriented quads with the normal of the first
2619:        pointing up at the second.

2621:          7---6
2622:         /|  /|
2623:        4---5 |
2624:        | 3-|-2
2625:        |/  |/
2626:        0---1

2628:        Faces are determined by the first 4 vertices (corners of faces) */
2629:     const PetscInt faceSizeQuadHex = 9;
2630:     PetscInt       sortedIndices[9], i, iFace;
2631:     PetscBool      found                         = PETSC_FALSE;
2632:     PetscInt       faceVerticesQuadHexSorted[54] = {
2633:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2634:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2635:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2636:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2637:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2638:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2639:     };
2640:     PetscInt       faceVerticesQuadHex[54] = {
2641:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2642:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2643:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2644:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2645:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2646:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2647:     };

2649:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2650:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2651:     for (iFace = 0; iFace < 6; ++iFace) {
2652:       const PetscInt ii = iFace*faceSizeQuadHex;
2653:       PetscInt       fVertex, cVertex;

2655:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2656:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2657:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2658:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2659:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2660:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2661:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2662:               faceVertices[fVertex] = origVertices[cVertex];
2663:               break;
2664:             }
2665:           }
2666:         }
2667:         found = PETSC_TRUE;
2668:         break;
2669:       }
2670:     }
2671:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2672:     if (posOriented) *posOriented = PETSC_TRUE;
2673:     return(0);
2674:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2675:   if (!posOrient) {
2676:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2677:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2678:   } else {
2679:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2680:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2681:   }
2682:   if (posOriented) *posOriented = posOrient;
2683:   return(0);
2684: }

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

2690:   Not collective

2692:   Input Parameters:
2693: + dm           - The original mesh
2694: . cell         - The cell mesh point
2695: . faceSize     - The number of vertices on the face
2696: . face         - The face vertices
2697: . numCorners   - The number of vertices on the cell
2698: . indices      - Local numbering of face vertices in cell cone
2699: - origVertices - Original face vertices

2701:   Output Parameter:
2702: + faceVertices - The face vertices properly oriented
2703: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2705:   Level: developer

2707: .seealso: DMPlexGetCone()
2708: @*/
2709: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2710: {
2711:   const PetscInt *cone = NULL;
2712:   PetscInt        coneSize, v, f, v2;
2713:   PetscInt        oppositeVertex = -1;
2714:   PetscErrorCode  ierr;

2717:   DMPlexGetConeSize(dm, cell, &coneSize);
2718:   DMPlexGetCone(dm, cell, &cone);
2719:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2720:     PetscBool found = PETSC_FALSE;

2722:     for (f = 0; f < faceSize; ++f) {
2723:       if (face[f] == cone[v]) {
2724:         found = PETSC_TRUE; break;
2725:       }
2726:     }
2727:     if (found) {
2728:       indices[v2]      = v;
2729:       origVertices[v2] = cone[v];
2730:       ++v2;
2731:     } else {
2732:       oppositeVertex = v;
2733:     }
2734:   }
2735:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2736:   return(0);
2737: }

2739: /*
2740:   DMPlexInsertFace_Internal - Puts a face into the mesh

2742:   Not collective

2744:   Input Parameters:
2745:   + dm              - The DMPlex
2746:   . numFaceVertex   - The number of vertices in the face
2747:   . faceVertices    - The vertices in the face for dm
2748:   . subfaceVertices - The vertices in the face for subdm
2749:   . numCorners      - The number of vertices in the cell
2750:   . cell            - A cell in dm containing the face
2751:   . subcell         - A cell in subdm containing the face
2752:   . firstFace       - First face in the mesh
2753:   - newFacePoint    - Next face in the mesh

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

2758:   Level: developer
2759: */
2760: 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)
2761: {
2762:   MPI_Comm        comm;
2763:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2764:   const PetscInt *faces;
2765:   PetscInt        numFaces, coneSize;
2766:   PetscErrorCode  ierr;

2769:   PetscObjectGetComm((PetscObject)dm,&comm);
2770:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2771:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2772: #if 0
2773:   /* Cannot use this because support() has not been constructed yet */
2774:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2775: #else
2776:   {
2777:     PetscInt f;

2779:     numFaces = 0;
2780:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2781:     for (f = firstFace; f < *newFacePoint; ++f) {
2782:       PetscInt dof, off, d;

2784:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2785:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2786:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2787:       for (d = 0; d < dof; ++d) {
2788:         const PetscInt p = submesh->cones[off+d];
2789:         PetscInt       v;

2791:         for (v = 0; v < numFaceVertices; ++v) {
2792:           if (subfaceVertices[v] == p) break;
2793:         }
2794:         if (v == numFaceVertices) break;
2795:       }
2796:       if (d == dof) {
2797:         numFaces               = 1;
2798:         ((PetscInt*) faces)[0] = f;
2799:       }
2800:     }
2801:   }
2802: #endif
2803:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2804:   else if (numFaces == 1) {
2805:     /* Add the other cell neighbor for this face */
2806:     DMPlexSetCone(subdm, subcell, faces);
2807:   } else {
2808:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2809:     PetscBool posOriented;

2811:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2812:     origVertices        = &orientedVertices[numFaceVertices];
2813:     indices             = &orientedVertices[numFaceVertices*2];
2814:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2815:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2816:     /* TODO: I know that routine should return a permutation, not the indices */
2817:     for (v = 0; v < numFaceVertices; ++v) {
2818:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2819:       for (ov = 0; ov < numFaceVertices; ++ov) {
2820:         if (orientedVertices[ov] == vertex) {
2821:           orientedSubVertices[ov] = subvertex;
2822:           break;
2823:         }
2824:       }
2825:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2826:     }
2827:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2828:     DMPlexSetCone(subdm, subcell, newFacePoint);
2829:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2830:     ++(*newFacePoint);
2831:   }
2832: #if 0
2833:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2834: #else
2835:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2836: #endif
2837:   return(0);
2838: }

2840: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2841: {
2842:   MPI_Comm        comm;
2843:   DMLabel         subpointMap;
2844:   IS              subvertexIS,  subcellIS;
2845:   const PetscInt *subVertices, *subCells;
2846:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2847:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2848:   PetscInt        vStart, vEnd, c, f;
2849:   PetscErrorCode  ierr;

2852:   PetscObjectGetComm((PetscObject)dm,&comm);
2853:   /* Create subpointMap which marks the submesh */
2854:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2855:   DMPlexSetSubpointMap(subdm, subpointMap);
2856:   DMLabelDestroy(&subpointMap);
2857:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2858:   /* Setup chart */
2859:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2860:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2861:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2862:   DMPlexSetVTKCellHeight(subdm, 1);
2863:   /* Set cone sizes */
2864:   firstSubVertex = numSubCells;
2865:   firstSubFace   = numSubCells+numSubVertices;
2866:   newFacePoint   = firstSubFace;
2867:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2868:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2869:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2870:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2871:   for (c = 0; c < numSubCells; ++c) {
2872:     DMPlexSetConeSize(subdm, c, 1);
2873:   }
2874:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2875:     DMPlexSetConeSize(subdm, f, nFV);
2876:   }
2877:   DMSetUp(subdm);
2878:   /* Create face cones */
2879:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2880:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2881:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2882:   for (c = 0; c < numSubCells; ++c) {
2883:     const PetscInt cell    = subCells[c];
2884:     const PetscInt subcell = c;
2885:     PetscInt      *closure = NULL;
2886:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2888:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2889:     for (cl = 0; cl < closureSize*2; cl += 2) {
2890:       const PetscInt point = closure[cl];
2891:       PetscInt       subVertex;

2893:       if ((point >= vStart) && (point < vEnd)) {
2894:         ++numCorners;
2895:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2896:         if (subVertex >= 0) {
2897:           closure[faceSize] = point;
2898:           subface[faceSize] = firstSubVertex+subVertex;
2899:           ++faceSize;
2900:         }
2901:       }
2902:     }
2903:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2904:     if (faceSize == nFV) {
2905:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2906:     }
2907:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2908:   }
2909:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2910:   DMPlexSymmetrize(subdm);
2911:   DMPlexStratify(subdm);
2912:   /* Build coordinates */
2913:   {
2914:     PetscSection coordSection, subCoordSection;
2915:     Vec          coordinates, subCoordinates;
2916:     PetscScalar *coords, *subCoords;
2917:     PetscInt     numComp, coordSize, v;
2918:     const char  *name;

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

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

2951:         PetscSectionGetDof(coordSection, vertex, &dof);
2952:         PetscSectionGetOffset(coordSection, vertex, &off);
2953:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2954:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2955:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2956:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2957:       }
2958:       VecRestoreArray(coordinates,    &coords);
2959:       VecRestoreArray(subCoordinates, &subCoords);
2960:     }
2961:     DMSetCoordinatesLocal(subdm, subCoordinates);
2962:     VecDestroy(&subCoordinates);
2963:   }
2964:   /* Cleanup */
2965:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2966:   ISDestroy(&subvertexIS);
2967:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2968:   ISDestroy(&subcellIS);
2969:   return(0);
2970: }

2972: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2973: {
2974:   PetscInt       subPoint;

2977:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2978:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2979: }

2981: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2982: {
2983:   MPI_Comm         comm;
2984:   DMLabel          subpointMap;
2985:   IS              *subpointIS;
2986:   const PetscInt **subpoints;
2987:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2988:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2989:   PetscMPIInt      rank;
2990:   PetscErrorCode   ierr;

2993:   PetscObjectGetComm((PetscObject)dm,&comm);
2994:   MPI_Comm_rank(comm, &rank);
2995:   /* Create subpointMap which marks the submesh */
2996:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2997:   DMPlexSetSubpointMap(subdm, subpointMap);
2998:   if (cellHeight) {
2999:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
3000:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
3001:   } else {
3002:     DMLabel         depth;
3003:     IS              pointIS;
3004:     const PetscInt *points;
3005:     PetscInt        numPoints;

3007:     DMPlexGetDepthLabel(dm, &depth);
3008:     DMLabelGetStratumSize(label, value, &numPoints);
3009:     DMLabelGetStratumIS(label, value, &pointIS);
3010:     ISGetIndices(pointIS, &points);
3011:     for (p = 0; p < numPoints; ++p) {
3012:       PetscInt *closure = NULL;
3013:       PetscInt  closureSize, c, pdim;

3015:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3016:       for (c = 0; c < closureSize*2; c += 2) {
3017:         DMLabelGetValue(depth, closure[c], &pdim);
3018:         DMLabelSetValue(subpointMap, closure[c], pdim);
3019:       }
3020:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3021:     }
3022:     ISRestoreIndices(pointIS, &points);
3023:     ISDestroy(&pointIS);
3024:   }
3025:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
3026:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3027:   cMax = (cMax < 0) ? cEnd : cMax;
3028:   /* Setup chart */
3029:   DMGetDimension(dm, &dim);
3030:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3031:   for (d = 0; d <= dim; ++d) {
3032:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3033:     totSubPoints += numSubPoints[d];
3034:   }
3035:   DMPlexSetChart(subdm, 0, totSubPoints);
3036:   DMPlexSetVTKCellHeight(subdm, cellHeight);
3037:   /* Set cone sizes */
3038:   firstSubPoint[dim] = 0;
3039:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
3040:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
3041:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3042:   for (d = 0; d <= dim; ++d) {
3043:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3044:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3045:   }
3046:   for (d = 0; d <= dim; ++d) {
3047:     for (p = 0; p < numSubPoints[d]; ++p) {
3048:       const PetscInt  point    = subpoints[d][p];
3049:       const PetscInt  subpoint = firstSubPoint[d] + p;
3050:       const PetscInt *cone;
3051:       PetscInt        coneSize, coneSizeNew, c, val;

3053:       DMPlexGetConeSize(dm, point, &coneSize);
3054:       DMPlexSetConeSize(subdm, subpoint, coneSize);
3055:       if (cellHeight && (d == dim)) {
3056:         DMPlexGetCone(dm, point, &cone);
3057:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3058:           DMLabelGetValue(subpointMap, cone[c], &val);
3059:           if (val >= 0) coneSizeNew++;
3060:         }
3061:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3062:       }
3063:     }
3064:   }
3065:   DMLabelDestroy(&subpointMap);
3066:   DMSetUp(subdm);
3067:   /* Set cones */
3068:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3069:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3070:   for (d = 0; d <= dim; ++d) {
3071:     for (p = 0; p < numSubPoints[d]; ++p) {
3072:       const PetscInt  point    = subpoints[d][p];
3073:       const PetscInt  subpoint = firstSubPoint[d] + p;
3074:       const PetscInt *cone, *ornt;
3075:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

3077:       if (d == dim-1) {
3078:         const PetscInt *support, *cone, *ornt;
3079:         PetscInt        supportSize, coneSize, s, subc;

3081:         DMPlexGetSupport(dm, point, &support);
3082:         DMPlexGetSupportSize(dm, point, &supportSize);
3083:         for (s = 0; s < supportSize; ++s) {
3084:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
3085:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3086:           if (subc >= 0) {
3087:             const PetscInt ccell = subpoints[d+1][subc];

3089:             DMPlexGetCone(dm, ccell, &cone);
3090:             DMPlexGetConeSize(dm, ccell, &coneSize);
3091:             DMPlexGetConeOrientation(dm, ccell, &ornt);
3092:             for (c = 0; c < coneSize; ++c) {
3093:               if (cone[c] == point) {
3094:                 fornt = ornt[c];
3095:                 break;
3096:               }
3097:             }
3098:             break;
3099:           }
3100:         }
3101:       }
3102:       DMPlexGetConeSize(dm, point, &coneSize);
3103:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3104:       DMPlexGetCone(dm, point, &cone);
3105:       DMPlexGetConeOrientation(dm, point, &ornt);
3106:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3107:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3108:         if (subc >= 0) {
3109:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3110:           orntNew[coneSizeNew] = ornt[c];
3111:           ++coneSizeNew;
3112:         }
3113:       }
3114:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3115:       if (fornt < 0) {
3116:         /* This should be replaced by a call to DMPlexReverseCell() */
3117: #if 0
3118:         DMPlexReverseCell(subdm, subpoint);
3119: #else
3120:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3121:           PetscInt faceSize, tmp;

3123:           tmp        = coneNew[c];
3124:           coneNew[c] = coneNew[coneSizeNew-1-c];
3125:           coneNew[coneSizeNew-1-c] = tmp;
3126:           DMPlexGetConeSize(dm, cone[c], &faceSize);
3127:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3128:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3129:           orntNew[coneSizeNew-1-c] = tmp;
3130:         }
3131:       }
3132:       DMPlexSetCone(subdm, subpoint, coneNew);
3133:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3134: #endif
3135:     }
3136:   }
3137:   PetscFree2(coneNew,orntNew);
3138:   DMPlexSymmetrize(subdm);
3139:   DMPlexStratify(subdm);
3140:   /* Build coordinates */
3141:   {
3142:     PetscSection coordSection, subCoordSection;
3143:     Vec          coordinates, subCoordinates;
3144:     PetscScalar *coords, *subCoords;
3145:     PetscInt     cdim, numComp, coordSize;
3146:     const char  *name;

3148:     DMGetCoordinateDim(dm, &cdim);
3149:     DMGetCoordinateSection(dm, &coordSection);
3150:     DMGetCoordinatesLocal(dm, &coordinates);
3151:     DMGetCoordinateSection(subdm, &subCoordSection);
3152:     PetscSectionSetNumFields(subCoordSection, 1);
3153:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3154:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3155:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3156:     for (v = 0; v < numSubPoints[0]; ++v) {
3157:       const PetscInt vertex    = subpoints[0][v];
3158:       const PetscInt subvertex = firstSubPoint[0]+v;
3159:       PetscInt       dof;

3161:       PetscSectionGetDof(coordSection, vertex, &dof);
3162:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3163:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3164:     }
3165:     PetscSectionSetUp(subCoordSection);
3166:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3167:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3168:     PetscObjectGetName((PetscObject)coordinates,&name);
3169:     PetscObjectSetName((PetscObject)subCoordinates,name);
3170:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3171:     VecSetBlockSize(subCoordinates, cdim);
3172:     VecSetType(subCoordinates,VECSTANDARD);
3173:     VecGetArray(coordinates,    &coords);
3174:     VecGetArray(subCoordinates, &subCoords);
3175:     for (v = 0; v < numSubPoints[0]; ++v) {
3176:       const PetscInt vertex    = subpoints[0][v];
3177:       const PetscInt subvertex = firstSubPoint[0]+v;
3178:       PetscInt dof, off, sdof, soff, d;

3180:       PetscSectionGetDof(coordSection, vertex, &dof);
3181:       PetscSectionGetOffset(coordSection, vertex, &off);
3182:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3183:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3184:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3185:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3186:     }
3187:     VecRestoreArray(coordinates,    &coords);
3188:     VecRestoreArray(subCoordinates, &subCoords);
3189:     DMSetCoordinatesLocal(subdm, subCoordinates);
3190:     VecDestroy(&subCoordinates);
3191:   }
3192:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3193:   {
3194:     PetscSF            sfPoint, sfPointSub;
3195:     IS                 subpIS;
3196:     const PetscSFNode *remotePoints;
3197:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3198:     const PetscInt    *localPoints, *subpoints;
3199:     PetscInt          *slocalPoints;
3200:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3201:     PetscMPIInt        rank;

3203:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3204:     DMGetPointSF(dm, &sfPoint);
3205:     DMGetPointSF(subdm, &sfPointSub);
3206:     DMPlexGetChart(dm, &pStart, &pEnd);
3207:     DMPlexGetChart(subdm, NULL, &numSubroots);
3208:     DMPlexCreateSubpointIS(subdm, &subpIS);
3209:     if (subpIS) {
3210:       ISGetIndices(subpIS, &subpoints);
3211:       ISGetLocalSize(subpIS, &numSubpoints);
3212:     }
3213:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3214:     if (numRoots >= 0) {
3215:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3216:       for (p = 0; p < pEnd-pStart; ++p) {
3217:         newLocalPoints[p].rank  = -2;
3218:         newLocalPoints[p].index = -2;
3219:       }
3220:       /* Set subleaves */
3221:       for (l = 0; l < numLeaves; ++l) {
3222:         const PetscInt point    = localPoints[l];
3223:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3225:         if (subpoint < 0) continue;
3226:         newLocalPoints[point-pStart].rank  = rank;
3227:         newLocalPoints[point-pStart].index = subpoint;
3228:         ++numSubleaves;
3229:       }
3230:       /* Must put in owned subpoints */
3231:       for (p = pStart; p < pEnd; ++p) {
3232:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3234:         if (subpoint < 0) {
3235:           newOwners[p-pStart].rank  = -3;
3236:           newOwners[p-pStart].index = -3;
3237:         } else {
3238:           newOwners[p-pStart].rank  = rank;
3239:           newOwners[p-pStart].index = subpoint;
3240:         }
3241:       }
3242:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3243:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3244:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3245:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3246:       PetscMalloc1(numSubleaves, &slocalPoints);
3247:       PetscMalloc1(numSubleaves, &sremotePoints);
3248:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3249:         const PetscInt point    = localPoints[l];
3250:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3252:         if (subpoint < 0) continue;
3253:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3254:         slocalPoints[sl]        = subpoint;
3255:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3256:         sremotePoints[sl].index = newLocalPoints[point].index;
3257:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3258:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3259:         ++sl;
3260:       }
3261:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3262:       PetscFree2(newLocalPoints,newOwners);
3263:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3264:     }
3265:     if (subpIS) {
3266:       ISRestoreIndices(subpIS, &subpoints);
3267:       ISDestroy(&subpIS);
3268:     }
3269:   }
3270:   /* Cleanup */
3271:   for (d = 0; d <= dim; ++d) {
3272:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3273:     ISDestroy(&subpointIS[d]);
3274:   }
3275:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3276:   return(0);
3277: }

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

3284:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3285:   return(0);
3286: }

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

3291:   Input Parameters:
3292: + dm           - The original mesh
3293: . vertexLabel  - The DMLabel marking points contained in the surface
3294: . value        - The label value to use
3295: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

3297:   Output Parameter:
3298: . subdm - The surface mesh

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

3302:   Level: developer

3304: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3305: @*/
3306: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3307: {
3308:   DMPlexInterpolatedFlag interpolated;
3309:   PetscInt       dim, cdim;

3315:   DMGetDimension(dm, &dim);
3316:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3317:   DMSetType(*subdm, DMPLEX);
3318:   DMSetDimension(*subdm, dim-1);
3319:   DMGetCoordinateDim(dm, &cdim);
3320:   DMSetCoordinateDim(*subdm, cdim);
3321:   DMPlexIsInterpolated(dm, &interpolated);
3322:   if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
3323:   if (interpolated) {
3324:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3325:   } else {
3326:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3327:   }
3328:   return(0);
3329: }

3331: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3332: {
3333:   MPI_Comm        comm;
3334:   DMLabel         subpointMap;
3335:   IS              subvertexIS;
3336:   const PetscInt *subVertices;
3337:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3338:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3339:   PetscInt        cMax, c, f;
3340:   PetscErrorCode  ierr;

3343:   PetscObjectGetComm((PetscObject)dm, &comm);
3344:   /* Create subpointMap which marks the submesh */
3345:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3346:   DMPlexSetSubpointMap(subdm, subpointMap);
3347:   DMLabelDestroy(&subpointMap);
3348:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3349:   /* Setup chart */
3350:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3351:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3352:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3353:   DMPlexSetVTKCellHeight(subdm, 1);
3354:   /* Set cone sizes */
3355:   firstSubVertex = numSubCells;
3356:   firstSubFace   = numSubCells+numSubVertices;
3357:   newFacePoint   = firstSubFace;
3358:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3359:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3360:   for (c = 0; c < numSubCells; ++c) {
3361:     DMPlexSetConeSize(subdm, c, 1);
3362:   }
3363:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3364:     DMPlexSetConeSize(subdm, f, nFV);
3365:   }
3366:   DMSetUp(subdm);
3367:   /* Create face cones */
3368:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3369:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3370:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3371:   for (c = 0; c < numSubCells; ++c) {
3372:     const PetscInt  cell    = subCells[c];
3373:     const PetscInt  subcell = c;
3374:     const PetscInt *cone, *cells;
3375:     PetscInt        numCells, subVertex, p, v;

3377:     if (cell < cMax) continue;
3378:     DMPlexGetCone(dm, cell, &cone);
3379:     for (v = 0; v < nFV; ++v) {
3380:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3381:       subface[v] = firstSubVertex+subVertex;
3382:     }
3383:     DMPlexSetCone(subdm, newFacePoint, subface);
3384:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3385:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3386:     /* Not true in parallel
3387:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3388:     for (p = 0; p < numCells; ++p) {
3389:       PetscInt negsubcell;

3391:       if (cells[p] >= cMax) continue;
3392:       /* I know this is a crap search */
3393:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3394:         if (subCells[negsubcell] == cells[p]) break;
3395:       }
3396:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3397:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3398:     }
3399:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3400:     ++newFacePoint;
3401:   }
3402:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3403:   DMPlexSymmetrize(subdm);
3404:   DMPlexStratify(subdm);
3405:   /* Build coordinates */
3406:   {
3407:     PetscSection coordSection, subCoordSection;
3408:     Vec          coordinates, subCoordinates;
3409:     PetscScalar *coords, *subCoords;
3410:     PetscInt     cdim, numComp, coordSize, v;
3411:     const char  *name;

3413:     DMGetCoordinateDim(dm, &cdim);
3414:     DMGetCoordinateSection(dm, &coordSection);
3415:     DMGetCoordinatesLocal(dm, &coordinates);
3416:     DMGetCoordinateSection(subdm, &subCoordSection);
3417:     PetscSectionSetNumFields(subCoordSection, 1);
3418:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3419:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3420:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3421:     for (v = 0; v < numSubVertices; ++v) {
3422:       const PetscInt vertex    = subVertices[v];
3423:       const PetscInt subvertex = firstSubVertex+v;
3424:       PetscInt       dof;

3426:       PetscSectionGetDof(coordSection, vertex, &dof);
3427:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3428:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3429:     }
3430:     PetscSectionSetUp(subCoordSection);
3431:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3432:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3433:     PetscObjectGetName((PetscObject)coordinates,&name);
3434:     PetscObjectSetName((PetscObject)subCoordinates,name);
3435:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3436:     VecSetBlockSize(subCoordinates, cdim);
3437:     VecSetType(subCoordinates,VECSTANDARD);
3438:     VecGetArray(coordinates,    &coords);
3439:     VecGetArray(subCoordinates, &subCoords);
3440:     for (v = 0; v < numSubVertices; ++v) {
3441:       const PetscInt vertex    = subVertices[v];
3442:       const PetscInt subvertex = firstSubVertex+v;
3443:       PetscInt       dof, off, sdof, soff, d;

3445:       PetscSectionGetDof(coordSection, vertex, &dof);
3446:       PetscSectionGetOffset(coordSection, vertex, &off);
3447:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3448:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3449:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3450:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3451:     }
3452:     VecRestoreArray(coordinates,    &coords);
3453:     VecRestoreArray(subCoordinates, &subCoords);
3454:     DMSetCoordinatesLocal(subdm, subCoordinates);
3455:     VecDestroy(&subCoordinates);
3456:   }
3457:   /* Build SF */
3458:   CHKMEMQ;
3459:   {
3460:     PetscSF            sfPoint, sfPointSub;
3461:     const PetscSFNode *remotePoints;
3462:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3463:     const PetscInt    *localPoints;
3464:     PetscInt          *slocalPoints;
3465:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3466:     PetscMPIInt        rank;

3468:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3469:     DMGetPointSF(dm, &sfPoint);
3470:     DMGetPointSF(subdm, &sfPointSub);
3471:     DMPlexGetChart(dm, &pStart, &pEnd);
3472:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3473:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3474:     if (numRoots >= 0) {
3475:       /* Only vertices should be shared */
3476:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3477:       for (p = 0; p < pEnd-pStart; ++p) {
3478:         newLocalPoints[p].rank  = -2;
3479:         newLocalPoints[p].index = -2;
3480:       }
3481:       /* Set subleaves */
3482:       for (l = 0; l < numLeaves; ++l) {
3483:         const PetscInt point    = localPoints[l];
3484:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3486:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3487:         if (subPoint < 0) continue;
3488:         newLocalPoints[point-pStart].rank  = rank;
3489:         newLocalPoints[point-pStart].index = subPoint;
3490:         ++numSubLeaves;
3491:       }
3492:       /* Must put in owned subpoints */
3493:       for (p = pStart; p < pEnd; ++p) {
3494:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3496:         if (subPoint < 0) {
3497:           newOwners[p-pStart].rank  = -3;
3498:           newOwners[p-pStart].index = -3;
3499:         } else {
3500:           newOwners[p-pStart].rank  = rank;
3501:           newOwners[p-pStart].index = subPoint;
3502:         }
3503:       }
3504:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3505:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3506:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3507:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3508:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3509:       PetscMalloc1(numSubLeaves, &sremotePoints);
3510:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3511:         const PetscInt point    = localPoints[l];
3512:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3514:         if (subPoint < 0) continue;
3515:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3516:         slocalPoints[sl]        = subPoint;
3517:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3518:         sremotePoints[sl].index = newLocalPoints[point].index;
3519:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3520:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3521:         ++sl;
3522:       }
3523:       PetscFree2(newLocalPoints,newOwners);
3524:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3525:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3526:     }
3527:   }
3528:   CHKMEMQ;
3529:   /* Cleanup */
3530:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3531:   ISDestroy(&subvertexIS);
3532:   PetscFree(subCells);
3533:   return(0);
3534: }

3536: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3537: {
3538:   DMLabel        label = NULL;

3542:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3543:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3544:   return(0);
3545: }

3547: /*@C
3548:   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.

3550:   Input Parameters:
3551: + dm          - The original mesh
3552: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3553: . label       - A label name, or NULL
3554: - value  - A label value

3556:   Output Parameter:
3557: . subdm - The surface mesh

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

3561:   Level: developer

3563: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3564: @*/
3565: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3566: {
3567:   PetscInt       dim, cdim, depth;

3573:   DMGetDimension(dm, &dim);
3574:   DMPlexGetDepth(dm, &depth);
3575:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3576:   DMSetType(*subdm, DMPLEX);
3577:   DMSetDimension(*subdm, dim-1);
3578:   DMGetCoordinateDim(dm, &cdim);
3579:   DMSetCoordinateDim(*subdm, cdim);
3580:   if (depth == dim) {
3581:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3582:   } else {
3583:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3584:   }
3585:   return(0);
3586: }

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

3591:   Input Parameters:
3592: + dm        - The original mesh
3593: . cellLabel - The DMLabel marking cells contained in the new mesh
3594: - value     - The label value to use

3596:   Output Parameter:
3597: . subdm - The new mesh

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

3601:   Level: developer

3603: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3604: @*/
3605: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3606: {
3607:   PetscInt       dim;

3613:   DMGetDimension(dm, &dim);
3614:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3615:   DMSetType(*subdm, DMPLEX);
3616:   DMSetDimension(*subdm, dim);
3617:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3618:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3619:   return(0);
3620: }

3622: /*@
3623:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3625:   Input Parameter:
3626: . dm - The submesh DM

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

3631:   Level: developer

3633: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3634: @*/
3635: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3636: {
3640:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3641:   return(0);
3642: }

3644: /*@
3645:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3647:   Input Parameters:
3648: + dm - The submesh DM
3649: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3653:   Level: developer

3655: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3656: @*/
3657: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3658: {
3659:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3660:   DMLabel        tmp;

3665:   tmp  = mesh->subpointMap;
3666:   mesh->subpointMap = subpointMap;
3667:   PetscObjectReference((PetscObject) mesh->subpointMap);
3668:   DMLabelDestroy(&tmp);
3669:   return(0);
3670: }

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

3675:   Input Parameter:
3676: . dm - The submesh DM

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

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

3683:   Level: developer

3685: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3686: @*/
3687: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3688: {
3689:   MPI_Comm        comm;
3690:   DMLabel         subpointMap;
3691:   IS              is;
3692:   const PetscInt *opoints;
3693:   PetscInt       *points, *depths;
3694:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3695:   PetscErrorCode  ierr;

3700:   PetscObjectGetComm((PetscObject)dm,&comm);
3701:   *subpointIS = NULL;
3702:   DMPlexGetSubpointMap(dm, &subpointMap);
3703:   DMPlexGetDepth(dm, &depth);
3704:   if (subpointMap && depth >= 0) {
3705:     DMPlexGetChart(dm, &pStart, &pEnd);
3706:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3707:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3708:     depths[0] = depth;
3709:     depths[1] = 0;
3710:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3711:     PetscMalloc1(pEnd, &points);
3712:     for(d = 0, off = 0; d <= depth; ++d) {
3713:       const PetscInt dep = depths[d];

3715:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3716:       DMLabelGetStratumSize(subpointMap, dep, &n);
3717:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3718:         if (n != depEnd-depStart) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
3719:       } else {
3720:         if (!n) {
3721:           if (d == 0) {
3722:             /* Missing cells */
3723:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3724:           } else {
3725:             /* Missing faces */
3726:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3727:           }
3728:         }
3729:       }
3730:       if (n) {
3731:         DMLabelGetStratumIS(subpointMap, dep, &is);
3732:         ISGetIndices(is, &opoints);
3733:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3734:         ISRestoreIndices(is, &opoints);
3735:         ISDestroy(&is);
3736:       }
3737:     }
3738:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3739:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3740:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3741:   }
3742:   return(0);
3743: }

3745: /*@
3746:   DMGetEnclosureRelation - Get the relationship between dmA and dmB

3748:   Input Parameters:
3749: + dmA - The first DM
3750: - dmB - The second DM

3752:   Output Parameter:
3753: . rel - The relation of dmA to dmB

3755:   Level: intermediate

3757: .seealso: DMPlexGetEnclosurePoint()
3758: @*/
3759: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
3760: {
3761:   DM             plexA, plexB, sdm;
3762:   DMLabel        spmap;
3763:   PetscInt       pStartA, pEndA, pStartB, pEndB, NpA, NpB;

3768:   *rel = DM_ENC_NONE;
3769:   if (!dmA || !dmB) return(0);
3772:   if (dmA == dmB) {*rel = DM_ENC_EQUALITY; return(0);}
3773:   DMConvert(dmA, DMPLEX, &plexA);
3774:   DMConvert(dmB, DMPLEX, &plexB);
3775:   DMPlexGetChart(plexA, &pStartA, &pEndA);
3776:   DMPlexGetChart(plexB, &pStartB, &pEndB);
3777:   /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
3778:     - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
3779:   if ((pStartA == pStartB) && (pEndA == pEndB)) {
3780:     *rel = DM_ENC_EQUALITY;
3781:     goto end;
3782:   }
3783:   NpA = pEndA - pStartA;
3784:   NpB = pEndB - pStartB;
3785:   if (NpA == NpB) goto end;
3786:   sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
3787:   DMPlexGetSubpointMap(sdm, &spmap);
3788:   if (!spmap) goto end;
3789:   /* TODO Check the space mapped to by subpointMap is same size as dm */
3790:   if (NpA > NpB) {
3791:     *rel = DM_ENC_SUPERMESH;
3792:   } else {
3793:     *rel = DM_ENC_SUBMESH;
3794:   }
3795:   end:
3796:   DMDestroy(&plexA);
3797:   DMDestroy(&plexB);
3798:   return(0);
3799: }

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

3804:   Input Parameters:
3805: + dmA   - The first DM
3806: . dmB   - The second DM
3807: . etype - The type of enclosure relation that dmA has to dmB
3808: - pB    - A point of dmB

3810:   Output Parameter:
3811: . pA    - The corresponding point of dmA

3813:   Level: intermediate

3815: .seealso: DMGetEnclosureRelation()
3816: @*/
3817: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
3818: {
3819:   DM              sdm;
3820:   DMLabel         spmap;
3821:   IS              subpointIS;
3822:   const PetscInt *subpoints;
3823:   PetscInt        numSubpoints;
3824:   PetscErrorCode  ierr;

3827:   /* TODO Cache the IS, making it look like an index */
3828:   switch (etype) {
3829:     case DM_ENC_SUPERMESH:
3830:     sdm  = dmB;
3831:     DMPlexGetSubpointMap(sdm, &spmap);
3832:     DMPlexCreateSubpointIS(sdm, &subpointIS);
3833:     ISGetIndices(subpointIS, &subpoints);
3834:     *pA  = subpoints[pB];
3835:     ISRestoreIndices(subpointIS, &subpoints);
3836:     ISDestroy(&subpointIS);
3837:     break;
3838:     case DM_ENC_SUBMESH:
3839:     sdm  = dmA;
3840:     DMPlexGetSubpointMap(sdm, &spmap);
3841:     DMPlexCreateSubpointIS(sdm, &subpointIS);
3842:     ISGetLocalSize(subpointIS, &numSubpoints);
3843:     ISGetIndices(subpointIS, &subpoints);
3844:     PetscFindInt(pB, numSubpoints, subpoints, pA);
3845:     if (*pA < 0) {
3846:       DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
3847:       DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
3848:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB);
3849:     }
3850:     ISRestoreIndices(subpointIS, &subpoints);
3851:     ISDestroy(&subpointIS);
3852:     break;
3853:     case DM_ENC_EQUALITY:
3854:     case DM_ENC_NONE:
3855:     *pA = pB;break;
3856:     case DM_ENC_UNKNOWN:
3857:     {
3858:       DMEnclosureType enc;

3860:       DMGetEnclosureRelation(dmA, dmB, &enc);
3861:       DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
3862:     }
3863:     break;
3864:     default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype);
3865:   }
3866:   return(0);
3867: }