Actual source code: plexsubmesh.c

petsc-master 2019-09-16
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:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
615:   for (f = fStart; f < fEnd; ++f) {
616:     PetscInt numCells;

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

625:       DMPlexGetSupport(dmNew, f, &cells);
626:       DMLabelGetValue(vtkLabel, cells[0], &vA);
627:       DMLabelGetValue(vtkLabel, cells[1], &vB);
628:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
629:     }
630:   }
631:   return(0);
632: }

634: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
635: {
636:   DM             refTree;
637:   PetscSection   pSec;
638:   PetscInt       *parents, *childIDs;

642:   DMPlexGetReferenceTree(dm,&refTree);
643:   DMPlexSetReferenceTree(dmNew,refTree);
644:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
645:   if (pSec) {
646:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
647:     PetscInt *childIDsShifted;
648:     PetscSection pSecShifted;

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

665:       PetscSectionGetDof(pSec,p,&dof);
666:       PetscSectionSetDof(pSecShifted,pNew,dof);
667:     }
668:     PetscSectionSetUp(pSecShifted);
669:     for (p = pStart; p < pEnd; p++) {
670:       PetscInt dof;
671:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

673:       PetscSectionGetDof(pSec,p,&dof);
674:       if (dof) {
675:         PetscInt off, offNew;

677:         PetscSectionGetOffset(pSec,p,&off);
678:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
679:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
680:         childIDsShifted[offNew] = childIDs[off];
681:       }
682:     }
683:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
684:     PetscFree2(parentsShifted,childIDsShifted);
685:     PetscSectionDestroy(&pSecShifted);
686:   }
687:   return(0);
688: }

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

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

