Actual source code: plexsubmesh.c

petsc-master 2016-02-06
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>    /*I      "petscdmplex.h"    I*/
  2: #include <petsc/private/dmlabelimpl.h>   /*I      "petscdmlabel.h"   I*/
  3: #include <petscsf.h>

  7: PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label)
  8: {
  9:   PetscInt       fStart, fEnd, f;

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

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

 26: /*@
 27:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 29:   Not Collective

 31:   Input Parameter:
 32: . dm - The original DM

 34:   Output Parameter:
 35: . label - The DMLabel marking boundary faces with value 1

 37:   Level: developer

 39: .seealso: DMLabelCreate(), DMCreateLabel()
 40: @*/
 41: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
 42: {

 46:   DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
 47:   return(0);
 48: }

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

 55:   Input Parameters:
 56: + dm - The DM
 57: - label - A DMLabel marking the surface points

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

 62:   Level: developer

 64: .seealso: DMPlexLabelCohesiveComplete()
 65: @*/
 66: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
 67: {
 68:   IS              valueIS;
 69:   const PetscInt *values;
 70:   PetscInt        numValues, v;
 71:   PetscErrorCode  ierr;

 74:   DMLabelGetNumValues(label, &numValues);
 75:   DMLabelGetValueIS(label, &valueIS);
 76:   ISGetIndices(valueIS, &values);
 77:   for (v = 0; v < numValues; ++v) {
 78:     IS              pointIS;
 79:     const PetscInt *points;
 80:     PetscInt        numPoints, p;

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

 89:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 90:       for (c = 0; c < closureSize*2; c += 2) {
 91:         DMLabelSetValue(label, closure[c], values[v]);
 92:       }
 93:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 94:     }
 95:     ISRestoreIndices(pointIS, &points);
 96:     ISDestroy(&pointIS);
 97:   }
 98:   ISRestoreIndices(valueIS, &values);
 99:   ISDestroy(&valueIS);
100:   return(0);
101: }

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

108:   Input Parameters:
109: + dm - The DM
110: - label - A DMLabel marking the surface points

112:   Output Parameter:
113: . label - A DMLabel incorporating cells

115:   Level: developer

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

119: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
120: @*/
121: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
122: {
123:   IS              valueIS;
124:   const PetscInt *values;
125:   PetscInt        numValues, v, cStart, cEnd, cEndInterior;
126:   PetscErrorCode  ierr;

129:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
130:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
131:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
132:   DMLabelGetNumValues(label, &numValues);
133:   DMLabelGetValueIS(label, &valueIS);
134:   ISGetIndices(valueIS, &values);
135:   for (v = 0; v < numValues; ++v) {
136:     IS              pointIS;
137:     const PetscInt *points;
138:     PetscInt        numPoints, p;

140:     DMLabelGetStratumSize(label, values[v], &numPoints);
141:     DMLabelGetStratumIS(label, values[v], &pointIS);
142:     ISGetIndices(pointIS, &points);
143:     for (p = 0; p < numPoints; ++p) {
144:       PetscInt *closure = NULL;
145:       PetscInt  closureSize, point, cl;

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

164: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
165:  * index (skipping first, which is (0,0)) */
166: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
167: {
168:   PetscInt d, off = 0;

171:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
172:   for (d = 0; d < depth; d++) {
173:     PetscInt firstd = d;
174:     PetscInt firstStart = depthShift[2*d];
175:     PetscInt e;

177:     for (e = d+1; e <= depth; e++) {
178:       if (depthShift[2*e] < firstStart) {
179:         firstd = e;
180:         firstStart = depthShift[2*d];
181:       }
182:     }
183:     if (firstd != d) {
184:       PetscInt swap[2];

186:       e = firstd;
187:       swap[0] = depthShift[2*d];
188:       swap[1] = depthShift[2*d+1];
189:       depthShift[2*d]   = depthShift[2*e];
190:       depthShift[2*d+1] = depthShift[2*e+1];
191:       depthShift[2*e]   = swap[0];
192:       depthShift[2*e+1] = swap[1];
193:     }
194:   }
195:   /* convert (oldstart, added) to (oldstart, newstart) */
196:   for (d = 0; d <= depth; d++) {
197:     off += depthShift[2*d+1];
198:     depthShift[2*d+1] = depthShift[2*d] + off;
199:   }
200:   return(0);
201: }

205: /* depthShift is a list of (old, new) pairs */
206: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
207: {
208:   PetscInt d;
209:   PetscInt newOff = 0;

211:   for (d = 0; d <= depth; d++) {
212:     if (p < depthShift[2*d]) return p + newOff;
213:     else newOff = depthShift[2*d+1] - depthShift[2*d];
214:   }
215:   return p + newOff;
216: }

220: /* depthShift is a list of (old, new) pairs */
221: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
222: {
223:   PetscInt d;
224:   PetscInt newOff = 0;

226:   for (d = 0; d <= depth; d++) {
227:     if (p < depthShift[2*d+1]) return p + newOff;
228:     else newOff = depthShift[2*d] - depthShift[2*d+1];
229:   }
230:   return p + newOff;
231: }

235: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
236: {
237:   PetscInt       depth = 0, d, pStart, pEnd, p;

241:   DMPlexGetDepth(dm, &depth);
242:   if (depth < 0) return(0);
243:   /* Step 1: Expand chart */
244:   DMPlexGetChart(dm, &pStart, &pEnd);
245:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
246:   DMPlexSetChart(dmNew, pStart, pEnd);
247:   /* Step 2: Set cone and support sizes */
248:   for (d = 0; d <= depth; ++d) {
249:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
250:     for (p = pStart; p < pEnd; ++p) {
251:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
252:       PetscInt size;

254:       DMPlexGetConeSize(dm, p, &size);
255:       DMPlexSetConeSize(dmNew, newp, size);
256:       DMPlexGetSupportSize(dm, p, &size);
257:       DMPlexSetSupportSize(dmNew, newp, size);
258:     }
259:   }
260:   return(0);
261: }

265: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
266: {
267:   PetscInt      *newpoints;
268:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

272:   DMPlexGetDepth(dm, &depth);
273:   if (depth < 0) return(0);
274:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
275:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
276:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
277:   /* Step 5: Set cones and supports */
278:   DMPlexGetChart(dm, &pStart, &pEnd);
279:   for (p = pStart; p < pEnd; ++p) {
280:     const PetscInt *points = NULL, *orientations = NULL;
281:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

283:     DMPlexGetConeSize(dm, p, &size);
284:     DMPlexGetCone(dm, p, &points);
285:     DMPlexGetConeOrientation(dm, p, &orientations);
286:     for (i = 0; i < size; ++i) {
287:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
288:     }
289:     DMPlexSetCone(dmNew, newp, newpoints);
290:     DMPlexSetConeOrientation(dmNew, newp, orientations);
291:     DMPlexGetSupportSize(dm, p, &size);
292:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
293:     DMPlexGetSupport(dm, p, &points);
294:     for (i = 0; i < size; ++i) {
295:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
296:     }
297:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
298:     DMPlexSetSupport(dmNew, newp, newpoints);
299:   }
300:   PetscFree(newpoints);
301:   return(0);
302: }

306: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
307: {
308:   PetscSection   coordSection, newCoordSection;
309:   Vec            coordinates, newCoordinates;
310:   PetscScalar   *coords, *newCoords;
311:   PetscInt       coordSize;
312:   PetscInt       dim, depth = 0, vStart, vEnd, vStartNew, vEndNew, v;

316:   DMGetDimension(dm, &dim);
317:   DMPlexGetDepth(dm, &depth);
318:   /* Step 8: Convert coordinates */
319:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
320:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
321:   DMGetCoordinateSection(dm, &coordSection);
322:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
323:   PetscSectionSetNumFields(newCoordSection, 1);
324:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
325:   PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);
326:   for (v = vStartNew; v < vEndNew; ++v) {
327:     PetscSectionSetDof(newCoordSection, v, dim);
328:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
329:   }
330:   PetscSectionSetUp(newCoordSection);
331:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
332:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
333:   VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);
334:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
335:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
336:   VecSetBlockSize(newCoordinates, dim);
337:   VecSetType(newCoordinates,VECSTANDARD);
338:   DMSetCoordinatesLocal(dmNew, newCoordinates);
339:   DMGetCoordinatesLocal(dm, &coordinates);
340:   VecGetArray(coordinates, &coords);
341:   VecGetArray(newCoordinates, &newCoords);
342:   for (v = vStart; v < vEnd; ++v) {
343:     PetscInt dof, off, noff, d;

345:     PetscSectionGetDof(coordSection, v, &dof);
346:     PetscSectionGetOffset(coordSection, v, &off);
347:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
348:     for (d = 0; d < dof; ++d) {
349:       newCoords[noff+d] = coords[off+d];
350:     }
351:   }
352:   VecRestoreArray(coordinates, &coords);
353:   VecRestoreArray(newCoordinates, &newCoords);
354:   VecDestroy(&newCoordinates);
355:   PetscSectionDestroy(&newCoordSection);
356:   return(0);
357: }

361: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
362: {
363:   PetscInt           depth = 0;
364:   PetscSF            sfPoint, sfPointNew;
365:   const PetscSFNode *remotePoints;
366:   PetscSFNode       *gremotePoints;
367:   const PetscInt    *localPoints;
368:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
369:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
370:   PetscErrorCode     ierr;

373:   DMPlexGetDepth(dm, &depth);
374:   /* Step 9: Convert pointSF */
375:   DMGetPointSF(dm, &sfPoint);
376:   DMGetPointSF(dmNew, &sfPointNew);
377:   DMPlexGetChart(dm, &pStart, &pEnd);
378:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
379:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
380:   if (numRoots >= 0) {
381:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
382:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
383:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
384:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
385:     PetscMalloc1(numLeaves,    &glocalPoints);
386:     PetscMalloc1(numLeaves, &gremotePoints);
387:     for (l = 0; l < numLeaves; ++l) {
388:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
389:       gremotePoints[l].rank  = remotePoints[l].rank;
390:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
391:     }
392:     PetscFree2(newLocation,newRemoteLocation);
393:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
394:   }
395:   return(0);
396: }

400: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
401: {
402:   PetscSF            sfPoint;
403:   DMLabel            vtkLabel, ghostLabel;
404:   const PetscSFNode *leafRemote;
405:   const PetscInt    *leafLocal;
406:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
407:   PetscMPIInt        rank;
408:   PetscErrorCode     ierr;

411:   DMPlexGetDepth(dm, &depth);
412:   /* Step 10: Convert labels */
413:   DMGetNumLabels(dm, &numLabels);
414:   for (l = 0; l < numLabels; ++l) {
415:     DMLabel         label, newlabel;
416:     const char     *lname;
417:     PetscBool       isDepth;
418:     IS              valueIS;
419:     const PetscInt *values;
420:     PetscInt        numValues, val;

422:     DMGetLabelName(dm, l, &lname);
423:     PetscStrcmp(lname, "depth", &isDepth);
424:     if (isDepth) continue;
425:     DMCreateLabel(dmNew, lname);
426:     DMGetLabel(dm, lname, &label);
427:     DMGetLabel(dmNew, lname, &newlabel);
428:     DMLabelGetValueIS(label, &valueIS);
429:     ISGetLocalSize(valueIS, &numValues);
430:     ISGetIndices(valueIS, &values);
431:     for (val = 0; val < numValues; ++val) {
432:       IS              pointIS;
433:       const PetscInt *points;
434:       PetscInt        numPoints, p;

436:       DMLabelGetStratumIS(label, values[val], &pointIS);
437:       ISGetLocalSize(pointIS, &numPoints);
438:       ISGetIndices(pointIS, &points);
439:       for (p = 0; p < numPoints; ++p) {
440:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

442:         DMLabelSetValue(newlabel, newpoint, values[val]);
443:       }
444:       ISRestoreIndices(pointIS, &points);
445:       ISDestroy(&pointIS);
446:     }
447:     ISRestoreIndices(valueIS, &values);
448:     ISDestroy(&valueIS);
449:   }
450:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
451:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
452:   DMGetPointSF(dm, &sfPoint);
453:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
454:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
455:   DMCreateLabel(dmNew, "vtk");
456:   DMCreateLabel(dmNew, "ghost");
457:   DMGetLabel(dmNew, "vtk", &vtkLabel);
458:   DMGetLabel(dmNew, "ghost", &ghostLabel);
459:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
460:     for (; c < leafLocal[l] && c < cEnd; ++c) {
461:       DMLabelSetValue(vtkLabel, c, 1);
462:     }
463:     if (leafLocal[l] >= cEnd) break;
464:     if (leafRemote[l].rank == rank) {
465:       DMLabelSetValue(vtkLabel, c, 1);
466:     } else {
467:       DMLabelSetValue(ghostLabel, c, 2);
468:     }
469:   }
470:   for (; c < cEnd; ++c) {
471:     DMLabelSetValue(vtkLabel, c, 1);
472:   }
473:   if (0) {
474:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
475:   }
476:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
477:   for (f = fStart; f < fEnd; ++f) {
478:     PetscInt numCells;

480:     DMPlexGetSupportSize(dmNew, f, &numCells);
481:     if (numCells < 2) {
482:       DMLabelSetValue(ghostLabel, f, 1);
483:     } else {
484:       const PetscInt *cells = NULL;
485:       PetscInt        vA, vB;

487:       DMPlexGetSupport(dmNew, f, &cells);
488:       DMLabelGetValue(vtkLabel, cells[0], &vA);
489:       DMLabelGetValue(vtkLabel, cells[1], &vB);
490:       if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
491:     }
492:   }
493:   if (0) {
494:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
495:   }
496:   return(0);
497: }

