Actual source code: plexsubmesh.c

petsc-3.9.2 2018-05-20
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;

352:   DMPlexGetDepth(dm, &depth);
353:   if (depth < 0) return(0);
354:   /* Step 1: Expand chart */
355:   DMPlexGetChart(dm, &pStart, &pEnd);
356:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
357:   DMPlexSetChart(dmNew, pStart, pEnd);
358:   /* Step 2: Set cone and support sizes */
359:   for (d = 0; d <= depth; ++d) {
360:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
361:     for (p = pStart; p < pEnd; ++p) {
362:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
363:       PetscInt size;

365:       DMPlexGetConeSize(dm, p, &size);
366:       DMPlexSetConeSize(dmNew, newp, size);
367:       DMPlexGetSupportSize(dm, p, &size);
368:       DMPlexSetSupportSize(dmNew, newp, size);
369:     }
370:   }
371:   return(0);
372: }

374: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
375: {
376:   PetscInt      *newpoints;
377:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

381:   DMPlexGetDepth(dm, &depth);
382:   if (depth < 0) return(0);
383:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
384:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
385:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
386:   /* Step 5: Set cones and supports */
387:   DMPlexGetChart(dm, &pStart, &pEnd);
388:   for (p = pStart; p < pEnd; ++p) {
389:     const PetscInt *points = NULL, *orientations = NULL;
390:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

392:     DMPlexGetConeSize(dm, p, &size);
393:     DMPlexGetCone(dm, p, &points);
394:     DMPlexGetConeOrientation(dm, p, &orientations);
395:     for (i = 0; i < size; ++i) {
396:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
397:     }
398:     DMPlexSetCone(dmNew, newp, newpoints);
399:     DMPlexSetConeOrientation(dmNew, newp, orientations);
400:     DMPlexGetSupportSize(dm, p, &size);
401:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
402:     DMPlexGetSupport(dm, p, &points);
403:     for (i = 0; i < size; ++i) {
404:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
405:     }
406:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
407:     DMPlexSetSupport(dmNew, newp, newpoints);
408:   }
409:   PetscFree(newpoints);
410:   return(0);
411: }

413: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
414: {
415:   PetscSection   coordSection, newCoordSection;
416:   Vec            coordinates, newCoordinates;
417:   PetscScalar   *coords, *newCoords;
418:   PetscInt       coordSize, sStart, sEnd;
419:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
420:   PetscBool      hasCells;

424:   DMGetCoordinateDim(dm, &dim);
425:   DMSetCoordinateDim(dmNew, dim);
426:   DMPlexGetDepth(dm, &depth);
427:   /* Step 8: Convert coordinates */
428:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
429:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
430:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
431:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
432:   DMGetCoordinateSection(dm, &coordSection);
433:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
434:   PetscSectionSetNumFields(newCoordSection, 1);
435:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
436:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
437:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
438:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
439:   if (hasCells) {
440:     for (c = cStart; c < cEnd; ++c) {
441:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

443:       PetscSectionGetDof(coordSection, c, &dof);
444:       PetscSectionSetDof(newCoordSection, cNew, dof);
445:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
446:     }
447:   }
448:   for (v = vStartNew; v < vEndNew; ++v) {
449:     PetscSectionSetDof(newCoordSection, v, dim);
450:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
451:   }
452:   PetscSectionSetUp(newCoordSection);
453:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
454:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
455:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
456:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
457:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
458:   VecSetBlockSize(newCoordinates, dim);
459:   VecSetType(newCoordinates,VECSTANDARD);
460:   DMSetCoordinatesLocal(dmNew, newCoordinates);
461:   DMGetCoordinatesLocal(dm, &coordinates);
462:   VecGetArray(coordinates, &coords);
463:   VecGetArray(newCoordinates, &newCoords);
464:   if (hasCells) {
465:     for (c = cStart; c < cEnd; ++c) {
466:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

468:       PetscSectionGetDof(coordSection, c, &dof);
469:       PetscSectionGetOffset(coordSection, c, &off);
470:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
471:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
472:     }
473:   }
474:   for (v = vStart; v < vEnd; ++v) {
475:     PetscInt dof, off, noff, d;

477:     PetscSectionGetDof(coordSection, v, &dof);
478:     PetscSectionGetOffset(coordSection, v, &off);
479:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
480:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
481:   }
482:   VecRestoreArray(coordinates, &coords);
483:   VecRestoreArray(newCoordinates, &newCoords);
484:   VecDestroy(&newCoordinates);
485:   PetscSectionDestroy(&newCoordSection);
486:   return(0);
487: }

489: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
490: {
491:   PetscInt           depth = 0;
492:   PetscSF            sfPoint, sfPointNew;
493:   const PetscSFNode *remotePoints;
494:   PetscSFNode       *gremotePoints;
495:   const PetscInt    *localPoints;
496:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
497:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
498:   PetscErrorCode     ierr;

501:   DMPlexGetDepth(dm, &depth);
502:   /* Step 9: Convert pointSF */
503:   DMGetPointSF(dm, &sfPoint);
504:   DMGetPointSF(dmNew, &sfPointNew);
505:   DMPlexGetChart(dm, &pStart, &pEnd);
506:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
507:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
508:   if (numRoots >= 0) {
509:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
510:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
511:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
512:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
513:     PetscMalloc1(numLeaves,    &glocalPoints);
514:     PetscMalloc1(numLeaves, &gremotePoints);
515:     for (l = 0; l < numLeaves; ++l) {
516:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
517:       gremotePoints[l].rank  = remotePoints[l].rank;
518:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
519:     }
520:     PetscFree2(newLocation,newRemoteLocation);
521:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
522:   }
523:   return(0);
524: }

