Actual source code: plexsubmesh.c

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

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

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

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

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

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

 42:   Not Collective

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

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

 51:   Level: developer

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

 60:   DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
 61:   return(0);
 62: }

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

 73:   DMLabelGetNumValues(label, &numValues);
 74:   DMLabelGetValueIS(label, &valueIS);
 75:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 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  q = points[p];
 87:       PetscInt *closure = NULL;
 88:       PetscInt  closureSize, c;

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

113:     DMGetPointSF(dm, &sfPoint);
114:     /* Pull point contributions from remote leaves into local roots */
115:     DMLabelGather(label, sfPoint, &lblLeaves);
116:     DMLabelGetValueIS(lblLeaves, &valueIS);
117:     ISGetLocalSize(valueIS, &numValues);
118:     ISGetIndices(valueIS, &values);
119:     for (v = 0; v < numValues; ++v) {
120:       const PetscInt value = values[v];
121: 
122:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
123:       DMLabelInsertIS(label, pointIS, value);
124:       ISDestroy(&pointIS);
125:     }
126:     ISRestoreIndices(valueIS, &values);
127:     ISDestroy(&valueIS);
128:     DMLabelDestroy(&lblLeaves);
129:     /* Push point contributions from roots into remote leaves */
130:     DMLabelDistribute(label, sfPoint, &lblRoots);
131:     DMLabelGetValueIS(lblRoots, &valueIS);
132:     ISGetLocalSize(valueIS, &numValues);
133:     ISGetIndices(valueIS, &values);
134:     for (v = 0; v < numValues; ++v) {
135:       const PetscInt value = values[v];
136: 
137:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
138:       DMLabelInsertIS(label, pointIS, value);
139:       ISDestroy(&pointIS);
140:     }
141:     ISRestoreIndices(valueIS, &values);
142:     ISDestroy(&valueIS);
143:     DMLabelDestroy(&lblRoots);
144:   }
145:   return(0);
146: }

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

151:   Input Parameters:
152: + dm - The DM
153: - label - A DMLabel marking the surface points

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

158:   Level: developer

160: .seealso: DMPlexLabelCohesiveComplete()
161: @*/
162: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
163: {

167:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
168:   return(0);
169: }

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

174:   Input Parameters:
175: + dm - The DM
176: - label - A DMLabel marking the surface points

178:   Output Parameter:
179: . label - A DMLabel incorporating cells

181:   Level: developer

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

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

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

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

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

228: /*@
229:   DMPlexLabelClearCells - Remove cells from a label

231:   Input Parameters:
232: + dm - The DM
233: - label - A DMLabel marking surface points and their adjacent cells

235:   Output Parameter:
236: . label - A DMLabel without cells

238:   Level: developer

240:   Note: This undoes DMPlexLabelAddCells()

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

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

263:     DMLabelGetStratumSize(label, values[v], &numPoints);
264:     DMLabelGetStratumIS(label, values[v], &pointIS);
265:     ISGetIndices(pointIS, &points);
266:     for (p = 0; p < numPoints; ++p) {
267:       PetscInt point = points[p];

269:       if (point >= cStart && point < cEnd) {
270:         DMLabelClearValue(label,point,values[v]);
271:       }
272:     }
273:     ISRestoreIndices(pointIS, &points);
274:     ISDestroy(&pointIS);
275:   }
276:   ISRestoreIndices(valueIS, &values);
277:   ISDestroy(&valueIS);
278:   return(0);
279: }

281: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
282:  * index (skipping first, which is (0,0)) */
283: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
284: {
285:   PetscInt d, off = 0;

288:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
289:   for (d = 0; d < depth; d++) {
290:     PetscInt firstd = d;
291:     PetscInt firstStart = depthShift[2*d];
292:     PetscInt e;

294:     for (e = d+1; e <= depth; e++) {
295:       if (depthShift[2*e] < firstStart) {
296:         firstd = e;
297:         firstStart = depthShift[2*d];
298:       }
299:     }
300:     if (firstd != d) {
301:       PetscInt swap[2];

303:       e = firstd;
304:       swap[0] = depthShift[2*d];
305:       swap[1] = depthShift[2*d+1];
306:       depthShift[2*d]   = depthShift[2*e];
307:       depthShift[2*d+1] = depthShift[2*e+1];
308:       depthShift[2*e]   = swap[0];
309:       depthShift[2*e+1] = swap[1];
310:     }
311:   }
312:   /* convert (oldstart, added) to (oldstart, newstart) */
313:   for (d = 0; d <= depth; d++) {
314:     off += depthShift[2*d+1];
315:     depthShift[2*d+1] = depthShift[2*d] + off;
316:   }
317:   return(0);
318: }

320: /* depthShift is a list of (old, new) pairs */
321: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
322: {
323:   PetscInt d;
324:   PetscInt newOff = 0;

326:   for (d = 0; d <= depth; d++) {
327:     if (p < depthShift[2*d]) return p + newOff;
328:     else newOff = depthShift[2*d+1] - depthShift[2*d];
329:   }
330:   return p + newOff;
331: }

333: /* depthShift is a list of (old, new) pairs */
334: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
335: {
336:   PetscInt d;
337:   PetscInt newOff = 0;

339:   for (d = 0; d <= depth; d++) {
340:     if (p < depthShift[2*d+1]) return p + newOff;
341:     else newOff = depthShift[2*d] - depthShift[2*d+1];
342:   }
343:   return p + newOff;
344: }

346: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
347: {
348:   PetscInt       depth = 0, d, pStart, pEnd, p;
349:   DMLabel        depthLabel;

353:   DMPlexGetDepth(dm, &depth);
354:   if (depth < 0) return(0);
355:   /* Step 1: Expand chart */
356:   DMPlexGetChart(dm, &pStart, &pEnd);
357:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
358:   DMPlexSetChart(dmNew, pStart, pEnd);
359:   DMCreateLabel(dmNew,"depth");
360:   DMPlexGetDepthLabel(dmNew,&depthLabel);
361:   /* Step 2: Set cone and support sizes */
362:   for (d = 0; d <= depth; ++d) {
363:     PetscInt pStartNew, pEndNew;
364:     IS pIS;

366:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
367:     pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
368:     pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
369:     ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
370:     DMLabelSetStratumIS(depthLabel, d, pIS);
371:     ISDestroy(&pIS);
372:     for (p = pStart; p < pEnd; ++p) {
373:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
374:       PetscInt size;

376:       DMPlexGetConeSize(dm, p, &size);
377:       DMPlexSetConeSize(dmNew, newp, size);
378:       DMPlexGetSupportSize(dm, p, &size);
379:       DMPlexSetSupportSize(dmNew, newp, size);
380:     }
381:   }
382:   return(0);
383: }

385: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
386: {
387:   PetscInt      *newpoints;
388:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

392:   DMPlexGetDepth(dm, &depth);
393:   if (depth < 0) return(0);
394:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
395:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
396:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
397:   /* Step 5: Set cones and supports */
398:   DMPlexGetChart(dm, &pStart, &pEnd);
399:   for (p = pStart; p < pEnd; ++p) {
400:     const PetscInt *points = NULL, *orientations = NULL;
401:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

403:     DMPlexGetConeSize(dm, p, &size);
404:     DMPlexGetCone(dm, p, &points);
405:     DMPlexGetConeOrientation(dm, p, &orientations);
406:     for (i = 0; i < size; ++i) {
407:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
408:     }
409:     DMPlexSetCone(dmNew, newp, newpoints);
410:     DMPlexSetConeOrientation(dmNew, newp, orientations);
411:     DMPlexGetSupportSize(dm, p, &size);
412:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
413:     DMPlexGetSupport(dm, p, &points);
414:     for (i = 0; i < size; ++i) {
415:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
416:     }
417:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
418:     DMPlexSetSupport(dmNew, newp, newpoints);
419:   }
420:   PetscFree(newpoints);
421:   return(0);
422: }

424: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
425: {
426:   PetscSection   coordSection, newCoordSection;
427:   Vec            coordinates, newCoordinates;
428:   PetscScalar   *coords, *newCoords;
429:   PetscInt       coordSize, sStart, sEnd;
430:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
431:   PetscBool      hasCells;

435:   DMGetCoordinateDim(dm, &dim);
436:   DMSetCoordinateDim(dmNew, dim);
437:   DMPlexGetDepth(dm, &depth);
438:   /* Step 8: Convert coordinates */
439:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
440:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
441:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
442:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
443:   DMGetCoordinateSection(dm, &coordSection);
444:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
445:   PetscSectionSetNumFields(newCoordSection, 1);
446:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
447:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
448:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
449:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
450:   if (hasCells) {
451:     for (c = cStart; c < cEnd; ++c) {
452:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

454:       PetscSectionGetDof(coordSection, c, &dof);
455:       PetscSectionSetDof(newCoordSection, cNew, dof);
456:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
457:     }
458:   }
459:   for (v = vStartNew; v < vEndNew; ++v) {
460:     PetscSectionSetDof(newCoordSection, v, dim);
461:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
462:   }
463:   PetscSectionSetUp(newCoordSection);
464:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
465:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
466:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
467:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
468:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
469:   VecSetBlockSize(newCoordinates, dim);
470:   VecSetType(newCoordinates,VECSTANDARD);
471:   DMSetCoordinatesLocal(dmNew, newCoordinates);
472:   DMGetCoordinatesLocal(dm, &coordinates);
473:   VecGetArray(coordinates, &coords);
474:   VecGetArray(newCoordinates, &newCoords);
475:   if (hasCells) {
476:     for (c = cStart; c < cEnd; ++c) {
477:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

479:       PetscSectionGetDof(coordSection, c, &dof);
480:       PetscSectionGetOffset(coordSection, c, &off);
481:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
482:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
483:     }
484:   }
485:   for (v = vStart; v < vEnd; ++v) {
486:     PetscInt dof, off, noff, d;

488:     PetscSectionGetDof(coordSection, v, &dof);
489:     PetscSectionGetOffset(coordSection, v, &off);
490:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
491:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
492:   }
493:   VecRestoreArray(coordinates, &coords);
494:   VecRestoreArray(newCoordinates, &newCoords);
495:   VecDestroy(&newCoordinates);
496:   PetscSectionDestroy(&newCoordSection);
497:   return(0);
498: }

500: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
501: {
502:   PetscInt           depth = 0;
503:   PetscSF            sfPoint, sfPointNew;
504:   const PetscSFNode *remotePoints;
505:   PetscSFNode       *gremotePoints;
506:   const PetscInt    *localPoints;
507:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
508:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
509:   PetscErrorCode     ierr;

512:   DMPlexGetDepth(dm, &depth);
513:   /* Step 9: Convert pointSF */
514:   DMGetPointSF(dm, &sfPoint);
515:   DMGetPointSF(dmNew, &sfPointNew);
516:   DMPlexGetChart(dm, &pStart, &pEnd);
517:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
518:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
519:   if (numRoots >= 0) {
520:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
521:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
522:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
523:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
524:     PetscMalloc1(numLeaves,    &glocalPoints);
525:     PetscMalloc1(numLeaves, &gremotePoints);
526:     for (l = 0; l < numLeaves; ++l) {
527:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
528:       gremotePoints[l].rank  = remotePoints[l].rank;
529:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
530:     }
531:     PetscFree2(newLocation,newRemoteLocation);
532:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
533:   }
534:   return(0);
535: }

537: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
538: {
539:   PetscSF            sfPoint;
540:   DMLabel            vtkLabel, ghostLabel;
541:   const PetscSFNode *leafRemote;
542:   const PetscInt    *leafLocal;
543:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
544:   PetscMPIInt        rank;
545:   PetscErrorCode     ierr;

548:   DMPlexGetDepth(dm, &depth);
549:   /* Step 10: Convert labels */
550:   DMGetNumLabels(dm, &numLabels);
551:   for (l = 0; l < numLabels; ++l) {
552:     DMLabel         label, newlabel;
553:     const char     *lname;
554:     PetscBool       isDepth, isDim;
555:     IS              valueIS;
556:     const PetscInt *values;
557:     PetscInt        numValues, val;

559:     DMGetLabelName(dm, l, &lname);
560:     PetscStrcmp(lname, "depth", &isDepth);
561:     if (isDepth) continue;
562:     PetscStrcmp(lname, "dim", &isDim);
563:     if (isDim) continue;
564:     DMCreateLabel(dmNew, lname);
565:     DMGetLabel(dm, lname, &label);
566:     DMGetLabel(dmNew, lname, &newlabel);
567:     DMLabelGetDefaultValue(label,&val);
568:     DMLabelSetDefaultValue(newlabel,val);
569:     DMLabelGetValueIS(label, &valueIS);
570:     ISGetLocalSize(valueIS, &numValues);
571:     ISGetIndices(valueIS, &values);
572:     for (val = 0; val < numValues; ++val) {
573:       IS              pointIS;
574:       const PetscInt *points;
575:       PetscInt        numPoints, p;

577:       DMLabelGetStratumIS(label, values[val], &pointIS);
578:       ISGetLocalSize(pointIS, &numPoints);
579:       ISGetIndices(pointIS, &points);
580:       for (p = 0; p < numPoints; ++p) {
581:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

583:         DMLabelSetValue(newlabel, newpoint, values[val]);
584:       }
585:       ISRestoreIndices(pointIS, &points);
586:       ISDestroy(&pointIS);
587:     }
588:     ISRestoreIndices(valueIS, &values);
589:     ISDestroy(&valueIS);
590:   }
591:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
592:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
593:   DMGetPointSF(dm, &sfPoint);
594:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
595:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
596:   DMCreateLabel(dmNew, "vtk");
597:   DMCreateLabel(dmNew, "ghost");
598:   DMGetLabel(dmNew, "vtk", &vtkLabel);
599:   DMGetLabel(dmNew, "ghost", &ghostLabel);
600:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
601:     for (; c < leafLocal[l] && c < cEnd; ++c) {
602:       DMLabelSetValue(vtkLabel, c, 1);
603:     }
604:     if (leafLocal[l] >= cEnd) break;
605:     if (leafRemote[l].rank == rank) {
606:       DMLabelSetValue(vtkLabel, c, 1);
607:     } else {
608:       DMLabelSetValue(ghostLabel, c, 2);
609:     }
610:   }
611:   for (; c < cEnd; ++c) {
612:     DMLabelSetValue(vtkLabel, c, 1);
613:   }
614:   if (0) {
615:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
616:   }
617:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
618:   for (f = fStart; f < fEnd; ++f) {
619:     PetscInt numCells;

621:     DMPlexGetSupportSize(dmNew, f, &numCells);
622:     if (numCells < 2) {
623:       DMLabelSetValue(ghostLabel, f, 1);
624:     } else {
625:       const PetscInt *cells = NULL;
626:       PetscInt        vA, vB;

628:       DMPlexGetSupport(dmNew, f, &cells);
629:       DMLabelGetValue(vtkLabel, cells[0], &vA);
630:       DMLabelGetValue(vtkLabel, cells[1], &vB);
631:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
632:     }
633:   }
634:   if (0) {
635:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
636:   }
637:   return(0);
638: }