717:     DMLabelGetStratumIS(label, values[fs], &faceIS);
718:     ISGetLocalSize(faceIS, &numFaces);
719:     ISGetIndices(faceIS, &faces);
720:     for (f = 0; f < numFaces; ++f) {
721:       PetscInt numChildren;

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

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

755:     DMLabelGetStratumIS(label, values[fs], &faceIS);
756:     ISGetLocalSize(faceIS, &numFaces);
757:     ISGetIndices(faceIS, &faces);
758:     for (f = 0; f < numFaces; ++f) {
759:       PetscInt size, numChildren;

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

782:     DMLabelGetStratumIS(label, values[fs], &faceIS);
783:     ISGetLocalSize(faceIS, &numFaces);
784:     ISGetIndices(faceIS, &faces);
785:     for (f = 0; f < numFaces; ++f) {
786:       PetscInt newFace = faces[f] + Ng, numChildren;

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

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

816:   Collective on dm

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

822:   Output Parameters:
823: + numGhostCells - The number of ghost cells added to the DM
824: - dmGhosted - The new DM

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

828:   Level: developer

830: .seealso: DMCreate()
831: @*/
832: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
833: {
834:   DM             gdm;
835:   DMLabel        label;
836:   const char    *name = labelName ? labelName : "Face Sets";
837:   PetscInt       dim;
838:   PetscBool      useCone, useClosure;

845:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
846:   DMSetType(gdm, DMPLEX);
847:   DMGetDimension(dm, &dim);
848:   DMSetDimension(gdm, dim);
849:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
850:   DMSetBasicAdjacency(gdm, useCone,  useClosure);
851:   DMGetLabel(dm, name, &label);
852:   if (!label) {
853:     /* Get label for boundary faces */
854:     DMCreateLabel(dm, name);
855:     DMGetLabel(dm, name, &label);
856:     DMPlexMarkBoundaryFaces(dm, 1, label);
857:   }
858:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
859:   DMCopyBoundary(dm, gdm);
860:   DMCopyDisc(dm, gdm);
861:   gdm->setfromoptionscalled = dm->setfromoptionscalled;
862:   *dmGhosted = gdm;
863:   return(0);
864: }

866: /*
867:   We are adding three kinds of points here:
868:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
869:     Non-replicated: Points which exist on the fault, but are not replicated
870:     Hybrid:         Entirely new points, such as cohesive cells

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

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

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

965:       DMPlexGetConeSize(dm, oldp, &coneSize);
966:       DMPlexSetConeSize(sdm, splitp, coneSize);
967:       DMPlexGetSupportSize(dm, oldp, &supportSize);
968:       DMPlexSetSupportSize(sdm, splitp, supportSize);
969:       if (dep == depth-1) {
970:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

977:         DMPlexGetSupport(dm, oldp, &support);
978:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
979:           PetscInt val;

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

996:         DMPlexGetSupport(dm, oldp, &support);
997:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
998:           PetscInt val;

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

1022:       DMPlexGetConeSize(dm, oldp, &coneSize);
1023:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1024:       DMPlexGetSupport(dm, oldp, &support);
1025:       if (dep == 0) {
1026:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1039:           DMLabelGetValue(label, support[e], &val);
1040:           /* Split and unsplit edges produce hybrid faces */
1041:           if (val == 1) ++qf;
1042:           if (val == (shift2 + 1)) ++qf;
1043:         }
1044:         DMPlexSetSupportSize(sdm, hybedge, qf);
1045:       } else if (dep == dim-2) {
1046:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1047:         PetscInt       val;

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

1080:       DMPlexGetConeSize(dm, oldp, &coneSize);
1081:       DMPlexGetCone(dm, oldp, &cone);
1082:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1083:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1084:       DMPlexGetSupport(dm, oldp, &support);
1085:       if (dep == depth-1) {
1086:         PetscBool       hasUnsplit = PETSC_FALSE;
1087:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1088:         const PetscInt *supportF;

1090:         /* Split face:       copy in old face to new face to start */
1091:         DMPlexGetSupport(sdm, newp,  &supportF);
1092:         DMPlexSetSupport(sdm, splitp, supportF);
1093:         /* Split old face:   old vertices/edges in cone so no change */
1094:         /* Split new face:   new vertices/edges in cone */
1095:         for (q = 0; q < coneSize; ++q) {
1096:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1097:           if (v < 0) {
1098:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1099:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
1100:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1101:             hasUnsplit   = PETSC_TRUE;
1102:           } else {
1103:             coneNew[2+q] = v + pMaxNew[dep-1];
1104:             if (dep > 1) {
1105:               const PetscInt *econe;
1106:               PetscInt        econeSize, r, vs, vu;

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

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

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

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

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

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

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

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

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

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

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

1294:       DMPlexGetConeSize(dm, oldp, &coneSize);
1295:       DMPlexGetCone(dm, oldp, &cone);
1296:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1297:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1298:       DMPlexGetSupport(dm, oldp, &support);
1299:       if (dep == 0) {
1300:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

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

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

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

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

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

1401:       /* Negative edge: replace split vertex */
1402:       /* Negative cell: replace split face */
1403:       DMPlexGetConeSize(sdm, newp, &coneSize);
1404:       DMPlexGetCone(sdm, newp, &cone);
1405:       for (c = 0; c < coneSize; ++c) {
1406:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1407:         PetscInt       csplitp, cp, val;

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

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

1451:       if (splitLabel) {
1452:         const PetscInt val = 100 + dep;

1454:         DMLabelSetValue(splitLabel, newp,    val);
1455:         DMLabelSetValue(splitLabel, splitp, -val);
1456:       }
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: + splitLabel - The label containing the split points, or NULL if no output is desired
1509: - dmSplit - The new DM

1511:   Level: developer

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

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

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

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

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

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

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

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

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

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

1599:   Level: developer

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

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

1641:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1642: #if 0
1643:     if (supportSize != 2) {
1644:       const PetscInt *lp;
1645:       PetscInt        Nlp, pind;

1647:       /* Check that for a cell with a single support face, that face is in the SF */
1648:       /*   THis check only works for the remote side. We would need root side information */
1649:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1650:       PetscFindInt(points[p], Nlp, lp, &pind);
1651:       if (pind < 0) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports, and the face is not shared with another process", points[p], supportSize);
1652:     }
1653: #endif
1654:     DMPlexGetSupport(dm, points[p], &support);
1655:     for (s = 0; s < supportSize; ++s) {
1656:       const PetscInt *cone;
1657:       PetscInt        coneSize, c;
1658:       PetscBool       pos;

1660:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1661:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1662:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1663:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1664:       /* Put faces touching the fault in the label */
1665:       DMPlexGetConeSize(dm, support[s], &coneSize);
1666:       DMPlexGetCone(dm, support[s], &cone);
1667:       for (c = 0; c < coneSize; ++c) {
1668:         const PetscInt point = cone[c];

1670:         DMLabelGetValue(label, point, &val);
1671:         if (val == -1) {
1672:           PetscInt *closure = NULL;
1673:           PetscInt  closureSize, cl;

1675:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1676:           for (cl = 0; cl < closureSize*2; cl += 2) {
1677:             const PetscInt clp  = closure[cl];
1678:             PetscInt       bval = -1;

1680:             DMLabelGetValue(label, clp, &val);
1681:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1682:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1683:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1684:               break;
1685:             }
1686:           }
1687:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1688:         }
1689:       }
1690:     }
1691:   }
1692:   ISRestoreIndices(dimIS, &points);
1693:   ISDestroy(&dimIS);
1694:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1695:   ISDestroy(&subpointIS);
1696:   /* Mark boundary points as unsplit */
1697:   if (blabel) {
1698:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1699:     ISGetLocalSize(dimIS, &numPoints);
1700:     ISGetIndices(dimIS, &points);
1701:     for (p = 0; p < numPoints; ++p) {
1702:       const PetscInt point = points[p];
1703:       PetscInt       val, bval;

1705:       DMLabelGetValue(blabel, point, &bval);
1706:       if (bval >= 0) {
1707:         DMLabelGetValue(label, point, &val);
1708:         if ((val < 0) || (val > dim)) {
1709:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1710:           DMLabelClearValue(blabel, point, bval);
1711:         }
1712:       }
1713:     }
1714:     for (p = 0; p < numPoints; ++p) {
1715:       const PetscInt point = points[p];
1716:       PetscInt       val, bval;

1718:       DMLabelGetValue(blabel, point, &bval);
1719:       if (bval >= 0) {
1720:         const PetscInt *cone,    *support;
1721:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

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

1758:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1759:     ISDestroy(&crossEdgeIS);
1760:     ISDestroy(&vertIS);
1761:   }
1762:   if (!dimIS) {
1763:     PetscFree(pMax);
1764:     return(0);
1765:   }
1766:   ISGetLocalSize(dimIS, &numPoints);
1767:   ISGetIndices(dimIS, &points);
1768:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1769:     PetscInt *star = NULL;
1770:     PetscInt  starSize, s;
1771:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1773:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1774:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1775:     while (again) {
1776:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1777:       again = 0;
1778:       for (s = 0; s < starSize*2; s += 2) {
1779:         const PetscInt  point = star[s];
1780:         const PetscInt *cone;
1781:         PetscInt        coneSize, c;

1783:         if ((point < cStart) || (point >= cMax)) continue;
1784:         DMLabelGetValue(label, point, &val);
1785:         if (val != -1) continue;
1786:         again = again == 1 ? 1 : 2;
1787:         DMPlexGetConeSize(dm, point, &coneSize);
1788:         DMPlexGetCone(dm, point, &cone);
1789:         for (c = 0; c < coneSize; ++c) {
1790:           DMLabelGetValue(label, cone[c], &val);
1791:           if (val != -1) {
1792:             const PetscInt *ccone;
1793:             PetscInt        cconeSize, cc, side;

1795:             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);
1796:             if (val > 0) side =  1;
1797:             else         side = -1;
1798:             DMLabelSetValue(label, point, side*(shift+dim));
1799:             /* Mark cell faces which touch the fault */
1800:             DMPlexGetConeSize(dm, point, &cconeSize);
1801:             DMPlexGetCone(dm, point, &ccone);
1802:             for (cc = 0; cc < cconeSize; ++cc) {
1803:               PetscInt *closure = NULL;
1804:               PetscInt  closureSize, cl;

1806:               DMLabelGetValue(label, ccone[cc], &val);
1807:               if (val != -1) continue;
1808:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1809:               for (cl = 0; cl < closureSize*2; cl += 2) {
1810:                 const PetscInt clp = closure[cl];

1812:                 DMLabelGetValue(label, clp, &val);
1813:                 if (val == -1) continue;
1814:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1815:                 break;
1816:               }
1817:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1818:             }
1819:             again = 1;
1820:             break;
1821:           }
1822:         }
1823:       }
1824:     }
1825:     /* Classify the rest by cell membership */
1826:     for (s = 0; s < starSize*2; s += 2) {
1827:       const PetscInt point = star[s];

1829:       DMLabelGetValue(label, point, &val);
1830:       if (val == -1) {
1831:         PetscInt  *sstar = NULL;
1832:         PetscInt   sstarSize, ss;
1833:         PetscBool  marked = PETSC_FALSE;

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

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

1873:     DMPlexGetSupportSize(dm, point, &supportSize);
1874:     if (supportSize != 2) continue;
1875:     DMPlexGetSupport(dm, point, &support);
1876:     DMLabelGetValue(label, support[0], &valA);
1877:     DMLabelGetValue(label, support[1], &valB);
1878:     if ((valA == -1) || (valB == -1)) continue;
1879:     if (valA*valB > 0) continue;
1880:     /* Split the face */
1881:     DMLabelGetValue(label, point, &valA);
1882:     DMLabelClearValue(label, point, valA);
1883:     DMLabelSetValue(label, point, dim-1);
1884:     /* Label its closure:
1885:       unmarked: label as unsplit
1886:       incident: relabel as split
1887:       split:    do nothing
1888:     */
1889:     {
1890:       PetscInt *closure = NULL;
1891:       PetscInt  closureSize, cl;

1893:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1894:       for (cl = 0; cl < closureSize*2; cl += 2) {
1895:         DMLabelGetValue(label, closure[cl], &valA);
1896:         if (valA == -1) { /* Mark as unsplit */
1897:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1898:           DMLabelSetValue(label, closure[cl], shift2+dep);
1899:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1900:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1901:           DMLabelClearValue(label, closure[cl], valA);
1902:           DMLabelSetValue(label, closure[cl], dep);
1903:         }
1904:       }
1905:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1906:     }
1907:   }
1908:   ISRestoreIndices(dimIS, &points);
1909:   ISDestroy(&dimIS);
1910:   PetscFree(pMax);
1911:   return(0);
1912: }

1914: /* Check that no cell have all vertices on the fault */
1915: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
1916: {
1917:   IS              subpointIS;
1918:   const PetscInt *dmpoints;
1919:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;
1920:   PetscErrorCode  ierr;

1923:   if (!label) return(0);
1924:   DMLabelGetDefaultValue(label, &defaultValue);
1925:   DMPlexCreateSubpointIS(subdm, &subpointIS);
1926:   if (!subpointIS) return(0);
1927:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
1928:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1929:   ISGetIndices(subpointIS, &dmpoints);
1930:   for (c = cStart; c < cEnd; ++c) {
1931:     PetscBool invalidCell = PETSC_TRUE;
1932:     PetscInt *closure     = NULL;
1933:     PetscInt  closureSize, cl;

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

1939:       if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
1940:       DMLabelGetValue(label, closure[cl], &value);
1941:       if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
1942:     }
1943:     DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
1944:     if (invalidCell) {
1945:       ISRestoreIndices(subpointIS, &dmpoints);
1946:       ISDestroy(&subpointIS);
1947:       DMDestroy(&subdm);
1948:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
1949:     }
1950:   }
1951:   ISRestoreIndices(subpointIS, &dmpoints);
1952:   ISDestroy(&subpointIS);
1953:   return(0);
1954: }


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

1960:   Collective on dm

1962:   Input Parameters:
1963: + dm - The original DM
1964: . label - The label specifying the interface vertices
1965: - bdlabel - The optional label specifying the interface boundary vertices

1967:   Output Parameters:
1968: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
1969: . splitLabel - The label containing the split points, or NULL if no output is desired
1970: . dmInterface - The new interface DM, or NULL
1971: - dmHybrid - The new DM with cohesive cells

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

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

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

1986:   Level: developer

1988: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
1989: @*/
1990: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
1991: {
1992:   DM             idm;
1993:   DMLabel        subpointMap, hlabel, slabel = NULL;
1994:   PetscInt       dim;

2004:   DMGetDimension(dm, &dim);
2005:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2006:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
2007:   DMPlexOrient(idm);
2008:   DMPlexGetSubpointMap(idm, &subpointMap);
2009:   DMLabelDuplicate(subpointMap, &hlabel);
2010:   DMLabelClearStratum(hlabel, dim);
2011:   if (splitLabel) {
2012:     const char *name;
2013:     char        sname[PETSC_MAX_PATH_LEN];

2015:     PetscObjectGetName((PetscObject) hlabel, &name);
2016:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2017:     PetscStrcat(sname, " split");
2018:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2019:   }
2020:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2021:   if (dmInterface) {*dmInterface = idm;}
2022:   else             {DMDestroy(&idm);}
2023:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2024:   if (hybridLabel) *hybridLabel = hlabel;
2025:   else             {DMLabelDestroy(&hlabel);}
2026:   if (splitLabel)  *splitLabel  = slabel;
2027:   return(0);
2028: }

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

2032:      For any marked cell, the marked vertices constitute a single face
2033: */
2034: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2035: {
2036:   IS               subvertexIS = NULL;
2037:   const PetscInt  *subvertices;
2038:   PetscInt        *pStart, *pEnd, *pMax, pSize;
2039:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
2040:   PetscErrorCode   ierr;

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

2065:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2066:     for (s = 0; s < starSize*2; s += 2) {
2067:       const PetscInt point = star[s];
2068:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2069:     }
2070:     for (c = 0; c < numCells; ++c) {
2071:       const PetscInt cell    = star[c];
2072:       PetscInt      *closure = NULL;
2073:       PetscInt       closureSize, cl;
2074:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

2076:       DMLabelGetValue(subpointMap, cell, &cellLoc);
2077:       if (cellLoc == 2) continue;
2078:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2079:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2080:       for (cl = 0; cl < closureSize*2; cl += 2) {
2081:         const PetscInt point = closure[cl];
2082:         PetscInt       vertexLoc;

2084:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2085:           ++numCorners;
2086:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2087:           if (vertexLoc == value) closure[faceSize++] = point;
2088:         }
2089:       }
2090:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2091:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2092:       if (faceSize == *nFV) {
2093:         const PetscInt *cells = NULL;
2094:         PetscInt        numCells, nc;

2096:         ++(*numFaces);
2097:         for (cl = 0; cl < faceSize; ++cl) {
2098:           DMLabelSetValue(subpointMap, closure[cl], 0);
2099:         }
2100:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2101:         for (nc = 0; nc < numCells; ++nc) {
2102:           DMLabelSetValue(subpointMap, cells[nc], 2);
2103:         }
2104:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2105:       }
2106:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2107:     }
2108:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2109:   }
2110:   if (subvertexIS) {
2111:     ISRestoreIndices(subvertexIS, &subvertices);
2112:   }
2113:   ISDestroy(&subvertexIS);
2114:   PetscFree3(pStart,pEnd,pMax);
2115:   return(0);
2116: }

2118: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2119: {
2120:   IS               subvertexIS = NULL;
2121:   const PetscInt  *subvertices;
2122:   PetscInt        *pStart, *pEnd, *pMax;
2123:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2124:   PetscErrorCode   ierr;

2127:   DMGetDimension(dm, &dim);
2128:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2129:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2130:   for (d = 0; d <= dim; ++d) {
2131:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2132:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2133:   }
2134:   /* Loop over initial vertices and mark all faces in the collective star() */
2135:   if (vertexLabel) {
2136:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2137:     if (subvertexIS) {
2138:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2139:       ISGetIndices(subvertexIS, &subvertices);
2140:     }
2141:   }
2142:   for (v = 0; v < numSubVerticesInitial; ++v) {
2143:     const PetscInt vertex = subvertices[v];
2144:     PetscInt      *star   = NULL;
2145:     PetscInt       starSize, s, numFaces = 0, f;

2147:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2148:     for (s = 0; s < starSize*2; s += 2) {
2149:       const PetscInt point = star[s];
2150:       PetscInt       faceLoc;

2152:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2153:         if (markedFaces) {
2154:           DMLabelGetValue(vertexLabel, point, &faceLoc);
2155:           if (faceLoc < 0) continue;
2156:         }
2157:         star[numFaces++] = point;
2158:       }
2159:     }
2160:     for (f = 0; f < numFaces; ++f) {
2161:       const PetscInt face    = star[f];
2162:       PetscInt      *closure = NULL;
2163:       PetscInt       closureSize, c;
2164:       PetscInt       faceLoc;

2166:       DMLabelGetValue(subpointMap, face, &faceLoc);
2167:       if (faceLoc == dim-1) continue;
2168:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2169:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2170:       for (c = 0; c < closureSize*2; c += 2) {
2171:         const PetscInt point = closure[c];
2172:         PetscInt       vertexLoc;

2174:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2175:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2176:           if (vertexLoc != value) break;
2177:         }
2178:       }
2179:       if (c == closureSize*2) {
2180:         const PetscInt *support;
2181:         PetscInt        supportSize, s;

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

2186:           for (d = 0; d < dim; ++d) {
2187:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2188:               DMLabelSetValue(subpointMap, point, d);
2189:               break;
2190:             }
2191:           }
2192:         }
2193:         DMPlexGetSupportSize(dm, face, &supportSize);
2194:         DMPlexGetSupport(dm, face, &support);
2195:         for (s = 0; s < supportSize; ++s) {
2196:           DMLabelSetValue(subpointMap, support[s], dim);
2197:         }
2198:       }
2199:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2200:     }
2201:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2202:   }
2203:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2204:   ISDestroy(&subvertexIS);
2205:   PetscFree3(pStart,pEnd,pMax);
2206:   return(0);
2207: }