526: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
527: {
528:   PetscSF            sfPoint;
529:   DMLabel            vtkLabel, ghostLabel;
530:   const PetscSFNode *leafRemote;
531:   const PetscInt    *leafLocal;
532:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
533:   PetscMPIInt        rank;
534:   PetscErrorCode     ierr;

537:   DMPlexGetDepth(dm, &depth);
538:   /* Step 10: Convert labels */
539:   DMGetNumLabels(dm, &numLabels);
540:   for (l = 0; l < numLabels; ++l) {
541:     DMLabel         label, newlabel;
542:     const char     *lname;
543:     PetscBool       isDepth;
544:     IS              valueIS;
545:     const PetscInt *values;
546:     PetscInt        numValues, val;

548:     DMGetLabelName(dm, l, &lname);
549:     PetscStrcmp(lname, "depth", &isDepth);
550:     if (isDepth) continue;
551:     DMCreateLabel(dmNew, lname);
552:     DMGetLabel(dm, lname, &label);
553:     DMGetLabel(dmNew, lname, &newlabel);
554:     DMLabelGetDefaultValue(label,&val);
555:     DMLabelSetDefaultValue(newlabel,val);
556:     DMLabelGetValueIS(label, &valueIS);
557:     ISGetLocalSize(valueIS, &numValues);
558:     ISGetIndices(valueIS, &values);
559:     for (val = 0; val < numValues; ++val) {
560:       IS              pointIS;
561:       const PetscInt *points;
562:       PetscInt        numPoints, p;

564:       DMLabelGetStratumIS(label, values[val], &pointIS);
565:       ISGetLocalSize(pointIS, &numPoints);
566:       ISGetIndices(pointIS, &points);
567:       for (p = 0; p < numPoints; ++p) {
568:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

570:         DMLabelSetValue(newlabel, newpoint, values[val]);
571:       }
572:       ISRestoreIndices(pointIS, &points);
573:       ISDestroy(&pointIS);
574:     }
575:     ISRestoreIndices(valueIS, &values);
576:     ISDestroy(&valueIS);
577:   }
578:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
579:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
580:   DMGetPointSF(dm, &sfPoint);
581:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
582:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
583:   DMCreateLabel(dmNew, "vtk");
584:   DMCreateLabel(dmNew, "ghost");
585:   DMGetLabel(dmNew, "vtk", &vtkLabel);
586:   DMGetLabel(dmNew, "ghost", &ghostLabel);
587:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
588:     for (; c < leafLocal[l] && c < cEnd; ++c) {
589:       DMLabelSetValue(vtkLabel, c, 1);
590:     }
591:     if (leafLocal[l] >= cEnd) break;
592:     if (leafRemote[l].rank == rank) {
593:       DMLabelSetValue(vtkLabel, c, 1);
594:     } else {
595:       DMLabelSetValue(ghostLabel, c, 2);
596:     }
597:   }
598:   for (; c < cEnd; ++c) {
599:     DMLabelSetValue(vtkLabel, c, 1);
600:   }
601:   if (0) {
602:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
603:   }
604:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
605:   for (f = fStart; f < fEnd; ++f) {
606:     PetscInt numCells;

608:     DMPlexGetSupportSize(dmNew, f, &numCells);
609:     if (numCells < 2) {
610:       DMLabelSetValue(ghostLabel, f, 1);
611:     } else {
612:       const PetscInt *cells = NULL;
613:       PetscInt        vA, vB;

615:       DMPlexGetSupport(dmNew, f, &cells);
616:       DMLabelGetValue(vtkLabel, cells[0], &vA);
617:       DMLabelGetValue(vtkLabel, cells[1], &vB);
618:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
619:     }
620:   }
621:   if (0) {
622:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
623:   }
624:   return(0);
625: }

627: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
628: {
629:   DM             refTree;
630:   PetscSection   pSec;
631:   PetscInt       *parents, *childIDs;

635:   DMPlexGetReferenceTree(dm,&refTree);
636:   DMPlexSetReferenceTree(dmNew,refTree);
637:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
638:   if (pSec) {
639:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
640:     PetscInt *childIDsShifted;
641:     PetscSection pSecShifted;

643:     PetscSectionGetChart(pSec,&pStart,&pEnd);
644:     DMPlexGetDepth(dm,&depth);
645:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
646:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
647:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
648:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
649:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
650:     for (p = pStartShifted; p < pEndShifted; p++) {
651:       /* start off assuming no children */
652:       PetscSectionSetDof(pSecShifted,p,0);
653:     }
654:     for (p = pStart; p < pEnd; p++) {
655:       PetscInt dof;
656:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

658:       PetscSectionGetDof(pSec,p,&dof);
659:       PetscSectionSetDof(pSecShifted,pNew,dof);
660:     }
661:     PetscSectionSetUp(pSecShifted);
662:     for (p = pStart; p < pEnd; p++) {
663:       PetscInt dof;
664:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

666:       PetscSectionGetDof(pSec,p,&dof);
667:       if (dof) {
668:         PetscInt off, offNew;

670:         PetscSectionGetOffset(pSec,p,&off);
671:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
672:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
673:         childIDsShifted[offNew] = childIDs[off];
674:       }
675:     }
676:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
677:     PetscFree2(parentsShifted,childIDsShifted);
678:     PetscSectionDestroy(&pSecShifted);
679:   }
680:   return(0);
681: }