640: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
641: {
642:   DM             refTree;
643:   PetscSection   pSec;
644:   PetscInt       *parents, *childIDs;

648:   DMPlexGetReferenceTree(dm,&refTree);
649:   DMPlexSetReferenceTree(dmNew,refTree);
650:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
651:   if (pSec) {
652:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
653:     PetscInt *childIDsShifted;
654:     PetscSection pSecShifted;

656:     PetscSectionGetChart(pSec,&pStart,&pEnd);
657:     DMPlexGetDepth(dm,&depth);
658:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
659:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
660:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
661:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
662:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
663:     for (p = pStartShifted; p < pEndShifted; p++) {
664:       /* start off assuming no children */
665:       PetscSectionSetDof(pSecShifted,p,0);
666:     }
667:     for (p = pStart; p < pEnd; p++) {
668:       PetscInt dof;
669:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

671:       PetscSectionGetDof(pSec,p,&dof);
672:       PetscSectionSetDof(pSecShifted,pNew,dof);
673:     }
674:     PetscSectionSetUp(pSecShifted);
675:     for (p = pStart; p < pEnd; p++) {
676:       PetscInt dof;
677:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

679:       PetscSectionGetDof(pSec,p,&dof);
680:       if (dof) {
681:         PetscInt off, offNew;

683:         PetscSectionGetOffset(pSec,p,&off);
684:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
685:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
686:         childIDsShifted[offNew] = childIDs[off];
687:       }
688:     }
689:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
690:     PetscFree2(parentsShifted,childIDsShifted);
691:     PetscSectionDestroy(&pSecShifted);
692:   }
693:   return(0);
694: }

696: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
697: {
698:   PetscSF               sf;
699:   IS                    valueIS;
700:   const PetscInt       *values, *leaves;
701:   PetscInt             *depthShift;
702:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
703:   PetscBool             isper;
704:   const PetscReal      *maxCell, *L;
705:   const DMBoundaryType *bd;
706:   PetscErrorCode        ierr;

709:   DMGetPointSF(dm, &sf);
710:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
711:   nleaves = PetscMax(0, nleaves);
712:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
713:   /* Count ghost cells */
714:   DMLabelGetValueIS(label, &valueIS);
715:   ISGetLocalSize(valueIS, &numFS);
716:   ISGetIndices(valueIS, &values);
717:   Ng   = 0;
718:   for (fs = 0; fs < numFS; ++fs) {
719:     IS              faceIS;
720:     const PetscInt *faces;
721:     PetscInt        numFaces, f, numBdFaces = 0;

723:     DMLabelGetStratumIS(label, values[fs], &faceIS);
724:     ISGetLocalSize(faceIS, &numFaces);
725:     ISGetIndices(faceIS, &faces);
726:     for (f = 0; f < numFaces; ++f) {
727:       PetscInt numChildren;

729:       PetscFindInt(faces[f], nleaves, leaves, &loc);
730:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
731:       /* non-local and ancestors points don't get to register ghosts */
732:       if (loc >= 0 || numChildren) continue;
733:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
734:     }
735:     Ng += numBdFaces;
736:     ISDestroy(&faceIS);
737:   }
738:   DMPlexGetDepth(dm, &depth);
739:   PetscMalloc1(2*(depth+1), &depthShift);
740:   for (d = 0; d <= depth; d++) {
741:     PetscInt dEnd;

743:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
744:     depthShift[2*d]   = dEnd;
745:     depthShift[2*d+1] = 0;
746:   }
747:   if (depth >= 0) depthShift[2*depth+1] = Ng;
748:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
749:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
750:   /* Step 3: Set cone/support sizes for new points */
751:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
752:   DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
753:   for (c = cEnd; c < cEnd + Ng; ++c) {
754:     DMPlexSetConeSize(gdm, c, 1);
755:   }
756:   for (fs = 0; fs < numFS; ++fs) {
757:     IS              faceIS;
758:     const PetscInt *faces;
759:     PetscInt        numFaces, f;

761:     DMLabelGetStratumIS(label, values[fs], &faceIS);
762:     ISGetLocalSize(faceIS, &numFaces);
763:     ISGetIndices(faceIS, &faces);
764:     for (f = 0; f < numFaces; ++f) {
765:       PetscInt size, numChildren;

767:       PetscFindInt(faces[f], nleaves, leaves, &loc);
768:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
769:       if (loc >= 0 || numChildren) continue;
770:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
771:       DMPlexGetSupportSize(dm, faces[f], &size);
772:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
773:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
774:     }
775:     ISRestoreIndices(faceIS, &faces);
776:     ISDestroy(&faceIS);
777:   }
778:   /* Step 4: Setup ghosted DM */
779:   DMSetUp(gdm);
780:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
781:   /* Step 6: Set cones and supports for new points */
782:   ghostCell = cEnd;
783:   for (fs = 0; fs < numFS; ++fs) {
784:     IS              faceIS;
785:     const PetscInt *faces;
786:     PetscInt        numFaces, f;

788:     DMLabelGetStratumIS(label, values[fs], &faceIS);
789:     ISGetLocalSize(faceIS, &numFaces);
790:     ISGetIndices(faceIS, &faces);
791:     for (f = 0; f < numFaces; ++f) {
792:       PetscInt newFace = faces[f] + Ng, numChildren;

794:       PetscFindInt(faces[f], nleaves, leaves, &loc);
795:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
796:       if (loc >= 0 || numChildren) continue;
797:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
798:       DMPlexSetCone(gdm, ghostCell, &newFace);
799:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
800:       ++ghostCell;
801:     }
802:     ISRestoreIndices(faceIS, &faces);
803:     ISDestroy(&faceIS);
804:   }
805:   ISRestoreIndices(valueIS, &values);
806:   ISDestroy(&valueIS);
807:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
808:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
809:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
810:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
811:   PetscFree(depthShift);
812:   /* Step 7: Periodicity */
813:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
814:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
815:   if (numGhostCells) *numGhostCells = Ng;
816:   return(0);
817: }

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

822:   Collective on dm

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

828:   Output Parameters:
829: + numGhostCells - The number of ghost cells added to the DM
830: - dmGhosted - The new DM

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

834:   Level: developer

836: .seealso: DMCreate()
837: @*/
838: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
839: {
840:   DM             gdm;
841:   DMLabel        label;
842:   const char    *name = labelName ? labelName : "Face Sets";
843:   PetscInt       dim;
844:   PetscBool      flag;

851:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
852:   DMSetType(gdm, DMPLEX);
853:   DMGetDimension(dm, &dim);
854:   DMSetDimension(gdm, dim);
855:   DMPlexGetAdjacencyUseCone(dm, &flag);
856:   DMPlexSetAdjacencyUseCone(gdm, flag);
857:   DMPlexGetAdjacencyUseClosure(dm, &flag);
858:   DMPlexSetAdjacencyUseClosure(gdm, flag);
859:   DMGetLabel(dm, name, &label);
860:   if (!label) {
861:     /* Get label for boundary faces */
862:     DMCreateLabel(dm, name);
863:     DMGetLabel(dm, name, &label);
864:     DMPlexMarkBoundaryFaces(dm, 1, label);
865:   }
866:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
867:   DMCopyBoundary(dm, gdm);
868:   *dmGhosted = gdm;
869:   return(0);
870: }