2209: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2210: {
2211:   DMLabel         label = NULL;
2212:   const PetscInt *cone;
2213:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2214:   PetscErrorCode  ierr;

2217:   *numFaces = 0;
2218:   *nFV = 0;
2219:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2220:   *subCells = NULL;
2221:   DMGetDimension(dm, &dim);
2222:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2223:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2224:   if (cMax < 0) return(0);
2225:   if (label) {
2226:     for (c = cMax; c < cEnd; ++c) {
2227:       PetscInt val;

2229:       DMLabelGetValue(label, c, &val);
2230:       if (val == value) {
2231:         ++(*numFaces);
2232:         DMPlexGetConeSize(dm, c, &coneSize);
2233:       }
2234:     }
2235:   } else {
2236:     *numFaces = cEnd - cMax;
2237:     DMPlexGetConeSize(dm, cMax, &coneSize);
2238:   }
2239:   PetscMalloc1(*numFaces *2, subCells);
2240:   if (!(*numFaces)) return(0);
2241:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2242:   for (c = cMax; c < cEnd; ++c) {
2243:     const PetscInt *cells;
2244:     PetscInt        numCells;

2246:     if (label) {
2247:       PetscInt val;

2249:       DMLabelGetValue(label, c, &val);
2250:       if (val != value) continue;
2251:     }
2252:     DMPlexGetCone(dm, c, &cone);
2253:     for (p = 0; p < *nFV; ++p) {
2254:       DMLabelSetValue(subpointMap, cone[p], 0);
2255:     }
2256:     /* Negative face */
2257:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2258:     /* Not true in parallel
2259:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2260:     for (p = 0; p < numCells; ++p) {
2261:       DMLabelSetValue(subpointMap, cells[p], 2);
2262:       (*subCells)[subc++] = cells[p];
2263:     }
2264:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2265:     /* Positive face is not included */
2266:   }
2267:   return(0);
2268: }

2270: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2271: {
2272:   PetscInt      *pStart, *pEnd;
2273:   PetscInt       dim, cMax, cEnd, c, d;

2277:   DMGetDimension(dm, &dim);
2278:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2279:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2280:   if (cMax < 0) return(0);
2281:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2282:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2283:   for (c = cMax; c < cEnd; ++c) {
2284:     const PetscInt *cone;
2285:     PetscInt       *closure = NULL;
2286:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2288:     if (label) {
2289:       DMLabelGetValue(label, c, &val);
2290:       if (val != value) continue;
2291:     }
2292:     DMPlexGetConeSize(dm, c, &coneSize);
2293:     DMPlexGetCone(dm, c, &cone);
2294:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2295:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2296:     /* Negative face */
2297:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2298:     for (cl = 0; cl < closureSize*2; cl += 2) {
2299:       const PetscInt point = closure[cl];

2301:       for (d = 0; d <= dim; ++d) {
2302:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2303:           DMLabelSetValue(subpointMap, point, d);
2304:           break;
2305:         }
2306:       }
2307:     }
2308:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2309:     /* Cells -- positive face is not included */
2310:     for (cl = 0; cl < 1; ++cl) {
2311:       const PetscInt *support;
2312:       PetscInt        supportSize, s;

2314:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2315:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2316:       DMPlexGetSupport(dm, cone[cl], &support);
2317:       for (s = 0; s < supportSize; ++s) {
2318:         DMLabelSetValue(subpointMap, support[s], dim);
2319:       }
2320:     }
2321:   }
2322:   PetscFree2(pStart, pEnd);
2323:   return(0);
2324: }