683: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
684: {
685:   PetscSF               sf;
686:   IS                    valueIS;
687:   const PetscInt       *values, *leaves;
688:   PetscInt             *depthShift;
689:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
690:   PetscBool             isper;
691:   const PetscReal      *maxCell, *L;
692:   const DMBoundaryType *bd;
693:   PetscErrorCode        ierr;

696:   DMGetPointSF(dm, &sf);
697:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
698:   nleaves = PetscMax(0, nleaves);
699:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
700:   /* Count ghost cells */
701:   DMLabelGetValueIS(label, &valueIS);
702:   ISGetLocalSize(valueIS, &numFS);
703:   ISGetIndices(valueIS, &values);
704:   Ng   = 0;
705:   for (fs = 0; fs < numFS; ++fs) {
706:     IS              faceIS;
707:     const PetscInt *faces;
708:     PetscInt        numFaces, f, numBdFaces = 0;

710:     DMLabelGetStratumIS(label, values[fs], &faceIS);
711:     ISGetLocalSize(faceIS, &numFaces);
712:     ISGetIndices(faceIS, &faces);
713:     for (f = 0; f < numFaces; ++f) {
714:       PetscInt numChildren;

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

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

748:     DMLabelGetStratumIS(label, values[fs], &faceIS);
749:     ISGetLocalSize(faceIS, &numFaces);
750:     ISGetIndices(faceIS, &faces);
751:     for (f = 0; f < numFaces; ++f) {
752:       PetscInt size, numChildren;

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

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

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

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

811:   Collective on dm

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

817:   Output Parameters:
818: + numGhostCells - The number of ghost cells added to the DM
819: - dmGhosted - The new DM

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

823:   Level: developer

825: .seealso: DMCreate()
826: @*/
827: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
828: {
829:   DM             gdm;
830:   DMLabel        label;
831:   const char    *name = labelName ? labelName : "Face Sets";
832:   PetscInt       dim;
833:   PetscBool      flag;

840:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
841:   DMSetType(gdm, DMPLEX);
842:   DMGetDimension(dm, &dim);
843:   DMSetDimension(gdm, dim);
844:   DMPlexGetAdjacencyUseCone(dm, &flag);
845:   DMPlexSetAdjacencyUseCone(gdm, flag);
846:   DMPlexGetAdjacencyUseClosure(dm, &flag);
847:   DMPlexSetAdjacencyUseClosure(gdm, flag);
848:   DMGetLabel(dm, name, &label);
849:   if (!label) {
850:     /* Get label for boundary faces */
851:     DMCreateLabel(dm, name);
852:     DMGetLabel(dm, name, &label);
853:     DMPlexMarkBoundaryFaces(dm, 1, label);
854:   }
855:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
856:   DMCopyBoundary(dm, gdm);
857:   *dmGhosted = gdm;
858:   return(0);
859: }

861: /*
862:   We are adding three kinds of points here:
863:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
864:     Non-replicated: Points which exist on the fault, but are not replicated
865:     Hybrid:         Entirely new points, such as cohesive cells

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

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

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

960:       DMPlexGetConeSize(dm, oldp, &coneSize);
961:       DMPlexSetConeSize(sdm, splitp, coneSize);
962:       DMPlexGetSupportSize(dm, oldp, &supportSize);
963:       DMPlexSetSupportSize(sdm, splitp, supportSize);
964:       if (dep == depth-1) {
965:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

972:         DMPlexGetSupport(dm, oldp, &support);
973:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
974:           PetscInt val;

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

991:         DMPlexGetSupport(dm, oldp, &support);
992:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
993:           PetscInt val;

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

1017:       DMPlexGetConeSize(dm, oldp, &coneSize);
1018:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1019:       DMPlexGetSupport(dm, oldp, &support);
1020:       if (dep == 0) {
1021:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

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

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

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

1085:         /* Split face:       copy in old face to new face to start */
1086:         DMPlexGetSupport(sdm, newp,  &supportF);
1087:         DMPlexSetSupport(sdm, splitp, supportF);
1088:         /* Split old face:   old vertices/edges in cone so no change */
1089:         /* Split new face:   new vertices/edges in cone */
1090:         for (q = 0; q < coneSize; ++q) {
1091:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1092:           if (v < 0) {
1093:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1094:             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);
1095:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1096:             hasUnsplit   = PETSC_TRUE;
1097:           } else {
1098:             coneNew[2+q] = v + pMaxNew[dep-1];
1099:             if (dep > 1) {
1100:               const PetscInt *econe;
1101:               PetscInt        econeSize, r, vs, vu;

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

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

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

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

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

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

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

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

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

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

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

1289:       DMPlexGetConeSize(dm, oldp, &coneSize);
1290:       DMPlexGetCone(dm, oldp, &cone);
1291:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1292:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1293:       DMPlexGetSupport(dm, oldp, &support);
1294:       if (dep == 0) {
1295:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

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

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

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

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

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

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

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

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

1448:       for (l = 0; l < numLabels; ++l) {
1449:         DMLabel     mlabel;
1450:         const char *lname;
1451:         PetscInt    val;
1452:         PetscBool   isDepth;

1454:         DMGetLabelName(sdm, l, &lname);
1455:         PetscStrcmp(lname, "depth", &isDepth);
1456:         if (isDepth) continue;
1457:         DMGetLabel(sdm, lname, &mlabel);
1458:         DMLabelGetValue(mlabel, newp, &val);
1459:         if (val >= 0) {
1460:           DMLabelSetValue(mlabel, splitp, val);
1461:         }
1462:       }
1463:     }
1464:   }
1465:   for (sp = 0; sp < numSP; ++sp) {
1466:     const PetscInt dep = values[sp];

1468:     if ((dep < 0) || (dep > depth)) continue;
1469:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1470:     ISDestroy(&splitIS[dep]);
1471:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1472:     ISDestroy(&unsplitIS[dep]);
1473:   }
1474:   if (label) {
1475:     ISRestoreIndices(valueIS, &values);
1476:     ISDestroy(&valueIS);
1477:   }
1478:   for (d = 0; d <= depth; ++d) {
1479:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1480:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1481:   }
1482:   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);
1483:   PetscFree3(coneNew, coneONew, supportNew);
1484:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1485:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1486:   return(0);
1487: }

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

1492:   Collective on dm

1494:   Input Parameters:
1495: + dm - The original DM
1496: - label - The label specifying the boundary faces (this could be auto-generated)

1498:   Output Parameters:
1499: - dmSplit - The new DM

1501:   Level: developer

1503: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1504: @*/
1505: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1506: {
1507:   DM             sdm;
1508:   PetscInt       dim;

1514:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1515:   DMSetType(sdm, DMPLEX);
1516:   DMGetDimension(dm, &dim);
1517:   DMSetDimension(sdm, dim);
1518:   switch (dim) {
1519:   case 2:
1520:   case 3:
1521:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1522:     break;
1523:   default:
1524:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1525:   }
1526:   *dmSplit = sdm;
1527:   return(0);
1528: }

1530: /* Returns the side of the surface for a given cell with a face on the surface */
1531: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1532: {
1533:   const PetscInt *cone, *ornt;
1534:   PetscInt        dim, coneSize, c;
1535:   PetscErrorCode  ierr;

1538:   *pos = PETSC_TRUE;
1539:   DMGetDimension(dm, &dim);
1540:   DMPlexGetConeSize(dm, cell, &coneSize);
1541:   DMPlexGetCone(dm, cell, &cone);
1542:   DMPlexGetConeOrientation(dm, cell, &ornt);
1543:   for (c = 0; c < coneSize; ++c) {
1544:     if (cone[c] == face) {
1545:       PetscInt o = ornt[c];

1547:       if (subdm) {
1548:         const PetscInt *subcone, *subornt;
1549:         PetscInt        subpoint, subface, subconeSize, sc;

1551:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1552:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1553:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1554:         DMPlexGetCone(subdm, subpoint, &subcone);
1555:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1556:         for (sc = 0; sc < subconeSize; ++sc) {
1557:           if (subcone[sc] == subface) {
1558:             o = subornt[0];
1559:             break;
1560:           }
1561:         }
1562:         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);
1563:       }
1564:       if (o >= 0) *pos = PETSC_TRUE;
1565:       else        *pos = PETSC_FALSE;
1566:       break;
1567:     }
1568:   }
1569:   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);
1570:   return(0);
1571: }

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

1577:   Input Parameters:
1578: + dm     - The DM
1579: . label  - A DMLabel marking the surface
1580: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1581: . flip   - Flag to flip the submesh normal and replace points on the other side
1582: - subdm  - The subDM associated with the label, or NULL

1584:   Output Parameter:
1585: . label - A DMLabel marking all surface points

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

1589:   Level: developer

1591: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1592: @*/
1593: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1594: {
1595:   DMLabel         depthLabel;
1596:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1597:   const PetscInt *points, *subpoints;
1598:   const PetscInt  rev   = flip ? -1 : 1;
1599:   PetscInt       *pMax;
1600:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1601:   PetscErrorCode  ierr;

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

1631:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1632:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1633:     DMPlexGetSupport(dm, points[p], &support);
1634:     for (s = 0; s < supportSize; ++s) {
1635:       const PetscInt *cone;
1636:       PetscInt        coneSize, c;
1637:       PetscBool       pos;

1639:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1640:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1641:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1642:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1643:       /* Put faces touching the fault in the label */
1644:       DMPlexGetConeSize(dm, support[s], &coneSize);
1645:       DMPlexGetCone(dm, support[s], &cone);
1646:       for (c = 0; c < coneSize; ++c) {
1647:         const PetscInt point = cone[c];

1649:         DMLabelGetValue(label, point, &val);
1650:         if (val == -1) {
1651:           PetscInt *closure = NULL;
1652:           PetscInt  closureSize, cl;

1654:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1655:           for (cl = 0; cl < closureSize*2; cl += 2) {
1656:             const PetscInt clp  = closure[cl];
1657:             PetscInt       bval = -1;

1659:             DMLabelGetValue(label, clp, &val);
1660:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1661:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1662:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1663:               break;
1664:             }
1665:           }
1666:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1667:         }
1668:       }
1669:     }
1670:   }
1671:   ISRestoreIndices(dimIS, &points);
1672:   ISDestroy(&dimIS);
1673:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1674:   ISDestroy(&subpointIS);
1675:   /* Mark boundary points as unsplit */
1676:   if (blabel) {
1677:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1678:     ISGetLocalSize(dimIS, &numPoints);
1679:     ISGetIndices(dimIS, &points);
1680:     for (p = 0; p < numPoints; ++p) {
1681:       const PetscInt point = points[p];
1682:       PetscInt       val, bval;

1684:       DMLabelGetValue(blabel, point, &bval);
1685:       if (bval >= 0) {
1686:         DMLabelGetValue(label, point, &val);
1687:         if ((val < 0) || (val > dim)) {
1688:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1689:           DMLabelClearValue(blabel, point, bval);
1690:         }
1691:       }
1692:     }
1693:     for (p = 0; p < numPoints; ++p) {
1694:       const PetscInt point = points[p];
1695:       PetscInt       val, bval;

1697:       DMLabelGetValue(blabel, point, &bval);
1698:       if (bval >= 0) {
1699:         const PetscInt *cone,    *support;
1700:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

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

1737:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1738:     ISDestroy(&crossEdgeIS);
1739:     ISDestroy(&vertIS);
1740:   }
1741:   if (!dimIS) {
1742:     PetscFree(pMax);
1743:     return(0);
1744:   }
1745:   ISGetLocalSize(dimIS, &numPoints);
1746:   ISGetIndices(dimIS, &points);
1747:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1748:     PetscInt *star = NULL;
1749:     PetscInt  starSize, s;
1750:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1752:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1753:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1754:     while (again) {
1755:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1756:       again = 0;
1757:       for (s = 0; s < starSize*2; s += 2) {
1758:         const PetscInt  point = star[s];
1759:         const PetscInt *cone;
1760:         PetscInt        coneSize, c;

1762:         if ((point < cStart) || (point >= cMax)) continue;
1763:         DMLabelGetValue(label, point, &val);
1764:         if (val != -1) continue;
1765:         again = again == 1 ? 1 : 2;
1766:         DMPlexGetConeSize(dm, point, &coneSize);
1767:         DMPlexGetCone(dm, point, &cone);
1768:         for (c = 0; c < coneSize; ++c) {
1769:           DMLabelGetValue(label, cone[c], &val);
1770:           if (val != -1) {
1771:             const PetscInt *ccone;
1772:             PetscInt        cconeSize, cc, side;

1774:             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);
1775:             if (val > 0) side =  1;
1776:             else         side = -1;
1777:             DMLabelSetValue(label, point, side*(shift+dim));
1778:             /* Mark cell faces which touch the fault */
1779:             DMPlexGetConeSize(dm, point, &cconeSize);
1780:             DMPlexGetCone(dm, point, &ccone);
1781:             for (cc = 0; cc < cconeSize; ++cc) {
1782:               PetscInt *closure = NULL;
1783:               PetscInt  closureSize, cl;

1785:               DMLabelGetValue(label, ccone[cc], &val);
1786:               if (val != -1) continue;
1787:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1788:               for (cl = 0; cl < closureSize*2; cl += 2) {
1789:                 const PetscInt clp = closure[cl];

1791:                 DMLabelGetValue(label, clp, &val);
1792:                 if (val == -1) continue;
1793:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1794:                 break;
1795:               }
1796:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1797:             }
1798:             again = 1;
1799:             break;
1800:           }
1801:         }
1802:       }
1803:     }
1804:     /* Classify the rest by cell membership */
1805:     for (s = 0; s < starSize*2; s += 2) {
1806:       const PetscInt point = star[s];

1808:       DMLabelGetValue(label, point, &val);
1809:       if (val == -1) {
1810:         PetscInt  *sstar = NULL;
1811:         PetscInt   sstarSize, ss;
1812:         PetscBool  marked = PETSC_FALSE;

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

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

1852:     DMPlexGetSupportSize(dm, point, &supportSize);
1853:     if (supportSize != 2) continue;
1854:     DMPlexGetSupport(dm, point, &support);
1855:     DMLabelGetValue(label, support[0], &valA);
1856:     DMLabelGetValue(label, support[1], &valB);
1857:     if ((valA == -1) || (valB == -1)) continue;
1858:     if (valA*valB > 0) continue;
1859:     /* Split the face */
1860:     DMLabelGetValue(label, point, &valA);
1861:     DMLabelClearValue(label, point, valA);
1862:     DMLabelSetValue(label, point, dim-1);
1863:     /* Label its closure:
1864:       unmarked: label as unsplit
1865:       incident: relabel as split
1866:       split:    do nothing
1867:     */
1868:     {
1869:       PetscInt *closure = NULL;
1870:       PetscInt  closureSize, cl;

1872:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1873:       for (cl = 0; cl < closureSize*2; cl += 2) {
1874:         DMLabelGetValue(label, closure[cl], &valA);
1875:         if (valA == -1) { /* Mark as unsplit */
1876:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1877:           DMLabelSetValue(label, closure[cl], shift2+dep);
1878:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1879:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1880:           DMLabelClearValue(label, closure[cl], valA);
1881:           DMLabelSetValue(label, closure[cl], dep);
1882:         }
1883:       }
1884:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1885:     }
1886:   }
1887:   ISRestoreIndices(dimIS, &points);
1888:   ISDestroy(&dimIS);
1889:   PetscFree(pMax);
1890:   return(0);
1891: }

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

1896:   Collective on dm

1898:   Input Parameters:
1899: + dm - The original DM
1900: - labelName - The label specifying the interface vertices

1902:   Output Parameters:
1903: + hybridLabel - The label fully marking the interface
1904: - dmHybrid - The new DM

1906:   Level: developer

1908: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1909: @*/
1910: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1911: {
1912:   DM             idm;
1913:   DMLabel        subpointMap, hlabel;
1914:   PetscInt       dim;

1921:   DMGetDimension(dm, &dim);
1922:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1923:   DMPlexOrient(idm);
1924:   DMPlexGetSubpointMap(idm, &subpointMap);
1925:   DMLabelDuplicate(subpointMap, &hlabel);
1926:   DMLabelClearStratum(hlabel, dim);
1927:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1928:   DMDestroy(&idm);
1929:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1930:   if (hybridLabel) *hybridLabel = hlabel;
1931:   else             {DMLabelDestroy(&hlabel);}
1932:   return(0);
1933: }

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

1937:      For any marked cell, the marked vertices constitute a single face
1938: */
1939: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1940: {
1941:   IS               subvertexIS = NULL;
1942:   const PetscInt  *subvertices;
1943:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1944:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1945:   PetscErrorCode   ierr;

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

1970:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1971:     for (s = 0; s < starSize*2; s += 2) {
1972:       const PetscInt point = star[s];
1973:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1974:     }
1975:     for (c = 0; c < numCells; ++c) {
1976:       const PetscInt cell    = star[c];
1977:       PetscInt      *closure = NULL;
1978:       PetscInt       closureSize, cl;
1979:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1981:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1982:       if (cellLoc == 2) continue;
1983:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1984:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1985:       for (cl = 0; cl < closureSize*2; cl += 2) {
1986:         const PetscInt point = closure[cl];
1987:         PetscInt       vertexLoc;

1989:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1990:           ++numCorners;
1991:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1992:           if (vertexLoc == value) closure[faceSize++] = point;
1993:         }
1994:       }
1995:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1996:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1997:       if (faceSize == *nFV) {
1998:         const PetscInt *cells = NULL;
1999:         PetscInt        numCells, nc;

2001:         ++(*numFaces);
2002:         for (cl = 0; cl < faceSize; ++cl) {
2003:           DMLabelSetValue(subpointMap, closure[cl], 0);
2004:         }
2005:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2006:         for (nc = 0; nc < numCells; ++nc) {
2007:           DMLabelSetValue(subpointMap, cells[nc], 2);
2008:         }
2009:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2010:       }
2011:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2012:     }
2013:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2014:   }
2015:   if (subvertexIS) {
2016:     ISRestoreIndices(subvertexIS, &subvertices);
2017:   }
2018:   ISDestroy(&subvertexIS);
2019:   PetscFree3(pStart,pEnd,pMax);
2020:   return(0);
2021: }