501: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
502: {
503:   DM             refTree;
504:   PetscSection   pSec;
505:   PetscInt       *parents, *childIDs;

509:   DMPlexGetReferenceTree(dm,&refTree);
510:   DMPlexSetReferenceTree(dmNew,refTree);
511:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
512:   if (pSec) {
513:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
514:     PetscInt *childIDsShifted;
515:     PetscSection pSecShifted;

517:     PetscSectionGetChart(pSec,&pStart,&pEnd);
518:     DMPlexGetDepth(dm,&depth);
519:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
520:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
521:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
522:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
523:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
524:     for (p = pStartShifted; p < pEndShifted; p++) {
525:       /* start off assuming no children */
526:       PetscSectionSetDof(pSecShifted,p,0);
527:     }
528:     for (p = pStart; p < pEnd; p++) {
529:       PetscInt dof;
530:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

532:       PetscSectionGetDof(pSec,p,&dof);
533:       PetscSectionSetDof(pSecShifted,pNew,dof);
534:     }
535:     PetscSectionSetUp(pSecShifted);
536:     for (p = pStart; p < pEnd; p++) {
537:       PetscInt dof;
538:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

540:       PetscSectionGetDof(pSec,p,&dof);
541:       if (dof) {
542:         PetscInt off, offNew;

544:         PetscSectionGetOffset(pSec,p,&off);
545:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
546:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
547:         childIDsShifted[offNew] = childIDs[off];
548:       }
549:     }
550:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
551:     PetscFree2(parentsShifted,childIDsShifted);
552:     PetscSectionDestroy(&pSecShifted);
553:   }
554:   return(0);
555: }

559: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
560: {
561:   PetscSF         sf;
562:   IS              valueIS;
563:   const PetscInt *values, *leaves;
564:   PetscInt       *depthShift;
565:   PetscInt        d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
566:   PetscErrorCode  ierr;

569:   DMGetPointSF(dm, &sf);
570:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
571:   nleaves = PetscMax(0, nleaves);
572:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
573:   /* Count ghost cells */
574:   DMLabelGetValueIS(label, &valueIS);
575:   ISGetLocalSize(valueIS, &numFS);
576:   ISGetIndices(valueIS, &values);
577:   Ng   = 0;
578:   for (fs = 0; fs < numFS; ++fs) {
579:     IS              faceIS;
580:     const PetscInt *faces;
581:     PetscInt        numFaces, f, numBdFaces = 0;

583:     DMLabelGetStratumIS(label, values[fs], &faceIS);
584:     ISGetLocalSize(faceIS, &numFaces);
585:     ISGetIndices(faceIS, &faces);
586:     for (f = 0; f < numFaces; ++f) {
587:       PetscInt numChildren;

589:       PetscFindInt(faces[f], nleaves, leaves, &loc);
590:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
591:       /* non-local and ancestors points don't get to register ghosts */
592:       if (loc >= 0 || numChildren) continue;
593:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
594:     }
595:     Ng += numBdFaces;
596:     ISDestroy(&faceIS);
597:   }
598:   DMPlexGetDepth(dm, &depth);
599:   PetscMalloc1(2*(depth+1), &depthShift);
600:   for (d = 0; d <= depth; d++) {
601:     PetscInt dEnd;

603:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
604:     depthShift[2*d]   = dEnd;
605:     depthShift[2*d+1] = 0;
606:   }
607:   if (depth >= 0) depthShift[2*depth+1] = Ng;
608:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
609:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
610:   /* Step 3: Set cone/support sizes for new points */
611:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
612:   DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
613:   for (c = cEnd; c < cEnd + Ng; ++c) {
614:     DMPlexSetConeSize(gdm, c, 1);
615:   }
616:   for (fs = 0; fs < numFS; ++fs) {
617:     IS              faceIS;
618:     const PetscInt *faces;
619:     PetscInt        numFaces, f;

621:     DMLabelGetStratumIS(label, values[fs], &faceIS);
622:     ISGetLocalSize(faceIS, &numFaces);
623:     ISGetIndices(faceIS, &faces);
624:     for (f = 0; f < numFaces; ++f) {
625:       PetscInt size, numChildren;

627:       PetscFindInt(faces[f], nleaves, leaves, &loc);
628:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
629:       if (loc >= 0 || numChildren) continue;
630:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
631:       DMPlexGetSupportSize(dm, faces[f], &size);
632:       if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
633:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
634:     }
635:     ISRestoreIndices(faceIS, &faces);
636:     ISDestroy(&faceIS);
637:   }
638:   /* Step 4: Setup ghosted DM */
639:   DMSetUp(gdm);
640:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
641:   /* Step 6: Set cones and supports for new points */
642:   ghostCell = cEnd;
643:   for (fs = 0; fs < numFS; ++fs) {
644:     IS              faceIS;
645:     const PetscInt *faces;
646:     PetscInt        numFaces, f;

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

654:       PetscFindInt(faces[f], nleaves, leaves, &loc);
655:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
656:       if (loc >= 0 || numChildren) continue;
657:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
658:       DMPlexSetCone(gdm, ghostCell, &newFace);
659:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
660:       ++ghostCell;
661:     }
662:     ISRestoreIndices(faceIS, &faces);
663:     ISDestroy(&faceIS);
664:   }
665:   ISRestoreIndices(valueIS, &values);
666:   ISDestroy(&valueIS);
667:   /* Step 7: Stratify */
668:   DMPlexStratify(gdm);
669:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
670:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
671:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
672:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
673:   PetscFree(depthShift);
674:   /* Step 7: Periodicity */
675:   if (dm->maxCell) {
676:     const PetscReal *maxCell, *L;
677:     const DMBoundaryType *bd;
678:     DMGetPeriodicity(dm,  &maxCell, &L, &bd);
679:     DMSetPeriodicity(gdm,  maxCell,  L,  bd);
680:   }
681:   if (numGhostCells) *numGhostCells = Ng;
682:   return(0);
683: }

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

690:   Collective on dm

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

696:   Output Parameters:
697: + numGhostCells - The number of ghost cells added to the DM
698: - dmGhosted - The new DM

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

702:   Level: developer

704: .seealso: DMCreate()
705: @*/
706: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
707: {
708:   DM             gdm;
709:   DMLabel        label;
710:   const char    *name = labelName ? labelName : "Face Sets";
711:   PetscInt       dim;
712:   PetscBool      flag;

719:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
720:   DMSetType(gdm, DMPLEX);
721:   DMGetDimension(dm, &dim);
722:   DMSetDimension(gdm, dim);
723:   DMPlexGetAdjacencyUseCone(dm, &flag);
724:   DMPlexSetAdjacencyUseCone(gdm, flag);
725:   DMPlexGetAdjacencyUseClosure(dm, &flag);
726:   DMPlexSetAdjacencyUseClosure(gdm, flag);
727:   DMGetLabel(dm, name, &label);
728:   if (!label) {
729:     /* Get label for boundary faces */
730:     DMCreateLabel(dm, name);
731:     DMGetLabel(dm, name, &label);
732:     DMPlexMarkBoundaryFaces(dm, label);
733:   }
734:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
735:   DMPlexCopyBoundary(dm, gdm);
736:   *dmGhosted = gdm;
737:   return(0);
738: }