872: /*
873:   We are adding three kinds of points here:
874:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
875:     Non-replicated: Points which exist on the fault, but are not replicated
876:     Hybrid:         Entirely new points, such as cohesive cells

878:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
879:   each depth so that the new split/hybrid points can be inserted as a block.
880: */
881: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
882: {
883:   MPI_Comm         comm;
884:   IS               valueIS;
885:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
886:   const PetscInt  *values;          /* List of depths for which we have replicated points */
887:   IS              *splitIS;
888:   IS              *unsplitIS;
889:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
890:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
891:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
892:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
893:   const PetscInt **splitPoints;        /* Replicated points for each depth */
894:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
895:   PetscSection     coordSection;
896:   Vec              coordinates;
897:   PetscScalar     *coords;
898:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
899:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
900:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
901:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
902:   PetscInt        *coneNew, *coneONew, *supportNew;
903:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
904:   PetscErrorCode   ierr;

907:   PetscObjectGetComm((PetscObject)dm,&comm);
908:   DMGetDimension(dm, &dim);
909:   DMPlexGetDepth(dm, &depth);
910:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
911:   /* Count split points and add cohesive cells */
912:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
913:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
914:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
915:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
916:   for (d = 0; d <= depth; ++d) {
917:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
918:     depthEnd[d]           = pMaxNew[d];
919:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
920:     numSplitPoints[d]     = 0;
921:     numUnsplitPoints[d]   = 0;
922:     numHybridPoints[d]    = 0;
923:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
924:     splitPoints[d]        = NULL;
925:     unsplitPoints[d]      = NULL;
926:     splitIS[d]            = NULL;
927:     unsplitIS[d]          = NULL;
928:     /* we are shifting the existing hybrid points with the stratum behind them, so
929:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
930:     depthShift[2*d]       = depthMax[d];
931:     depthShift[2*d+1]     = 0;
932:   }
933:   if (label) {
934:     DMLabelGetValueIS(label, &valueIS);
935:     ISGetLocalSize(valueIS, &numSP);
936:     ISGetIndices(valueIS, &values);
937:   }
938:   for (sp = 0; sp < numSP; ++sp) {
939:     const PetscInt dep = values[sp];

941:     if ((dep < 0) || (dep > depth)) continue;
942:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
943:     if (splitIS[dep]) {
944:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
945:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
946:     }
947:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
948:     if (unsplitIS[dep]) {
949:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
950:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
951:     }
952:   }
953:   /* Calculate number of hybrid points */
954:   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   */
955:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
956:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
957:   /* the end of the points in this stratum that come before the new points:
958:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
959:    * added points */
960:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
961:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
962:   /* Step 3: Set cone/support sizes for new points */
963:   for (dep = 0; dep <= depth; ++dep) {
964:     for (p = 0; p < numSplitPoints[dep]; ++p) {
965:       const PetscInt  oldp   = splitPoints[dep][p];
966:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
967:       const PetscInt  splitp = p    + pMaxNew[dep];
968:       const PetscInt *support;
969:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

971:       DMPlexGetConeSize(dm, oldp, &coneSize);
972:       DMPlexSetConeSize(sdm, splitp, coneSize);
973:       DMPlexGetSupportSize(dm, oldp, &supportSize);
974:       DMPlexSetSupportSize(sdm, splitp, supportSize);
975:       if (dep == depth-1) {
976:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

983:         DMPlexGetSupport(dm, oldp, &support);
984:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
985:           PetscInt val;

987:           DMLabelGetValue(label, support[e], &val);
988:           if (val == 1) ++qf;
989:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
990:           if ((val == 1) || (val == -(shift + 1))) ++qp;
991:         }
992:         /* Split old vertex: Edges into original vertex and new cohesive edge */
993:         DMPlexSetSupportSize(sdm, newp, qn+1);
994:         /* Split new vertex: Edges into split vertex and new cohesive edge */
995:         DMPlexSetSupportSize(sdm, splitp, qp+1);
996:         /* Add hybrid edge */
997:         DMPlexSetConeSize(sdm, hybedge, 2);
998:         DMPlexSetSupportSize(sdm, hybedge, qf);
999:       } else if (dep == dim-2) {
1000:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1002:         DMPlexGetSupport(dm, oldp, &support);
1003:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1004:           PetscInt val;

1006:           DMLabelGetValue(label, support[e], &val);
1007:           if (val == dim-1) ++qf;
1008:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
1009:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1010:         }
1011:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1012:         DMPlexSetSupportSize(sdm, newp, qn+1);
1013:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1014:         DMPlexSetSupportSize(sdm, splitp, qp+1);
1015:         /* Add hybrid face */
1016:         DMPlexSetConeSize(sdm, hybface, 4);
1017:         DMPlexSetSupportSize(sdm, hybface, qf);
1018:       }
1019:     }
1020:   }
1021:   for (dep = 0; dep <= depth; ++dep) {
1022:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1023:       const PetscInt  oldp   = unsplitPoints[dep][p];
1024:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1025:       const PetscInt *support;
1026:       PetscInt        coneSize, supportSize, qf, e, s;

1028:       DMPlexGetConeSize(dm, oldp, &coneSize);
1029:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1030:       DMPlexGetSupport(dm, oldp, &support);
1031:       if (dep == 0) {
1032:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1034:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1035:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1036:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1037:           if (e >= 0) ++qf;
1038:         }
1039:         DMPlexSetSupportSize(sdm, newp, qf+2);
1040:         /* Add hybrid edge */
1041:         DMPlexSetConeSize(sdm, hybedge, 2);
1042:         for (e = 0, qf = 0; e < supportSize; ++e) {
1043:           PetscInt val;

1045:           DMLabelGetValue(label, support[e], &val);
1046:           /* Split and unsplit edges produce hybrid faces */
1047:           if (val == 1) ++qf;
1048:           if (val == (shift2 + 1)) ++qf;
1049:         }
1050:         DMPlexSetSupportSize(sdm, hybedge, qf);
1051:       } else if (dep == dim-2) {
1052:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1053:         PetscInt       val;

1055:         for (e = 0, qf = 0; e < supportSize; ++e) {
1056:           DMLabelGetValue(label, support[e], &val);
1057:           if (val == dim-1) qf += 2;
1058:           else              ++qf;
1059:         }
1060:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1061:         DMPlexSetSupportSize(sdm, newp, qf+2);
1062:         /* Add hybrid face */
1063:         for (e = 0, qf = 0; e < supportSize; ++e) {
1064:           DMLabelGetValue(label, support[e], &val);
1065:           if (val == dim-1) ++qf;
1066:         }
1067:         DMPlexSetConeSize(sdm, hybface, 4);
1068:         DMPlexSetSupportSize(sdm, hybface, qf);
1069:       }
1070:     }
1071:   }
1072:   /* Step 4: Setup split DM */
1073:   DMSetUp(sdm);
1074:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1075:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1076:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1077:   /* Step 6: Set cones and supports for new points */
1078:   for (dep = 0; dep <= depth; ++dep) {
1079:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1080:       const PetscInt  oldp   = splitPoints[dep][p];
1081:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1082:       const PetscInt  splitp = p    + pMaxNew[dep];
1083:       const PetscInt *cone, *support, *ornt;
1084:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1086:       DMPlexGetConeSize(dm, oldp, &coneSize);
1087:       DMPlexGetCone(dm, oldp, &cone);
1088:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1089:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1090:       DMPlexGetSupport(dm, oldp, &support);
1091:       if (dep == depth-1) {
1092:         PetscBool       hasUnsplit = PETSC_FALSE;
1093:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1094:         const PetscInt *supportF;

1096:         /* Split face:       copy in old face to new face to start */
1097:         DMPlexGetSupport(sdm, newp,  &supportF);
1098:         DMPlexSetSupport(sdm, splitp, supportF);
1099:         /* Split old face:   old vertices/edges in cone so no change */
1100:         /* Split new face:   new vertices/edges in cone */
1101:         for (q = 0; q < coneSize; ++q) {
1102:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1103:           if (v < 0) {
1104:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1105:             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);
1106:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1107:             hasUnsplit   = PETSC_TRUE;
1108:           } else {
1109:             coneNew[2+q] = v + pMaxNew[dep-1];
1110:             if (dep > 1) {
1111:               const PetscInt *econe;
1112:               PetscInt        econeSize, r, vs, vu;

1114:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1115:               DMPlexGetCone(dm, cone[q], &econe);
1116:               for (r = 0; r < econeSize; ++r) {
1117:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
1118:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1119:                 if (vs >= 0) continue;
1120:                 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);
1121:                 hasUnsplit   = PETSC_TRUE;
1122:               }
1123:             }
1124:           }
1125:         }
1126:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1127:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1128:         /* Face support */
1129:         for (s = 0; s < supportSize; ++s) {
1130:           PetscInt val;

1132:           DMLabelGetValue(label, support[s], &val);
1133:           if (val < 0) {
1134:             /* Split old face:   Replace negative side cell with cohesive cell */
1135:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1136:           } else {
1137:             /* Split new face:   Replace positive side cell with cohesive cell */
1138:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1139:             /* Get orientation for cohesive face */
1140:             {
1141:               const PetscInt *ncone, *nconeO;
1142:               PetscInt        nconeSize, nc;

1144:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1145:               DMPlexGetCone(dm, support[s], &ncone);
1146:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1147:               for (nc = 0; nc < nconeSize; ++nc) {
1148:                 if (ncone[nc] == oldp) {
1149:                   coneONew[0] = nconeO[nc];
1150:                   break;
1151:                 }
1152:               }
1153:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1154:             }
1155:           }
1156:         }
1157:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1158:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1159:         coneNew[1]  = splitp;
1160:         coneONew[1] = coneONew[0];
1161:         for (q = 0; q < coneSize; ++q) {
1162:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1163:           if (v < 0) {
1164:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1165:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1166:             coneONew[2+q] = 0;
1167:           } else {
1168:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1169:           }
1170:           coneONew[2+q] = 0;
1171:         }
1172:         DMPlexSetCone(sdm, hybcell, coneNew);
1173:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1174:         /* Label the hybrid cells on the boundary of the split */
1175:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1176:       } else if (dep == 0) {
1177:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1183:           DMLabelGetValue(label, support[e], &val);
1184:           if ((val == 1) || (val == (shift + 1))) {
1185:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1186:           }
1187:         }
1188:         supportNew[qn] = hybedge;
1189:         DMPlexSetSupport(sdm, newp, supportNew);
1190:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1191:         for (e = 0, qp = 0; e < supportSize; ++e) {
1192:           PetscInt val, edge;

1194:           DMLabelGetValue(label, support[e], &val);
1195:           if (val == 1) {
1196:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1197:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1198:             supportNew[qp++] = edge + pMaxNew[dep+1];
1199:           } else if (val == -(shift + 1)) {
1200:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1201:           }
1202:         }
1203:         supportNew[qp] = hybedge;
1204:         DMPlexSetSupport(sdm, splitp, supportNew);
1205:         /* Hybrid edge:    Old and new split vertex */
1206:         coneNew[0] = newp;
1207:         coneNew[1] = splitp;
1208:         DMPlexSetCone(sdm, hybedge, coneNew);
1209:         for (e = 0, qf = 0; e < supportSize; ++e) {
1210:           PetscInt val, edge;

1212:           DMLabelGetValue(label, support[e], &val);
1213:           if (val == 1) {
1214:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1215:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1216:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1217:           }
1218:         }
1219:         DMPlexSetSupport(sdm, hybedge, supportNew);
1220:       } else if (dep == dim-2) {
1221:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1223:         /* Split old edge:   old vertices in cone so no change */
1224:         /* Split new edge:   new vertices in cone */
1225:         for (q = 0; q < coneSize; ++q) {
1226:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1227:           if (v < 0) {
1228:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1229:             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);
1230:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1231:           } else {
1232:             coneNew[q] = v + pMaxNew[dep-1];
1233:           }
1234:         }
1235:         DMPlexSetCone(sdm, splitp, coneNew);
1236:         /* Split old edge: Faces in positive side cells and old split faces */
1237:         for (e = 0, q = 0; e < supportSize; ++e) {
1238:           PetscInt val;

1240:           DMLabelGetValue(label, support[e], &val);
1241:           if (val == dim-1) {
1242:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1243:           } else if (val == (shift + dim-1)) {
1244:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1245:           }
1246:         }
1247:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1248:         DMPlexSetSupport(sdm, newp, supportNew);
1249:         /* Split new edge: Faces in negative side cells and new split faces */
1250:         for (e = 0, q = 0; e < supportSize; ++e) {
1251:           PetscInt val, face;

1253:           DMLabelGetValue(label, support[e], &val);
1254:           if (val == dim-1) {
1255:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1256:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1257:             supportNew[q++] = face + pMaxNew[dep+1];
1258:           } else if (val == -(shift + dim-1)) {
1259:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1260:           }
1261:         }
1262:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1263:         DMPlexSetSupport(sdm, splitp, supportNew);
1264:         /* Hybrid face */
1265:         coneNew[0] = newp;
1266:         coneNew[1] = splitp;
1267:         for (v = 0; v < coneSize; ++v) {
1268:           PetscInt vertex;
1269:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1270:           if (vertex < 0) {
1271:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1272:             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);
1273:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1274:           } else {
1275:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1276:           }
1277:         }
1278:         DMPlexSetCone(sdm, hybface, coneNew);
1279:         for (e = 0, qf = 0; e < supportSize; ++e) {
1280:           PetscInt val, face;

1282:           DMLabelGetValue(label, support[e], &val);
1283:           if (val == dim-1) {
1284:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1285:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1286:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1287:           }
1288:         }
1289:         DMPlexSetSupport(sdm, hybface, supportNew);
1290:       }
1291:     }
1292:   }
1293:   for (dep = 0; dep <= depth; ++dep) {
1294:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1295:       const PetscInt  oldp   = unsplitPoints[dep][p];
1296:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1297:       const PetscInt *cone, *support, *ornt;
1298:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1300:       DMPlexGetConeSize(dm, oldp, &coneSize);
1301:       DMPlexGetCone(dm, oldp, &cone);
1302:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1303:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1304:       DMPlexGetSupport(dm, oldp, &support);
1305:       if (dep == 0) {
1306:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1308:         /* Unsplit vertex */
1309:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1310:         for (s = 0, q = 0; s < supportSize; ++s) {
1311:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1312:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1313:           if (e >= 0) {
1314:             supportNew[q++] = e + pMaxNew[dep+1];
1315:           }
1316:         }
1317:         supportNew[q++] = hybedge;
1318:         supportNew[q++] = hybedge;
1319:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1320:         DMPlexSetSupport(sdm, newp, supportNew);
1321:         /* Hybrid edge */
1322:         coneNew[0] = newp;
1323:         coneNew[1] = newp;
1324:         DMPlexSetCone(sdm, hybedge, coneNew);
1325:         for (e = 0, qf = 0; e < supportSize; ++e) {
1326:           PetscInt val, edge;

1328:           DMLabelGetValue(label, support[e], &val);
1329:           if (val == 1) {
1330:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1331:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1332:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1333:           } else if  (val ==  (shift2 + 1)) {
1334:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1335:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1336:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1337:           }
1338:         }
1339:         DMPlexSetSupport(sdm, hybedge, supportNew);
1340:       } else if (dep == dim-2) {
1341:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1347:           DMLabelGetValue(label, support[f], &val);
1348:           if (val == dim-1) {
1349:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1350:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1351:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1352:             supportNew[qf++] = face + pMaxNew[dep+1];
1353:           } else {
1354:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1355:           }
1356:         }
1357:         supportNew[qf++] = hybface;
1358:         supportNew[qf++] = hybface;
1359:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1360:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1361:         DMPlexSetSupport(sdm, newp, supportNew);
1362:         /* Add hybrid face */
1363:         coneNew[0] = newp;
1364:         coneNew[1] = newp;
1365:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1366:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1367:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1368:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1369:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1370:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1371:         DMPlexSetCone(sdm, hybface, coneNew);
1372:         for (f = 0, qf = 0; f < supportSize; ++f) {
1373:           PetscInt val, face;

1375:           DMLabelGetValue(label, support[f], &val);
1376:           if (val == dim-1) {
1377:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1378:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1379:           }
1380:         }
1381:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1382:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1383:         DMPlexSetSupport(sdm, hybface, supportNew);
1384:       }
1385:     }
1386:   }
1387:   /* Step 6b: Replace split points in negative side cones */
1388:   for (sp = 0; sp < numSP; ++sp) {
1389:     PetscInt        dep = values[sp];
1390:     IS              pIS;
1391:     PetscInt        numPoints;
1392:     const PetscInt *points;

1394:     if (dep >= 0) continue;
1395:     DMLabelGetStratumIS(label, dep, &pIS);
1396:     if (!pIS) continue;
1397:     dep  = -dep - shift;
1398:     ISGetLocalSize(pIS, &numPoints);
1399:     ISGetIndices(pIS, &points);
1400:     for (p = 0; p < numPoints; ++p) {
1401:       const PetscInt  oldp = points[p];
1402:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1403:       const PetscInt *cone;
1404:       PetscInt        coneSize, c;
1405:       /* PetscBool       replaced = PETSC_FALSE; */

1407:       /* Negative edge: replace split vertex */
1408:       /* Negative cell: replace split face */
1409:       DMPlexGetConeSize(sdm, newp, &coneSize);
1410:       DMPlexGetCone(sdm, newp, &cone);
1411:       for (c = 0; c < coneSize; ++c) {
1412:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1413:         PetscInt       csplitp, cp, val;

1415:         DMLabelGetValue(label, coldp, &val);
1416:         if (val == dep-1) {
1417:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1418:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1419:           csplitp  = pMaxNew[dep-1] + cp;
1420:           DMPlexInsertCone(sdm, newp, c, csplitp);
1421:           /* replaced = PETSC_TRUE; */
1422:         }
1423:       }
1424:       /* Cells with only a vertex or edge on the submesh have no replacement */
1425:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1426:     }
1427:     ISRestoreIndices(pIS, &points);
1428:     ISDestroy(&pIS);
1429:   }
1430:   /* Step 7: Coordinates */
1431:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1432:   DMGetCoordinateSection(sdm, &coordSection);
1433:   DMGetCoordinatesLocal(sdm, &coordinates);
1434:   VecGetArray(coordinates, &coords);
1435:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1436:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1437:     const PetscInt splitp = pMaxNew[0] + v;
1438:     PetscInt       dof, off, soff, d;