2023: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2024: {
2025:   IS               subvertexIS = NULL;
2026:   const PetscInt  *subvertices;
2027:   PetscInt        *pStart, *pEnd, *pMax;
2028:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2029:   PetscErrorCode   ierr;

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

2052:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2053:     for (s = 0; s < starSize*2; s += 2) {
2054:       const PetscInt point = star[s];
2055:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2056:     }
2057:     for (f = 0; f < numFaces; ++f) {
2058:       const PetscInt face    = star[f];
2059:       PetscInt      *closure = NULL;
2060:       PetscInt       closureSize, c;
2061:       PetscInt       faceLoc;

2063:       DMLabelGetValue(subpointMap, face, &faceLoc);
2064:       if (faceLoc == dim-1) continue;
2065:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2066:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2067:       for (c = 0; c < closureSize*2; c += 2) {
2068:         const PetscInt point = closure[c];
2069:         PetscInt       vertexLoc;

2071:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2072:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2073:           if (vertexLoc != value) break;
2074:         }
2075:       }
2076:       if (c == closureSize*2) {
2077:         const PetscInt *support;
2078:         PetscInt        supportSize, s;

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

2083:           for (d = 0; d < dim; ++d) {
2084:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2085:               DMLabelSetValue(subpointMap, point, d);
2086:               break;
2087:             }
2088:           }
2089:         }
2090:         DMPlexGetSupportSize(dm, face, &supportSize);
2091:         DMPlexGetSupport(dm, face, &support);
2092:         for (s = 0; s < supportSize; ++s) {
2093:           DMLabelSetValue(subpointMap, support[s], dim);
2094:         }
2095:       }
2096:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2097:     }
2098:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2099:   }
2100:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2101:   ISDestroy(&subvertexIS);
2102:   PetscFree3(pStart,pEnd,pMax);
2103:   return(0);
2104: }