742: /*
743:   We are adding three kinds of points here:
744:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
745:     Non-replicated: Points which exist on the fault, but are not replicated
746:     Hybrid:         Entirely new points, such as cohesive cells

748:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
749:   each depth so that the new split/hybrid points can be inserted as a block.
750: */
751: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
752: {
753:   MPI_Comm         comm;
754:   IS               valueIS;
755:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
756:   const PetscInt  *values;          /* List of depths for which we have replicated points */
757:   IS              *splitIS;
758:   IS              *unsplitIS;
759:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
760:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
761:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
762:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
763:   const PetscInt **splitPoints;        /* Replicated points for each depth */
764:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
765:   PetscSection     coordSection;
766:   Vec              coordinates;
767:   PetscScalar     *coords;
768:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
769:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
770:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
771:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
772:   PetscInt        *coneNew, *coneONew, *supportNew;
773:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
774:   PetscErrorCode   ierr;

777:   PetscObjectGetComm((PetscObject)dm,&comm);
778:   DMGetDimension(dm, &dim);
779:   DMPlexGetDepth(dm, &depth);
780:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
781:   /* Count split points and add cohesive cells */
782:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
783:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
784:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
785:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
786:   for (d = 0; d <= depth; ++d) {
787:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
788:     depthEnd[d]           = pMaxNew[d];
789:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
790:     numSplitPoints[d]     = 0;
791:     numUnsplitPoints[d]   = 0;
792:     numHybridPoints[d]    = 0;
793:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
794:     splitPoints[d]        = NULL;
795:     unsplitPoints[d]      = NULL;
796:     splitIS[d]            = NULL;
797:     unsplitIS[d]          = NULL;
798:     /* we are shifting the existing hybrid points with the stratum behind them, so
799:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
800:     depthShift[2*d]       = depthMax[d];
801:     depthShift[2*d+1]     = 0;
802:   }
803:   if (label) {
804:     DMLabelGetValueIS(label, &valueIS);
805:     ISGetLocalSize(valueIS, &numSP);
806:     ISGetIndices(valueIS, &values);
807:   }
808:   for (sp = 0; sp < numSP; ++sp) {
809:     const PetscInt dep = values[sp];

811:     if ((dep < 0) || (dep > depth)) continue;
812:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
813:     if (splitIS[dep]) {
814:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
815:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
816:     }
817:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
818:     if (unsplitIS[dep]) {
819:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
820:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
821:     }
822:   }
823:   /* Calculate number of hybrid points */
824:   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   */
825:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
826:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
827:   /* the end of the points in this stratum that come before the new points:
828:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
829:    * added points */
830:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
831:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
832:   /* Step 3: Set cone/support sizes for new points */
833:   for (dep = 0; dep <= depth; ++dep) {
834:     for (p = 0; p < numSplitPoints[dep]; ++p) {
835:       const PetscInt  oldp   = splitPoints[dep][p];
836:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
837:       const PetscInt  splitp = p    + pMaxNew[dep];
838:       const PetscInt *support;
839:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

841:       DMPlexGetConeSize(dm, oldp, &coneSize);
842:       DMPlexSetConeSize(sdm, splitp, coneSize);
843:       DMPlexGetSupportSize(dm, oldp, &supportSize);
844:       DMPlexSetSupportSize(sdm, splitp, supportSize);
845:       if (dep == depth-1) {
846:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

853:         DMPlexGetSupport(dm, oldp, &support);
854:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
855:           PetscInt val;

857:           DMLabelGetValue(label, support[e], &val);
858:           if (val == 1) ++qf;
859:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
860:           if ((val == 1) || (val == -(shift + 1))) ++qp;
861:         }
862:         /* Split old vertex: Edges into original vertex and new cohesive edge */
863:         DMPlexSetSupportSize(sdm, newp, qn+1);
864:         /* Split new vertex: Edges into split vertex and new cohesive edge */
865:         DMPlexSetSupportSize(sdm, splitp, qp+1);
866:         /* Add hybrid edge */
867:         DMPlexSetConeSize(sdm, hybedge, 2);
868:         DMPlexSetSupportSize(sdm, hybedge, qf);
869:       } else if (dep == dim-2) {
870:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

872:         DMPlexGetSupport(dm, oldp, &support);
873:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
874:           PetscInt val;

876:           DMLabelGetValue(label, support[e], &val);
877:           if (val == dim-1) ++qf;
878:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
879:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
880:         }
881:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
882:         DMPlexSetSupportSize(sdm, newp, qn+1);
883:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
884:         DMPlexSetSupportSize(sdm, splitp, qp+1);
885:         /* Add hybrid face */
886:         DMPlexSetConeSize(sdm, hybface, 4);
887:         DMPlexSetSupportSize(sdm, hybface, qf);
888:       }
889:     }
890:   }
891:   for (dep = 0; dep <= depth; ++dep) {
892:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
893:       const PetscInt  oldp   = unsplitPoints[dep][p];
894:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
895:       const PetscInt *support;
896:       PetscInt        coneSize, supportSize, qf, e, s;

898:       DMPlexGetConeSize(dm, oldp, &coneSize);
899:       DMPlexGetSupportSize(dm, oldp, &supportSize);
900:       DMPlexGetSupport(dm, oldp, &support);
901:       if (dep == 0) {
902:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

904:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
905:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
906:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
907:           if (e >= 0) ++qf;
908:         }
909:         DMPlexSetSupportSize(sdm, newp, qf+2);
910:         /* Add hybrid edge */
911:         DMPlexSetConeSize(sdm, hybedge, 2);
912:         for (e = 0, qf = 0; e < supportSize; ++e) {
913:           PetscInt val;

915:           DMLabelGetValue(label, support[e], &val);
916:           /* Split and unsplit edges produce hybrid faces */
917:           if (val == 1) ++qf;
918:           if (val == (shift2 + 1)) ++qf;
919:         }
920:         DMPlexSetSupportSize(sdm, hybedge, qf);
921:       } else if (dep == dim-2) {
922:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
923:         PetscInt       val;

925:         for (e = 0, qf = 0; e < supportSize; ++e) {
926:           DMLabelGetValue(label, support[e], &val);
927:           if (val == dim-1) qf += 2;
928:           else              ++qf;
929:         }
930:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
931:         DMPlexSetSupportSize(sdm, newp, qf+2);
932:         /* Add hybrid face */
933:         for (e = 0, qf = 0; e < supportSize; ++e) {
934:           DMLabelGetValue(label, support[e], &val);
935:           if (val == dim-1) ++qf;
936:         }
937:         DMPlexSetConeSize(sdm, hybface, 4);
938:         DMPlexSetSupportSize(sdm, hybface, qf);
939:       }
940:     }
941:   }
942:   /* Step 4: Setup split DM */
943:   DMSetUp(sdm);
944:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
945:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
946:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
947:   /* Step 6: Set cones and supports for new points */
948:   for (dep = 0; dep <= depth; ++dep) {
949:     for (p = 0; p < numSplitPoints[dep]; ++p) {
950:       const PetscInt  oldp   = splitPoints[dep][p];
951:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
952:       const PetscInt  splitp = p    + pMaxNew[dep];
953:       const PetscInt *cone, *support, *ornt;
954:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

956:       DMPlexGetConeSize(dm, oldp, &coneSize);
957:       DMPlexGetCone(dm, oldp, &cone);
958:       DMPlexGetConeOrientation(dm, oldp, &ornt);
959:       DMPlexGetSupportSize(dm, oldp, &supportSize);
960:       DMPlexGetSupport(dm, oldp, &support);
961:       if (dep == depth-1) {
962:         PetscBool       hasUnsplit = PETSC_FALSE;
963:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
964:         const PetscInt *supportF;

966:         /* Split face:       copy in old face to new face to start */
967:         DMPlexGetSupport(sdm, newp,  &supportF);
968:         DMPlexSetSupport(sdm, splitp, supportF);
969:         /* Split old face:   old vertices/edges in cone so no change */
970:         /* Split new face:   new vertices/edges in cone */
971:         for (q = 0; q < coneSize; ++q) {
972:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
973:           if (v < 0) {
974:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
975:             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);
976:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
977:             hasUnsplit   = PETSC_TRUE;
978:           } else {
979:             coneNew[2+q] = v + pMaxNew[dep-1];
980:             if (dep > 1) {
981:               const PetscInt *econe;
982:               PetscInt        econeSize, r, vs, vu;

984:               DMPlexGetConeSize(dm, cone[q], &econeSize);
985:               DMPlexGetCone(dm, cone[q], &econe);
986:               for (r = 0; r < econeSize; ++r) {
987:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
988:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
989:                 if (vs >= 0) continue;
990:                 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);
991:                 hasUnsplit   = PETSC_TRUE;
992:               }
993:             }
994:           }
995:         }
996:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
997:         DMPlexSetConeOrientation(sdm, splitp, ornt);
998:         /* Face support */
999:         for (s = 0; s < supportSize; ++s) {
1000:           PetscInt val;

1002:           DMLabelGetValue(label, support[s], &val);
1003:           if (val < 0) {
1004:             /* Split old face:   Replace negative side cell with cohesive cell */
1005:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1006:           } else {
1007:             /* Split new face:   Replace positive side cell with cohesive cell */
1008:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1009:             /* Get orientation for cohesive face */
1010:             {
1011:               const PetscInt *ncone, *nconeO;
1012:               PetscInt        nconeSize, nc;

1014:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1015:               DMPlexGetCone(dm, support[s], &ncone);
1016:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1017:               for (nc = 0; nc < nconeSize; ++nc) {
1018:                 if (ncone[nc] == oldp) {
1019:                   coneONew[0] = nconeO[nc];
1020:                   break;
1021:                 }
1022:               }
1023:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1024:             }
1025:           }
1026:         }
1027:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1028:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1029:         coneNew[1]  = splitp;
1030:         coneONew[1] = coneONew[0];
1031:         for (q = 0; q < coneSize; ++q) {
1032:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1033:           if (v < 0) {
1034:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1035:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1036:             coneONew[2+q] = 0;
1037:           } else {
1038:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1039:           }
1040:           coneONew[2+q] = 0;
1041:         }
1042:         DMPlexSetCone(sdm, hybcell, coneNew);
1043:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1044:         /* Label the hybrid cells on the boundary of the split */
1045:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1046:       } else if (dep == 0) {
1047:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1053:           DMLabelGetValue(label, support[e], &val);
1054:           if ((val == 1) || (val == (shift + 1))) {
1055:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1056:           }
1057:         }
1058:         supportNew[qn] = hybedge;
1059:         DMPlexSetSupport(sdm, newp, supportNew);
1060:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1061:         for (e = 0, qp = 0; e < supportSize; ++e) {
1062:           PetscInt val, edge;

1064:           DMLabelGetValue(label, support[e], &val);
1065:           if (val == 1) {
1066:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1067:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1068:             supportNew[qp++] = edge + pMaxNew[dep+1];
1069:           } else if (val == -(shift + 1)) {
1070:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1071:           }
1072:         }
1073:         supportNew[qp] = hybedge;
1074:         DMPlexSetSupport(sdm, splitp, supportNew);
1075:         /* Hybrid edge:    Old and new split vertex */
1076:         coneNew[0] = newp;
1077:         coneNew[1] = splitp;
1078:         DMPlexSetCone(sdm, hybedge, coneNew);
1079:         for (e = 0, qf = 0; e < supportSize; ++e) {
1080:           PetscInt val, edge;

1082:           DMLabelGetValue(label, support[e], &val);
1083:           if (val == 1) {
1084:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1085:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1086:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1087:           }
1088:         }
1089:         DMPlexSetSupport(sdm, hybedge, supportNew);
1090:       } else if (dep == dim-2) {
1091:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1093:         /* Split old edge:   old vertices in cone so no change */
1094:         /* Split new edge:   new vertices in cone */
1095:         for (q = 0; q < coneSize; ++q) {
1096:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1097:           if (v < 0) {
1098:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1099:             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);
1100:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1101:           } else {
1102:             coneNew[q] = v + pMaxNew[dep-1];
1103:           }
1104:         }
1105:         DMPlexSetCone(sdm, splitp, coneNew);
1106:         /* Split old edge: Faces in positive side cells and old split faces */
1107:         for (e = 0, q = 0; e < supportSize; ++e) {
1108:           PetscInt val;

1110:           DMLabelGetValue(label, support[e], &val);
1111:           if (val == dim-1) {
1112:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1113:           } else if (val == (shift + dim-1)) {
1114:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1115:           }
1116:         }
1117:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1118:         DMPlexSetSupport(sdm, newp, supportNew);
1119:         /* Split new edge: Faces in negative side cells and new split faces */
1120:         for (e = 0, q = 0; e < supportSize; ++e) {
1121:           PetscInt val, face;

1123:           DMLabelGetValue(label, support[e], &val);
1124:           if (val == dim-1) {
1125:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1126:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1127:             supportNew[q++] = face + pMaxNew[dep+1];
1128:           } else if (val == -(shift + dim-1)) {
1129:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1130:           }
1131:         }
1132:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1133:         DMPlexSetSupport(sdm, splitp, supportNew);
1134:         /* Hybrid face */
1135:         coneNew[0] = newp;
1136:         coneNew[1] = splitp;
1137:         for (v = 0; v < coneSize; ++v) {
1138:           PetscInt vertex;
1139:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1140:           if (vertex < 0) {
1141:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1142:             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);
1143:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1144:           } else {
1145:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1146:           }
1147:         }
1148:         DMPlexSetCone(sdm, hybface, coneNew);
1149:         for (e = 0, qf = 0; e < supportSize; ++e) {
1150:           PetscInt val, face;

1152:           DMLabelGetValue(label, support[e], &val);
1153:           if (val == dim-1) {
1154:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1155:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1156:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1157:           }
1158:         }
1159:         DMPlexSetSupport(sdm, hybface, supportNew);
1160:       }
1161:     }
1162:   }
1163:   for (dep = 0; dep <= depth; ++dep) {
1164:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1165:       const PetscInt  oldp   = unsplitPoints[dep][p];
1166:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1167:       const PetscInt *cone, *support, *ornt;
1168:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1170:       DMPlexGetConeSize(dm, oldp, &coneSize);
1171:       DMPlexGetCone(dm, oldp, &cone);
1172:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1173:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1174:       DMPlexGetSupport(dm, oldp, &support);
1175:       if (dep == 0) {
1176:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1178:         /* Unsplit vertex */
1179:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1180:         for (s = 0, q = 0; s < supportSize; ++s) {
1181:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1182:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1183:           if (e >= 0) {
1184:             supportNew[q++] = e + pMaxNew[dep+1];
1185:           }
1186:         }
1187:         supportNew[q++] = hybedge;
1188:         supportNew[q++] = hybedge;
1189:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1190:         DMPlexSetSupport(sdm, newp, supportNew);
1191:         /* Hybrid edge */
1192:         coneNew[0] = newp;
1193:         coneNew[1] = newp;
1194:         DMPlexSetCone(sdm, hybedge, coneNew);
1195:         for (e = 0, qf = 0; e < supportSize; ++e) {
1196:           PetscInt val, edge;

1198:           DMLabelGetValue(label, support[e], &val);
1199:           if (val == 1) {
1200:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1201:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1202:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1203:           } else if  (val ==  (shift2 + 1)) {
1204:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1205:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1206:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1207:           }
1208:         }
1209:         DMPlexSetSupport(sdm, hybedge, supportNew);
1210:       } else if (dep == dim-2) {
1211:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1217:           DMLabelGetValue(label, support[f], &val);
1218:           if (val == dim-1) {
1219:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1220:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1221:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1222:             supportNew[qf++] = face + pMaxNew[dep+1];
1223:           } else {
1224:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1225:           }
1226:         }
1227:         supportNew[qf++] = hybface;
1228:         supportNew[qf++] = hybface;
1229:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1230:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1231:         DMPlexSetSupport(sdm, newp, supportNew);
1232:         /* Add hybrid face */
1233:         coneNew[0] = newp;
1234:         coneNew[1] = newp;
1235:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1236:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1237:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1238:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1239:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1240:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1241:         DMPlexSetCone(sdm, hybface, coneNew);
1242:         for (f = 0, qf = 0; f < supportSize; ++f) {
1243:           PetscInt val, face;

1245:           DMLabelGetValue(label, support[f], &val);
1246:           if (val == dim-1) {
1247:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1248:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1249:           }
1250:         }
1251:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1252:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1253:         DMPlexSetSupport(sdm, hybface, supportNew);
1254:       }
1255:     }
1256:   }
1257:   /* Step 6b: Replace split points in negative side cones */
1258:   for (sp = 0; sp < numSP; ++sp) {
1259:     PetscInt        dep = values[sp];
1260:     IS              pIS;
1261:     PetscInt        numPoints;
1262:     const PetscInt *points;

1264:     if (dep >= 0) continue;
1265:     DMLabelGetStratumIS(label, dep, &pIS);
1266:     if (!pIS) continue;
1267:     dep  = -dep - shift;
1268:     ISGetLocalSize(pIS, &numPoints);
1269:     ISGetIndices(pIS, &points);
1270:     for (p = 0; p < numPoints; ++p) {
1271:       const PetscInt  oldp = points[p];
1272:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1273:       const PetscInt *cone;
1274:       PetscInt        coneSize, c;
1275:       /* PetscBool       replaced = PETSC_FALSE; */

1277:       /* Negative edge: replace split vertex */
1278:       /* Negative cell: replace split face */
1279:       DMPlexGetConeSize(sdm, newp, &coneSize);
1280:       DMPlexGetCone(sdm, newp, &cone);
1281:       for (c = 0; c < coneSize; ++c) {
1282:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1283:         PetscInt       csplitp, cp, val;

1285:         DMLabelGetValue(label, coldp, &val);
1286:         if (val == dep-1) {
1287:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1288:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1289:           csplitp  = pMaxNew[dep-1] + cp;
1290:           DMPlexInsertCone(sdm, newp, c, csplitp);
1291:           /* replaced = PETSC_TRUE; */
1292:         }
1293:       }
1294:       /* Cells with only a vertex or edge on the submesh have no replacement */
1295:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1296:     }
1297:     ISRestoreIndices(pIS, &points);
1298:     ISDestroy(&pIS);
1299:   }
1300:   /* Step 7: Stratify */
1301:   DMPlexStratify(sdm);
1302:   /* Step 8: Coordinates */
1303:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1304:   DMGetCoordinateSection(sdm, &coordSection);
1305:   DMGetCoordinatesLocal(sdm, &coordinates);
1306:   VecGetArray(coordinates, &coords);
1307:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1308:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1309:     const PetscInt splitp = pMaxNew[0] + v;
1310:     PetscInt       dof, off, soff, d;