1440:     PetscSectionGetDof(coordSection, newp, &dof);
1441:     PetscSectionGetOffset(coordSection, newp, &off);
1442:     PetscSectionGetOffset(coordSection, splitp, &soff);
1443:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1444:   }
1445:   VecRestoreArray(coordinates, &coords);
1446:   /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1447:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1448:   /* Step 9: Labels */
1449:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1450:   DMGetNumLabels(sdm, &numLabels);
1451:   for (dep = 0; dep <= depth; ++dep) {
1452:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1453:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1454:       const PetscInt splitp = pMaxNew[dep] + p;
1455:       PetscInt       l;

1457:       for (l = 0; l < numLabels; ++l) {
1458:         DMLabel     mlabel;
1459:         const char *lname;
1460:         PetscInt    val;
1461:         PetscBool   isDepth;

1463:         DMGetLabelName(sdm, l, &lname);
1464:         PetscStrcmp(lname, "depth", &isDepth);
1465:         if (isDepth) continue;
1466:         DMGetLabel(sdm, lname, &mlabel);
1467:         DMLabelGetValue(mlabel, newp, &val);
1468:         if (val >= 0) {
1469:           DMLabelSetValue(mlabel, splitp, val);
1470:         }
1471:       }
1472:     }
1473:   }
1474:   for (sp = 0; sp < numSP; ++sp) {
1475:     const PetscInt dep = values[sp];

1477:     if ((dep < 0) || (dep > depth)) continue;
1478:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1479:     ISDestroy(&splitIS[dep]);
1480:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1481:     ISDestroy(&unsplitIS[dep]);
1482:   }
1483:   if (label) {
1484:     ISRestoreIndices(valueIS, &values);
1485:     ISDestroy(&valueIS);
1486:   }
1487:   for (d = 0; d <= depth; ++d) {
1488:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1489:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1490:   }
1491:   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);
1492:   PetscFree3(coneNew, coneONew, supportNew);
1493:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1494:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1495:   return(0);
1496: }

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

1501:   Collective on dm

1503:   Input Parameters:
1504: + dm - The original DM
1505: - label - The label specifying the boundary faces (this could be auto-generated)

1507:   Output Parameters:
1508: - dmSplit - The new DM

1510:   Level: developer

1512: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1513: @*/
1514: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1515: {
1516:   DM             sdm;
1517:   PetscInt       dim;

1523:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1524:   DMSetType(sdm, DMPLEX);
1525:   DMGetDimension(dm, &dim);
1526:   DMSetDimension(sdm, dim);
1527:   switch (dim) {
1528:   case 2:
1529:   case 3:
1530:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1531:     break;
1532:   default:
1533:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1534:   }
1535:   *dmSplit = sdm;
1536:   return(0);
1537: }

1539: /* Returns the side of the surface for a given cell with a face on the surface */
1540: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1541: {
1542:   const PetscInt *cone, *ornt;
1543:   PetscInt        dim, coneSize, c;
1544:   PetscErrorCode  ierr;

1547:   *pos = PETSC_TRUE;
1548:   DMGetDimension(dm, &dim);
1549:   DMPlexGetConeSize(dm, cell, &coneSize);
1550:   DMPlexGetCone(dm, cell, &cone);
1551:   DMPlexGetConeOrientation(dm, cell, &ornt);
1552:   for (c = 0; c < coneSize; ++c) {
1553:     if (cone[c] == face) {
1554:       PetscInt o = ornt[c];

1556:       if (subdm) {
1557:         const PetscInt *subcone, *subornt;
1558:         PetscInt        subpoint, subface, subconeSize, sc;

1560:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1561:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1562:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1563:         DMPlexGetCone(subdm, subpoint, &subcone);
1564:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1565:         for (sc = 0; sc < subconeSize; ++sc) {
1566:           if (subcone[sc] == subface) {
1567:             o = subornt[0];
1568:             break;
1569:           }
1570:         }
1571:         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);
1572:       }
1573:       if (o >= 0) *pos = PETSC_TRUE;
1574:       else        *pos = PETSC_FALSE;
1575:       break;
1576:     }
1577:   }
1578:   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);
1579:   return(0);
1580: }

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

1586:   Input Parameters:
1587: + dm     - The DM
1588: . label  - A DMLabel marking the surface
1589: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1590: . flip   - Flag to flip the submesh normal and replace points on the other side
1591: - subdm  - The subDM associated with the label, or NULL

1593:   Output Parameter:
1594: . label - A DMLabel marking all surface points

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

1598:   Level: developer

1600: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1601: @*/
1602: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1603: {
1604:   DMLabel         depthLabel;
1605:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1606:   const PetscInt *points, *subpoints;
1607:   const PetscInt  rev   = flip ? -1 : 1;
1608:   PetscInt       *pMax;
1609:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1610:   PetscErrorCode  ierr;

1613:   DMPlexGetDepth(dm, &depth);
1614:   DMGetDimension(dm, &dim);
1615:   pSize = PetscMax(depth, dim) + 1;
1616:   PetscMalloc1(pSize,&pMax);
1617:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1618:   DMPlexGetDepthLabel(dm, &depthLabel);
1619:   DMGetDimension(dm, &dim);
1620:   if (subdm) {
1621:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1622:     if (subpointIS) {
1623:       ISGetLocalSize(subpointIS, &numSubpoints);
1624:       ISGetIndices(subpointIS, &subpoints);
1625:     }
1626:   }
1627:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1628:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1629:   if (!dimIS) {
1630:     PetscFree(pMax);
1631:     ISDestroy(&subpointIS);
1632:     return(0);
1633:   }
1634:   ISGetLocalSize(dimIS, &numPoints);
1635:   ISGetIndices(dimIS, &points);
1636:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1637:     const PetscInt *support;
1638:     PetscInt        supportSize, s;

1640:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1641:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1642:     DMPlexGetSupport(dm, points[p], &support);
1643:     for (s = 0; s < supportSize; ++s) {
1644:       const PetscInt *cone;
1645:       PetscInt        coneSize, c;
1646:       PetscBool       pos;

1648:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1649:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1650:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1651:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1652:       /* Put faces touching the fault in the label */
1653:       DMPlexGetConeSize(dm, support[s], &coneSize);
1654:       DMPlexGetCone(dm, support[s], &cone);
1655:       for (c = 0; c < coneSize; ++c) {
1656:         const PetscInt point = cone[c];

1658:         DMLabelGetValue(label, point, &val);
1659:         if (val == -1) {
1660:           PetscInt *closure = NULL;
1661:           PetscInt  closureSize, cl;

1663:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1664:           for (cl = 0; cl < closureSize*2; cl += 2) {
1665:             const PetscInt clp  = closure[cl];
1666:             PetscInt       bval = -1;

1668:             DMLabelGetValue(label, clp, &val);
1669:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1670:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1671:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1672:               break;
1673:             }
1674:           }
1675:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1676:         }
1677:       }
1678:     }
1679:   }
1680:   ISRestoreIndices(dimIS, &points);
1681:   ISDestroy(&dimIS);
1682:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1683:   ISDestroy(&subpointIS);
1684:   /* Mark boundary points as unsplit */
1685:   if (blabel) {
1686:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1687:     ISGetLocalSize(dimIS, &numPoints);
1688:     ISGetIndices(dimIS, &points);
1689:     for (p = 0; p < numPoints; ++p) {
1690:       const PetscInt point = points[p];
1691:       PetscInt       val, bval;

1693:       DMLabelGetValue(blabel, point, &bval);
1694:       if (bval >= 0) {
1695:         DMLabelGetValue(label, point, &val);
1696:         if ((val < 0) || (val > dim)) {
1697:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1698:           DMLabelClearValue(blabel, point, bval);
1699:         }
1700:       }
1701:     }
1702:     for (p = 0; p < numPoints; ++p) {
1703:       const PetscInt point = points[p];
1704:       PetscInt       val, bval;

1706:       DMLabelGetValue(blabel, point, &bval);
1707:       if (bval >= 0) {
1708:         const PetscInt *cone,    *support;
1709:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1711:         /* Mark as unsplit */
1712:         DMLabelGetValue(label, point, &val);
1713:         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);
1714:         DMLabelClearValue(label, point, val);
1715:         DMLabelSetValue(label, point, shift2+val);
1716:         /* Check for cross-edge
1717:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1718:         if (val != 0) continue;
1719:         DMPlexGetSupport(dm, point, &support);
1720:         DMPlexGetSupportSize(dm, point, &supportSize);
1721:         for (s = 0; s < supportSize; ++s) {
1722:           DMPlexGetCone(dm, support[s], &cone);
1723:           DMPlexGetConeSize(dm, support[s], &coneSize);
1724:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1725:           DMLabelGetValue(blabel, cone[0], &valA);
1726:           DMLabelGetValue(blabel, cone[1], &valB);
1727:           DMLabelGetValue(blabel, support[s], &valE);
1728:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1729:         }
1730:       }
1731:     }
1732:     ISRestoreIndices(dimIS, &points);
1733:     ISDestroy(&dimIS);
1734:   }
1735:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1736:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1737:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1738:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1739:   cMax = cMax < 0 ? cEnd : cMax;
1740:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1741:   DMLabelGetStratumIS(label, 0, &dimIS);
1742:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1743:   if (dimIS && crossEdgeIS) {
1744:     IS vertIS = dimIS;

1746:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1747:     ISDestroy(&crossEdgeIS);
1748:     ISDestroy(&vertIS);
1749:   }
1750:   if (!dimIS) {
1751:     PetscFree(pMax);
1752:     return(0);
1753:   }
1754:   ISGetLocalSize(dimIS, &numPoints);
1755:   ISGetIndices(dimIS, &points);
1756:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1757:     PetscInt *star = NULL;
1758:     PetscInt  starSize, s;
1759:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1761:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1762:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1763:     while (again) {
1764:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1765:       again = 0;
1766:       for (s = 0; s < starSize*2; s += 2) {
1767:         const PetscInt  point = star[s];
1768:         const PetscInt *cone;
1769:         PetscInt        coneSize, c;

1771:         if ((point < cStart) || (point >= cMax)) continue;
1772:         DMLabelGetValue(label, point, &val);
1773:         if (val != -1) continue;
1774:         again = again == 1 ? 1 : 2;
1775:         DMPlexGetConeSize(dm, point, &coneSize);
1776:         DMPlexGetCone(dm, point, &cone);
1777:         for (c = 0; c < coneSize; ++c) {
1778:           DMLabelGetValue(label, cone[c], &val);
1779:           if (val != -1) {
1780:             const PetscInt *ccone;
1781:             PetscInt        cconeSize, cc, side;

1783:             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);
1784:             if (val > 0) side =  1;
1785:             else         side = -1;
1786:             DMLabelSetValue(label, point, side*(shift+dim));
1787:             /* Mark cell faces which touch the fault */
1788:             DMPlexGetConeSize(dm, point, &cconeSize);
1789:             DMPlexGetCone(dm, point, &ccone);
1790:             for (cc = 0; cc < cconeSize; ++cc) {
1791:               PetscInt *closure = NULL;
1792:               PetscInt  closureSize, cl;

1794:               DMLabelGetValue(label, ccone[cc], &val);
1795:               if (val != -1) continue;
1796:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1797:               for (cl = 0; cl < closureSize*2; cl += 2) {
1798:                 const PetscInt clp = closure[cl];

1800:                 DMLabelGetValue(label, clp, &val);
1801:                 if (val == -1) continue;
1802:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1803:                 break;
1804:               }
1805:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1806:             }
1807:             again = 1;
1808:             break;
1809:           }
1810:         }
1811:       }
1812:     }
1813:     /* Classify the rest by cell membership */
1814:     for (s = 0; s < starSize*2; s += 2) {
1815:       const PetscInt point = star[s];

1817:       DMLabelGetValue(label, point, &val);
1818:       if (val == -1) {
1819:         PetscInt  *sstar = NULL;
1820:         PetscInt   sstarSize, ss;
1821:         PetscBool  marked = PETSC_FALSE;

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

1827:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1828:           DMLabelGetValue(label, spoint, &val);
1829:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1830:           DMLabelGetValue(depthLabel, point, &dep);
1831:           if (val > 0) {
1832:             DMLabelSetValue(label, point,   shift+dep);
1833:           } else {
1834:             DMLabelSetValue(label, point, -(shift+dep));
1835:           }
1836:           marked = PETSC_TRUE;
1837:           break;
1838:         }
1839:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1840:         DMLabelGetValue(depthLabel, point, &dep);
1841:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1842:       }
1843:     }
1844:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1845:   }
1846:   ISRestoreIndices(dimIS, &points);
1847:   ISDestroy(&dimIS);
1848:   /* If any faces touching the fault divide cells on either side, split them */
1849:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1850:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1851:   ISExpand(facePosIS, faceNegIS, &dimIS);
1852:   ISDestroy(&facePosIS);
1853:   ISDestroy(&faceNegIS);
1854:   ISGetLocalSize(dimIS, &numPoints);
1855:   ISGetIndices(dimIS, &points);
1856:   for (p = 0; p < numPoints; ++p) {
1857:     const PetscInt  point = points[p];
1858:     const PetscInt *support;
1859:     PetscInt        supportSize, valA, valB;

1861:     DMPlexGetSupportSize(dm, point, &supportSize);
1862:     if (supportSize != 2) continue;
1863:     DMPlexGetSupport(dm, point, &support);
1864:     DMLabelGetValue(label, support[0], &valA);
1865:     DMLabelGetValue(label, support[1], &valB);
1866:     if ((valA == -1) || (valB == -1)) continue;
1867:     if (valA*valB > 0) continue;
1868:     /* Split the face */
1869:     DMLabelGetValue(label, point, &valA);
1870:     DMLabelClearValue(label, point, valA);
1871:     DMLabelSetValue(label, point, dim-1);
1872:     /* Label its closure:
1873:       unmarked: label as unsplit
1874:       incident: relabel as split
1875:       split:    do nothing
1876:     */
1877:     {
1878:       PetscInt *closure = NULL;
1879:       PetscInt  closureSize, cl;

1881:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1882:       for (cl = 0; cl < closureSize*2; cl += 2) {
1883:         DMLabelGetValue(label, closure[cl], &valA);
1884:         if (valA == -1) { /* Mark as unsplit */
1885:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1886:           DMLabelSetValue(label, closure[cl], shift2+dep);
1887:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1888:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1889:           DMLabelClearValue(label, closure[cl], valA);
1890:           DMLabelSetValue(label, closure[cl], dep);
1891:         }
1892:       }
1893:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1894:     }
1895:   }
1896:   ISRestoreIndices(dimIS, &points);
1897:   ISDestroy(&dimIS);
1898:   PetscFree(pMax);
1899:   return(0);
1900: }

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