2326: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2327: {
2328:   MPI_Comm       comm;
2329:   PetscBool      posOrient = PETSC_FALSE;
2330:   const PetscInt debug     = 0;
2331:   PetscInt       cellDim, faceSize, f;

2335:   PetscObjectGetComm((PetscObject)dm,&comm);
2336:   DMGetDimension(dm, &cellDim);
2337:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2339:   if (cellDim == 1 && numCorners == 2) {
2340:     /* Triangle */
2341:     faceSize  = numCorners-1;
2342:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2343:   } else if (cellDim == 2 && numCorners == 3) {
2344:     /* Triangle */
2345:     faceSize  = numCorners-1;
2346:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2347:   } else if (cellDim == 3 && numCorners == 4) {
2348:     /* Tetrahedron */
2349:     faceSize  = numCorners-1;
2350:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2351:   } else if (cellDim == 1 && numCorners == 3) {
2352:     /* Quadratic line */
2353:     faceSize  = 1;
2354:     posOrient = PETSC_TRUE;
2355:   } else if (cellDim == 2 && numCorners == 4) {
2356:     /* Quads */
2357:     faceSize = 2;
2358:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2359:       posOrient = PETSC_TRUE;
2360:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2361:       posOrient = PETSC_TRUE;
2362:     } else {
2363:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2364:         posOrient = PETSC_FALSE;
2365:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2366:     }
2367:   } else if (cellDim == 2 && numCorners == 6) {
2368:     /* Quadratic triangle (I hate this) */
2369:     /* Edges are determined by the first 2 vertices (corners of edges) */
2370:     const PetscInt faceSizeTri = 3;
2371:     PetscInt       sortedIndices[3], i, iFace;
2372:     PetscBool      found                    = PETSC_FALSE;
2373:     PetscInt       faceVerticesTriSorted[9] = {
2374:       0, 3,  4, /* bottom */
2375:       1, 4,  5, /* right */
2376:       2, 3,  5, /* left */
2377:     };
2378:     PetscInt       faceVerticesTri[9] = {
2379:       0, 3,  4, /* bottom */
2380:       1, 4,  5, /* right */
2381:       2, 5,  3, /* left */
2382:     };

2384:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2385:     PetscSortInt(faceSizeTri, sortedIndices);
2386:     for (iFace = 0; iFace < 3; ++iFace) {
2387:       const PetscInt ii = iFace*faceSizeTri;
2388:       PetscInt       fVertex, cVertex;

2390:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2391:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2392:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2393:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2394:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2395:               faceVertices[fVertex] = origVertices[cVertex];
2396:               break;
2397:             }
2398:           }
2399:         }
2400:         found = PETSC_TRUE;
2401:         break;
2402:       }
2403:     }
2404:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2405:     if (posOriented) *posOriented = PETSC_TRUE;
2406:     return(0);
2407:   } else if (cellDim == 2 && numCorners == 9) {
2408:     /* Quadratic quad (I hate this) */
2409:     /* Edges are determined by the first 2 vertices (corners of edges) */
2410:     const PetscInt faceSizeQuad = 3;
2411:     PetscInt       sortedIndices[3], i, iFace;
2412:     PetscBool      found                      = PETSC_FALSE;
2413:     PetscInt       faceVerticesQuadSorted[12] = {
2414:       0, 1,  4, /* bottom */
2415:       1, 2,  5, /* right */
2416:       2, 3,  6, /* top */
2417:       0, 3,  7, /* left */
2418:     };
2419:     PetscInt       faceVerticesQuad[12] = {
2420:       0, 1,  4, /* bottom */
2421:       1, 2,  5, /* right */
2422:       2, 3,  6, /* top */
2423:       3, 0,  7, /* left */
2424:     };

2426:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2427:     PetscSortInt(faceSizeQuad, sortedIndices);
2428:     for (iFace = 0; iFace < 4; ++iFace) {
2429:       const PetscInt ii = iFace*faceSizeQuad;
2430:       PetscInt       fVertex, cVertex;

2432:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2433:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2434:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2435:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2436:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2437:               faceVertices[fVertex] = origVertices[cVertex];
2438:               break;
2439:             }
2440:           }
2441:         }
2442:         found = PETSC_TRUE;
2443:         break;
2444:       }
2445:     }
2446:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2447:     if (posOriented) *posOriented = PETSC_TRUE;
2448:     return(0);
2449:   } else if (cellDim == 3 && numCorners == 8) {
2450:     /* Hexes
2451:        A hex is two oriented quads with the normal of the first
2452:        pointing up at the second.

2454:           7---6
2455:          /|  /|
2456:         4---5 |
2457:         | 1-|-2
2458:         |/  |/
2459:         0---3

2461:         Faces are determined by the first 4 vertices (corners of faces) */
2462:     const PetscInt faceSizeHex = 4;
2463:     PetscInt       sortedIndices[4], i, iFace;
2464:     PetscBool      found                     = PETSC_FALSE;
2465:     PetscInt       faceVerticesHexSorted[24] = {
2466:       0, 1, 2, 3,  /* bottom */
2467:       4, 5, 6, 7,  /* top */
2468:       0, 3, 4, 5,  /* front */
2469:       2, 3, 5, 6,  /* right */
2470:       1, 2, 6, 7,  /* back */
2471:       0, 1, 4, 7,  /* left */
2472:     };
2473:     PetscInt       faceVerticesHex[24] = {
2474:       1, 2, 3, 0,  /* bottom */
2475:       4, 5, 6, 7,  /* top */
2476:       0, 3, 5, 4,  /* front */
2477:       3, 2, 6, 5,  /* right */
2478:       2, 1, 7, 6,  /* back */
2479:       1, 0, 4, 7,  /* left */
2480:     };

2482:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2483:     PetscSortInt(faceSizeHex, sortedIndices);
2484:     for (iFace = 0; iFace < 6; ++iFace) {
2485:       const PetscInt ii = iFace*faceSizeHex;
2486:       PetscInt       fVertex, cVertex;

2488:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2489:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2490:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2491:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2492:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2493:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2494:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2495:               faceVertices[fVertex] = origVertices[cVertex];
2496:               break;
2497:             }
2498:           }
2499:         }
2500:         found = PETSC_TRUE;
2501:         break;
2502:       }
2503:     }
2504:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2505:     if (posOriented) *posOriented = PETSC_TRUE;
2506:     return(0);
2507:   } else if (cellDim == 3 && numCorners == 10) {
2508:     /* Quadratic tet */
2509:     /* Faces are determined by the first 3 vertices (corners of faces) */
2510:     const PetscInt faceSizeTet = 6;
2511:     PetscInt       sortedIndices[6], i, iFace;
2512:     PetscBool      found                     = PETSC_FALSE;
2513:     PetscInt       faceVerticesTetSorted[24] = {
2514:       0, 1, 2,  6, 7, 8, /* bottom */
2515:       0, 3, 4,  6, 7, 9,  /* front */
2516:       1, 4, 5,  7, 8, 9,  /* right */
2517:       2, 3, 5,  6, 8, 9,  /* left */
2518:     };
2519:     PetscInt       faceVerticesTet[24] = {
2520:       0, 1, 2,  6, 7, 8, /* bottom */
2521:       0, 4, 3,  6, 7, 9,  /* front */
2522:       1, 5, 4,  7, 8, 9,  /* right */
2523:       2, 3, 5,  8, 6, 9,  /* left */
2524:     };

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

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

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

2563:        Faces are determined by the first 4 vertices (corners of faces) */
2564:     const PetscInt faceSizeQuadHex = 9;
2565:     PetscInt       sortedIndices[9], i, iFace;
2566:     PetscBool      found                         = PETSC_FALSE;
2567:     PetscInt       faceVerticesQuadHexSorted[54] = {
2568:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2569:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2570:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2571:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2572:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2573:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2574:     };
2575:     PetscInt       faceVerticesQuadHex[54] = {
2576:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2577:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2578:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2579:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2580:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2581:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2582:     };

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

2590:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2591:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2592:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2593:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2594:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2595:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2596:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2597:               faceVertices[fVertex] = origVertices[cVertex];
2598:               break;
2599:             }
2600:           }
2601:         }
2602:         found = PETSC_TRUE;
2603:         break;
2604:       }
2605:     }
2606:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2607:     if (posOriented) *posOriented = PETSC_TRUE;
2608:     return(0);
2609:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2610:   if (!posOrient) {
2611:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2612:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2613:   } else {
2614:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2615:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2616:   }
2617:   if (posOriented) *posOriented = posOrient;
2618:   return(0);
2619: }

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