1312:     PetscSectionGetDof(coordSection, newp, &dof);
1313:     PetscSectionGetOffset(coordSection, newp, &off);
1314:     PetscSectionGetOffset(coordSection, splitp, &soff);
1315:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1316:   }
1317:   VecRestoreArray(coordinates, &coords);
1318:   /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1319:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1320:   /* Step 10: Labels */
1321:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1322:   DMGetNumLabels(sdm, &numLabels);
1323:   for (dep = 0; dep <= depth; ++dep) {
1324:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1325:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1326:       const PetscInt splitp = pMaxNew[dep] + p;
1327:       PetscInt       l;

1329:       for (l = 0; l < numLabels; ++l) {
1330:         DMLabel     mlabel;
1331:         const char *lname;
1332:         PetscInt    val;
1333:         PetscBool   isDepth;

1335:         DMGetLabelName(sdm, l, &lname);
1336:         PetscStrcmp(lname, "depth", &isDepth);
1337:         if (isDepth) continue;
1338:         DMGetLabel(sdm, lname, &mlabel);
1339:         DMLabelGetValue(mlabel, newp, &val);
1340:         if (val >= 0) {
1341:           DMLabelSetValue(mlabel, splitp, val);
1342: #if 0
1343:           /* Do not put cohesive edges into the label */
1344:           if (dep == 0) {
1345:             const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1346:             DMLabelSetValue(mlabel, cedge, val);
1347:           } else if (dep == dim-2) {
1348:             const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1349:             DMLabelSetValue(mlabel, cface, val);
1350:           }
1351:           /* Do not put cohesive faces into the label */
1352: #endif
1353:         }
1354:       }
1355:     }
1356:   }
1357:   for (sp = 0; sp < numSP; ++sp) {
1358:     const PetscInt dep = values[sp];

1360:     if ((dep < 0) || (dep > depth)) continue;
1361:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1362:     ISDestroy(&splitIS[dep]);
1363:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1364:     ISDestroy(&unsplitIS[dep]);
1365:   }
1366:   if (label) {
1367:     ISRestoreIndices(valueIS, &values);
1368:     ISDestroy(&valueIS);
1369:   }
1370:   for (d = 0; d <= depth; ++d) {
1371:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1372:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1373:   }
1374:   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);
1375:   PetscFree3(coneNew, coneONew, supportNew);
1376:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1377:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1378:   return(0);
1379: }

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

1386:   Collective on dm

1388:   Input Parameters:
1389: + dm - The original DM
1390: - label - The label specifying the boundary faces (this could be auto-generated)

1392:   Output Parameters:
1393: - dmSplit - The new DM

1395:   Level: developer

1397: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1398: @*/
1399: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1400: {
1401:   DM             sdm;
1402:   PetscInt       dim;

1408:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1409:   DMSetType(sdm, DMPLEX);
1410:   DMGetDimension(dm, &dim);
1411:   DMSetDimension(sdm, dim);
1412:   switch (dim) {
1413:   case 2:
1414:   case 3:
1415:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1416:     break;
1417:   default:
1418:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1419:   }
1420:   *dmSplit = sdm;
1421:   return(0);
1422: }

1426: /* Returns the side of the surface for a given cell with a face on the surface */
1427: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1428: {
1429:   const PetscInt *cone, *ornt;
1430:   PetscInt        dim, coneSize, c;
1431:   PetscErrorCode  ierr;

1434:   *pos = PETSC_TRUE;
1435:   DMGetDimension(dm, &dim);
1436:   DMPlexGetConeSize(dm, cell, &coneSize);
1437:   DMPlexGetCone(dm, cell, &cone);
1438:   DMPlexGetConeOrientation(dm, cell, &ornt);
1439:   for (c = 0; c < coneSize; ++c) {
1440:     if (cone[c] == face) {
1441:       PetscInt o = ornt[c];

1443:       if (subdm) {
1444:         const PetscInt *subcone, *subornt;
1445:         PetscInt        subpoint, subface, subconeSize, sc;

1447:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1448:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1449:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1450:         DMPlexGetCone(subdm, subpoint, &subcone);
1451:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1452:         for (sc = 0; sc < subconeSize; ++sc) {
1453:           if (subcone[sc] == subface) {
1454:             o = subornt[0];
1455:             break;
1456:           }
1457:         }
1458:         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);
1459:       }
1460:       if (o >= 0) *pos = PETSC_TRUE;
1461:       else        *pos = PETSC_FALSE;
1462:       break;
1463:     }
1464:   }
1465:   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);
1466:   return(0);
1467: }

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

1475:   Input Parameters:
1476: + dm     - The DM
1477: . label  - A DMLabel marking the surface
1478: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1479: . flip   - Flag to flip the submesh normal and replace points on the other side
1480: - subdm  - The subDM associated with the label, or NULL

1482:   Output Parameter:
1483: . label - A DMLabel marking all surface points

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

1487:   Level: developer

1489: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1490: @*/
1491: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1492: {
1493:   DMLabel         depthLabel;
1494:   IS              dimIS, subpointIS, facePosIS, faceNegIS, crossEdgeIS = NULL;
1495:   const PetscInt *points, *subpoints;
1496:   const PetscInt  rev   = flip ? -1 : 1;
1497:   PetscInt       *pMax;
1498:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1499:   PetscErrorCode  ierr;

1502:   DMPlexGetDepth(dm, &depth);
1503:   DMGetDimension(dm, &dim);
1504:   pSize = PetscMax(depth, dim) + 1;
1505:   PetscMalloc1(pSize,&pMax);
1506:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1507:   DMPlexGetDepthLabel(dm, &depthLabel);
1508:   DMGetDimension(dm, &dim);
1509:   if (subdm) {
1510:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1511:     if (subpointIS) {
1512:       ISGetLocalSize(subpointIS, &numSubpoints);
1513:       ISGetIndices(subpointIS, &subpoints);
1514:     }
1515:   }
1516:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1517:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1518:   if (!dimIS) {
1519:     PetscFree(pMax);
1520:     return(0);
1521:   }
1522:   ISGetLocalSize(dimIS, &numPoints);
1523:   ISGetIndices(dimIS, &points);
1524:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1525:     const PetscInt *support;
1526:     PetscInt        supportSize, s;

1528:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1529:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1530:     DMPlexGetSupport(dm, points[p], &support);
1531:     for (s = 0; s < supportSize; ++s) {
1532:       const PetscInt *cone;
1533:       PetscInt        coneSize, c;
1534:       PetscBool       pos;

1536:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1537:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1538:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1539:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1540:       /* Put faces touching the fault in the label */
1541:       DMPlexGetConeSize(dm, support[s], &coneSize);
1542:       DMPlexGetCone(dm, support[s], &cone);
1543:       for (c = 0; c < coneSize; ++c) {
1544:         const PetscInt point = cone[c];

1546:         DMLabelGetValue(label, point, &val);
1547:         if (val == -1) {
1548:           PetscInt *closure = NULL;
1549:           PetscInt  closureSize, cl;

1551:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1552:           for (cl = 0; cl < closureSize*2; cl += 2) {
1553:             const PetscInt clp  = closure[cl];
1554:             PetscInt       bval = -1;

1556:             DMLabelGetValue(label, clp, &val);
1557:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1558:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1559:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1560:               break;
1561:             }
1562:           }
1563:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1564:         }
1565:       }
1566:     }
1567:   }
1568:   ISRestoreIndices(dimIS, &points);
1569:   ISDestroy(&dimIS);
1570:   if (subdm) {
1571:     if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1572:     ISDestroy(&subpointIS);
1573:   }
1574:   /* Mark boundary points as unsplit */
1575:   if (blabel) {
1576:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1577:     ISGetLocalSize(dimIS, &numPoints);
1578:     ISGetIndices(dimIS, &points);
1579:     for (p = 0; p < numPoints; ++p) {
1580:       const PetscInt point = points[p];
1581:       PetscInt       val, bval;

1583:       DMLabelGetValue(blabel, point, &bval);
1584:       if (bval >= 0) {
1585:         DMLabelGetValue(label, point, &val);
1586:         if ((val < 0) || (val > dim)) {
1587:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1588:           DMLabelClearValue(blabel, point, bval);
1589:         }
1590:       }
1591:     }
1592:     for (p = 0; p < numPoints; ++p) {
1593:       const PetscInt point = points[p];
1594:       PetscInt       val, bval;

1596:       DMLabelGetValue(blabel, point, &bval);
1597:       if (bval >= 0) {
1598:         const PetscInt *cone,    *support;
1599:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1601:         /* Mark as unsplit */
1602:         DMLabelGetValue(label, point, &val);
1603:         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);
1604:         DMLabelClearValue(label, point, val);
1605:         DMLabelSetValue(label, point, shift2+val);
1606:         /* Check for cross-edge
1607:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1608:         if (val != 0) continue;
1609:         DMPlexGetSupport(dm, point, &support);
1610:         DMPlexGetSupportSize(dm, point, &supportSize);
1611:         for (s = 0; s < supportSize; ++s) {
1612:           DMPlexGetCone(dm, support[s], &cone);
1613:           DMPlexGetConeSize(dm, support[s], &coneSize);
1614:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1615:           DMLabelGetValue(blabel, cone[0], &valA);
1616:           DMLabelGetValue(blabel, cone[1], &valB);
1617:           DMLabelGetValue(blabel, support[s], &valE);
1618:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1619:         }
1620:       }
1621:     }
1622:     ISRestoreIndices(dimIS, &points);
1623:     ISDestroy(&dimIS);
1624:   }
1625:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1626:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1627:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1628:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1629:   cMax = cMax < 0 ? cEnd : cMax;
1630:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1631:   DMLabelGetStratumIS(label, 0, &dimIS);
1632:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1633:   if (dimIS && crossEdgeIS) {
1634:     IS vertIS = dimIS;

1636:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1637:     ISDestroy(&crossEdgeIS);
1638:     ISDestroy(&vertIS);
1639:   }
1640:   if (!dimIS) {
1641:     PetscFree(pMax);
1642:     return(0);
1643:   }
1644:   ISGetLocalSize(dimIS, &numPoints);
1645:   ISGetIndices(dimIS, &points);
1646:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1647:     PetscInt *star = NULL;
1648:     PetscInt  starSize, s;
1649:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1651:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1652:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1653:     while (again) {
1654:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1655:       again = 0;
1656:       for (s = 0; s < starSize*2; s += 2) {
1657:         const PetscInt  point = star[s];
1658:         const PetscInt *cone;
1659:         PetscInt        coneSize, c;

1661:         if ((point < cStart) || (point >= cMax)) continue;
1662:         DMLabelGetValue(label, point, &val);
1663:         if (val != -1) continue;
1664:         again = again == 1 ? 1 : 2;
1665:         DMPlexGetConeSize(dm, point, &coneSize);
1666:         DMPlexGetCone(dm, point, &cone);
1667:         for (c = 0; c < coneSize; ++c) {
1668:           DMLabelGetValue(label, cone[c], &val);
1669:           if (val != -1) {
1670:             const PetscInt *ccone;
1671:             PetscInt        cconeSize, cc, side;

1673:             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);
1674:             if (val > 0) side =  1;
1675:             else         side = -1;
1676:             DMLabelSetValue(label, point, side*(shift+dim));
1677:             /* Mark cell faces which touch the fault */
1678:             DMPlexGetConeSize(dm, point, &cconeSize);
1679:             DMPlexGetCone(dm, point, &ccone);
1680:             for (cc = 0; cc < cconeSize; ++cc) {
1681:               PetscInt *closure = NULL;
1682:               PetscInt  closureSize, cl;

1684:               DMLabelGetValue(label, ccone[cc], &val);
1685:               if (val != -1) continue;
1686:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1687:               for (cl = 0; cl < closureSize*2; cl += 2) {
1688:                 const PetscInt clp = closure[cl];

1690:                 DMLabelGetValue(label, clp, &val);
1691:                 if (val == -1) continue;
1692:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1693:                 break;
1694:               }
1695:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1696:             }
1697:             again = 1;
1698:             break;
1699:           }
1700:         }
1701:       }
1702:     }
1703:     /* Classify the rest by cell membership */
1704:     for (s = 0; s < starSize*2; s += 2) {
1705:       const PetscInt point = star[s];

1707:       DMLabelGetValue(label, point, &val);
1708:       if (val == -1) {
1709:         PetscInt  *sstar = NULL;
1710:         PetscInt   sstarSize, ss;
1711:         PetscBool  marked = PETSC_FALSE;

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

1717:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1718:           DMLabelGetValue(label, spoint, &val);
1719:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1720:           DMLabelGetValue(depthLabel, point, &dep);
1721:           if (val > 0) {
1722:             DMLabelSetValue(label, point,   shift+dep);
1723:           } else {
1724:             DMLabelSetValue(label, point, -(shift+dep));
1725:           }
1726:           marked = PETSC_TRUE;
1727:           break;
1728:         }
1729:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1730:         DMLabelGetValue(depthLabel, point, &dep);
1731:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1732:       }
1733:     }
1734:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1735:   }
1736:   ISRestoreIndices(dimIS, &points);
1737:   ISDestroy(&dimIS);
1738:   /* If any faces touching the fault divide cells on either side, split them */
1739:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1740:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1741:   ISExpand(facePosIS, faceNegIS, &dimIS);
1742:   ISDestroy(&facePosIS);
1743:   ISDestroy(&faceNegIS);
1744:   ISGetLocalSize(dimIS, &numPoints);
1745:   ISGetIndices(dimIS, &points);
1746:   for (p = 0; p < numPoints; ++p) {
1747:     const PetscInt  point = points[p];
1748:     const PetscInt *support;
1749:     PetscInt        supportSize, valA, valB;

1751:     DMPlexGetSupportSize(dm, point, &supportSize);
1752:     if (supportSize != 2) continue;
1753:     DMPlexGetSupport(dm, point, &support);
1754:     DMLabelGetValue(label, support[0], &valA);
1755:     DMLabelGetValue(label, support[1], &valB);
1756:     if ((valA == -1) || (valB == -1)) continue;
1757:     if (valA*valB > 0) continue;
1758:     /* Split the face */
1759:     DMLabelGetValue(label, point, &valA);
1760:     DMLabelClearValue(label, point, valA);
1761:     DMLabelSetValue(label, point, dim-1);
1762:     /* Label its closure:
1763:       unmarked: label as unsplit
1764:       incident: relabel as split
1765:       split:    do nothing
1766:     */
1767:     {
1768:       PetscInt *closure = NULL;
1769:       PetscInt  closureSize, cl;

1771:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1772:       for (cl = 0; cl < closureSize*2; cl += 2) {
1773:         DMLabelGetValue(label, closure[cl], &valA);
1774:         if (valA == -1) { /* Mark as unsplit */
1775:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1776:           DMLabelSetValue(label, closure[cl], shift2+dep);
1777:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1778:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1779:           DMLabelClearValue(label, closure[cl], valA);
1780:           DMLabelSetValue(label, closure[cl], dep);
1781:         }
1782:       }
1783:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1784:     }
1785:   }
1786:   ISRestoreIndices(dimIS, &points);
1787:   ISDestroy(&dimIS);
1788:   PetscFree(pMax);
1789:   return(0);
1790: }

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