1905:   Collective on dm

1907:   Input Parameters:
1908: + dm - The original DM
1909: - labelName - The label specifying the interface vertices

1911:   Output Parameters:
1912: + hybridLabel - The label fully marking the interface
1913: - dmHybrid - The new DM

1915:   Level: developer

1917: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1918: @*/
1919: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1920: {
1921:   DM             idm;
1922:   DMLabel        subpointMap, hlabel;
1923:   PetscInt       dim;

1930:   DMGetDimension(dm, &dim);
1931:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1932:   DMPlexOrient(idm);
1933:   DMPlexGetSubpointMap(idm, &subpointMap);
1934:   DMLabelDuplicate(subpointMap, &hlabel);
1935:   DMLabelClearStratum(hlabel, dim);
1936:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1937:   DMDestroy(&idm);
1938:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1939:   if (hybridLabel) *hybridLabel = hlabel;
1940:   else             {DMLabelDestroy(&hlabel);}
1941:   return(0);
1942: }

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

1946:      For any marked cell, the marked vertices constitute a single face
1947: */
1948: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1949: {
1950:   IS               subvertexIS = NULL;
1951:   const PetscInt  *subvertices;
1952:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1953:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1954:   PetscErrorCode   ierr;

1957:   *numFaces = 0;
1958:   *nFV      = 0;
1959:   DMPlexGetDepth(dm, &depth);
1960:   DMGetDimension(dm, &dim);
1961:   pSize = PetscMax(depth, dim) + 1;
1962:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1963:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1964:   for (d = 0; d <= depth; ++d) {
1965:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1966:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1967:   }
1968:   /* Loop over initial vertices and mark all faces in the collective star() */
1969:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1970:   if (subvertexIS) {
1971:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1972:     ISGetIndices(subvertexIS, &subvertices);
1973:   }
1974:   for (v = 0; v < numSubVerticesInitial; ++v) {
1975:     const PetscInt vertex = subvertices[v];
1976:     PetscInt      *star   = NULL;
1977:     PetscInt       starSize, s, numCells = 0, c;

1979:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1980:     for (s = 0; s < starSize*2; s += 2) {
1981:       const PetscInt point = star[s];
1982:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1983:     }
1984:     for (c = 0; c < numCells; ++c) {
1985:       const PetscInt cell    = star[c];
1986:       PetscInt      *closure = NULL;
1987:       PetscInt       closureSize, cl;
1988:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1990:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1991:       if (cellLoc == 2) continue;
1992:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1993:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1994:       for (cl = 0; cl < closureSize*2; cl += 2) {
1995:         const PetscInt point = closure[cl];
1996:         PetscInt       vertexLoc;

1998:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1999:           ++numCorners;
2000:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2001:           if (vertexLoc == value) closure[faceSize++] = point;
2002:         }
2003:       }
2004:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2005:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2006:       if (faceSize == *nFV) {
2007:         const PetscInt *cells = NULL;
2008:         PetscInt        numCells, nc;

2010:         ++(*numFaces);
2011:         for (cl = 0; cl < faceSize; ++cl) {
2012:           DMLabelSetValue(subpointMap, closure[cl], 0);
2013:         }
2014:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2015:         for (nc = 0; nc < numCells; ++nc) {
2016:           DMLabelSetValue(subpointMap, cells[nc], 2);
2017:         }
2018:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2019:       }
2020:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2021:     }
2022:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2023:   }
2024:   if (subvertexIS) {
2025:     ISRestoreIndices(subvertexIS, &subvertices);
2026:   }
2027:   ISDestroy(&subvertexIS);
2028:   PetscFree3(pStart,pEnd,pMax);
2029:   return(0);
2030: }

2032: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2033: {
2034:   IS               subvertexIS = NULL;
2035:   const PetscInt  *subvertices;
2036:   PetscInt        *pStart, *pEnd, *pMax;
2037:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2038:   PetscErrorCode   ierr;

2041:   DMGetDimension(dm, &dim);
2042:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2043:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2044:   for (d = 0; d <= dim; ++d) {
2045:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2046:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2047:   }
2048:   /* Loop over initial vertices and mark all faces in the collective star() */
2049:   if (vertexLabel) {
2050:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2051:     if (subvertexIS) {
2052:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2053:       ISGetIndices(subvertexIS, &subvertices);
2054:     }
2055:   }
2056:   for (v = 0; v < numSubVerticesInitial; ++v) {
2057:     const PetscInt vertex = subvertices[v];
2058:     PetscInt      *star   = NULL;
2059:     PetscInt       starSize, s, numFaces = 0, f;

2061:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2062:     for (s = 0; s < starSize*2; s += 2) {
2063:       const PetscInt point = star[s];
2064:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2065:     }
2066:     for (f = 0; f < numFaces; ++f) {
2067:       const PetscInt face    = star[f];
2068:       PetscInt      *closure = NULL;
2069:       PetscInt       closureSize, c;
2070:       PetscInt       faceLoc;

2072:       DMLabelGetValue(subpointMap, face, &faceLoc);
2073:       if (faceLoc == dim-1) continue;
2074:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2075:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2076:       for (c = 0; c < closureSize*2; c += 2) {
2077:         const PetscInt point = closure[c];
2078:         PetscInt       vertexLoc;

2080:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2081:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2082:           if (vertexLoc != value) break;
2083:         }
2084:       }
2085:       if (c == closureSize*2) {
2086:         const PetscInt *support;
2087:         PetscInt        supportSize, s;

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

2092:           for (d = 0; d < dim; ++d) {
2093:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2094:               DMLabelSetValue(subpointMap, point, d);
2095:               break;
2096:             }
2097:           }
2098:         }
2099:         DMPlexGetSupportSize(dm, face, &supportSize);
2100:         DMPlexGetSupport(dm, face, &support);
2101:         for (s = 0; s < supportSize; ++s) {
2102:           DMLabelSetValue(subpointMap, support[s], dim);
2103:         }
2104:       }
2105:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2106:     }
2107:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2108:   }
2109:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2110:   ISDestroy(&subvertexIS);
2111:   PetscFree3(pStart,pEnd,pMax);
2112:   return(0);
2113: }

2115: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2116: {
2117:   DMLabel         label = NULL;
2118:   const PetscInt *cone;
2119:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2120:   PetscErrorCode  ierr;

2123:   *numFaces = 0;
2124:   *nFV = 0;
2125:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2126:   *subCells = NULL;
2127:   DMGetDimension(dm, &dim);
2128:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2129:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2130:   if (cMax < 0) return(0);
2131:   if (label) {
2132:     for (c = cMax; c < cEnd; ++c) {
2133:       PetscInt val;

2135:       DMLabelGetValue(label, c, &val);
2136:       if (val == value) {
2137:         ++(*numFaces);
2138:         DMPlexGetConeSize(dm, c, &coneSize);
2139:       }
2140:     }
2141:   } else {
2142:     *numFaces = cEnd - cMax;
2143:     DMPlexGetConeSize(dm, cMax, &coneSize);
2144:   }
2145:   PetscMalloc1(*numFaces *2, subCells);
2146:   if (!(*numFaces)) return(0);
2147:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2148:   for (c = cMax; c < cEnd; ++c) {
2149:     const PetscInt *cells;
2150:     PetscInt        numCells;

2152:     if (label) {
2153:       PetscInt val;

2155:       DMLabelGetValue(label, c, &val);
2156:       if (val != value) continue;
2157:     }
2158:     DMPlexGetCone(dm, c, &cone);
2159:     for (p = 0; p < *nFV; ++p) {
2160:       DMLabelSetValue(subpointMap, cone[p], 0);
2161:     }
2162:     /* Negative face */
2163:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2164:     /* Not true in parallel
2165:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2166:     for (p = 0; p < numCells; ++p) {
2167:       DMLabelSetValue(subpointMap, cells[p], 2);
2168:       (*subCells)[subc++] = cells[p];
2169:     }
2170:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2171:     /* Positive face is not included */
2172:   }
2173:   return(0);
2174: }

2176: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2177: {
2178:   PetscInt      *pStart, *pEnd;
2179:   PetscInt       dim, cMax, cEnd, c, d;

2183:   DMGetDimension(dm, &dim);
2184:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2185:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2186:   if (cMax < 0) return(0);
2187:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2188:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2189:   for (c = cMax; c < cEnd; ++c) {
2190:     const PetscInt *cone;
2191:     PetscInt       *closure = NULL;
2192:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2194:     if (label) {
2195:       DMLabelGetValue(label, c, &val);
2196:       if (val != value) continue;
2197:     }
2198:     DMPlexGetConeSize(dm, c, &coneSize);
2199:     DMPlexGetCone(dm, c, &cone);
2200:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2201:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2202:     /* Negative face */
2203:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2204:     for (cl = 0; cl < closureSize*2; cl += 2) {
2205:       const PetscInt point = closure[cl];

2207:       for (d = 0; d <= dim; ++d) {
2208:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2209:           DMLabelSetValue(subpointMap, point, d);
2210:           break;
2211:         }
2212:       }
2213:     }
2214:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2215:     /* Cells -- positive face is not included */
2216:     for (cl = 0; cl < 1; ++cl) {
2217:       const PetscInt *support;
2218:       PetscInt        supportSize, s;

2220:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2221:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2222:       DMPlexGetSupport(dm, cone[cl], &support);
2223:       for (s = 0; s < supportSize; ++s) {
2224:         DMLabelSetValue(subpointMap, support[s], dim);
2225:       }
2226:     }
2227:   }
2228:   PetscFree2(pStart, pEnd);
2229:   return(0);
2230: }