2625:   Not collective

2627:   Input Parameters:
2628: + dm           - The original mesh
2629: . cell         - The cell mesh point
2630: . faceSize     - The number of vertices on the face
2631: . face         - The face vertices
2632: . numCorners   - The number of vertices on the cell
2633: . indices      - Local numbering of face vertices in cell cone
2634: - origVertices - Original face vertices

2636:   Output Parameter:
2637: + faceVertices - The face vertices properly oriented
2638: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2640:   Level: developer

2642: .seealso: DMPlexGetCone()
2643: @*/
2644: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2645: {
2646:   const PetscInt *cone = NULL;
2647:   PetscInt        coneSize, v, f, v2;
2648:   PetscInt        oppositeVertex = -1;
2649:   PetscErrorCode  ierr;

2652:   DMPlexGetConeSize(dm, cell, &coneSize);
2653:   DMPlexGetCone(dm, cell, &cone);
2654:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2655:     PetscBool found = PETSC_FALSE;

2657:     for (f = 0; f < faceSize; ++f) {
2658:       if (face[f] == cone[v]) {
2659:         found = PETSC_TRUE; break;
2660:       }
2661:     }
2662:     if (found) {
2663:       indices[v2]      = v;
2664:       origVertices[v2] = cone[v];
2665:       ++v2;
2666:     } else {
2667:       oppositeVertex = v;
2668:     }
2669:   }
2670:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2671:   return(0);
2672: }

2674: /*
2675:   DMPlexInsertFace_Internal - Puts a face into the mesh

2677:   Not collective

2679:   Input Parameters:
2680:   + dm              - The DMPlex
2681:   . numFaceVertex   - The number of vertices in the face
2682:   . faceVertices    - The vertices in the face for dm
2683:   . subfaceVertices - The vertices in the face for subdm
2684:   . numCorners      - The number of vertices in the cell
2685:   . cell            - A cell in dm containing the face
2686:   . subcell         - A cell in subdm containing the face
2687:   . firstFace       - First face in the mesh
2688:   - newFacePoint    - Next face in the mesh

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

2693:   Level: developer
2694: */
2695: 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)
2696: {
2697:   MPI_Comm        comm;
2698:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2699:   const PetscInt *faces;
2700:   PetscInt        numFaces, coneSize;
2701:   PetscErrorCode  ierr;

2704:   PetscObjectGetComm((PetscObject)dm,&comm);
2705:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2706:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2707: #if 0
2708:   /* Cannot use this because support() has not been constructed yet */
2709:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2710: #else
2711:   {
2712:     PetscInt f;

2714:     numFaces = 0;
2715:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2716:     for (f = firstFace; f < *newFacePoint; ++f) {
2717:       PetscInt dof, off, d;

2719:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2720:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2721:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2722:       for (d = 0; d < dof; ++d) {
2723:         const PetscInt p = submesh->cones[off+d];
2724:         PetscInt       v;

2726:         for (v = 0; v < numFaceVertices; ++v) {
2727:           if (subfaceVertices[v] == p) break;
2728:         }
2729:         if (v == numFaceVertices) break;
2730:       }
2731:       if (d == dof) {
2732:         numFaces               = 1;
2733:         ((PetscInt*) faces)[0] = f;
2734:       }
2735:     }
2736:   }
2737: #endif
2738:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2739:   else if (numFaces == 1) {
2740:     /* Add the other cell neighbor for this face */
2741:     DMPlexSetCone(subdm, subcell, faces);
2742:   } else {
2743:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2744:     PetscBool posOriented;

2746:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2747:     origVertices        = &orientedVertices[numFaceVertices];
2748:     indices             = &orientedVertices[numFaceVertices*2];
2749:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2750:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2751:     /* TODO: I know that routine should return a permutation, not the indices */
2752:     for (v = 0; v < numFaceVertices; ++v) {
2753:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2754:       for (ov = 0; ov < numFaceVertices; ++ov) {
2755:         if (orientedVertices[ov] == vertex) {
2756:           orientedSubVertices[ov] = subvertex;
2757:           break;
2758:         }
2759:       }
2760:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2761:     }
2762:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2763:     DMPlexSetCone(subdm, subcell, newFacePoint);
2764:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2765:     ++(*newFacePoint);
2766:   }
2767: #if 0
2768:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2769: #else
2770:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2771: #endif
2772:   return(0);
2773: }

2775: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2776: {
2777:   MPI_Comm        comm;
2778:   DMLabel         subpointMap;
2779:   IS              subvertexIS,  subcellIS;
2780:   const PetscInt *subVertices, *subCells;
2781:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2782:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2783:   PetscInt        vStart, vEnd, c, f;
2784:   PetscErrorCode  ierr;

2787:   PetscObjectGetComm((PetscObject)dm,&comm);
2788:   /* Create subpointMap which marks the submesh */
2789:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2790:   DMPlexSetSubpointMap(subdm, subpointMap);
2791:   DMLabelDestroy(&subpointMap);
2792:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2793:   /* Setup chart */
2794:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2795:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2796:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2797:   DMPlexSetVTKCellHeight(subdm, 1);
2798:   /* Set cone sizes */
2799:   firstSubVertex = numSubCells;
2800:   firstSubFace   = numSubCells+numSubVertices;
2801:   newFacePoint   = firstSubFace;
2802:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2803:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2804:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2805:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2806:   for (c = 0; c < numSubCells; ++c) {
2807:     DMPlexSetConeSize(subdm, c, 1);
2808:   }
2809:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2810:     DMPlexSetConeSize(subdm, f, nFV);
2811:   }
2812:   DMSetUp(subdm);
2813:   /* Create face cones */
2814:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2815:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2816:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2817:   for (c = 0; c < numSubCells; ++c) {
2818:     const PetscInt cell    = subCells[c];
2819:     const PetscInt subcell = c;
2820:     PetscInt      *closure = NULL;
2821:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2823:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2824:     for (cl = 0; cl < closureSize*2; cl += 2) {
2825:       const PetscInt point = closure[cl];
2826:       PetscInt       subVertex;

2828:       if ((point >= vStart) && (point < vEnd)) {
2829:         ++numCorners;
2830:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2831:         if (subVertex >= 0) {
2832:           closure[faceSize] = point;
2833:           subface[faceSize] = firstSubVertex+subVertex;
2834:           ++faceSize;
2835:         }
2836:       }
2837:     }
2838:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2839:     if (faceSize == nFV) {
2840:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2841:     }
2842:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2843:   }
2844:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2845:   DMPlexSymmetrize(subdm);
2846:   DMPlexStratify(subdm);
2847:   /* Build coordinates */
2848:   {
2849:     PetscSection coordSection, subCoordSection;
2850:     Vec          coordinates, subCoordinates;
2851:     PetscScalar *coords, *subCoords;
2852:     PetscInt     numComp, coordSize, v;
2853:     const char  *name;

2855:     DMGetCoordinateSection(dm, &coordSection);
2856:     DMGetCoordinatesLocal(dm, &coordinates);
2857:     DMGetCoordinateSection(subdm, &subCoordSection);
2858:     PetscSectionSetNumFields(subCoordSection, 1);
2859:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2860:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2861:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2862:     for (v = 0; v < numSubVertices; ++v) {
2863:       const PetscInt vertex    = subVertices[v];
2864:       const PetscInt subvertex = firstSubVertex+v;
2865:       PetscInt       dof;

2867:       PetscSectionGetDof(coordSection, vertex, &dof);
2868:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2869:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2870:     }
2871:     PetscSectionSetUp(subCoordSection);
2872:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2873:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2874:     PetscObjectGetName((PetscObject)coordinates,&name);
2875:     PetscObjectSetName((PetscObject)subCoordinates,name);
2876:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2877:     VecSetType(subCoordinates,VECSTANDARD);
2878:     if (coordSize) {
2879:       VecGetArray(coordinates,    &coords);
2880:       VecGetArray(subCoordinates, &subCoords);
2881:       for (v = 0; v < numSubVertices; ++v) {
2882:         const PetscInt vertex    = subVertices[v];
2883:         const PetscInt subvertex = firstSubVertex+v;
2884:         PetscInt       dof, off, sdof, soff, d;

2886:         PetscSectionGetDof(coordSection, vertex, &dof);
2887:         PetscSectionGetOffset(coordSection, vertex, &off);
2888:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2889:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2890:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2891:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2892:       }
2893:       VecRestoreArray(coordinates,    &coords);
2894:       VecRestoreArray(subCoordinates, &subCoords);
2895:     }
2896:     DMSetCoordinatesLocal(subdm, subCoordinates);
2897:     VecDestroy(&subCoordinates);
2898:   }
2899:   /* Cleanup */
2900:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2901:   ISDestroy(&subvertexIS);
2902:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2903:   ISDestroy(&subcellIS);
2904:   return(0);
2905: }