1797:   Collective on dm

1799:   Input Parameters:
1800: + dm - The original DM
1801: - labelName - The label specifying the interface vertices

1803:   Output Parameters:
1804: + hybridLabel - The label fully marking the interface
1805: - dmHybrid - The new DM

1807:   Level: developer

1809: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1810: @*/
1811: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1812: {
1813:   DM             idm;
1814:   DMLabel        subpointMap, hlabel;
1815:   PetscInt       dim;

1822:   DMGetDimension(dm, &dim);
1823:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1824:   DMPlexOrient(idm);
1825:   DMPlexGetSubpointMap(idm, &subpointMap);
1826:   DMLabelDuplicate(subpointMap, &hlabel);
1827:   DMLabelClearStratum(hlabel, dim);
1828:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1829:   DMDestroy(&idm);
1830:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1831:   if (hybridLabel) *hybridLabel = hlabel;
1832:   else             {DMLabelDestroy(&hlabel);}
1833:   return(0);
1834: }

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

1840:      For any marked cell, the marked vertices constitute a single face
1841: */
1842: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1843: {
1844:   IS               subvertexIS = NULL;
1845:   const PetscInt  *subvertices;
1846:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1847:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1848:   PetscErrorCode   ierr;

1851:   *numFaces = 0;
1852:   *nFV      = 0;
1853:   DMPlexGetDepth(dm, &depth);
1854:   DMGetDimension(dm, &dim);
1855:   pSize = PetscMax(depth, dim) + 1;
1856:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1857:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1858:   for (d = 0; d <= depth; ++d) {
1859:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1860:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1861:   }
1862:   /* Loop over initial vertices and mark all faces in the collective star() */
1863:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1864:   if (subvertexIS) {
1865:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1866:     ISGetIndices(subvertexIS, &subvertices);
1867:   }
1868:   for (v = 0; v < numSubVerticesInitial; ++v) {
1869:     const PetscInt vertex = subvertices[v];
1870:     PetscInt      *star   = NULL;
1871:     PetscInt       starSize, s, numCells = 0, c;

1873:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1874:     for (s = 0; s < starSize*2; s += 2) {
1875:       const PetscInt point = star[s];
1876:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1877:     }
1878:     for (c = 0; c < numCells; ++c) {
1879:       const PetscInt cell    = star[c];
1880:       PetscInt      *closure = NULL;
1881:       PetscInt       closureSize, cl;
1882:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1884:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1885:       if (cellLoc == 2) continue;
1886:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1887:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1888:       for (cl = 0; cl < closureSize*2; cl += 2) {
1889:         const PetscInt point = closure[cl];
1890:         PetscInt       vertexLoc;

1892:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1893:           ++numCorners;
1894:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1895:           if (vertexLoc == value) closure[faceSize++] = point;
1896:         }
1897:       }
1898:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1899:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1900:       if (faceSize == *nFV) {
1901:         const PetscInt *cells = NULL;
1902:         PetscInt        numCells, nc;

1904:         ++(*numFaces);
1905:         for (cl = 0; cl < faceSize; ++cl) {
1906:           DMLabelSetValue(subpointMap, closure[cl], 0);
1907:         }
1908:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1909:         for (nc = 0; nc < numCells; ++nc) {
1910:           DMLabelSetValue(subpointMap, cells[nc], 2);
1911:         }
1912:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1913:       }
1914:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1915:     }
1916:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1917:   }
1918:   if (subvertexIS) {
1919:     ISRestoreIndices(subvertexIS, &subvertices);
1920:   }
1921:   ISDestroy(&subvertexIS);
1922:   PetscFree3(pStart,pEnd,pMax);
1923:   return(0);
1924: }

1928: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1929: {
1930:   IS               subvertexIS = NULL;
1931:   const PetscInt  *subvertices;
1932:   PetscInt        *pStart, *pEnd, *pMax;
1933:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
1934:   PetscErrorCode   ierr;

1937:   DMGetDimension(dm, &dim);
1938:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1939:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1940:   for (d = 0; d <= dim; ++d) {
1941:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1942:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1943:   }
1944:   /* Loop over initial vertices and mark all faces in the collective star() */
1945:   if (vertexLabel) {
1946:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1947:     if (subvertexIS) {
1948:       ISGetSize(subvertexIS, &numSubVerticesInitial);
1949:       ISGetIndices(subvertexIS, &subvertices);
1950:     }
1951:   }
1952:   for (v = 0; v < numSubVerticesInitial; ++v) {
1953:     const PetscInt vertex = subvertices[v];
1954:     PetscInt      *star   = NULL;
1955:     PetscInt       starSize, s, numFaces = 0, f;

1957:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1958:     for (s = 0; s < starSize*2; s += 2) {
1959:       const PetscInt point = star[s];
1960:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1961:     }
1962:     for (f = 0; f < numFaces; ++f) {
1963:       const PetscInt face    = star[f];
1964:       PetscInt      *closure = NULL;
1965:       PetscInt       closureSize, c;
1966:       PetscInt       faceLoc;

1968:       DMLabelGetValue(subpointMap, face, &faceLoc);
1969:       if (faceLoc == dim-1) continue;
1970:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1971:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1972:       for (c = 0; c < closureSize*2; c += 2) {
1973:         const PetscInt point = closure[c];
1974:         PetscInt       vertexLoc;

1976:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1977:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1978:           if (vertexLoc != value) break;
1979:         }
1980:       }
1981:       if (c == closureSize*2) {
1982:         const PetscInt *support;
1983:         PetscInt        supportSize, s;

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

1988:           for (d = 0; d < dim; ++d) {
1989:             if ((point >= pStart[d]) && (point < pEnd[d])) {
1990:               DMLabelSetValue(subpointMap, point, d);
1991:               break;
1992:             }
1993:           }
1994:         }
1995:         DMPlexGetSupportSize(dm, face, &supportSize);
1996:         DMPlexGetSupport(dm, face, &support);
1997:         for (s = 0; s < supportSize; ++s) {
1998:           DMLabelSetValue(subpointMap, support[s], dim);
1999:         }
2000:       }
2001:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2002:     }
2003:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2004:   }
2005:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2006:   ISDestroy(&subvertexIS);
2007:   PetscFree3(pStart,pEnd,pMax);
2008:   return(0);
2009: }

2013: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2014: {
2015:   DMLabel         label = NULL;
2016:   const PetscInt *cone;
2017:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize;
2018:   PetscErrorCode  ierr;

2021:   *numFaces = 0;
2022:   *nFV = 0;
2023:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2024:   *subCells = NULL;
2025:   DMGetDimension(dm, &dim);
2026:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2027:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2028:   if (cMax < 0) return(0);
2029:   if (label) {
2030:     for (c = cMax; c < cEnd; ++c) {
2031:       PetscInt val;

2033:       DMLabelGetValue(label, c, &val);
2034:       if (val == value) {
2035:         ++(*numFaces);
2036:         DMPlexGetConeSize(dm, c, &coneSize);
2037:       }
2038:     }
2039:   } else {
2040:     *numFaces = cEnd - cMax;
2041:     DMPlexGetConeSize(dm, cMax, &coneSize);
2042:   }
2043:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2044:   PetscMalloc1(*numFaces *2, subCells);
2045:   for (c = cMax; c < cEnd; ++c) {
2046:     const PetscInt *cells;
2047:     PetscInt        numCells;

2049:     if (label) {
2050:       PetscInt val;

2052:       DMLabelGetValue(label, c, &val);
2053:       if (val != value) continue;
2054:     }
2055:     DMPlexGetCone(dm, c, &cone);
2056:     for (p = 0; p < *nFV; ++p) {
2057:       DMLabelSetValue(subpointMap, cone[p], 0);
2058:     }
2059:     /* Negative face */
2060:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2061:     /* Not true in parallel
2062:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2063:     for (p = 0; p < numCells; ++p) {
2064:       DMLabelSetValue(subpointMap, cells[p], 2);
2065:       (*subCells)[subc++] = cells[p];
2066:     }
2067:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2068:     /* Positive face is not included */
2069:   }
2070:   return(0);
2071: }

2075: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2076: {
2077:   PetscInt      *pStart, *pEnd;
2078:   PetscInt       dim, cMax, cEnd, c, d;

2082:   DMGetDimension(dm, &dim);
2083:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2084:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2085:   if (cMax < 0) return(0);
2086:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2087:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2088:   for (c = cMax; c < cEnd; ++c) {
2089:     const PetscInt *cone;
2090:     PetscInt       *closure = NULL;
2091:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2093:     if (label) {
2094:       DMLabelGetValue(label, c, &val);
2095:       if (val != value) continue;
2096:     }
2097:     DMPlexGetConeSize(dm, c, &coneSize);
2098:     DMPlexGetCone(dm, c, &cone);
2099:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2100:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2101:     /* Negative face */
2102:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2103:     for (cl = 0; cl < closureSize*2; cl += 2) {
2104:       const PetscInt point = closure[cl];

2106:       for (d = 0; d <= dim; ++d) {
2107:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2108:           DMLabelSetValue(subpointMap, point, d);
2109:           break;
2110:         }
2111:       }
2112:     }
2113:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2114:     /* Cells -- positive face is not included */
2115:     for (cl = 0; cl < 1; ++cl) {
2116:       const PetscInt *support;
2117:       PetscInt        supportSize, s;

2119:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2120:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2121:       DMPlexGetSupport(dm, cone[cl], &support);
2122:       for (s = 0; s < supportSize; ++s) {
2123:         DMLabelSetValue(subpointMap, support[s], dim);
2124:       }
2125:     }
2126:   }
2127:   PetscFree2(pStart, pEnd);
2128:   return(0);
2129: }