2106: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2107: {
2108:   DMLabel         label = NULL;
2109:   const PetscInt *cone;
2110:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2111:   PetscErrorCode  ierr;

2114:   *numFaces = 0;
2115:   *nFV = 0;
2116:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2117:   *subCells = NULL;
2118:   DMGetDimension(dm, &dim);
2119:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2120:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2121:   if (cMax < 0) return(0);
2122:   if (label) {
2123:     for (c = cMax; c < cEnd; ++c) {
2124:       PetscInt val;

2126:       DMLabelGetValue(label, c, &val);
2127:       if (val == value) {
2128:         ++(*numFaces);
2129:         DMPlexGetConeSize(dm, c, &coneSize);
2130:       }
2131:     }
2132:   } else {
2133:     *numFaces = cEnd - cMax;
2134:     DMPlexGetConeSize(dm, cMax, &coneSize);
2135:   }
2136:   PetscMalloc1(*numFaces *2, subCells);
2137:   if (!(*numFaces)) return(0);
2138:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2139:   for (c = cMax; c < cEnd; ++c) {
2140:     const PetscInt *cells;
2141:     PetscInt        numCells;

2143:     if (label) {
2144:       PetscInt val;

2146:       DMLabelGetValue(label, c, &val);
2147:       if (val != value) continue;
2148:     }
2149:     DMPlexGetCone(dm, c, &cone);
2150:     for (p = 0; p < *nFV; ++p) {
2151:       DMLabelSetValue(subpointMap, cone[p], 0);
2152:     }
2153:     /* Negative face */
2154:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2155:     /* Not true in parallel
2156:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2157:     for (p = 0; p < numCells; ++p) {
2158:       DMLabelSetValue(subpointMap, cells[p], 2);
2159:       (*subCells)[subc++] = cells[p];
2160:     }
2161:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2162:     /* Positive face is not included */
2163:   }
2164:   return(0);
2165: }

2167: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2168: {
2169:   PetscInt      *pStart, *pEnd;
2170:   PetscInt       dim, cMax, cEnd, c, d;

2174:   DMGetDimension(dm, &dim);
2175:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2176:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2177:   if (cMax < 0) return(0);
2178:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2179:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2180:   for (c = cMax; c < cEnd; ++c) {
2181:     const PetscInt *cone;
2182:     PetscInt       *closure = NULL;
2183:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2185:     if (label) {
2186:       DMLabelGetValue(label, c, &val);
2187:       if (val != value) continue;
2188:     }
2189:     DMPlexGetConeSize(dm, c, &coneSize);
2190:     DMPlexGetCone(dm, c, &cone);
2191:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2192:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2193:     /* Negative face */
2194:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2195:     for (cl = 0; cl < closureSize*2; cl += 2) {
2196:       const PetscInt point = closure[cl];

2198:       for (d = 0; d <= dim; ++d) {
2199:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2200:           DMLabelSetValue(subpointMap, point, d);
2201:           break;
2202:         }
2203:       }
2204:     }
2205:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2206:     /* Cells -- positive face is not included */
2207:     for (cl = 0; cl < 1; ++cl) {
2208:       const PetscInt *support;
2209:       PetscInt        supportSize, s;

2211:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2212:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2213:       DMPlexGetSupport(dm, cone[cl], &support);
2214:       for (s = 0; s < supportSize; ++s) {
2215:         DMLabelSetValue(subpointMap, support[s], dim);
2216:       }
2217:     }
2218:   }
2219:   PetscFree2(pStart, pEnd);
2220:   return(0);
2221: }