2907: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2908: {
2909:   PetscInt       subPoint;

2912:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2913:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2914: }

2916: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2917: {
2918:   MPI_Comm         comm;
2919:   DMLabel          subpointMap;
2920:   IS              *subpointIS;
2921:   const PetscInt **subpoints;
2922:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2923:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2924:   PetscMPIInt      rank;
2925:   PetscErrorCode   ierr;

2928:   PetscObjectGetComm((PetscObject)dm,&comm);
2929:   MPI_Comm_rank(comm, &rank);
2930:   /* Create subpointMap which marks the submesh */
2931:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2932:   DMPlexSetSubpointMap(subdm, subpointMap);
2933:   if (cellHeight) {
2934:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2935:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
2936:   } else {
2937:     DMLabel         depth;
2938:     IS              pointIS;
2939:     const PetscInt *points;
2940:     PetscInt        numPoints;

2942:     DMPlexGetDepthLabel(dm, &depth);
2943:     DMLabelGetStratumSize(label, value, &numPoints);
2944:     DMLabelGetStratumIS(label, value, &pointIS);
2945:     ISGetIndices(pointIS, &points);
2946:     for (p = 0; p < numPoints; ++p) {
2947:       PetscInt *closure = NULL;
2948:       PetscInt  closureSize, c, pdim;

2950:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2951:       for (c = 0; c < closureSize*2; c += 2) {
2952:         DMLabelGetValue(depth, closure[c], &pdim);
2953:         DMLabelSetValue(subpointMap, closure[c], pdim);
2954:       }
2955:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2956:     }
2957:     ISRestoreIndices(pointIS, &points);
2958:     ISDestroy(&pointIS);
2959:   }
2960:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2961:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2962:   cMax = (cMax < 0) ? cEnd : cMax;
2963:   /* Setup chart */
2964:   DMGetDimension(dm, &dim);
2965:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2966:   for (d = 0; d <= dim; ++d) {
2967:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2968:     totSubPoints += numSubPoints[d];
2969:   }
2970:   DMPlexSetChart(subdm, 0, totSubPoints);
2971:   DMPlexSetVTKCellHeight(subdm, cellHeight);
2972:   /* Set cone sizes */
2973:   firstSubPoint[dim] = 0;
2974:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2975:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2976:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2977:   for (d = 0; d <= dim; ++d) {
2978:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2979:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2980:   }
2981:   for (d = 0; d <= dim; ++d) {
2982:     for (p = 0; p < numSubPoints[d]; ++p) {
2983:       const PetscInt  point    = subpoints[d][p];
2984:       const PetscInt  subpoint = firstSubPoint[d] + p;
2985:       const PetscInt *cone;
2986:       PetscInt        coneSize, coneSizeNew, c, val;

2988:       DMPlexGetConeSize(dm, point, &coneSize);
2989:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2990:       if (cellHeight && (d == dim)) {
2991:         DMPlexGetCone(dm, point, &cone);
2992:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2993:           DMLabelGetValue(subpointMap, cone[c], &val);
2994:           if (val >= 0) coneSizeNew++;
2995:         }
2996:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2997:       }
2998:     }
2999:   }
3000:   DMLabelDestroy(&subpointMap);
3001:   DMSetUp(subdm);
3002:   /* Set cones */
3003:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3004:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3005:   for (d = 0; d <= dim; ++d) {
3006:     for (p = 0; p < numSubPoints[d]; ++p) {
3007:       const PetscInt  point    = subpoints[d][p];
3008:       const PetscInt  subpoint = firstSubPoint[d] + p;
3009:       const PetscInt *cone, *ornt;
3010:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

3012:       if (d == dim-1) {
3013:         const PetscInt *support, *cone, *ornt;
3014:         PetscInt        supportSize, coneSize, s, subc;

3016:         DMPlexGetSupport(dm, point, &support);
3017:         DMPlexGetSupportSize(dm, point, &supportSize);
3018:         for (s = 0; s < supportSize; ++s) {
3019:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
3020:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3021:           if (subc >= 0) {
3022:             const PetscInt ccell = subpoints[d+1][subc];

3024:             DMPlexGetCone(dm, ccell, &cone);
3025:             DMPlexGetConeSize(dm, ccell, &coneSize);
3026:             DMPlexGetConeOrientation(dm, ccell, &ornt);
3027:             for (c = 0; c < coneSize; ++c) {
3028:               if (cone[c] == point) {
3029:                 fornt = ornt[c];
3030:                 break;
3031:               }
3032:             }
3033:             break;
3034:           }
3035:         }
3036:       }
3037:       DMPlexGetConeSize(dm, point, &coneSize);
3038:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3039:       DMPlexGetCone(dm, point, &cone);
3040:       DMPlexGetConeOrientation(dm, point, &ornt);
3041:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3042:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3043:         if (subc >= 0) {
3044:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3045:           orntNew[coneSizeNew] = ornt[c];
3046:           ++coneSizeNew;
3047:         }
3048:       }
3049:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3050:       if (fornt < 0) {
3051:         /* This should be replaced by a call to DMPlexReverseCell() */
3052: #if 0
3053:         DMPlexReverseCell(subdm, subpoint);
3054: #else
3055:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3056:           PetscInt faceSize, tmp;

3058:           tmp        = coneNew[c];
3059:           coneNew[c] = coneNew[coneSizeNew-1-c];
3060:           coneNew[coneSizeNew-1-c] = tmp;
3061:           DMPlexGetConeSize(dm, cone[c], &faceSize);
3062:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3063:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3064:           orntNew[coneSizeNew-1-c] = tmp;
3065:         }
3066:       }
3067:       DMPlexSetCone(subdm, subpoint, coneNew);
3068:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3069: #endif
3070:     }
3071:   }
3072:   PetscFree2(coneNew,orntNew);
3073:   DMPlexSymmetrize(subdm);
3074:   DMPlexStratify(subdm);
3075:   /* Build coordinates */
3076:   {
3077:     PetscSection coordSection, subCoordSection;
3078:     Vec          coordinates, subCoordinates;
3079:     PetscScalar *coords, *subCoords;
3080:     PetscInt     cdim, numComp, coordSize;
3081:     const char  *name;

3083:     DMGetCoordinateDim(dm, &cdim);
3084:     DMGetCoordinateSection(dm, &coordSection);
3085:     DMGetCoordinatesLocal(dm, &coordinates);
3086:     DMGetCoordinateSection(subdm, &subCoordSection);
3087:     PetscSectionSetNumFields(subCoordSection, 1);
3088:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3089:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3090:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3091:     for (v = 0; v < numSubPoints[0]; ++v) {
3092:       const PetscInt vertex    = subpoints[0][v];
3093:       const PetscInt subvertex = firstSubPoint[0]+v;
3094:       PetscInt       dof;

3096:       PetscSectionGetDof(coordSection, vertex, &dof);
3097:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3098:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3099:     }
3100:     PetscSectionSetUp(subCoordSection);
3101:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3102:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3103:     PetscObjectGetName((PetscObject)coordinates,&name);
3104:     PetscObjectSetName((PetscObject)subCoordinates,name);
3105:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3106:     VecSetBlockSize(subCoordinates, cdim);
3107:     VecSetType(subCoordinates,VECSTANDARD);
3108:     VecGetArray(coordinates,    &coords);
3109:     VecGetArray(subCoordinates, &subCoords);
3110:     for (v = 0; v < numSubPoints[0]; ++v) {
3111:       const PetscInt vertex    = subpoints[0][v];
3112:       const PetscInt subvertex = firstSubPoint[0]+v;
3113:       PetscInt dof, off, sdof, soff, d;

3115:       PetscSectionGetDof(coordSection, vertex, &dof);
3116:       PetscSectionGetOffset(coordSection, vertex, &off);
3117:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3118:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3119:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3120:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3121:     }
3122:     VecRestoreArray(coordinates,    &coords);
3123:     VecRestoreArray(subCoordinates, &subCoords);
3124:     DMSetCoordinatesLocal(subdm, subCoordinates);
3125:     VecDestroy(&subCoordinates);
3126:   }
3127:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3128:   {
3129:     PetscSF            sfPoint, sfPointSub;
3130:     IS                 subpIS;
3131:     const PetscSFNode *remotePoints;
3132:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3133:     const PetscInt    *localPoints, *subpoints;
3134:     PetscInt          *slocalPoints;
3135:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3136:     PetscMPIInt        rank;

3138:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3139:     DMGetPointSF(dm, &sfPoint);
3140:     DMGetPointSF(subdm, &sfPointSub);
3141:     DMPlexGetChart(dm, &pStart, &pEnd);
3142:     DMPlexGetChart(subdm, NULL, &numSubroots);
3143:     DMPlexCreateSubpointIS(subdm, &subpIS);
3144:     if (subpIS) {
3145:       ISGetIndices(subpIS, &subpoints);
3146:       ISGetLocalSize(subpIS, &numSubpoints);
3147:     }
3148:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3149:     if (numRoots >= 0) {
3150:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3151:       for (p = 0; p < pEnd-pStart; ++p) {
3152:         newLocalPoints[p].rank  = -2;
3153:         newLocalPoints[p].index = -2;
3154:       }
3155:       /* Set subleaves */
3156:       for (l = 0; l < numLeaves; ++l) {
3157:         const PetscInt point    = localPoints[l];
3158:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3160:         if (subpoint < 0) continue;
3161:         newLocalPoints[point-pStart].rank  = rank;
3162:         newLocalPoints[point-pStart].index = subpoint;
3163:         ++numSubleaves;
3164:       }
3165:       /* Must put in owned subpoints */
3166:       for (p = pStart; p < pEnd; ++p) {
3167:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3169:         if (subpoint < 0) {
3170:           newOwners[p-pStart].rank  = -3;
3171:           newOwners[p-pStart].index = -3;
3172:         } else {
3173:           newOwners[p-pStart].rank  = rank;
3174:           newOwners[p-pStart].index = subpoint;
3175:         }
3176:       }
3177:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3178:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3179:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3180:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3181:       PetscMalloc1(numSubleaves, &slocalPoints);
3182:       PetscMalloc1(numSubleaves, &sremotePoints);
3183:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3184:         const PetscInt point    = localPoints[l];
3185:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3187:         if (subpoint < 0) continue;
3188:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3189:         slocalPoints[sl]        = subpoint;
3190:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3191:         sremotePoints[sl].index = newLocalPoints[point].index;
3192:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3193:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3194:         ++sl;
3195:       }
3196:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3197:       PetscFree2(newLocalPoints,newOwners);
3198:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3199:     }
3200:     if (subpIS) {
3201:       ISRestoreIndices(subpIS, &subpoints);
3202:       ISDestroy(&subpIS);
3203:     }
3204:   }
3205:   /* Cleanup */
3206:   for (d = 0; d <= dim; ++d) {
3207:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3208:     ISDestroy(&subpointIS[d]);
3209:   }
3210:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3211:   return(0);
3212: }

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