2133: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2134: {
2135:   MPI_Comm       comm;
2136:   PetscBool      posOrient = PETSC_FALSE;
2137:   const PetscInt debug     = 0;
2138:   PetscInt       cellDim, faceSize, f;

2142:   PetscObjectGetComm((PetscObject)dm,&comm);
2143:   DMGetDimension(dm, &cellDim);
2144:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2146:   if (cellDim == 1 && numCorners == 2) {
2147:     /* Triangle */
2148:     faceSize  = numCorners-1;
2149:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2150:   } else if (cellDim == 2 && numCorners == 3) {
2151:     /* Triangle */
2152:     faceSize  = numCorners-1;
2153:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2154:   } else if (cellDim == 3 && numCorners == 4) {
2155:     /* Tetrahedron */
2156:     faceSize  = numCorners-1;
2157:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2158:   } else if (cellDim == 1 && numCorners == 3) {
2159:     /* Quadratic line */
2160:     faceSize  = 1;
2161:     posOrient = PETSC_TRUE;
2162:   } else if (cellDim == 2 && numCorners == 4) {
2163:     /* Quads */
2164:     faceSize = 2;
2165:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2166:       posOrient = PETSC_TRUE;
2167:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2168:       posOrient = PETSC_TRUE;
2169:     } else {
2170:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2171:         posOrient = PETSC_FALSE;
2172:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2173:     }
2174:   } else if (cellDim == 2 && numCorners == 6) {
2175:     /* Quadratic triangle (I hate this) */
2176:     /* Edges are determined by the first 2 vertices (corners of edges) */
2177:     const PetscInt faceSizeTri = 3;
2178:     PetscInt       sortedIndices[3], i, iFace;
2179:     PetscBool      found                    = PETSC_FALSE;
2180:     PetscInt       faceVerticesTriSorted[9] = {
2181:       0, 3,  4, /* bottom */
2182:       1, 4,  5, /* right */
2183:       2, 3,  5, /* left */
2184:     };
2185:     PetscInt       faceVerticesTri[9] = {
2186:       0, 3,  4, /* bottom */
2187:       1, 4,  5, /* right */
2188:       2, 5,  3, /* left */
2189:     };

2191:     faceSize = faceSizeTri;
2192:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2193:     PetscSortInt(faceSizeTri, sortedIndices);
2194:     for (iFace = 0; iFace < 3; ++iFace) {
2195:       const PetscInt ii = iFace*faceSizeTri;
2196:       PetscInt       fVertex, cVertex;

2198:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2199:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2200:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2201:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2202:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2203:               faceVertices[fVertex] = origVertices[cVertex];
2204:               break;
2205:             }
2206:           }
2207:         }
2208:         found = PETSC_TRUE;
2209:         break;
2210:       }
2211:     }
2212:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2213:     if (posOriented) *posOriented = PETSC_TRUE;
2214:     return(0);
2215:   } else if (cellDim == 2 && numCorners == 9) {
2216:     /* Quadratic quad (I hate this) */
2217:     /* Edges are determined by the first 2 vertices (corners of edges) */
2218:     const PetscInt faceSizeQuad = 3;
2219:     PetscInt       sortedIndices[3], i, iFace;
2220:     PetscBool      found                      = PETSC_FALSE;
2221:     PetscInt       faceVerticesQuadSorted[12] = {
2222:       0, 1,  4, /* bottom */
2223:       1, 2,  5, /* right */
2224:       2, 3,  6, /* top */
2225:       0, 3,  7, /* left */
2226:     };
2227:     PetscInt       faceVerticesQuad[12] = {
2228:       0, 1,  4, /* bottom */
2229:       1, 2,  5, /* right */
2230:       2, 3,  6, /* top */
2231:       3, 0,  7, /* left */
2232:     };

2234:     faceSize = faceSizeQuad;
2235:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2236:     PetscSortInt(faceSizeQuad, sortedIndices);
2237:     for (iFace = 0; iFace < 4; ++iFace) {
2238:       const PetscInt ii = iFace*faceSizeQuad;
2239:       PetscInt       fVertex, cVertex;

2241:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2242:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2243:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2244:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2245:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2246:               faceVertices[fVertex] = origVertices[cVertex];
2247:               break;
2248:             }
2249:           }
2250:         }
2251:         found = PETSC_TRUE;
2252:         break;
2253:       }
2254:     }
2255:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2256:     if (posOriented) *posOriented = PETSC_TRUE;
2257:     return(0);
2258:   } else if (cellDim == 3 && numCorners == 8) {
2259:     /* Hexes
2260:        A hex is two oriented quads with the normal of the first
2261:        pointing up at the second.

2263:           7---6
2264:          /|  /|
2265:         4---5 |
2266:         | 1-|-2
2267:         |/  |/
2268:         0---3

2270:         Faces are determined by the first 4 vertices (corners of faces) */
2271:     const PetscInt faceSizeHex = 4;
2272:     PetscInt       sortedIndices[4], i, iFace;
2273:     PetscBool      found                     = PETSC_FALSE;
2274:     PetscInt       faceVerticesHexSorted[24] = {
2275:       0, 1, 2, 3,  /* bottom */
2276:       4, 5, 6, 7,  /* top */
2277:       0, 3, 4, 5,  /* front */
2278:       2, 3, 5, 6,  /* right */
2279:       1, 2, 6, 7,  /* back */
2280:       0, 1, 4, 7,  /* left */
2281:     };
2282:     PetscInt       faceVerticesHex[24] = {
2283:       1, 2, 3, 0,  /* bottom */
2284:       4, 5, 6, 7,  /* top */
2285:       0, 3, 5, 4,  /* front */
2286:       3, 2, 6, 5,  /* right */
2287:       2, 1, 7, 6,  /* back */
2288:       1, 0, 4, 7,  /* left */
2289:     };

2291:     faceSize = faceSizeHex;
2292:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2293:     PetscSortInt(faceSizeHex, sortedIndices);
2294:     for (iFace = 0; iFace < 6; ++iFace) {
2295:       const PetscInt ii = iFace*faceSizeHex;
2296:       PetscInt       fVertex, cVertex;

2298:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2299:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2300:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2301:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2302:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2303:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2304:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2305:               faceVertices[fVertex] = origVertices[cVertex];
2306:               break;
2307:             }
2308:           }
2309:         }
2310:         found = PETSC_TRUE;
2311:         break;
2312:       }
2313:     }
2314:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2315:     if (posOriented) *posOriented = PETSC_TRUE;
2316:     return(0);
2317:   } else if (cellDim == 3 && numCorners == 10) {
2318:     /* Quadratic tet */
2319:     /* Faces are determined by the first 3 vertices (corners of faces) */
2320:     const PetscInt faceSizeTet = 6;
2321:     PetscInt       sortedIndices[6], i, iFace;
2322:     PetscBool      found                     = PETSC_FALSE;
2323:     PetscInt       faceVerticesTetSorted[24] = {
2324:       0, 1, 2,  6, 7, 8, /* bottom */
2325:       0, 3, 4,  6, 7, 9,  /* front */
2326:       1, 4, 5,  7, 8, 9,  /* right */
2327:       2, 3, 5,  6, 8, 9,  /* left */
2328:     };
2329:     PetscInt       faceVerticesTet[24] = {
2330:       0, 1, 2,  6, 7, 8, /* bottom */
2331:       0, 4, 3,  6, 7, 9,  /* front */
2332:       1, 5, 4,  7, 8, 9,  /* right */
2333:       2, 3, 5,  8, 6, 9,  /* left */
2334:     };

2336:     faceSize = faceSizeTet;
2337:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2338:     PetscSortInt(faceSizeTet, sortedIndices);
2339:     for (iFace=0; iFace < 4; ++iFace) {
2340:       const PetscInt ii = iFace*faceSizeTet;
2341:       PetscInt       fVertex, cVertex;

2343:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2344:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2345:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2346:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2347:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2348:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2349:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2350:               faceVertices[fVertex] = origVertices[cVertex];
2351:               break;
2352:             }
2353:           }
2354:         }
2355:         found = PETSC_TRUE;
2356:         break;
2357:       }
2358:     }
2359:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2360:     if (posOriented) *posOriented = PETSC_TRUE;
2361:     return(0);
2362:   } else if (cellDim == 3 && numCorners == 27) {
2363:     /* Quadratic hexes (I hate this)
2364:        A hex is two oriented quads with the normal of the first
2365:        pointing up at the second.

2367:          7---6
2368:         /|  /|
2369:        4---5 |
2370:        | 3-|-2
2371:        |/  |/
2372:        0---1

2374:        Faces are determined by the first 4 vertices (corners of faces) */
2375:     const PetscInt faceSizeQuadHex = 9;
2376:     PetscInt       sortedIndices[9], i, iFace;
2377:     PetscBool      found                         = PETSC_FALSE;
2378:     PetscInt       faceVerticesQuadHexSorted[54] = {
2379:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2380:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2381:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2382:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2383:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2384:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2385:     };
2386:     PetscInt       faceVerticesQuadHex[54] = {
2387:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2388:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2389:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2390:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2391:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2392:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2393:     };

2395:     faceSize = faceSizeQuadHex;
2396:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2397:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2398:     for (iFace = 0; iFace < 6; ++iFace) {
2399:       const PetscInt ii = iFace*faceSizeQuadHex;
2400:       PetscInt       fVertex, cVertex;

2402:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2403:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2404:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2405:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2406:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2407:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2408:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2409:               faceVertices[fVertex] = origVertices[cVertex];
2410:               break;
2411:             }
2412:           }
2413:         }
2414:         found = PETSC_TRUE;
2415:         break;
2416:       }
2417:     }
2418:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2419:     if (posOriented) *posOriented = PETSC_TRUE;
2420:     return(0);
2421:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2422:   if (!posOrient) {
2423:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2424:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2425:   } else {
2426:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2427:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2428:   }
2429:   if (posOriented) *posOriented = posOrient;
2430:   return(0);
2431: }

2435: /*
2436:     Given a cell and a face, as a set of vertices,
2437:       return the oriented face, as a set of vertices, in faceVertices
2438:     The orientation is such that the face normal points out of the cell
2439: */
2440: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2441: {
2442:   const PetscInt *cone = NULL;
2443:   PetscInt        coneSize, v, f, v2;
2444:   PetscInt        oppositeVertex = -1;
2445:   PetscErrorCode  ierr;

2448:   DMPlexGetConeSize(dm, cell, &coneSize);
2449:   DMPlexGetCone(dm, cell, &cone);
2450:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2451:     PetscBool found = PETSC_FALSE;

2453:     for (f = 0; f < faceSize; ++f) {
2454:       if (face[f] == cone[v]) {
2455:         found = PETSC_TRUE; break;
2456:       }
2457:     }
2458:     if (found) {
2459:       indices[v2]      = v;
2460:       origVertices[v2] = cone[v];
2461:       ++v2;
2462:     } else {
2463:       oppositeVertex = v;
2464:     }
2465:   }
2466:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2467:   return(0);
2468: }

2472: /*
2473:   DMPlexInsertFace_Internal - Puts a face into the mesh

2475:   Not collective

2477:   Input Parameters:
2478:   + dm              - The DMPlex
2479:   . numFaceVertex   - The number of vertices in the face
2480:   . faceVertices    - The vertices in the face for dm
2481:   . subfaceVertices - The vertices in the face for subdm
2482:   . numCorners      - The number of vertices in the cell
2483:   . cell            - A cell in dm containing the face
2484:   . subcell         - A cell in subdm containing the face
2485:   . firstFace       - First face in the mesh
2486:   - newFacePoint    - Next face in the mesh

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

2491:   Level: developer
2492: */
2493: 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)
2494: {
2495:   MPI_Comm        comm;
2496:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2497:   const PetscInt *faces;
2498:   PetscInt        numFaces, coneSize;
2499:   PetscErrorCode  ierr;

2502:   PetscObjectGetComm((PetscObject)dm,&comm);
2503:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2504:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2505: #if 0
2506:   /* Cannot use this because support() has not been constructed yet */
2507:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2508: #else
2509:   {
2510:     PetscInt f;

2512:     numFaces = 0;
2513:     DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2514:     for (f = firstFace; f < *newFacePoint; ++f) {
2515:       PetscInt dof, off, d;

2517:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2518:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2519:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2520:       for (d = 0; d < dof; ++d) {
2521:         const PetscInt p = submesh->cones[off+d];
2522:         PetscInt       v;

2524:         for (v = 0; v < numFaceVertices; ++v) {
2525:           if (subfaceVertices[v] == p) break;
2526:         }
2527:         if (v == numFaceVertices) break;
2528:       }
2529:       if (d == dof) {
2530:         numFaces               = 1;
2531:         ((PetscInt*) faces)[0] = f;
2532:       }
2533:     }
2534:   }
2535: #endif
2536:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2537:   else if (numFaces == 1) {
2538:     /* Add the other cell neighbor for this face */
2539:     DMPlexSetCone(subdm, subcell, faces);
2540:   } else {
2541:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2542:     PetscBool posOriented;

2544:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2545:     origVertices        = &orientedVertices[numFaceVertices];
2546:     indices             = &orientedVertices[numFaceVertices*2];
2547:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2548:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2549:     /* TODO: I know that routine should return a permutation, not the indices */
2550:     for (v = 0; v < numFaceVertices; ++v) {
2551:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2552:       for (ov = 0; ov < numFaceVertices; ++ov) {
2553:         if (orientedVertices[ov] == vertex) {
2554:           orientedSubVertices[ov] = subvertex;
2555:           break;
2556:         }
2557:       }
2558:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2559:     }
2560:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2561:     DMPlexSetCone(subdm, subcell, newFacePoint);
2562:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2563:     ++(*newFacePoint);
2564:   }
2565: #if 0
2566:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2567: #else
2568:   DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2569: #endif
2570:   return(0);
2571: }