2223: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2224: {
2225:   MPI_Comm       comm;
2226:   PetscBool      posOrient = PETSC_FALSE;
2227:   const PetscInt debug     = 0;
2228:   PetscInt       cellDim, faceSize, f;

2232:   PetscObjectGetComm((PetscObject)dm,&comm);
2233:   DMGetDimension(dm, &cellDim);
2234:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

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

2281:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2282:     PetscSortInt(faceSizeTri, sortedIndices);
2283:     for (iFace = 0; iFace < 3; ++iFace) {
2284:       const PetscInt ii = iFace*faceSizeTri;
2285:       PetscInt       fVertex, cVertex;

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

2323:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2324:     PetscSortInt(faceSizeQuad, sortedIndices);
2325:     for (iFace = 0; iFace < 4; ++iFace) {
2326:       const PetscInt ii = iFace*faceSizeQuad;
2327:       PetscInt       fVertex, cVertex;

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

2351:           7---6
2352:          /|  /|
2353:         4---5 |
2354:         | 1-|-2
2355:         |/  |/
2356:         0---3

2358:         Faces are determined by the first 4 vertices (corners of faces) */
2359:     const PetscInt faceSizeHex = 4;
2360:     PetscInt       sortedIndices[4], i, iFace;
2361:     PetscBool      found                     = PETSC_FALSE;
2362:     PetscInt       faceVerticesHexSorted[24] = {
2363:       0, 1, 2, 3,  /* bottom */
2364:       4, 5, 6, 7,  /* top */
2365:       0, 3, 4, 5,  /* front */
2366:       2, 3, 5, 6,  /* right */
2367:       1, 2, 6, 7,  /* back */
2368:       0, 1, 4, 7,  /* left */
2369:     };
2370:     PetscInt       faceVerticesHex[24] = {
2371:       1, 2, 3, 0,  /* bottom */
2372:       4, 5, 6, 7,  /* top */
2373:       0, 3, 5, 4,  /* front */
2374:       3, 2, 6, 5,  /* right */
2375:       2, 1, 7, 6,  /* back */
2376:       1, 0, 4, 7,  /* left */
2377:     };

2379:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2380:     PetscSortInt(faceSizeHex, sortedIndices);
2381:     for (iFace = 0; iFace < 6; ++iFace) {
2382:       const PetscInt ii = iFace*faceSizeHex;
2383:       PetscInt       fVertex, cVertex;

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

2423:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2424:     PetscSortInt(faceSizeTet, sortedIndices);
2425:     for (iFace=0; iFace < 4; ++iFace) {
2426:       const PetscInt ii = iFace*faceSizeTet;
2427:       PetscInt       fVertex, cVertex;

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

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

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

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

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

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

2522:   Not collective

2524:   Input Parameters:
2525: + dm           - The original mesh
2526: . cell         - The cell mesh point
2527: . faceSize     - The number of vertices on the face
2528: . face         - The face vertices
2529: . numCorners   - The number of vertices on the cell
2530: . indices      - Local numbering of face vertices in cell cone
2531: - origVertices - Original face vertices

2533:   Output Parameter:
2534: + faceVertices - The face vertices properly oriented
2535: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2537:   Level: developer

2539: .seealso: DMPlexGetCone()
2540: @*/
2541: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2542: {
2543:   const PetscInt *cone = NULL;
2544:   PetscInt        coneSize, v, f, v2;
2545:   PetscInt        oppositeVertex = -1;
2546:   PetscErrorCode  ierr;

2549:   DMPlexGetConeSize(dm, cell, &coneSize);
2550:   DMPlexGetCone(dm, cell, &cone);
2551:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2552:     PetscBool found = PETSC_FALSE;

2554:     for (f = 0; f < faceSize; ++f) {
2555:       if (face[f] == cone[v]) {
2556:         found = PETSC_TRUE; break;
2557:       }
2558:     }
2559:     if (found) {
2560:       indices[v2]      = v;
2561:       origVertices[v2] = cone[v];
2562:       ++v2;
2563:     } else {
2564:       oppositeVertex = v;
2565:     }
2566:   }
2567:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2568:   return(0);
2569: }

2571: /*
2572:   DMPlexInsertFace_Internal - Puts a face into the mesh

2574:   Not collective

2576:   Input Parameters:
2577:   + dm              - The DMPlex
2578:   . numFaceVertex   - The number of vertices in the face
2579:   . faceVertices    - The vertices in the face for dm
2580:   . subfaceVertices - The vertices in the face for subdm
2581:   . numCorners      - The number of vertices in the cell
2582:   . cell            - A cell in dm containing the face
2583:   . subcell         - A cell in subdm containing the face
2584:   . firstFace       - First face in the mesh
2585:   - newFacePoint    - Next face in the mesh

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

2590:   Level: developer
2591: */
2592: 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)
2593: {
2594:   MPI_Comm        comm;
2595:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2596:   const PetscInt *faces;
2597:   PetscInt        numFaces, coneSize;
2598:   PetscErrorCode  ierr;

2601:   PetscObjectGetComm((PetscObject)dm,&comm);
2602:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2603:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2604: #if 0
2605:   /* Cannot use this because support() has not been constructed yet */
2606:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2607: #else
2608:   {
2609:     PetscInt f;

2611:     numFaces = 0;
2612:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2613:     for (f = firstFace; f < *newFacePoint; ++f) {
2614:       PetscInt dof, off, d;

2616:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2617:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2618:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2619:       for (d = 0; d < dof; ++d) {
2620:         const PetscInt p = submesh->cones[off+d];
2621:         PetscInt       v;

2623:         for (v = 0; v < numFaceVertices; ++v) {
2624:           if (subfaceVertices[v] == p) break;
2625:         }
2626:         if (v == numFaceVertices) break;
2627:       }
2628:       if (d == dof) {
2629:         numFaces               = 1;
2630:         ((PetscInt*) faces)[0] = f;
2631:       }
2632:     }
2633:   }
2634: #endif
2635:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2636:   else if (numFaces == 1) {
2637:     /* Add the other cell neighbor for this face */
2638:     DMPlexSetCone(subdm, subcell, faces);
2639:   } else {
2640:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2641:     PetscBool posOriented;

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

2672: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2673: {
2674:   MPI_Comm        comm;
2675:   DMLabel         subpointMap;
2676:   IS              subvertexIS,  subcellIS;
2677:   const PetscInt *subVertices, *subCells;
2678:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2679:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2680:   PetscInt        vStart, vEnd, c, f;
2681:   PetscErrorCode  ierr;

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

2720:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2721:     for (cl = 0; cl < closureSize*2; cl += 2) {
2722:       const PetscInt point = closure[cl];
2723:       PetscInt       subVertex;

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

2752:     DMGetCoordinateSection(dm, &coordSection);
2753:     DMGetCoordinatesLocal(dm, &coordinates);
2754:     DMGetCoordinateSection(subdm, &subCoordSection);
2755:     PetscSectionSetNumFields(subCoordSection, 1);
2756:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2757:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2758:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2759:     for (v = 0; v < numSubVertices; ++v) {
2760:       const PetscInt vertex    = subVertices[v];
2761:       const PetscInt subvertex = firstSubVertex+v;
2762:       PetscInt       dof;

2764:       PetscSectionGetDof(coordSection, vertex, &dof);
2765:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2766:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2767:     }
2768:     PetscSectionSetUp(subCoordSection);
2769:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2770:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2771:     PetscObjectGetName((PetscObject)coordinates,&name);
2772:     PetscObjectSetName((PetscObject)subCoordinates,name);
2773:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2774:     VecSetType(subCoordinates,VECSTANDARD);
2775:     if (coordSize) {
2776:       VecGetArray(coordinates,    &coords);
2777:       VecGetArray(subCoordinates, &subCoords);
2778:       for (v = 0; v < numSubVertices; ++v) {
2779:         const PetscInt vertex    = subVertices[v];
2780:         const PetscInt subvertex = firstSubVertex+v;
2781:         PetscInt       dof, off, sdof, soff, d;

2783:         PetscSectionGetDof(coordSection, vertex, &dof);
2784:         PetscSectionGetOffset(coordSection, vertex, &off);
2785:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2786:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2787:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2788:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2789:       }
2790:       VecRestoreArray(coordinates,    &coords);
2791:       VecRestoreArray(subCoordinates, &subCoords);
2792:     }
2793:     DMSetCoordinatesLocal(subdm, subCoordinates);
2794:     VecDestroy(&subCoordinates);
2795:   }
2796:   /* Cleanup */
2797:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2798:   ISDestroy(&subvertexIS);
2799:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2800:   ISDestroy(&subcellIS);
2801:   return(0);
2802: }

2804: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2805: {
2806:   PetscInt       subPoint;

2809:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2810:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2811: }

2813: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2814: {
2815:   MPI_Comm         comm;
2816:   DMLabel          subpointMap;
2817:   IS              *subpointIS;
2818:   const PetscInt **subpoints;
2819:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2820:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2821:   PetscMPIInt      rank;
2822:   PetscErrorCode   ierr;

2825:   PetscObjectGetComm((PetscObject)dm,&comm);
2826:   MPI_Comm_rank(comm, &rank);
2827:   /* Create subpointMap which marks the submesh */
2828:   DMLabelCreate("subpoint_map", &subpointMap);
2829:   DMPlexSetSubpointMap(subdm, subpointMap);
2830:   if (cellHeight) {
2831:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2832:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2833:   } else {
2834:     DMLabel         depth;
2835:     IS              pointIS;
2836:     const PetscInt *points;
2837:     PetscInt        numPoints;

2839:     DMPlexGetDepthLabel(dm, &depth);
2840:     DMLabelGetStratumSize(label, value, &numPoints);
2841:     DMLabelGetStratumIS(label, value, &pointIS);
2842:     ISGetIndices(pointIS, &points);
2843:     for (p = 0; p < numPoints; ++p) {
2844:       PetscInt *closure = NULL;
2845:       PetscInt  closureSize, c, pdim;

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

2886:       DMPlexGetConeSize(dm, point, &coneSize);
2887:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2888:       if (cellHeight && (d == dim)) {
2889:         DMPlexGetCone(dm, point, &cone);
2890:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2891:           DMLabelGetValue(subpointMap, cone[c], &val);
2892:           if (val >= 0) coneSizeNew++;
2893:         }
2894:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2895:       }
2896:     }
2897:   }
2898:   DMSetUp(subdm);
2899:   /* Set cones */
2900:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2901:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2902:   for (d = 0; d <= dim; ++d) {
2903:     for (p = 0; p < numSubPoints[d]; ++p) {
2904:       const PetscInt  point    = subpoints[d][p];
2905:       const PetscInt  subpoint = firstSubPoint[d] + p;
2906:       const PetscInt *cone, *ornt;
2907:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

2909:       if (d == dim-1) {
2910:         const PetscInt *support, *cone, *ornt;
2911:         PetscInt        supportSize, coneSize, s, subc;

2913:         DMPlexGetSupport(dm, point, &support);
2914:         DMPlexGetSupportSize(dm, point, &supportSize);
2915:         for (s = 0; s < supportSize; ++s) {
2916:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2917:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2918:           if (subc >= 0) {
2919:             const PetscInt ccell = subpoints[d+1][subc];

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

2955:           tmp        = coneNew[c];
2956:           coneNew[c] = coneNew[coneSizeNew-1-c];
2957:           coneNew[coneSizeNew-1-c] = tmp;
2958:           DMPlexGetConeSize(dm, cone[c], &faceSize);
2959:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2960:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2961:           orntNew[coneSizeNew-1-c] = tmp;
2962:         }
2963:       }
2964:       DMPlexSetCone(subdm, subpoint, coneNew);
2965:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2966: #endif
2967:     }
2968:   }
2969:   PetscFree2(coneNew,orntNew);
2970:   DMPlexSymmetrize(subdm);
2971:   DMPlexStratify(subdm);
2972:   /* Build coordinates */
2973:   {
2974:     PetscSection coordSection, subCoordSection;
2975:     Vec          coordinates, subCoordinates;
2976:     PetscScalar *coords, *subCoords;
2977:     PetscInt     cdim, numComp, coordSize;
2978:     const char  *name;

2980:     DMGetCoordinateDim(dm, &cdim);
2981:     DMGetCoordinateSection(dm, &coordSection);
2982:     DMGetCoordinatesLocal(dm, &coordinates);
2983:     DMGetCoordinateSection(subdm, &subCoordSection);
2984:     PetscSectionSetNumFields(subCoordSection, 1);
2985:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2986:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2987:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2988:     for (v = 0; v < numSubPoints[0]; ++v) {
2989:       const PetscInt vertex    = subpoints[0][v];
2990:       const PetscInt subvertex = firstSubPoint[0]+v;
2991:       PetscInt       dof;

2993:       PetscSectionGetDof(coordSection, vertex, &dof);
2994:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2995:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2996:     }
2997:     PetscSectionSetUp(subCoordSection);
2998:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2999:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3000:     PetscObjectGetName((PetscObject)coordinates,&name);
3001:     PetscObjectSetName((PetscObject)subCoordinates,name);
3002:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3003:     VecSetBlockSize(subCoordinates, cdim);
3004:     VecSetType(subCoordinates,VECSTANDARD);
3005:     VecGetArray(coordinates,    &coords);
3006:     VecGetArray(subCoordinates, &subCoords);
3007:     for (v = 0; v < numSubPoints[0]; ++v) {
3008:       const PetscInt vertex    = subpoints[0][v];
3009:       const PetscInt subvertex = firstSubPoint[0]+v;
3010:       PetscInt dof, off, sdof, soff, d;

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

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

3057:         if (subpoint < 0) continue;
3058:         newLocalPoints[point-pStart].rank  = rank;
3059:         newLocalPoints[point-pStart].index = subpoint;
3060:         ++numSubleaves;
3061:       }
3062:       /* Must put in owned subpoints */
3063:       for (p = pStart; p < pEnd; ++p) {
3064:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3066:         if (subpoint < 0) {
3067:           newOwners[p-pStart].rank  = -3;
3068:           newOwners[p-pStart].index = -3;
3069:         } else {
3070:           newOwners[p-pStart].rank  = rank;
3071:           newOwners[p-pStart].index = subpoint;
3072:         }
3073:       }
3074:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3075:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3076:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3077:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3078:       PetscMalloc1(numSubleaves, &slocalPoints);
3079:       PetscMalloc1(numSubleaves, &sremotePoints);
3080:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3081:         const PetscInt point    = localPoints[l];
3082:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

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

3111: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3112: {

3116:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3117:   return(0);
3118: }

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

3123:   Input Parameters:
3124: + dm           - The original mesh
3125: . vertexLabel  - The DMLabel marking vertices contained in the surface
3126: - value        - The label value to use

3128:   Output Parameter:
3129: . subdm - The surface mesh

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

3133:   Level: developer

3135: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3136: @*/
3137: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3138: {
3139:   PetscInt       dim, cdim, depth;

3145:   DMGetDimension(dm, &dim);
3146:   DMPlexGetDepth(dm, &depth);
3147:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3148:   DMSetType(*subdm, DMPLEX);
3149:   DMSetDimension(*subdm, dim-1);
3150:   DMGetCoordinateDim(dm, &cdim);
3151:   DMSetCoordinateDim(*subdm, cdim);
3152:   if (depth == dim) {
3153:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3154:   } else {
3155:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3156:   }
3157:   return(0);
3158: }

3160: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3161: {
3162:   MPI_Comm        comm;
3163:   DMLabel         subpointMap;
3164:   IS              subvertexIS;
3165:   const PetscInt *subVertices;
3166:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3167:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3168:   PetscInt        cMax, c, f;
3169:   PetscErrorCode  ierr;

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

3206:     if (cell < cMax) continue;
3207:     DMPlexGetCone(dm, cell, &cone);
3208:     for (v = 0; v < nFV; ++v) {
3209:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3210:       subface[v] = firstSubVertex+subVertex;
3211:     }
3212:     DMPlexSetCone(subdm, newFacePoint, subface);
3213:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3214:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3215:     /* Not true in parallel
3216:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3217:     for (p = 0; p < numCells; ++p) {
3218:       PetscInt negsubcell;

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

3242:     DMGetCoordinateDim(dm, &cdim);
3243:     DMGetCoordinateSection(dm, &coordSection);
3244:     DMGetCoordinatesLocal(dm, &coordinates);
3245:     DMGetCoordinateSection(subdm, &subCoordSection);
3246:     PetscSectionSetNumFields(subCoordSection, 1);
3247:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3248:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3249:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3250:     for (v = 0; v < numSubVertices; ++v) {
3251:       const PetscInt vertex    = subVertices[v];
3252:       const PetscInt subvertex = firstSubVertex+v;
3253:       PetscInt       dof;

3255:       PetscSectionGetDof(coordSection, vertex, &dof);
3256:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3257:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3258:     }
3259:     PetscSectionSetUp(subCoordSection);
3260:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3261:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3262:     PetscObjectGetName((PetscObject)coordinates,&name);
3263:     PetscObjectSetName((PetscObject)subCoordinates,name);
3264:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3265:     VecSetBlockSize(subCoordinates, cdim);
3266:     VecSetType(subCoordinates,VECSTANDARD);
3267:     VecGetArray(coordinates,    &coords);
3268:     VecGetArray(subCoordinates, &subCoords);
3269:     for (v = 0; v < numSubVertices; ++v) {
3270:       const PetscInt vertex    = subVertices[v];
3271:       const PetscInt subvertex = firstSubVertex+v;
3272:       PetscInt       dof, off, sdof, soff, d;

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

3297:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3298:     DMGetPointSF(dm, &sfPoint);
3299:     DMGetPointSF(subdm, &sfPointSub);
3300:     DMPlexGetChart(dm, &pStart, &pEnd);
3301:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3302:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3303:     if (numRoots >= 0) {
3304:       /* Only vertices should be shared */
3305:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3306:       for (p = 0; p < pEnd-pStart; ++p) {
3307:         newLocalPoints[p].rank  = -2;
3308:         newLocalPoints[p].index = -2;
3309:       }
3310:       /* Set subleaves */
3311:       for (l = 0; l < numLeaves; ++l) {
3312:         const PetscInt point    = localPoints[l];
3313:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3315:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3316:         if (subPoint < 0) continue;
3317:         newLocalPoints[point-pStart].rank  = rank;
3318:         newLocalPoints[point-pStart].index = subPoint;
3319:         ++numSubLeaves;
3320:       }
3321:       /* Must put in owned subpoints */
3322:       for (p = pStart; p < pEnd; ++p) {
3323:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3325:         if (subPoint < 0) {
3326:           newOwners[p-pStart].rank  = -3;
3327:           newOwners[p-pStart].index = -3;
3328:         } else {
3329:           newOwners[p-pStart].rank  = rank;
3330:           newOwners[p-pStart].index = subPoint;
3331:         }
3332:       }
3333:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3334:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3335:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3336:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3337:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3338:       PetscMalloc1(numSubLeaves, &sremotePoints);
3339:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3340:         const PetscInt point    = localPoints[l];
3341:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

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

3365: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3366: {
3367:   DMLabel        label = NULL;

3371:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3372:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3373:   return(0);
3374: }

3376: /*@C
3377:   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.

3379:   Input Parameters:
3380: + dm          - The original mesh
3381: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3382: . label       - A label name, or NULL
3383: - value  - A label value

3385:   Output Parameter:
3386: . subdm - The surface mesh

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

3390:   Level: developer

3392: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3393: @*/
3394: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3395: {
3396:   PetscInt       dim, cdim, depth;

3402:   DMGetDimension(dm, &dim);
3403:   DMPlexGetDepth(dm, &depth);
3404:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3405:   DMSetType(*subdm, DMPLEX);
3406:   DMSetDimension(*subdm, dim-1);
3407:   DMGetCoordinateDim(dm, &cdim);
3408:   DMSetCoordinateDim(*subdm, cdim);
3409:   if (depth == dim) {
3410:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3411:   } else {
3412:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3413:   }
3414:   return(0);
3415: }

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

3420:   Input Parameters:
3421: + dm        - The original mesh
3422: . cellLabel - The DMLabel marking cells contained in the new mesh
3423: - value     - The label value to use

3425:   Output Parameter:
3426: . subdm - The new mesh

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

3430:   Level: developer

3432: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3433: @*/
3434: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3435: {
3436:   PetscInt       dim;

3442:   DMGetDimension(dm, &dim);
3443:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3444:   DMSetType(*subdm, DMPLEX);
3445:   DMSetDimension(*subdm, dim);
3446:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3447:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3448:   return(0);
3449: }

3451: /*@
3452:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3454:   Input Parameter:
3455: . dm - The submesh DM

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

3460:   Level: developer

3462: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3463: @*/
3464: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3465: {
3469:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3470:   return(0);
3471: }

3473: /*@
3474:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3476:   Input Parameters:
3477: + dm - The submesh DM
3478: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3482:   Level: developer

3484: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3485: @*/
3486: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3487: {
3488:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3489:   DMLabel        tmp;

3494:   tmp  = mesh->subpointMap;
3495:   mesh->subpointMap = subpointMap;
3496:   ++mesh->subpointMap->refct;
3497:   DMLabelDestroy(&tmp);
3498:   return(0);
3499: }

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

3504:   Input Parameter:
3505: . dm - The submesh DM

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

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

3512:   Level: developer

3514: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3515: @*/
3516: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3517: {
3518:   MPI_Comm        comm;
3519:   DMLabel         subpointMap;
3520:   IS              is;
3521:   const PetscInt *opoints;
3522:   PetscInt       *points, *depths;
3523:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3524:   PetscErrorCode  ierr;

3529:   PetscObjectGetComm((PetscObject)dm,&comm);
3530:   *subpointIS = NULL;
3531:   DMPlexGetSubpointMap(dm, &subpointMap);
3532:   DMPlexGetDepth(dm, &depth);
3533:   if (subpointMap && depth >= 0) {
3534:     DMPlexGetChart(dm, &pStart, &pEnd);
3535:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3536:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3537:     depths[0] = depth;
3538:     depths[1] = 0;
3539:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3540:     PetscMalloc1(pEnd, &points);
3541:     for(d = 0, off = 0; d <= depth; ++d) {
3542:       const PetscInt dep = depths[d];

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