3219:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3220:   return(0);
3221: }

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

3226:   Input Parameters:
3227: + dm           - The original mesh
3228: . vertexLabel  - The DMLabel marking points contained in the surface
3229: . value        - The label value to use
3230: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

3232:   Output Parameter:
3233: . subdm - The surface mesh

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

3237:   Level: developer

3239: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3240: @*/
3241: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3242: {
3243:   PetscInt       dim, cdim, depth;

3249:   DMGetDimension(dm, &dim);
3250:   DMPlexGetDepth(dm, &depth);
3251:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3252:   DMSetType(*subdm, DMPLEX);
3253:   DMSetDimension(*subdm, dim-1);
3254:   DMGetCoordinateDim(dm, &cdim);
3255:   DMSetCoordinateDim(*subdm, cdim);
3256:   if (depth == dim) {
3257:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3258:   } else {
3259:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3260:   }
3261:   return(0);
3262: }

3264: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3265: {
3266:   MPI_Comm        comm;
3267:   DMLabel         subpointMap;
3268:   IS              subvertexIS;
3269:   const PetscInt *subVertices;
3270:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3271:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3272:   PetscInt        cMax, c, f;
3273:   PetscErrorCode  ierr;

3276:   PetscObjectGetComm((PetscObject)dm, &comm);
3277:   /* Create subpointMap which marks the submesh */
3278:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3279:   DMPlexSetSubpointMap(subdm, subpointMap);
3280:   DMLabelDestroy(&subpointMap);
3281:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3282:   /* Setup chart */
3283:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3284:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3285:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3286:   DMPlexSetVTKCellHeight(subdm, 1);
3287:   /* Set cone sizes */
3288:   firstSubVertex = numSubCells;
3289:   firstSubFace   = numSubCells+numSubVertices;
3290:   newFacePoint   = firstSubFace;
3291:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3292:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3293:   for (c = 0; c < numSubCells; ++c) {
3294:     DMPlexSetConeSize(subdm, c, 1);
3295:   }
3296:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3297:     DMPlexSetConeSize(subdm, f, nFV);
3298:   }
3299:   DMSetUp(subdm);
3300:   /* Create face cones */
3301:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3302:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3303:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3304:   for (c = 0; c < numSubCells; ++c) {
3305:     const PetscInt  cell    = subCells[c];
3306:     const PetscInt  subcell = c;
3307:     const PetscInt *cone, *cells;
3308:     PetscInt        numCells, subVertex, p, v;

3310:     if (cell < cMax) continue;
3311:     DMPlexGetCone(dm, cell, &cone);
3312:     for (v = 0; v < nFV; ++v) {
3313:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3314:       subface[v] = firstSubVertex+subVertex;
3315:     }
3316:     DMPlexSetCone(subdm, newFacePoint, subface);
3317:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3318:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3319:     /* Not true in parallel
3320:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3321:     for (p = 0; p < numCells; ++p) {
3322:       PetscInt negsubcell;

3324:       if (cells[p] >= cMax) continue;
3325:       /* I know this is a crap search */
3326:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3327:         if (subCells[negsubcell] == cells[p]) break;
3328:       }
3329:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3330:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3331:     }
3332:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3333:     ++newFacePoint;
3334:   }
3335:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3336:   DMPlexSymmetrize(subdm);
3337:   DMPlexStratify(subdm);
3338:   /* Build coordinates */
3339:   {
3340:     PetscSection coordSection, subCoordSection;
3341:     Vec          coordinates, subCoordinates;
3342:     PetscScalar *coords, *subCoords;
3343:     PetscInt     cdim, numComp, coordSize, v;
3344:     const char  *name;

3346:     DMGetCoordinateDim(dm, &cdim);
3347:     DMGetCoordinateSection(dm, &coordSection);
3348:     DMGetCoordinatesLocal(dm, &coordinates);
3349:     DMGetCoordinateSection(subdm, &subCoordSection);
3350:     PetscSectionSetNumFields(subCoordSection, 1);
3351:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3352:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3353:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3354:     for (v = 0; v < numSubVertices; ++v) {
3355:       const PetscInt vertex    = subVertices[v];
3356:       const PetscInt subvertex = firstSubVertex+v;
3357:       PetscInt       dof;

3359:       PetscSectionGetDof(coordSection, vertex, &dof);
3360:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3361:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3362:     }
3363:     PetscSectionSetUp(subCoordSection);
3364:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3365:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3366:     PetscObjectGetName((PetscObject)coordinates,&name);
3367:     PetscObjectSetName((PetscObject)subCoordinates,name);
3368:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3369:     VecSetBlockSize(subCoordinates, cdim);
3370:     VecSetType(subCoordinates,VECSTANDARD);
3371:     VecGetArray(coordinates,    &coords);
3372:     VecGetArray(subCoordinates, &subCoords);
3373:     for (v = 0; v < numSubVertices; ++v) {
3374:       const PetscInt vertex    = subVertices[v];
3375:       const PetscInt subvertex = firstSubVertex+v;
3376:       PetscInt       dof, off, sdof, soff, d;

3378:       PetscSectionGetDof(coordSection, vertex, &dof);
3379:       PetscSectionGetOffset(coordSection, vertex, &off);
3380:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3381:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3382:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3383:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3384:     }
3385:     VecRestoreArray(coordinates,    &coords);
3386:     VecRestoreArray(subCoordinates, &subCoords);
3387:     DMSetCoordinatesLocal(subdm, subCoordinates);
3388:     VecDestroy(&subCoordinates);
3389:   }
3390:   /* Build SF */
3391:   CHKMEMQ;
3392:   {
3393:     PetscSF            sfPoint, sfPointSub;
3394:     const PetscSFNode *remotePoints;
3395:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3396:     const PetscInt    *localPoints;
3397:     PetscInt          *slocalPoints;
3398:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3399:     PetscMPIInt        rank;

3401:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3402:     DMGetPointSF(dm, &sfPoint);
3403:     DMGetPointSF(subdm, &sfPointSub);
3404:     DMPlexGetChart(dm, &pStart, &pEnd);
3405:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3406:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3407:     if (numRoots >= 0) {
3408:       /* Only vertices should be shared */
3409:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3410:       for (p = 0; p < pEnd-pStart; ++p) {
3411:         newLocalPoints[p].rank  = -2;
3412:         newLocalPoints[p].index = -2;
3413:       }
3414:       /* Set subleaves */
3415:       for (l = 0; l < numLeaves; ++l) {
3416:         const PetscInt point    = localPoints[l];
3417:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3419:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3420:         if (subPoint < 0) continue;
3421:         newLocalPoints[point-pStart].rank  = rank;
3422:         newLocalPoints[point-pStart].index = subPoint;
3423:         ++numSubLeaves;
3424:       }
3425:       /* Must put in owned subpoints */
3426:       for (p = pStart; p < pEnd; ++p) {
3427:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3429:         if (subPoint < 0) {
3430:           newOwners[p-pStart].rank  = -3;
3431:           newOwners[p-pStart].index = -3;
3432:         } else {
3433:           newOwners[p-pStart].rank  = rank;
3434:           newOwners[p-pStart].index = subPoint;
3435:         }
3436:       }
3437:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3438:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3439:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3440:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3441:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3442:       PetscMalloc1(numSubLeaves, &sremotePoints);
3443:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3444:         const PetscInt point    = localPoints[l];
3445:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3447:         if (subPoint < 0) continue;
3448:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3449:         slocalPoints[sl]        = subPoint;
3450:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3451:         sremotePoints[sl].index = newLocalPoints[point].index;
3452:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3453:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3454:         ++sl;
3455:       }
3456:       PetscFree2(newLocalPoints,newOwners);
3457:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3458:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3459:     }
3460:   }
3461:   CHKMEMQ;
3462:   /* Cleanup */
3463:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3464:   ISDestroy(&subvertexIS);
3465:   PetscFree(subCells);
3466:   return(0);
3467: }