2575: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2576: {
2577:   MPI_Comm        comm;
2578:   DMLabel         subpointMap;
2579:   IS              subvertexIS,  subcellIS;
2580:   const PetscInt *subVertices, *subCells;
2581:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2582:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2583:   PetscInt        vStart, vEnd, c, f;
2584:   PetscErrorCode  ierr;

2587:   PetscObjectGetComm((PetscObject)dm,&comm);
2588:   /* Create subpointMap which marks the submesh */
2589:   DMLabelCreate("subpoint_map", &subpointMap);
2590:   DMPlexSetSubpointMap(subdm, subpointMap);
2591:   DMLabelDestroy(&subpointMap);
2592:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2593:   /* Setup chart */
2594:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2595:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2596:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2597:   DMPlexSetVTKCellHeight(subdm, 1);
2598:   /* Set cone sizes */
2599:   firstSubVertex = numSubCells;
2600:   firstSubFace   = numSubCells+numSubVertices;
2601:   newFacePoint   = firstSubFace;
2602:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2603:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2604:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2605:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2606:   for (c = 0; c < numSubCells; ++c) {
2607:     DMPlexSetConeSize(subdm, c, 1);
2608:   }
2609:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2610:     DMPlexSetConeSize(subdm, f, nFV);
2611:   }
2612:   DMSetUp(subdm);
2613:   /* Create face cones */
2614:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2615:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2616:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2617:   for (c = 0; c < numSubCells; ++c) {
2618:     const PetscInt cell    = subCells[c];
2619:     const PetscInt subcell = c;
2620:     PetscInt      *closure = NULL;
2621:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2623:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2624:     for (cl = 0; cl < closureSize*2; cl += 2) {
2625:       const PetscInt point = closure[cl];
2626:       PetscInt       subVertex;

2628:       if ((point >= vStart) && (point < vEnd)) {
2629:         ++numCorners;
2630:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2631:         if (subVertex >= 0) {
2632:           closure[faceSize] = point;
2633:           subface[faceSize] = firstSubVertex+subVertex;
2634:           ++faceSize;
2635:         }
2636:       }
2637:     }
2638:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2639:     if (faceSize == nFV) {
2640:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2641:     }
2642:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2643:   }
2644:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2645:   DMPlexSymmetrize(subdm);
2646:   DMPlexStratify(subdm);
2647:   /* Build coordinates */
2648:   {
2649:     PetscSection coordSection, subCoordSection;
2650:     Vec          coordinates, subCoordinates;
2651:     PetscScalar *coords, *subCoords;
2652:     PetscInt     numComp, coordSize, v;
2653:     const char  *name;

2655:     DMGetCoordinateSection(dm, &coordSection);
2656:     DMGetCoordinatesLocal(dm, &coordinates);
2657:     DMGetCoordinateSection(subdm, &subCoordSection);
2658:     PetscSectionSetNumFields(subCoordSection, 1);
2659:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2660:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2661:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2662:     for (v = 0; v < numSubVertices; ++v) {
2663:       const PetscInt vertex    = subVertices[v];
2664:       const PetscInt subvertex = firstSubVertex+v;
2665:       PetscInt       dof;

2667:       PetscSectionGetDof(coordSection, vertex, &dof);
2668:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2669:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2670:     }
2671:     PetscSectionSetUp(subCoordSection);
2672:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2673:     VecCreate(comm, &subCoordinates);
2674:     PetscObjectGetName((PetscObject)coordinates,&name);
2675:     PetscObjectSetName((PetscObject)subCoordinates,name);
2676:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2677:     VecSetType(subCoordinates,VECSTANDARD);
2678:     if (coordSize) {
2679:       VecGetArray(coordinates,    &coords);
2680:       VecGetArray(subCoordinates, &subCoords);
2681:       for (v = 0; v < numSubVertices; ++v) {
2682:         const PetscInt vertex    = subVertices[v];
2683:         const PetscInt subvertex = firstSubVertex+v;
2684:         PetscInt       dof, off, sdof, soff, d;

2686:         PetscSectionGetDof(coordSection, vertex, &dof);
2687:         PetscSectionGetOffset(coordSection, vertex, &off);
2688:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2689:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2690:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2691:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2692:       }
2693:       VecRestoreArray(coordinates,    &coords);
2694:       VecRestoreArray(subCoordinates, &subCoords);
2695:     }
2696:     DMSetCoordinatesLocal(subdm, subCoordinates);
2697:     VecDestroy(&subCoordinates);
2698:   }
2699:   /* Cleanup */
2700:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2701:   ISDestroy(&subvertexIS);
2702:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2703:   ISDestroy(&subcellIS);
2704:   return(0);
2705: }

2709: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2710: {
2711:   PetscInt       subPoint;

2714:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2715:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2716: }

2720: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2721: {
2722:   MPI_Comm         comm;
2723:   DMLabel          subpointMap;
2724:   IS              *subpointIS;
2725:   const PetscInt **subpoints;
2726:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2727:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2728:   PetscErrorCode   ierr;

2731:   PetscObjectGetComm((PetscObject)dm,&comm);
2732:   /* Create subpointMap which marks the submesh */
2733:   DMLabelCreate("subpoint_map", &subpointMap);
2734:   DMPlexSetSubpointMap(subdm, subpointMap);
2735:   if (cellHeight) {
2736:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2737:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2738:   } else {
2739:     DMLabel         depth;
2740:     IS              pointIS;
2741:     const PetscInt *points;
2742:     PetscInt        numPoints;

2744:     DMPlexGetDepthLabel(dm, &depth);
2745:     DMLabelGetStratumSize(label, value, &numPoints);
2746:     DMLabelGetStratumIS(label, value, &pointIS);
2747:     ISGetIndices(pointIS, &points);
2748:     for (p = 0; p < numPoints; ++p) {
2749:       PetscInt *closure = NULL;
2750:       PetscInt  closureSize, c, pdim;

2752:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2753:       for (c = 0; c < closureSize*2; c += 2) {
2754:         DMLabelGetValue(depth, closure[c], &pdim);
2755:         DMLabelSetValue(subpointMap, closure[c], pdim);
2756:       }
2757:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2758:     }
2759:     ISRestoreIndices(pointIS, &points);
2760:     ISDestroy(&pointIS);
2761:   }
2762:   DMLabelDestroy(&subpointMap);
2763:   /* Setup chart */
2764:   DMGetDimension(dm, &dim);
2765:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2766:   for (d = 0; d <= dim; ++d) {
2767:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2768:     totSubPoints += numSubPoints[d];
2769:   }
2770:   DMPlexSetChart(subdm, 0, totSubPoints);
2771:   DMPlexSetVTKCellHeight(subdm, cellHeight);
2772:   /* Set cone sizes */
2773:   firstSubPoint[dim] = 0;
2774:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2775:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2776:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2777:   for (d = 0; d <= dim; ++d) {
2778:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2779:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2780:   }
2781:   for (d = 0; d <= dim; ++d) {
2782:     for (p = 0; p < numSubPoints[d]; ++p) {
2783:       const PetscInt  point    = subpoints[d][p];
2784:       const PetscInt  subpoint = firstSubPoint[d] + p;
2785:       const PetscInt *cone;
2786:       PetscInt        coneSize, coneSizeNew, c, val;

2788:       DMPlexGetConeSize(dm, point, &coneSize);
2789:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2790:       if (cellHeight && (d == dim)) {
2791:         DMPlexGetCone(dm, point, &cone);
2792:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2793:           DMLabelGetValue(subpointMap, cone[c], &val);
2794:           if (val >= 0) coneSizeNew++;
2795:         }
2796:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2797:       }
2798:     }
2799:   }
2800:   DMSetUp(subdm);
2801:   /* Set cones */
2802:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2803:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2804:   for (d = 0; d <= dim; ++d) {
2805:     for (p = 0; p < numSubPoints[d]; ++p) {
2806:       const PetscInt  point    = subpoints[d][p];
2807:       const PetscInt  subpoint = firstSubPoint[d] + p;
2808:       const PetscInt *cone, *ornt;
2809:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

2811:       DMPlexGetConeSize(dm, point, &coneSize);
2812:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2813:       DMPlexGetCone(dm, point, &cone);
2814:       DMPlexGetConeOrientation(dm, point, &ornt);
2815:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2816:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2817:         if (subc >= 0) {
2818:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2819:           orntNew[coneSizeNew] = ornt[c];
2820:           ++coneSizeNew;
2821:         }
2822:       }
2823:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2824:       DMPlexSetCone(subdm, subpoint, coneNew);
2825:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2826:     }
2827:   }
2828:   PetscFree2(coneNew,orntNew);
2829:   DMPlexSymmetrize(subdm);
2830:   DMPlexStratify(subdm);
2831:   /* Build coordinates */
2832:   {
2833:     PetscSection coordSection, subCoordSection;
2834:     Vec          coordinates, subCoordinates;
2835:     PetscScalar *coords, *subCoords;
2836:     PetscInt     numComp, coordSize;
2837:     const char  *name;

2839:     DMGetCoordinateSection(dm, &coordSection);
2840:     DMGetCoordinatesLocal(dm, &coordinates);
2841:     DMGetCoordinateSection(subdm, &subCoordSection);
2842:     PetscSectionSetNumFields(subCoordSection, 1);
2843:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2844:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2845:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2846:     for (v = 0; v < numSubPoints[0]; ++v) {
2847:       const PetscInt vertex    = subpoints[0][v];
2848:       const PetscInt subvertex = firstSubPoint[0]+v;
2849:       PetscInt       dof;

2851:       PetscSectionGetDof(coordSection, vertex, &dof);
2852:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2853:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2854:     }
2855:     PetscSectionSetUp(subCoordSection);
2856:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2857:     VecCreate(comm, &subCoordinates);
2858:     PetscObjectGetName((PetscObject)coordinates,&name);
2859:     PetscObjectSetName((PetscObject)subCoordinates,name);
2860:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2861:     VecSetType(subCoordinates,VECSTANDARD);
2862:     VecGetArray(coordinates,    &coords);
2863:     VecGetArray(subCoordinates, &subCoords);
2864:     for (v = 0; v < numSubPoints[0]; ++v) {
2865:       const PetscInt vertex    = subpoints[0][v];
2866:       const PetscInt subvertex = firstSubPoint[0]+v;
2867:       PetscInt dof, off, sdof, soff, d;

2869:       PetscSectionGetDof(coordSection, vertex, &dof);
2870:       PetscSectionGetOffset(coordSection, vertex, &off);
2871:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2872:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2873:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2874:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2875:     }
2876:     VecRestoreArray(coordinates,    &coords);
2877:     VecRestoreArray(subCoordinates, &subCoords);
2878:     DMSetCoordinatesLocal(subdm, subCoordinates);
2879:     VecDestroy(&subCoordinates);
2880:   }
2881:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
2882:   {
2883:     PetscSF            sfPoint, sfPointSub;
2884:     IS                 subpIS;
2885:     const PetscSFNode *remotePoints;
2886:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
2887:     const PetscInt    *localPoints, *subpoints;
2888:     PetscInt          *slocalPoints;
2889:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
2890:     PetscMPIInt        rank;

2892:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2893:     DMGetPointSF(dm, &sfPoint);
2894:     DMGetPointSF(subdm, &sfPointSub);
2895:     DMPlexGetChart(dm, &pStart, &pEnd);
2896:     DMPlexGetChart(subdm, NULL, &numSubroots);
2897:     DMPlexCreateSubpointIS(subdm, &subpIS);
2898:     if (subpIS) {
2899:       ISGetIndices(subpIS, &subpoints);
2900:       ISGetLocalSize(subpIS, &numSubpoints);
2901:     }
2902:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2903:     if (numRoots >= 0) {
2904:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2905:       for (p = 0; p < pEnd-pStart; ++p) {
2906:         newLocalPoints[p].rank  = -2;
2907:         newLocalPoints[p].index = -2;
2908:       }
2909:       /* Set subleaves */
2910:       for (l = 0; l < numLeaves; ++l) {
2911:         const PetscInt point    = localPoints[l];
2912:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

2914:         if (subpoint < 0) continue;
2915:         newLocalPoints[point-pStart].rank  = rank;
2916:         newLocalPoints[point-pStart].index = subpoint;
2917:         ++numSubleaves;
2918:       }
2919:       /* Must put in owned subpoints */
2920:       for (p = pStart; p < pEnd; ++p) {
2921:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

2923:         if (subpoint < 0) {
2924:           newOwners[p-pStart].rank  = -3;
2925:           newOwners[p-pStart].index = -3;
2926:         } else {
2927:           newOwners[p-pStart].rank  = rank;
2928:           newOwners[p-pStart].index = subpoint;
2929:         }
2930:       }
2931:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2932:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2933:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2934:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2935:       PetscMalloc1(numSubleaves, &slocalPoints);
2936:       PetscMalloc1(numSubleaves, &sremotePoints);
2937:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
2938:         const PetscInt point    = localPoints[l];
2939:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

2941:         if (subpoint < 0) continue;
2942:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
2943:         slocalPoints[sl]        = subpoint;
2944:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
2945:         sremotePoints[sl].index = newLocalPoints[point].index;
2946:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
2947:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
2948:         ++sl;
2949:       }
2950:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
2951:       PetscFree2(newLocalPoints,newOwners);
2952:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
2953:     }
2954:     if (subpIS) {
2955:       ISRestoreIndices(subpIS, &subpoints);
2956:       ISDestroy(&subpIS);
2957:     }
2958:   }
2959:   /* Cleanup */
2960:   for (d = 0; d <= dim; ++d) {
2961:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
2962:     ISDestroy(&subpointIS[d]);
2963:   }
2964:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2965:   return(0);
2966: }