2232: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2233: {
2234:   MPI_Comm       comm;
2235:   PetscBool      posOrient = PETSC_FALSE;
2236:   const PetscInt debug     = 0;
2237:   PetscInt       cellDim, faceSize, f;

2241:   PetscObjectGetComm((PetscObject)dm,&comm);
2242:   DMGetDimension(dm, &cellDim);
2243:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2245:   if (cellDim == 1 && numCorners == 2) {
2246:     /* Triangle */
2247:     faceSize  = numCorners-1;
2248:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2249:   } else if (cellDim == 2 && numCorners == 3) {
2250:     /* Triangle */
2251:     faceSize  = numCorners-1;
2252:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2253:   } else if (cellDim == 3 && numCorners == 4) {
2254:     /* Tetrahedron */
2255:     faceSize  = numCorners-1;
2256:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2257:   } else if (cellDim == 1 && numCorners == 3) {
2258:     /* Quadratic line */
2259:     faceSize  = 1;
2260:     posOrient = PETSC_TRUE;
2261:   } else if (cellDim == 2 && numCorners == 4) {
2262:     /* Quads */
2263:     faceSize = 2;
2264:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2265:       posOrient = PETSC_TRUE;
2266:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2267:       posOrient = PETSC_TRUE;
2268:     } else {
2269:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2270:         posOrient = PETSC_FALSE;
2271:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2272:     }
2273:   } else if (cellDim == 2 && numCorners == 6) {
2274:     /* Quadratic triangle (I hate this) */
2275:     /* Edges are determined by the first 2 vertices (corners of edges) */
2276:     const PetscInt faceSizeTri = 3;
2277:     PetscInt       sortedIndices[3], i, iFace;
2278:     PetscBool      found                    = PETSC_FALSE;
2279:     PetscInt       faceVerticesTriSorted[9] = {
2280:       0, 3,  4, /* bottom */
2281:       1, 4,  5, /* right */
2282:       2, 3,  5, /* left */
2283:     };
2284:     PetscInt       faceVerticesTri[9] = {
2285:       0, 3,  4, /* bottom */
2286:       1, 4,  5, /* right */
2287:       2, 5,  3, /* left */
2288:     };

2290:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2291:     PetscSortInt(faceSizeTri, sortedIndices);
2292:     for (iFace = 0; iFace < 3; ++iFace) {
2293:       const PetscInt ii = iFace*faceSizeTri;
2294:       PetscInt       fVertex, cVertex;

2296:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2297:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2298:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2299:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2300:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2301:               faceVertices[fVertex] = origVertices[cVertex];
2302:               break;
2303:             }
2304:           }
2305:         }
2306:         found = PETSC_TRUE;
2307:         break;
2308:       }
2309:     }
2310:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2311:     if (posOriented) *posOriented = PETSC_TRUE;
2312:     return(0);
2313:   } else if (cellDim == 2 && numCorners == 9) {
2314:     /* Quadratic quad (I hate this) */
2315:     /* Edges are determined by the first 2 vertices (corners of edges) */
2316:     const PetscInt faceSizeQuad = 3;
2317:     PetscInt       sortedIndices[3], i, iFace;
2318:     PetscBool      found                      = PETSC_FALSE;
2319:     PetscInt       faceVerticesQuadSorted[12] = {
2320:       0, 1,  4, /* bottom */
2321:       1, 2,  5, /* right */
2322:       2, 3,  6, /* top */
2323:       0, 3,  7, /* left */
2324:     };
2325:     PetscInt       faceVerticesQuad[12] = {
2326:       0, 1,  4, /* bottom */
2327:       1, 2,  5, /* right */
2328:       2, 3,  6, /* top */
2329:       3, 0,  7, /* left */
2330:     };

2332:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2333:     PetscSortInt(faceSizeQuad, sortedIndices);
2334:     for (iFace = 0; iFace < 4; ++iFace) {
2335:       const PetscInt ii = iFace*faceSizeQuad;
2336:       PetscInt       fVertex, cVertex;

2338:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2339:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2340:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2341:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2342:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2343:               faceVertices[fVertex] = origVertices[cVertex];
2344:               break;
2345:             }
2346:           }
2347:         }
2348:         found = PETSC_TRUE;
2349:         break;
2350:       }
2351:     }
2352:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2353:     if (posOriented) *posOriented = PETSC_TRUE;
2354:     return(0);
2355:   } else if (cellDim == 3 && numCorners == 8) {
2356:     /* Hexes
2357:        A hex is two oriented quads with the normal of the first
2358:        pointing up at the second.

2360:           7---6
2361:          /|  /|
2362:         4---5 |
2363:         | 1-|-2
2364:         |/  |/
2365:         0---3

2367:         Faces are determined by the first 4 vertices (corners of faces) */
2368:     const PetscInt faceSizeHex = 4;
2369:     PetscInt       sortedIndices[4], i, iFace;
2370:     PetscBool      found                     = PETSC_FALSE;
2371:     PetscInt       faceVerticesHexSorted[24] = {
2372:       0, 1, 2, 3,  /* bottom */
2373:       4, 5, 6, 7,  /* top */
2374:       0, 3, 4, 5,  /* front */
2375:       2, 3, 5, 6,  /* right */
2376:       1, 2, 6, 7,  /* back */
2377:       0, 1, 4, 7,  /* left */
2378:     };
2379:     PetscInt       faceVerticesHex[24] = {
2380:       1, 2, 3, 0,  /* bottom */
2381:       4, 5, 6, 7,  /* top */
2382:       0, 3, 5, 4,  /* front */
2383:       3, 2, 6, 5,  /* right */
2384:       2, 1, 7, 6,  /* back */
2385:       1, 0, 4, 7,  /* left */
2386:     };

2388:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2389:     PetscSortInt(faceSizeHex, sortedIndices);
2390:     for (iFace = 0; iFace < 6; ++iFace) {
2391:       const PetscInt ii = iFace*faceSizeHex;
2392:       PetscInt       fVertex, cVertex;

2394:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2395:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2396:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2397:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2398:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2399:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2400:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2401:               faceVertices[fVertex] = origVertices[cVertex];
2402:               break;
2403:             }
2404:           }
2405:         }
2406:         found = PETSC_TRUE;
2407:         break;
2408:       }
2409:     }
2410:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2411:     if (posOriented) *posOriented = PETSC_TRUE;
2412:     return(0);
2413:   } else if (cellDim == 3 && numCorners == 10) {
2414:     /* Quadratic tet */
2415:     /* Faces are determined by the first 3 vertices (corners of faces) */
2416:     const PetscInt faceSizeTet = 6;
2417:     PetscInt       sortedIndices[6], i, iFace;
2418:     PetscBool      found                     = PETSC_FALSE;
2419:     PetscInt       faceVerticesTetSorted[24] = {
2420:       0, 1, 2,  6, 7, 8, /* bottom */
2421:       0, 3, 4,  6, 7, 9,  /* front */
2422:       1, 4, 5,  7, 8, 9,  /* right */
2423:       2, 3, 5,  6, 8, 9,  /* left */
2424:     };
2425:     PetscInt       faceVerticesTet[24] = {
2426:       0, 1, 2,  6, 7, 8, /* bottom */
2427:       0, 4, 3,  6, 7, 9,  /* front */
2428:       1, 5, 4,  7, 8, 9,  /* right */
2429:       2, 3, 5,  8, 6, 9,  /* left */
2430:     };

2432:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2433:     PetscSortInt(faceSizeTet, sortedIndices);
2434:     for (iFace=0; iFace < 4; ++iFace) {
2435:       const PetscInt ii = iFace*faceSizeTet;
2436:       PetscInt       fVertex, cVertex;

2438:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2439:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2440:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2441:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2442:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2443:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2444:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2445:               faceVertices[fVertex] = origVertices[cVertex];
2446:               break;
2447:             }
2448:           }
2449:         }
2450:         found = PETSC_TRUE;
2451:         break;
2452:       }
2453:     }
2454:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2455:     if (posOriented) *posOriented = PETSC_TRUE;
2456:     return(0);
2457:   } else if (cellDim == 3 && numCorners == 27) {
2458:     /* Quadratic hexes (I hate this)
2459:        A hex is two oriented quads with the normal of the first
2460:        pointing up at the second.

2462:          7---6
2463:         /|  /|
2464:        4---5 |
2465:        | 3-|-2
2466:        |/  |/
2467:        0---1

2469:        Faces are determined by the first 4 vertices (corners of faces) */
2470:     const PetscInt faceSizeQuadHex = 9;
2471:     PetscInt       sortedIndices[9], i, iFace;
2472:     PetscBool      found                         = PETSC_FALSE;
2473:     PetscInt       faceVerticesQuadHexSorted[54] = {
2474:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2475:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2476:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2477:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2478:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2479:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2480:     };
2481:     PetscInt       faceVerticesQuadHex[54] = {
2482:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2483:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2484:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2485:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2486:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2487:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2488:     };

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

2496:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2497:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2498:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2499:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2500:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2501:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2502:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2503:               faceVertices[fVertex] = origVertices[cVertex];
2504:               break;
2505:             }
2506:           }
2507:         }
2508:         found = PETSC_TRUE;
2509:         break;
2510:       }
2511:     }
2512:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2513:     if (posOriented) *posOriented = PETSC_TRUE;
2514:     return(0);
2515:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2516:   if (!posOrient) {
2517:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2518:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2519:   } else {
2520:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2521:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2522:   }
2523:   if (posOriented) *posOriented = posOrient;
2524:   return(0);
2525: }

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

2531:   Not collective

2533:   Input Parameters:
2534: + dm           - The original mesh
2535: . cell         - The cell mesh point
2536: . faceSize     - The number of vertices on the face
2537: . face         - The face vertices
2538: . numCorners   - The number of vertices on the cell
2539: . indices      - Local numbering of face vertices in cell cone
2540: - origVertices - Original face vertices

2542:   Output Parameter:
2543: + faceVertices - The face vertices properly oriented
2544: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2546:   Level: developer

2548: .seealso: DMPlexGetCone()
2549: @*/
2550: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2551: {
2552:   const PetscInt *cone = NULL;
2553:   PetscInt        coneSize, v, f, v2;
2554:   PetscInt        oppositeVertex = -1;
2555:   PetscErrorCode  ierr;

2558:   DMPlexGetConeSize(dm, cell, &coneSize);
2559:   DMPlexGetCone(dm, cell, &cone);
2560:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2561:     PetscBool found = PETSC_FALSE;

2563:     for (f = 0; f < faceSize; ++f) {
2564:       if (face[f] == cone[v]) {
2565:         found = PETSC_TRUE; break;
2566:       }
2567:     }
2568:     if (found) {
2569:       indices[v2]      = v;
2570:       origVertices[v2] = cone[v];
2571:       ++v2;
2572:     } else {
2573:       oppositeVertex = v;
2574:     }
2575:   }
2576:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2577:   return(0);
2578: }

2580: /*
2581:   DMPlexInsertFace_Internal - Puts a face into the mesh

2583:   Not collective

2585:   Input Parameters:
2586:   + dm              - The DMPlex
2587:   . numFaceVertex   - The number of vertices in the face
2588:   . faceVertices    - The vertices in the face for dm
2589:   . subfaceVertices - The vertices in the face for subdm
2590:   . numCorners      - The number of vertices in the cell
2591:   . cell            - A cell in dm containing the face
2592:   . subcell         - A cell in subdm containing the face
2593:   . firstFace       - First face in the mesh
2594:   - newFacePoint    - Next face in the mesh

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

2599:   Level: developer
2600: */
2601: 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)
2602: {
2603:   MPI_Comm        comm;
2604:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2605:   const PetscInt *faces;
2606:   PetscInt        numFaces, coneSize;
2607:   PetscErrorCode  ierr;

2610:   PetscObjectGetComm((PetscObject)dm,&comm);
2611:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2612:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2613: #if 0
2614:   /* Cannot use this because support() has not been constructed yet */
2615:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2616: #else
2617:   {
2618:     PetscInt f;

2620:     numFaces = 0;
2621:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2622:     for (f = firstFace; f < *newFacePoint; ++f) {
2623:       PetscInt dof, off, d;

2625:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2626:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2627:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2628:       for (d = 0; d < dof; ++d) {
2629:         const PetscInt p = submesh->cones[off+d];
2630:         PetscInt       v;

2632:         for (v = 0; v < numFaceVertices; ++v) {
2633:           if (subfaceVertices[v] == p) break;
2634:         }
2635:         if (v == numFaceVertices) break;
2636:       }
2637:       if (d == dof) {
2638:         numFaces               = 1;
2639:         ((PetscInt*) faces)[0] = f;
2640:       }
2641:     }
2642:   }
2643: #endif
2644:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2645:   else if (numFaces == 1) {
2646:     /* Add the other cell neighbor for this face */
2647:     DMPlexSetCone(subdm, subcell, faces);
2648:   } else {
2649:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2650:     PetscBool posOriented;

2652:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2653:     origVertices        = &orientedVertices[numFaceVertices];
2654:     indices             = &orientedVertices[numFaceVertices*2];
2655:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2656:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2657:     /* TODO: I know that routine should return a permutation, not the indices */
2658:     for (v = 0; v < numFaceVertices; ++v) {
2659:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2660:       for (ov = 0; ov < numFaceVertices; ++ov) {
2661:         if (orientedVertices[ov] == vertex) {
2662:           orientedSubVertices[ov] = subvertex;
2663:           break;
2664:         }
2665:       }
2666:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2667:     }
2668:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2669:     DMPlexSetCone(subdm, subcell, newFacePoint);
2670:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2671:     ++(*newFacePoint);
2672:   }
2673: #if 0
2674:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2675: #else
2676:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2677: #endif
2678:   return(0);
2679: }