3469: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3470: {
3471:   DMLabel        label = NULL;

3475:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3476:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3477:   return(0);
3478: }

3480: /*@C
3481:   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.

3483:   Input Parameters:
3484: + dm          - The original mesh
3485: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3486: . label       - A label name, or NULL
3487: - value  - A label value

3489:   Output Parameter:
3490: . subdm - The surface mesh

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

3494:   Level: developer

3496: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3497: @*/
3498: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3499: {
3500:   PetscInt       dim, cdim, depth;

3506:   DMGetDimension(dm, &dim);
3507:   DMPlexGetDepth(dm, &depth);
3508:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3509:   DMSetType(*subdm, DMPLEX);
3510:   DMSetDimension(*subdm, dim-1);
3511:   DMGetCoordinateDim(dm, &cdim);
3512:   DMSetCoordinateDim(*subdm, cdim);
3513:   if (depth == dim) {
3514:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3515:   } else {
3516:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3517:   }
3518:   return(0);
3519: }

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

3524:   Input Parameters:
3525: + dm        - The original mesh
3526: . cellLabel - The DMLabel marking cells contained in the new mesh
3527: - value     - The label value to use

3529:   Output Parameter:
3530: . subdm - The new mesh

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

3534:   Level: developer

3536: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3537: @*/
3538: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3539: {
3540:   PetscInt       dim;

3546:   DMGetDimension(dm, &dim);
3547:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3548:   DMSetType(*subdm, DMPLEX);
3549:   DMSetDimension(*subdm, dim);
3550:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3551:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3552:   return(0);
3553: }

3555: /*@
3556:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3558:   Input Parameter:
3559: . dm - The submesh DM

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

3564:   Level: developer

3566: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3567: @*/
3568: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3569: {
3573:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3574:   return(0);
3575: }

3577: /*@
3578:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3580:   Input Parameters:
3581: + dm - The submesh DM
3582: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3586:   Level: developer

3588: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3589: @*/
3590: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3591: {
3592:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3593:   DMLabel        tmp;

3598:   tmp  = mesh->subpointMap;
3599:   mesh->subpointMap = subpointMap;
3600:   PetscObjectReference((PetscObject) mesh->subpointMap);
3601:   DMLabelDestroy(&tmp);
3602:   return(0);
3603: }

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

3608:   Input Parameter:
3609: . dm - The submesh DM

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

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

3616:   Level: developer

3618: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3619: @*/
3620: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3621: {
3622:   MPI_Comm        comm;
3623:   DMLabel         subpointMap;
3624:   IS              is;
3625:   const PetscInt *opoints;
3626:   PetscInt       *points, *depths;
3627:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3628:   PetscErrorCode  ierr;

3633:   PetscObjectGetComm((PetscObject)dm,&comm);
3634:   *subpointIS = NULL;
3635:   DMPlexGetSubpointMap(dm, &subpointMap);
3636:   DMPlexGetDepth(dm, &depth);
3637:   if (subpointMap && depth >= 0) {
3638:     DMPlexGetChart(dm, &pStart, &pEnd);
3639:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3640:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3641:     depths[0] = depth;
3642:     depths[1] = 0;
3643:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3644:     PetscMalloc1(pEnd, &points);
3645:     for(d = 0, off = 0; d <= depth; ++d) {
3646:       const PetscInt dep = depths[d];

3648:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3649:       DMLabelGetStratumSize(subpointMap, dep, &n);
3650:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3651:         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);
3652:       } else {
3653:         if (!n) {
3654:           if (d == 0) {
3655:             /* Missing cells */
3656:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3657:           } else {
3658:             /* Missing faces */
3659:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3660:           }
3661:         }
3662:       }
3663:       if (n) {
3664:         DMLabelGetStratumIS(subpointMap, dep, &is);
3665:         ISGetIndices(is, &opoints);
3666:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3667:         ISRestoreIndices(is, &opoints);
3668:         ISDestroy(&is);
3669:       }
3670:     }
3671:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3672:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3673:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3674:   }
3675:   return(0);
3676: }

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

3682:   Note collective

3684:   Input Parameters:
3685: + dm - The submesh DM
3686: - p  - The point in the original, from which the submesh was created

3688:   Output Parameter:
3689: . subp - The point in the submesh

3691:   Level: developer

3693: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3694: @*/
3695: PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp)
3696: {
3697:   DMLabel        spmap;

3701:   *subp = p;
3702:   DMPlexGetSubpointMap(dm, &spmap);
3703:   if (spmap) {
3704:     IS              subpointIS;
3705:     const PetscInt *subpoints;
3706:     PetscInt        numSubpoints;

3708:     /* TODO Cache the IS, making it look like an index */
3709:     DMPlexCreateSubpointIS(dm, &subpointIS);
3710:     ISGetLocalSize(subpointIS, &numSubpoints);
3711:     ISGetIndices(subpointIS, &subpoints);
3712:     PetscFindInt(p, numSubpoints, subpoints, subp);
3713:     if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p);
3714:     ISRestoreIndices(subpointIS, &subpoints);
3715:     ISDestroy(&subpointIS);
3716:   }
3717:   return(0);
3718: }

3720: /*@
3721:   DMPlexGetAuxiliaryPoint - For a given point in the DM, return the matching point in the auxiliary DM.

3723:   Note collective

3725:   Input Parameters:
3726: + dm    - The DM
3727: . dmAux - The related auxiliary DM
3728: - p     - The point in the original DM

3730:   Output Parameter:
3731: . subp - The point in the auxiliary DM

3733:   Notes: If the DM is a submesh, we assume the dmAux is as well and just return the point. If only dmAux is a submesh,
3734:   then we map the point back to the original space.

3736:   Level: developer

3738: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3739: @*/
3740: PetscErrorCode DMPlexGetAuxiliaryPoint(DM dm, DM dmAux, PetscInt p, PetscInt *subp)
3741: {
3742:   DMLabel        spmap;

3746:   *subp = p;
3747:   /* If dm is a submesh, do not get subpoint */
3748:   DMPlexGetSubpointMap(dm, &spmap);
3749:   if (dmAux && !spmap) {
3750:     PetscInt h;

3752:     DMPlexGetVTKCellHeight(dmAux, &h);
3753:     DMPlexGetSubpointMap(dmAux, &spmap);
3754:     if (spmap && !h) {DMLabelGetValue(spmap, p, subp);}
3755:     else             {DMPlexGetSubpoint(dmAux, p, subp);}
3756:   }
3757:   return(0);
3758: }