2970: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2971: {

2975:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
2976:   return(0);
2977: }

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

2984:   Input Parameters:
2985: + dm           - The original mesh
2986: . vertexLabel  - The DMLabel marking vertices contained in the surface
2987: - value        - The label value to use

2989:   Output Parameter:
2990: . subdm - The surface mesh

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

2994:   Level: developer

2996: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
2997: @*/
2998: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
2999: {
3000:   PetscInt       dim, depth;

3006:   DMGetDimension(dm, &dim);
3007:   DMPlexGetDepth(dm, &depth);
3008:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3009:   DMSetType(*subdm, DMPLEX);
3010:   DMSetDimension(*subdm, dim-1);
3011:   if (depth == dim) {
3012:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3013:   } else {
3014:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3015:   }
3016:   return(0);
3017: }

3021: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3022: {
3023:   MPI_Comm        comm;
3024:   DMLabel         subpointMap;
3025:   IS              subvertexIS;
3026:   const PetscInt *subVertices;
3027:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3028:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3029:   PetscInt        cMax, c, f;
3030:   PetscErrorCode  ierr;

3033:   PetscObjectGetComm((PetscObject)dm, &comm);
3034:   /* Create subpointMap which marks the submesh */
3035:   DMLabelCreate("subpoint_map", &subpointMap);
3036:   DMPlexSetSubpointMap(subdm, subpointMap);
3037:   DMLabelDestroy(&subpointMap);
3038:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3039:   /* Setup chart */
3040:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3041:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3042:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3043:   DMPlexSetVTKCellHeight(subdm, 1);
3044:   /* Set cone sizes */
3045:   firstSubVertex = numSubCells;
3046:   firstSubFace   = numSubCells+numSubVertices;
3047:   newFacePoint   = firstSubFace;
3048:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3049:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3050:   for (c = 0; c < numSubCells; ++c) {
3051:     DMPlexSetConeSize(subdm, c, 1);
3052:   }
3053:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3054:     DMPlexSetConeSize(subdm, f, nFV);
3055:   }
3056:   DMSetUp(subdm);
3057:   /* Create face cones */
3058:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3059:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3060:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3061:   for (c = 0; c < numSubCells; ++c) {
3062:     const PetscInt  cell    = subCells[c];
3063:     const PetscInt  subcell = c;
3064:     const PetscInt *cone, *cells;
3065:     PetscInt        numCells, subVertex, p, v;

3067:     if (cell < cMax) continue;
3068:     DMPlexGetCone(dm, cell, &cone);
3069:     for (v = 0; v < nFV; ++v) {
3070:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3071:       subface[v] = firstSubVertex+subVertex;
3072:     }
3073:     DMPlexSetCone(subdm, newFacePoint, subface);
3074:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3075:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3076:     /* Not true in parallel
3077:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3078:     for (p = 0; p < numCells; ++p) {
3079:       PetscInt negsubcell;

3081:       if (cells[p] >= cMax) continue;
3082:       /* I know this is a crap search */
3083:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3084:         if (subCells[negsubcell] == cells[p]) break;
3085:       }
3086:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3087:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3088:     }
3089:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3090:     ++newFacePoint;
3091:   }
3092:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3093:   DMPlexSymmetrize(subdm);
3094:   DMPlexStratify(subdm);
3095:   /* Build coordinates */
3096:   {
3097:     PetscSection coordSection, subCoordSection;
3098:     Vec          coordinates, subCoordinates;
3099:     PetscScalar *coords, *subCoords;
3100:     PetscInt     numComp, coordSize, v;
3101:     const char  *name;

3103:     DMGetCoordinateSection(dm, &coordSection);
3104:     DMGetCoordinatesLocal(dm, &coordinates);
3105:     DMGetCoordinateSection(subdm, &subCoordSection);
3106:     PetscSectionSetNumFields(subCoordSection, 1);
3107:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3108:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3109:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3110:     for (v = 0; v < numSubVertices; ++v) {
3111:       const PetscInt vertex    = subVertices[v];
3112:       const PetscInt subvertex = firstSubVertex+v;
3113:       PetscInt       dof;

3115:       PetscSectionGetDof(coordSection, vertex, &dof);
3116:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3117:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3118:     }
3119:     PetscSectionSetUp(subCoordSection);
3120:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3121:     VecCreate(comm, &subCoordinates);
3122:     PetscObjectGetName((PetscObject)coordinates,&name);
3123:     PetscObjectSetName((PetscObject)subCoordinates,name);
3124:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3125:     VecSetType(subCoordinates,VECSTANDARD);
3126:     VecGetArray(coordinates,    &coords);
3127:     VecGetArray(subCoordinates, &subCoords);
3128:     for (v = 0; v < numSubVertices; ++v) {
3129:       const PetscInt vertex    = subVertices[v];
3130:       const PetscInt subvertex = firstSubVertex+v;
3131:       PetscInt       dof, off, sdof, soff, d;

3133:       PetscSectionGetDof(coordSection, vertex, &dof);
3134:       PetscSectionGetOffset(coordSection, vertex, &off);
3135:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3136:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3137:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3138:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3139:     }
3140:     VecRestoreArray(coordinates,    &coords);
3141:     VecRestoreArray(subCoordinates, &subCoords);
3142:     DMSetCoordinatesLocal(subdm, subCoordinates);
3143:     VecDestroy(&subCoordinates);
3144:   }
3145:   /* Build SF */
3146:   CHKMEMQ;
3147:   {
3148:     PetscSF            sfPoint, sfPointSub;
3149:     const PetscSFNode *remotePoints;
3150:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3151:     const PetscInt    *localPoints;
3152:     PetscInt          *slocalPoints;
3153:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3154:     PetscMPIInt        rank;

3156:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3157:     DMGetPointSF(dm, &sfPoint);
3158:     DMGetPointSF(subdm, &sfPointSub);
3159:     DMPlexGetChart(dm, &pStart, &pEnd);
3160:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3161:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3162:     if (numRoots >= 0) {
3163:       /* Only vertices should be shared */
3164:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3165:       for (p = 0; p < pEnd-pStart; ++p) {
3166:         newLocalPoints[p].rank  = -2;
3167:         newLocalPoints[p].index = -2;
3168:       }
3169:       /* Set subleaves */
3170:       for (l = 0; l < numLeaves; ++l) {
3171:         const PetscInt point    = localPoints[l];
3172:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3174:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3175:         if (subPoint < 0) continue;
3176:         newLocalPoints[point-pStart].rank  = rank;
3177:         newLocalPoints[point-pStart].index = subPoint;
3178:         ++numSubLeaves;
3179:       }
3180:       /* Must put in owned subpoints */
3181:       for (p = pStart; p < pEnd; ++p) {
3182:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3184:         if (subPoint < 0) {
3185:           newOwners[p-pStart].rank  = -3;
3186:           newOwners[p-pStart].index = -3;
3187:         } else {
3188:           newOwners[p-pStart].rank  = rank;
3189:           newOwners[p-pStart].index = subPoint;
3190:         }
3191:       }
3192:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3193:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3194:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3195:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3196:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3197:       PetscMalloc1(numSubLeaves, &sremotePoints);
3198:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3199:         const PetscInt point    = localPoints[l];
3200:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3202:         if (subPoint < 0) continue;
3203:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3204:         slocalPoints[sl]        = subPoint;
3205:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3206:         sremotePoints[sl].index = newLocalPoints[point].index;
3207:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3208:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3209:         ++sl;
3210:       }
3211:       PetscFree2(newLocalPoints,newOwners);
3212:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3213:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3214:     }
3215:   }
3216:   CHKMEMQ;
3217:   /* Cleanup */
3218:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3219:   ISDestroy(&subvertexIS);
3220:   PetscFree(subCells);
3221:   return(0);
3222: }

3226: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3227: {
3228:   DMLabel        label = NULL;

3232:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3233:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3234:   return(0);
3235: }

3239: /*
3240:   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.

3242:   Input Parameters:
3243: + dm          - The original mesh
3244: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3245: . label       - A label name, or NULL
3246: - value  - A label value

3248:   Output Parameter:
3249: . subdm - The surface mesh

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

3253:   Level: developer

3255: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3256: */
3257: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3258: {
3259:   PetscInt       dim, depth;

3265:   DMGetDimension(dm, &dim);
3266:   DMPlexGetDepth(dm, &depth);
3267:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3268:   DMSetType(*subdm, DMPLEX);
3269:   DMSetDimension(*subdm, dim-1);
3270:   if (depth == dim) {
3271:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3272:   } else {
3273:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3274:   }
3275:   return(0);
3276: }

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

3283:   Input Parameters:
3284: + dm        - The original mesh
3285: . cellLabel - The DMLabel marking cells contained in the new mesh
3286: - value     - The label value to use

3288:   Output Parameter:
3289: . subdm - The new mesh

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

3293:   Level: developer

3295: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3296: @*/
3297: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3298: {
3299:   PetscInt       dim;

3305:   DMGetDimension(dm, &dim);
3306:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3307:   DMSetType(*subdm, DMPLEX);
3308:   DMSetDimension(*subdm, dim);
3309:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3310:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3311:   return(0);
3312: }

3316: /*@
3317:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3319:   Input Parameter:
3320: . dm - The submesh DM

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

3325:   Level: developer

3327: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3328: @*/
3329: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3330: {
3334:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3335:   return(0);
3336: }

3340: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3341: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3342: {
3343:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3344:   DMLabel        tmp;

3349:   tmp  = mesh->subpointMap;
3350:   mesh->subpointMap = subpointMap;
3351:   ++mesh->subpointMap->refct;
3352:   DMLabelDestroy(&tmp);
3353:   return(0);
3354: }

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

3361:   Input Parameter:
3362: . dm - The submesh DM

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

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

3369:   Level: developer

3371: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3372: @*/
3373: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3374: {
3375:   MPI_Comm        comm;
3376:   DMLabel         subpointMap;
3377:   IS              is;
3378:   const PetscInt *opoints;
3379:   PetscInt       *points, *depths;
3380:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3381:   PetscErrorCode  ierr;

3386:   PetscObjectGetComm((PetscObject)dm,&comm);
3387:   *subpointIS = NULL;
3388:   DMPlexGetSubpointMap(dm, &subpointMap);
3389:   DMPlexGetDepth(dm, &depth);
3390:   if (subpointMap && depth >= 0) {
3391:     DMPlexGetChart(dm, &pStart, &pEnd);
3392:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3393:     DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3394:     depths[0] = depth;
3395:     depths[1] = 0;
3396:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3397:     PetscMalloc1(pEnd, &points);
3398:     for(d = 0, off = 0; d <= depth; ++d) {
3399:       const PetscInt dep = depths[d];

3401:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3402:       DMLabelGetStratumSize(subpointMap, dep, &n);
3403:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3404:         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);
3405:       } else {
3406:         if (!n) {
3407:           if (d == 0) {
3408:             /* Missing cells */
3409:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3410:           } else {
3411:             /* Missing faces */
3412:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3413:           }
3414:         }
3415:       }
3416:       if (n) {
3417:         DMLabelGetStratumIS(subpointMap, dep, &is);
3418:         ISGetIndices(is, &opoints);
3419:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3420:         ISRestoreIndices(is, &opoints);
3421:         ISDestroy(&is);
3422:       }
3423:     }
3424:     DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3425:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3426:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3427:   }
3428:   return(0);
3429: }