2681: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2682: {
2683:   MPI_Comm        comm;
2684:   DMLabel         subpointMap;
2685:   IS              subvertexIS,  subcellIS;
2686:   const PetscInt *subVertices, *subCells;
2687:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2688:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2689:   PetscInt        vStart, vEnd, c, f;
2690:   PetscErrorCode  ierr;

2693:   PetscObjectGetComm((PetscObject)dm,&comm);
2694:   /* Create subpointMap which marks the submesh */
2695:   DMLabelCreate("subpoint_map", &subpointMap);
2696:   DMPlexSetSubpointMap(subdm, subpointMap);
2697:   DMLabelDestroy(&subpointMap);
2698:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2699:   /* Setup chart */
2700:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2701:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2702:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2703:   DMPlexSetVTKCellHeight(subdm, 1);
2704:   /* Set cone sizes */
2705:   firstSubVertex = numSubCells;
2706:   firstSubFace   = numSubCells+numSubVertices;
2707:   newFacePoint   = firstSubFace;
2708:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2709:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2710:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2711:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2712:   for (c = 0; c < numSubCells; ++c) {
2713:     DMPlexSetConeSize(subdm, c, 1);
2714:   }
2715:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2716:     DMPlexSetConeSize(subdm, f, nFV);
2717:   }
2718:   DMSetUp(subdm);
2719:   /* Create face cones */
2720:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2721:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2722:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2723:   for (c = 0; c < numSubCells; ++c) {
2724:     const PetscInt cell    = subCells[c];
2725:     const PetscInt subcell = c;
2726:     PetscInt      *closure = NULL;
2727:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2729:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2730:     for (cl = 0; cl < closureSize*2; cl += 2) {
2731:       const PetscInt point = closure[cl];
2732:       PetscInt       subVertex;

2734:       if ((point >= vStart) && (point < vEnd)) {
2735:         ++numCorners;
2736:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2737:         if (subVertex >= 0) {
2738:           closure[faceSize] = point;
2739:           subface[faceSize] = firstSubVertex+subVertex;
2740:           ++faceSize;
2741:         }
2742:       }
2743:     }
2744:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2745:     if (faceSize == nFV) {
2746:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2747:     }
2748:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2749:   }
2750:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2751:   DMPlexSymmetrize(subdm);
2752:   DMPlexStratify(subdm);
2753:   /* Build coordinates */
2754:   {
2755:     PetscSection coordSection, subCoordSection;
2756:     Vec          coordinates, subCoordinates;
2757:     PetscScalar *coords, *subCoords;
2758:     PetscInt     numComp, coordSize, v;
2759:     const char  *name;

2761:     DMGetCoordinateSection(dm, &coordSection);
2762:     DMGetCoordinatesLocal(dm, &coordinates);
2763:     DMGetCoordinateSection(subdm, &subCoordSection);
2764:     PetscSectionSetNumFields(subCoordSection, 1);
2765:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2766:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2767:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2768:     for (v = 0; v < numSubVertices; ++v) {
2769:       const PetscInt vertex    = subVertices[v];
2770:       const PetscInt subvertex = firstSubVertex+v;
2771:       PetscInt       dof;

2773:       PetscSectionGetDof(coordSection, vertex, &dof);
2774:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2775:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2776:     }
2777:     PetscSectionSetUp(subCoordSection);
2778:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2779:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2780:     PetscObjectGetName((PetscObject)coordinates,&name);
2781:     PetscObjectSetName((PetscObject)subCoordinates,name);
2782:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2783:     VecSetType(subCoordinates,VECSTANDARD);
2784:     if (coordSize) {
2785:       VecGetArray(coordinates,    &coords);
2786:       VecGetArray(subCoordinates, &subCoords);
2787:       for (v = 0; v < numSubVertices; ++v) {
2788:         const PetscInt vertex    = subVertices[v];
2789:         const PetscInt subvertex = firstSubVertex+v;
2790:         PetscInt       dof, off, sdof, soff, d;

2792:         PetscSectionGetDof(coordSection, vertex, &dof);
2793:         PetscSectionGetOffset(coordSection, vertex, &off);
2794:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2795:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2796:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2797:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2798:       }
2799:       VecRestoreArray(coordinates,    &coords);
2800:       VecRestoreArray(subCoordinates, &subCoords);
2801:     }
2802:     DMSetCoordinatesLocal(subdm, subCoordinates);
2803:     VecDestroy(&subCoordinates);
2804:   }
2805:   /* Cleanup */
2806:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2807:   ISDestroy(&subvertexIS);
2808:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2809:   ISDestroy(&subcellIS);
2810:   return(0);
2811: }

2813: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2814: {
2815:   PetscInt       subPoint;

2818:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2819:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2820: }

2822: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2823: {
2824:   MPI_Comm         comm;
2825:   DMLabel          subpointMap;
2826:   IS              *subpointIS;
2827:   const PetscInt **subpoints;
2828:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2829:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2830:   PetscMPIInt      rank;
2831:   PetscErrorCode   ierr;

2834:   PetscObjectGetComm((PetscObject)dm,&comm);
2835:   MPI_Comm_rank(comm, &rank);
2836:   /* Create subpointMap which marks the submesh */
2837:   DMLabelCreate("subpoint_map", &subpointMap);
2838:   DMPlexSetSubpointMap(subdm, subpointMap);
2839:   if (cellHeight) {
2840:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2841:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2842:   } else {
2843:     DMLabel         depth;
2844:     IS              pointIS;
2845:     const PetscInt *points;
2846:     PetscInt        numPoints;

2848:     DMPlexGetDepthLabel(dm, &depth);
2849:     DMLabelGetStratumSize(label, value, &numPoints);
2850:     DMLabelGetStratumIS(label, value, &pointIS);
2851:     ISGetIndices(pointIS, &points);
2852:     for (p = 0; p < numPoints; ++p) {
2853:       PetscInt *closure = NULL;
2854:       PetscInt  closureSize, c, pdim;

2856:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2857:       for (c = 0; c < closureSize*2; c += 2) {
2858:         DMLabelGetValue(depth, closure[c], &pdim);
2859:         DMLabelSetValue(subpointMap, closure[c], pdim);
2860:       }
2861:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2862:     }
2863:     ISRestoreIndices(pointIS, &points);
2864:     ISDestroy(&pointIS);
2865:   }
2866:   DMLabelDestroy(&subpointMap);
2867:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2868:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2869:   cMax = (cMax < 0) ? cEnd : cMax;
2870:   /* Setup chart */
2871:   DMGetDimension(dm, &dim);
2872:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2873:   for (d = 0; d <= dim; ++d) {
2874:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2875:     totSubPoints += numSubPoints[d];
2876:   }
2877:   DMPlexSetChart(subdm, 0, totSubPoints);
2878:   DMPlexSetVTKCellHeight(subdm, cellHeight);
2879:   /* Set cone sizes */
2880:   firstSubPoint[dim] = 0;
2881:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2882:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2883:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2884:   for (d = 0; d <= dim; ++d) {
2885:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2886:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2887:   }
2888:   for (d = 0; d <= dim; ++d) {
2889:     for (p = 0; p < numSubPoints[d]; ++p) {
2890:       const PetscInt  point    = subpoints[d][p];
2891:       const PetscInt  subpoint = firstSubPoint[d] + p;
2892:       const PetscInt *cone;
2893:       PetscInt        coneSize, coneSizeNew, c, val;

2895:       DMPlexGetConeSize(dm, point, &coneSize);
2896:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2897:       if (cellHeight && (d == dim)) {
2898:         DMPlexGetCone(dm, point, &cone);
2899:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2900:           DMLabelGetValue(subpointMap, cone[c], &val);
2901:           if (val >= 0) coneSizeNew++;
2902:         }
2903:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2904:       }
2905:     }
2906:   }
2907:   DMSetUp(subdm);
2908:   /* Set cones */
2909:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2910:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2911:   for (d = 0; d <= dim; ++d) {
2912:     for (p = 0; p < numSubPoints[d]; ++p) {
2913:       const PetscInt  point    = subpoints[d][p];
2914:       const PetscInt  subpoint = firstSubPoint[d] + p;
2915:       const PetscInt *cone, *ornt;
2916:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

2918:       if (d == dim-1) {
2919:         const PetscInt *support, *cone, *ornt;
2920:         PetscInt        supportSize, coneSize, s, subc;

2922:         DMPlexGetSupport(dm, point, &support);
2923:         DMPlexGetSupportSize(dm, point, &supportSize);
2924:         for (s = 0; s < supportSize; ++s) {
2925:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2926:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2927:           if (subc >= 0) {
2928:             const PetscInt ccell = subpoints[d+1][subc];

2930:             DMPlexGetCone(dm, ccell, &cone);
2931:             DMPlexGetConeSize(dm, ccell, &coneSize);
2932:             DMPlexGetConeOrientation(dm, ccell, &ornt);
2933:             for (c = 0; c < coneSize; ++c) {
2934:               if (cone[c] == point) {
2935:                 fornt = ornt[c];
2936:                 break;
2937:               }
2938:             }
2939:             break;
2940:           }
2941:         }
2942:       }
2943:       DMPlexGetConeSize(dm, point, &coneSize);
2944:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2945:       DMPlexGetCone(dm, point, &cone);
2946:       DMPlexGetConeOrientation(dm, point, &ornt);
2947:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2948:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2949:         if (subc >= 0) {
2950:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2951:           orntNew[coneSizeNew] = ornt[c];
2952:           ++coneSizeNew;
2953:         }
2954:       }
2955:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2956:       if (fornt < 0) {
2957:         /* This should be replaced by a call to DMPlexReverseCell() */
2958: #if 0
2959:         DMPlexReverseCell(subdm, subpoint);
2960: #else
2961:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
2962:           PetscInt faceSize, tmp;

2964:           tmp        = coneNew[c];
2965:           coneNew[c] = coneNew[coneSizeNew-1-c];
2966:           coneNew[coneSizeNew-1-c] = tmp;
2967:           DMPlexGetConeSize(dm, cone[c], &faceSize);
2968:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2969:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2970:           orntNew[coneSizeNew-1-c] = tmp;
2971:         }
2972:       }
2973:       DMPlexSetCone(subdm, subpoint, coneNew);
2974:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2975: #endif
2976:     }
2977:   }
2978:   PetscFree2(coneNew,orntNew);
2979:   DMPlexSymmetrize(subdm);
2980:   DMPlexStratify(subdm);
2981:   /* Build coordinates */
2982:   {
2983:     PetscSection coordSection, subCoordSection;
2984:     Vec          coordinates, subCoordinates;
2985:     PetscScalar *coords, *subCoords;
2986:     PetscInt     cdim, numComp, coordSize;
2987:     const char  *name;

2989:     DMGetCoordinateDim(dm, &cdim);
2990:     DMGetCoordinateSection(dm, &coordSection);
2991:     DMGetCoordinatesLocal(dm, &coordinates);
2992:     DMGetCoordinateSection(subdm, &subCoordSection);
2993:     PetscSectionSetNumFields(subCoordSection, 1);
2994:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2995:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2996:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2997:     for (v = 0; v < numSubPoints[0]; ++v) {
2998:       const PetscInt vertex    = subpoints[0][v];
2999:       const PetscInt subvertex = firstSubPoint[0]+v;
3000:       PetscInt       dof;

3002:       PetscSectionGetDof(coordSection, vertex, &dof);
3003:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3004:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3005:     }
3006:     PetscSectionSetUp(subCoordSection);
3007:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3008:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3009:     PetscObjectGetName((PetscObject)coordinates,&name);
3010:     PetscObjectSetName((PetscObject)subCoordinates,name);
3011:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3012:     VecSetBlockSize(subCoordinates, cdim);
3013:     VecSetType(subCoordinates,VECSTANDARD);
3014:     VecGetArray(coordinates,    &coords);
3015:     VecGetArray(subCoordinates, &subCoords);
3016:     for (v = 0; v < numSubPoints[0]; ++v) {
3017:       const PetscInt vertex    = subpoints[0][v];
3018:       const PetscInt subvertex = firstSubPoint[0]+v;
3019:       PetscInt dof, off, sdof, soff, d;

3021:       PetscSectionGetDof(coordSection, vertex, &dof);
3022:       PetscSectionGetOffset(coordSection, vertex, &off);
3023:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3024:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3025:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3026:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3027:     }
3028:     VecRestoreArray(coordinates,    &coords);
3029:     VecRestoreArray(subCoordinates, &subCoords);
3030:     DMSetCoordinatesLocal(subdm, subCoordinates);
3031:     VecDestroy(&subCoordinates);
3032:   }
3033:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3034:   {
3035:     PetscSF            sfPoint, sfPointSub;
3036:     IS                 subpIS;
3037:     const PetscSFNode *remotePoints;
3038:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3039:     const PetscInt    *localPoints, *subpoints;
3040:     PetscInt          *slocalPoints;
3041:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3042:     PetscMPIInt        rank;

3044:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3045:     DMGetPointSF(dm, &sfPoint);
3046:     DMGetPointSF(subdm, &sfPointSub);
3047:     DMPlexGetChart(dm, &pStart, &pEnd);
3048:     DMPlexGetChart(subdm, NULL, &numSubroots);
3049:     DMPlexCreateSubpointIS(subdm, &subpIS);
3050:     if (subpIS) {
3051:       ISGetIndices(subpIS, &subpoints);
3052:       ISGetLocalSize(subpIS, &numSubpoints);
3053:     }
3054:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3055:     if (numRoots >= 0) {
3056:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3057:       for (p = 0; p < pEnd-pStart; ++p) {
3058:         newLocalPoints[p].rank  = -2;
3059:         newLocalPoints[p].index = -2;
3060:       }
3061:       /* Set subleaves */
3062:       for (l = 0; l < numLeaves; ++l) {
3063:         const PetscInt point    = localPoints[l];
3064:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3066:         if (subpoint < 0) continue;
3067:         newLocalPoints[point-pStart].rank  = rank;
3068:         newLocalPoints[point-pStart].index = subpoint;
3069:         ++numSubleaves;
3070:       }
3071:       /* Must put in owned subpoints */
3072:       for (p = pStart; p < pEnd; ++p) {
3073:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3075:         if (subpoint < 0) {
3076:           newOwners[p-pStart].rank  = -3;
3077:           newOwners[p-pStart].index = -3;
3078:         } else {
3079:           newOwners[p-pStart].rank  = rank;
3080:           newOwners[p-pStart].index = subpoint;
3081:         }
3082:       }
3083:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3084:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3085:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3086:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3087:       PetscMalloc1(numSubleaves, &slocalPoints);
3088:       PetscMalloc1(numSubleaves, &sremotePoints);
3089:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3090:         const PetscInt point    = localPoints[l];
3091:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3093:         if (subpoint < 0) continue;
3094:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3095:         slocalPoints[sl]        = subpoint;
3096:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3097:         sremotePoints[sl].index = newLocalPoints[point].index;
3098:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3099:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3100:         ++sl;
3101:       }
3102:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3103:       PetscFree2(newLocalPoints,newOwners);
3104:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3105:     }
3106:     if (subpIS) {
3107:       ISRestoreIndices(subpIS, &subpoints);
3108:       ISDestroy(&subpIS);
3109:     }
3110:   }
3111:   /* Cleanup */
3112:   for (d = 0; d <= dim; ++d) {
3113:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3114:     ISDestroy(&subpointIS[d]);
3115:   }
3116:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3117:   return(0);
3118: }

3120: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3121: {

3125:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3126:   return(0);
3127: }

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

3132:   Input Parameters:
3133: + dm           - The original mesh
3134: . vertexLabel  - The DMLabel marking vertices contained in the surface
3135: - value        - The label value to use

3137:   Output Parameter:
3138: . subdm - The surface mesh

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

3142:   Level: developer

3144: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3145: @*/
3146: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3147: {
3148:   PetscInt       dim, cdim, depth;

3154:   DMGetDimension(dm, &dim);
3155:   DMPlexGetDepth(dm, &depth);
3156:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3157:   DMSetType(*subdm, DMPLEX);
3158:   DMSetDimension(*subdm, dim-1);
3159:   DMGetCoordinateDim(dm, &cdim);
3160:   DMSetCoordinateDim(*subdm, cdim);
3161:   if (depth == dim) {
3162:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3163:   } else {
3164:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3165:   }
3166:   return(0);
3167: }

3169: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3170: {
3171:   MPI_Comm        comm;
3172:   DMLabel         subpointMap;
3173:   IS              subvertexIS;
3174:   const PetscInt *subVertices;
3175:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3176:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3177:   PetscInt        cMax, c, f;
3178:   PetscErrorCode  ierr;

3181:   PetscObjectGetComm((PetscObject)dm, &comm);
3182:   /* Create subpointMap which marks the submesh */
3183:   DMLabelCreate("subpoint_map", &subpointMap);
3184:   DMPlexSetSubpointMap(subdm, subpointMap);
3185:   DMLabelDestroy(&subpointMap);
3186:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3187:   /* Setup chart */
3188:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3189:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3190:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3191:   DMPlexSetVTKCellHeight(subdm, 1);
3192:   /* Set cone sizes */
3193:   firstSubVertex = numSubCells;
3194:   firstSubFace   = numSubCells+numSubVertices;
3195:   newFacePoint   = firstSubFace;
3196:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3197:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3198:   for (c = 0; c < numSubCells; ++c) {
3199:     DMPlexSetConeSize(subdm, c, 1);
3200:   }
3201:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3202:     DMPlexSetConeSize(subdm, f, nFV);
3203:   }
3204:   DMSetUp(subdm);
3205:   /* Create face cones */
3206:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3207:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3208:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3209:   for (c = 0; c < numSubCells; ++c) {
3210:     const PetscInt  cell    = subCells[c];
3211:     const PetscInt  subcell = c;
3212:     const PetscInt *cone, *cells;
3213:     PetscInt        numCells, subVertex, p, v;

3215:     if (cell < cMax) continue;
3216:     DMPlexGetCone(dm, cell, &cone);
3217:     for (v = 0; v < nFV; ++v) {
3218:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3219:       subface[v] = firstSubVertex+subVertex;
3220:     }
3221:     DMPlexSetCone(subdm, newFacePoint, subface);
3222:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3223:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3224:     /* Not true in parallel
3225:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3226:     for (p = 0; p < numCells; ++p) {
3227:       PetscInt negsubcell;

3229:       if (cells[p] >= cMax) continue;
3230:       /* I know this is a crap search */
3231:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3232:         if (subCells[negsubcell] == cells[p]) break;
3233:       }
3234:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3235:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3236:     }
3237:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3238:     ++newFacePoint;
3239:   }
3240:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3241:   DMPlexSymmetrize(subdm);
3242:   DMPlexStratify(subdm);
3243:   /* Build coordinates */
3244:   {
3245:     PetscSection coordSection, subCoordSection;
3246:     Vec          coordinates, subCoordinates;
3247:     PetscScalar *coords, *subCoords;
3248:     PetscInt     cdim, numComp, coordSize, v;
3249:     const char  *name;

3251:     DMGetCoordinateDim(dm, &cdim);
3252:     DMGetCoordinateSection(dm, &coordSection);
3253:     DMGetCoordinatesLocal(dm, &coordinates);
3254:     DMGetCoordinateSection(subdm, &subCoordSection);
3255:     PetscSectionSetNumFields(subCoordSection, 1);
3256:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3257:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3258:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3259:     for (v = 0; v < numSubVertices; ++v) {
3260:       const PetscInt vertex    = subVertices[v];
3261:       const PetscInt subvertex = firstSubVertex+v;
3262:       PetscInt       dof;

3264:       PetscSectionGetDof(coordSection, vertex, &dof);
3265:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3266:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3267:     }
3268:     PetscSectionSetUp(subCoordSection);
3269:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3270:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3271:     PetscObjectGetName((PetscObject)coordinates,&name);
3272:     PetscObjectSetName((PetscObject)subCoordinates,name);
3273:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3274:     VecSetBlockSize(subCoordinates, cdim);
3275:     VecSetType(subCoordinates,VECSTANDARD);
3276:     VecGetArray(coordinates,    &coords);
3277:     VecGetArray(subCoordinates, &subCoords);
3278:     for (v = 0; v < numSubVertices; ++v) {
3279:       const PetscInt vertex    = subVertices[v];
3280:       const PetscInt subvertex = firstSubVertex+v;
3281:       PetscInt       dof, off, sdof, soff, d;

3283:       PetscSectionGetDof(coordSection, vertex, &dof);
3284:       PetscSectionGetOffset(coordSection, vertex, &off);
3285:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3286:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3287:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3288:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3289:     }
3290:     VecRestoreArray(coordinates,    &coords);
3291:     VecRestoreArray(subCoordinates, &subCoords);
3292:     DMSetCoordinatesLocal(subdm, subCoordinates);
3293:     VecDestroy(&subCoordinates);
3294:   }
3295:   /* Build SF */
3296:   CHKMEMQ;
3297:   {
3298:     PetscSF            sfPoint, sfPointSub;
3299:     const PetscSFNode *remotePoints;
3300:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3301:     const PetscInt    *localPoints;
3302:     PetscInt          *slocalPoints;
3303:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3304:     PetscMPIInt        rank;

3306:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3307:     DMGetPointSF(dm, &sfPoint);
3308:     DMGetPointSF(subdm, &sfPointSub);
3309:     DMPlexGetChart(dm, &pStart, &pEnd);
3310:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3311:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3312:     if (numRoots >= 0) {
3313:       /* Only vertices should be shared */
3314:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3315:       for (p = 0; p < pEnd-pStart; ++p) {
3316:         newLocalPoints[p].rank  = -2;
3317:         newLocalPoints[p].index = -2;
3318:       }
3319:       /* Set subleaves */
3320:       for (l = 0; l < numLeaves; ++l) {
3321:         const PetscInt point    = localPoints[l];
3322:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3324:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3325:         if (subPoint < 0) continue;
3326:         newLocalPoints[point-pStart].rank  = rank;
3327:         newLocalPoints[point-pStart].index = subPoint;
3328:         ++numSubLeaves;
3329:       }
3330:       /* Must put in owned subpoints */
3331:       for (p = pStart; p < pEnd; ++p) {
3332:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3334:         if (subPoint < 0) {
3335:           newOwners[p-pStart].rank  = -3;
3336:           newOwners[p-pStart].index = -3;
3337:         } else {
3338:           newOwners[p-pStart].rank  = rank;
3339:           newOwners[p-pStart].index = subPoint;
3340:         }
3341:       }
3342:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3343:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3344:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3345:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3346:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3347:       PetscMalloc1(numSubLeaves, &sremotePoints);
3348:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3349:         const PetscInt point    = localPoints[l];
3350:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3352:         if (subPoint < 0) continue;
3353:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3354:         slocalPoints[sl]        = subPoint;
3355:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3356:         sremotePoints[sl].index = newLocalPoints[point].index;
3357:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3358:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3359:         ++sl;
3360:       }
3361:       PetscFree2(newLocalPoints,newOwners);
3362:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3363:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3364:     }
3365:   }
3366:   CHKMEMQ;
3367:   /* Cleanup */
3368:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3369:   ISDestroy(&subvertexIS);
3370:   PetscFree(subCells);
3371:   return(0);
3372: }

3374: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3375: {
3376:   DMLabel        label = NULL;

3380:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3381:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3382:   return(0);
3383: }

3385: /*@C
3386:   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.

3388:   Input Parameters:
3389: + dm          - The original mesh
3390: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3391: . label       - A label name, or NULL
3392: - value  - A label value

3394:   Output Parameter:
3395: . subdm - The surface mesh

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

3399:   Level: developer

3401: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3402: @*/
3403: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3404: {
3405:   PetscInt       dim, cdim, depth;

3411:   DMGetDimension(dm, &dim);
3412:   DMPlexGetDepth(dm, &depth);
3413:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3414:   DMSetType(*subdm, DMPLEX);
3415:   DMSetDimension(*subdm, dim-1);
3416:   DMGetCoordinateDim(dm, &cdim);
3417:   DMSetCoordinateDim(*subdm, cdim);
3418:   if (depth == dim) {
3419:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3420:   } else {
3421:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3422:   }
3423:   return(0);
3424: }

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

3429:   Input Parameters:
3430: + dm        - The original mesh
3431: . cellLabel - The DMLabel marking cells contained in the new mesh
3432: - value     - The label value to use

3434:   Output Parameter:
3435: . subdm - The new mesh

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

3439:   Level: developer

3441: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3442: @*/
3443: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3444: {
3445:   PetscInt       dim;

3451:   DMGetDimension(dm, &dim);
3452:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3453:   DMSetType(*subdm, DMPLEX);
3454:   DMSetDimension(*subdm, dim);
3455:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3456:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3457:   return(0);
3458: }

3460: /*@
3461:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3463:   Input Parameter:
3464: . dm - The submesh DM

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

3469:   Level: developer

3471: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3472: @*/
3473: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3474: {
3478:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3479:   return(0);
3480: }

3482: /*@
3483:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3485:   Input Parameters:
3486: + dm - The submesh DM
3487: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3491:   Level: developer

3493: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3494: @*/
3495: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3496: {
3497:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3498:   DMLabel        tmp;

3503:   tmp  = mesh->subpointMap;
3504:   mesh->subpointMap = subpointMap;
3505:   ++mesh->subpointMap->refct;
3506:   DMLabelDestroy(&tmp);
3507:   return(0);
3508: }

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

3513:   Input Parameter:
3514: . dm - The submesh DM

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

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

3521:   Level: developer

3523: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3524: @*/
3525: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3526: {
3527:   MPI_Comm        comm;
3528:   DMLabel         subpointMap;
3529:   IS              is;
3530:   const PetscInt *opoints;
3531:   PetscInt       *points, *depths;
3532:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3533:   PetscErrorCode  ierr;

3538:   PetscObjectGetComm((PetscObject)dm,&comm);
3539:   *subpointIS = NULL;
3540:   DMPlexGetSubpointMap(dm, &subpointMap);
3541:   DMPlexGetDepth(dm, &depth);
3542:   if (subpointMap && depth >= 0) {
3543:     DMPlexGetChart(dm, &pStart, &pEnd);
3544:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3545:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3546:     depths[0] = depth;
3547:     depths[1] = 0;
3548:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3549:     PetscMalloc1(pEnd, &points);
3550:     for(d = 0, off = 0; d <= depth; ++d) {
3551:       const PetscInt dep = depths[d];

3553:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3554:       DMLabelGetStratumSize(subpointMap, dep, &n);
3555:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3556:         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);
3557:       } else {
3558:         if (!n) {
3559:           if (d == 0) {
3560:             /* Missing cells */
3561:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3562:           } else {
3563:             /* Missing faces */
3564:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3565:           }
3566:         }
3567:       }
3568:       if (n) {
3569:         DMLabelGetStratumIS(subpointMap, dep, &is);
3570:         ISGetIndices(is, &opoints);
3571:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3572:         ISRestoreIndices(is, &opoints);
3573:         ISDestroy(&is);
3574:       }
3575:     }
3576:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3577:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3578:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3579:   }
3580:   return(0);
3581: }

3583: /*@
3584:   DMPlexGetSubpoint - Return the subpoint corresponding to a point in the original mesh. If the DM
3585:                       is not a submesh, just return the input point.

3587:   Note collective

3589:   Input Parameters:
3590: + dm - The submesh DM
3591: - p  - The point in the original, from which the submesh was created

3593:   Output Parameter:
3594: . subp - The point in the submesh

3596:   Level: developer

3598: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3599: @*/
3600: PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp)
3601: {
3602:   DMLabel        spmap;

3605:   *subp = p;
3606:   DMPlexGetSubpointMap(dm, &spmap);
3607:   if (spmap) {
3608:     IS              subpointIS;
3609:     const PetscInt *subpoints;
3610:     PetscInt        numSubpoints;

3612:     /* TODO Cache the IS, making it look like an index */
3613:     DMPlexCreateSubpointIS(dm, &subpointIS);
3614:     ISGetLocalSize(subpointIS, &numSubpoints);
3615:     ISGetIndices(subpointIS, &subpoints);
3616:     PetscFindInt(p, numSubpoints, subpoints, subp);
3617:     if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p);
3618:     ISRestoreIndices(subpointIS, &subpoints);
3619:     ISDestroy(&subpointIS);
3620:   }
3621:   return(0);
3622: }