Actual source code: plexsubmesh.c

petsc-master 2017-12-14
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 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) {DMLabelSetValue(label, f, 1);}
 18:   }
 19:   return(0);
 20: }

 22: /*@
 23:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 25:   Not Collective

 27:   Input Parameter:
 28: . dm - The original DM

 30:   Output Parameter:
 31: . label - The DMLabel marking boundary faces with value 1

 33:   Level: developer

 35: .seealso: DMLabelCreate(), DMCreateLabel()
 36: @*/
 37: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
 38: {

 42:   DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
 43:   return(0);
 44: }

 46: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
 47: {
 48:   IS              valueIS;
 49:   PetscSF         sfPoint;
 50:   const PetscInt *values;
 51:   PetscInt        numValues, v, cStart, cEnd, nroots;
 52:   PetscErrorCode  ierr;

 55:   DMLabelGetNumValues(label, &numValues);
 56:   DMLabelGetValueIS(label, &valueIS);
 57:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 58:   ISGetIndices(valueIS, &values);
 59:   for (v = 0; v < numValues; ++v) {
 60:     IS              pointIS;
 61:     const PetscInt *points;
 62:     PetscInt        numPoints, p;

 64:     DMLabelGetStratumSize(label, values[v], &numPoints);
 65:     DMLabelGetStratumIS(label, values[v], &pointIS);
 66:     ISGetIndices(pointIS, &points);
 67:     for (p = 0; p < numPoints; ++p) {
 68:       PetscInt  q = points[p];
 69:       PetscInt *closure = NULL;
 70:       PetscInt  closureSize, c;

 72:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
 73:         continue;
 74:       }
 75:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
 76:       for (c = 0; c < closureSize*2; c += 2) {
 77:         DMLabelSetValue(label, closure[c], values[v]);
 78:       }
 79:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
 80:     }
 81:     ISRestoreIndices(pointIS, &points);
 82:     ISDestroy(&pointIS);
 83:   }
 84:   ISRestoreIndices(valueIS, &values);
 85:   ISDestroy(&valueIS);
 86:   DMGetPointSF(dm, &sfPoint);
 87:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
 88:   if (nroots >= 0) {
 89:     DMLabel         lblRoots, lblLeaves;
 90:     IS              valueIS, pointIS;
 91:     const PetscInt *values;
 92:     PetscInt        numValues, v;
 93:     PetscErrorCode  ierr;

 95:     DMGetPointSF(dm, &sfPoint);
 96:     /* Pull point contributions from remote leaves into local roots */
 97:     DMLabelGather(label, sfPoint, &lblLeaves);
 98:     DMLabelGetValueIS(lblLeaves, &valueIS);
 99:     ISGetLocalSize(valueIS, &numValues);
100:     ISGetIndices(valueIS, &values);
101:     for (v = 0; v < numValues; ++v) {
102:       const PetscInt value = values[v];
103: 
104:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
105:       DMLabelInsertIS(label, pointIS, value);
106:       ISDestroy(&pointIS);
107:     }
108:     ISRestoreIndices(valueIS, &values);
109:     ISDestroy(&valueIS);
110:     DMLabelDestroy(&lblLeaves);
111:     /* Push point contributions from roots into remote leaves */
112:     DMLabelDistribute(label, sfPoint, &lblRoots);
113:     DMLabelGetValueIS(lblRoots, &valueIS);
114:     ISGetLocalSize(valueIS, &numValues);
115:     ISGetIndices(valueIS, &values);
116:     for (v = 0; v < numValues; ++v) {
117:       const PetscInt value = values[v];
118: 
119:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
120:       DMLabelInsertIS(label, pointIS, value);
121:       ISDestroy(&pointIS);
122:     }
123:     ISRestoreIndices(valueIS, &values);
124:     ISDestroy(&valueIS);
125:     DMLabelDestroy(&lblRoots);
126:   }
127:   return(0);
128: }

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

133:   Input Parameters:
134: + dm - The DM
135: - label - A DMLabel marking the surface points

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

140:   Level: developer

142: .seealso: DMPlexLabelCohesiveComplete()
143: @*/
144: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
145: {

149:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
150:   return(0);
151: }

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

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

160:   Output Parameter:
161: . label - A DMLabel incorporating cells

163:   Level: developer

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

167: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
168: @*/
169: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
170: {
171:   IS              valueIS;
172:   const PetscInt *values;
173:   PetscInt        numValues, v, cStart, cEnd, cEndInterior;
174:   PetscErrorCode  ierr;

177:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
178:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
179:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
180:   DMLabelGetNumValues(label, &numValues);
181:   DMLabelGetValueIS(label, &valueIS);
182:   ISGetIndices(valueIS, &values);
183:   for (v = 0; v < numValues; ++v) {
184:     IS              pointIS;
185:     const PetscInt *points;
186:     PetscInt        numPoints, p;

188:     DMLabelGetStratumSize(label, values[v], &numPoints);
189:     DMLabelGetStratumIS(label, values[v], &pointIS);
190:     ISGetIndices(pointIS, &points);
191:     for (p = 0; p < numPoints; ++p) {
192:       PetscInt *closure = NULL;
193:       PetscInt  closureSize, point, cl;

195:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
196:       for (cl = closureSize-1; cl > 0; --cl) {
197:         point = closure[cl*2];
198:         if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
199:       }
200:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
201:     }
202:     ISRestoreIndices(pointIS, &points);
203:     ISDestroy(&pointIS);
204:   }
205:   ISRestoreIndices(valueIS, &values);
206:   ISDestroy(&valueIS);
207:   return(0);
208: }

210: /*@
211:   DMPlexLabelClearCells - Remove cells from a label

213:   Input Parameters:
214: + dm - The DM
215: - label - A DMLabel marking surface points and their adjacent cells

217:   Output Parameter:
218: . label - A DMLabel without cells

220:   Level: developer

222:   Note: This undoes DMPlexLabelAddCells()

224: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
225: @*/
226: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
227: {
228:   IS              valueIS;
229:   const PetscInt *values;
230:   PetscInt        numValues, v, cStart, cEnd, cEndInterior;
231:   PetscErrorCode  ierr;

234:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
235:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
236:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
237:   DMLabelGetNumValues(label, &numValues);
238:   DMLabelGetValueIS(label, &valueIS);
239:   ISGetIndices(valueIS, &values);
240:   for (v = 0; v < numValues; ++v) {
241:     IS              pointIS;
242:     const PetscInt *points;
243:     PetscInt        numPoints, p;

245:     DMLabelGetStratumSize(label, values[v], &numPoints);
246:     DMLabelGetStratumIS(label, values[v], &pointIS);
247:     ISGetIndices(pointIS, &points);
248:     for (p = 0; p < numPoints; ++p) {
249:       PetscInt point = points[p];

251:       if (point >= cStart && point < cEnd) {
252:         DMLabelClearValue(label,point,values[v]);
253:       }
254:     }
255:     ISRestoreIndices(pointIS, &points);
256:     ISDestroy(&pointIS);
257:   }
258:   ISRestoreIndices(valueIS, &values);
259:   ISDestroy(&valueIS);
260:   return(0);
261: }

263: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
264:  * index (skipping first, which is (0,0)) */
265: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
266: {
267:   PetscInt d, off = 0;

270:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
271:   for (d = 0; d < depth; d++) {
272:     PetscInt firstd = d;
273:     PetscInt firstStart = depthShift[2*d];
274:     PetscInt e;

276:     for (e = d+1; e <= depth; e++) {
277:       if (depthShift[2*e] < firstStart) {
278:         firstd = e;
279:         firstStart = depthShift[2*d];
280:       }
281:     }
282:     if (firstd != d) {
283:       PetscInt swap[2];

285:       e = firstd;
286:       swap[0] = depthShift[2*d];
287:       swap[1] = depthShift[2*d+1];
288:       depthShift[2*d]   = depthShift[2*e];
289:       depthShift[2*d+1] = depthShift[2*e+1];
290:       depthShift[2*e]   = swap[0];
291:       depthShift[2*e+1] = swap[1];
292:     }
293:   }
294:   /* convert (oldstart, added) to (oldstart, newstart) */
295:   for (d = 0; d <= depth; d++) {
296:     off += depthShift[2*d+1];
297:     depthShift[2*d+1] = depthShift[2*d] + off;
298:   }
299:   return(0);
300: }

302: /* depthShift is a list of (old, new) pairs */
303: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
304: {
305:   PetscInt d;
306:   PetscInt newOff = 0;

308:   for (d = 0; d <= depth; d++) {
309:     if (p < depthShift[2*d]) return p + newOff;
310:     else newOff = depthShift[2*d+1] - depthShift[2*d];
311:   }
312:   return p + newOff;
313: }

315: /* depthShift is a list of (old, new) pairs */
316: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
317: {
318:   PetscInt d;
319:   PetscInt newOff = 0;

321:   for (d = 0; d <= depth; d++) {
322:     if (p < depthShift[2*d+1]) return p + newOff;
323:     else newOff = depthShift[2*d] - depthShift[2*d+1];
324:   }
325:   return p + newOff;
326: }

328: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
329: {
330:   PetscInt       depth = 0, d, pStart, pEnd, p;

334:   DMPlexGetDepth(dm, &depth);
335:   if (depth < 0) return(0);
336:   /* Step 1: Expand chart */
337:   DMPlexGetChart(dm, &pStart, &pEnd);
338:   pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
339:   DMPlexSetChart(dmNew, pStart, pEnd);
340:   /* Step 2: Set cone and support sizes */
341:   for (d = 0; d <= depth; ++d) {
342:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
343:     for (p = pStart; p < pEnd; ++p) {
344:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
345:       PetscInt size;

347:       DMPlexGetConeSize(dm, p, &size);
348:       DMPlexSetConeSize(dmNew, newp, size);
349:       DMPlexGetSupportSize(dm, p, &size);
350:       DMPlexSetSupportSize(dmNew, newp, size);
351:     }
352:   }
353:   return(0);
354: }

356: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
357: {
358:   PetscInt      *newpoints;
359:   PetscInt       depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

363:   DMPlexGetDepth(dm, &depth);
364:   if (depth < 0) return(0);
365:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
366:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
367:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
368:   /* Step 5: Set cones and supports */
369:   DMPlexGetChart(dm, &pStart, &pEnd);
370:   for (p = pStart; p < pEnd; ++p) {
371:     const PetscInt *points = NULL, *orientations = NULL;
372:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

374:     DMPlexGetConeSize(dm, p, &size);
375:     DMPlexGetCone(dm, p, &points);
376:     DMPlexGetConeOrientation(dm, p, &orientations);
377:     for (i = 0; i < size; ++i) {
378:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
379:     }
380:     DMPlexSetCone(dmNew, newp, newpoints);
381:     DMPlexSetConeOrientation(dmNew, newp, orientations);
382:     DMPlexGetSupportSize(dm, p, &size);
383:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
384:     DMPlexGetSupport(dm, p, &points);
385:     for (i = 0; i < size; ++i) {
386:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
387:     }
388:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
389:     DMPlexSetSupport(dmNew, newp, newpoints);
390:   }
391:   PetscFree(newpoints);
392:   return(0);
393: }

395: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
396: {
397:   PetscSection   coordSection, newCoordSection;
398:   Vec            coordinates, newCoordinates;
399:   PetscScalar   *coords, *newCoords;
400:   PetscInt       coordSize, sStart, sEnd;
401:   PetscInt       dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
402:   PetscBool      hasCells;

406:   DMGetCoordinateDim(dm, &dim);
407:   DMSetCoordinateDim(dmNew, dim);
408:   DMPlexGetDepth(dm, &depth);
409:   /* Step 8: Convert coordinates */
410:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
411:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
412:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
413:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
414:   DMGetCoordinateSection(dm, &coordSection);
415:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
416:   PetscSectionSetNumFields(newCoordSection, 1);
417:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
418:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
419:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
420:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
421:   if (hasCells) {
422:     for (c = cStart; c < cEnd; ++c) {
423:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

425:       PetscSectionGetDof(coordSection, c, &dof);
426:       PetscSectionSetDof(newCoordSection, cNew, dof);
427:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
428:     }
429:   }
430:   for (v = vStartNew; v < vEndNew; ++v) {
431:     PetscSectionSetDof(newCoordSection, v, dim);
432:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
433:   }
434:   PetscSectionSetUp(newCoordSection);
435:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
436:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
437:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
438:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
439:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
440:   VecSetBlockSize(newCoordinates, dim);
441:   VecSetType(newCoordinates,VECSTANDARD);
442:   DMSetCoordinatesLocal(dmNew, newCoordinates);
443:   DMGetCoordinatesLocal(dm, &coordinates);
444:   VecGetArray(coordinates, &coords);
445:   VecGetArray(newCoordinates, &newCoords);
446:   if (hasCells) {
447:     for (c = cStart; c < cEnd; ++c) {
448:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

450:       PetscSectionGetDof(coordSection, c, &dof);
451:       PetscSectionGetOffset(coordSection, c, &off);
452:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
453:       for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
454:     }
455:   }
456:   for (v = vStart; v < vEnd; ++v) {
457:     PetscInt dof, off, noff, d;

459:     PetscSectionGetDof(coordSection, v, &dof);
460:     PetscSectionGetOffset(coordSection, v, &off);
461:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
462:     for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
463:   }
464:   VecRestoreArray(coordinates, &coords);
465:   VecRestoreArray(newCoordinates, &newCoords);
466:   VecDestroy(&newCoordinates);
467:   PetscSectionDestroy(&newCoordSection);
468:   return(0);
469: }

471: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
472: {
473:   PetscInt           depth = 0;
474:   PetscSF            sfPoint, sfPointNew;
475:   const PetscSFNode *remotePoints;
476:   PetscSFNode       *gremotePoints;
477:   const PetscInt    *localPoints;
478:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
479:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
480:   PetscErrorCode     ierr;

483:   DMPlexGetDepth(dm, &depth);
484:   /* Step 9: Convert pointSF */
485:   DMGetPointSF(dm, &sfPoint);
486:   DMGetPointSF(dmNew, &sfPointNew);
487:   DMPlexGetChart(dm, &pStart, &pEnd);
488:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
489:   totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
490:   if (numRoots >= 0) {
491:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
492:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
493:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
494:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
495:     PetscMalloc1(numLeaves,    &glocalPoints);
496:     PetscMalloc1(numLeaves, &gremotePoints);
497:     for (l = 0; l < numLeaves; ++l) {
498:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
499:       gremotePoints[l].rank  = remotePoints[l].rank;
500:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
501:     }
502:     PetscFree2(newLocation,newRemoteLocation);
503:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
504:   }
505:   return(0);
506: }

508: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
509: {
510:   PetscSF            sfPoint;
511:   DMLabel            vtkLabel, ghostLabel;
512:   const PetscSFNode *leafRemote;
513:   const PetscInt    *leafLocal;
514:   PetscInt           depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
515:   PetscMPIInt        rank;
516:   PetscErrorCode     ierr;

519:   DMPlexGetDepth(dm, &depth);
520:   /* Step 10: Convert labels */
521:   DMGetNumLabels(dm, &numLabels);
522:   for (l = 0; l < numLabels; ++l) {
523:     DMLabel         label, newlabel;
524:     const char     *lname;
525:     PetscBool       isDepth;
526:     IS              valueIS;
527:     const PetscInt *values;
528:     PetscInt        numValues, val;

530:     DMGetLabelName(dm, l, &lname);
531:     PetscStrcmp(lname, "depth", &isDepth);
532:     if (isDepth) continue;
533:     DMCreateLabel(dmNew, lname);
534:     DMGetLabel(dm, lname, &label);
535:     DMGetLabel(dmNew, lname, &newlabel);
536:     DMLabelGetDefaultValue(label,&val);
537:     DMLabelSetDefaultValue(newlabel,val);
538:     DMLabelGetValueIS(label, &valueIS);
539:     ISGetLocalSize(valueIS, &numValues);
540:     ISGetIndices(valueIS, &values);
541:     for (val = 0; val < numValues; ++val) {
542:       IS              pointIS;
543:       const PetscInt *points;
544:       PetscInt        numPoints, p;

546:       DMLabelGetStratumIS(label, values[val], &pointIS);
547:       ISGetLocalSize(pointIS, &numPoints);
548:       ISGetIndices(pointIS, &points);
549:       for (p = 0; p < numPoints; ++p) {
550:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

552:         DMLabelSetValue(newlabel, newpoint, values[val]);
553:       }
554:       ISRestoreIndices(pointIS, &points);
555:       ISDestroy(&pointIS);
556:     }
557:     ISRestoreIndices(valueIS, &values);
558:     ISDestroy(&valueIS);
559:   }
560:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
561:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
562:   DMGetPointSF(dm, &sfPoint);
563:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
564:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
565:   DMCreateLabel(dmNew, "vtk");
566:   DMCreateLabel(dmNew, "ghost");
567:   DMGetLabel(dmNew, "vtk", &vtkLabel);
568:   DMGetLabel(dmNew, "ghost", &ghostLabel);
569:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
570:     for (; c < leafLocal[l] && c < cEnd; ++c) {
571:       DMLabelSetValue(vtkLabel, c, 1);
572:     }
573:     if (leafLocal[l] >= cEnd) break;
574:     if (leafRemote[l].rank == rank) {
575:       DMLabelSetValue(vtkLabel, c, 1);
576:     } else {
577:       DMLabelSetValue(ghostLabel, c, 2);
578:     }
579:   }
580:   for (; c < cEnd; ++c) {
581:     DMLabelSetValue(vtkLabel, c, 1);
582:   }
583:   if (0) {
584:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
585:   }
586:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
587:   for (f = fStart; f < fEnd; ++f) {
588:     PetscInt numCells;

590:     DMPlexGetSupportSize(dmNew, f, &numCells);
591:     if (numCells < 2) {
592:       DMLabelSetValue(ghostLabel, f, 1);
593:     } else {
594:       const PetscInt *cells = NULL;
595:       PetscInt        vA, vB;

597:       DMPlexGetSupport(dmNew, f, &cells);
598:       DMLabelGetValue(vtkLabel, cells[0], &vA);
599:       DMLabelGetValue(vtkLabel, cells[1], &vB);
600:       if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
601:     }
602:   }
603:   if (0) {
604:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
605:   }
606:   return(0);
607: }

609: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
610: {
611:   DM             refTree;
612:   PetscSection   pSec;
613:   PetscInt       *parents, *childIDs;

617:   DMPlexGetReferenceTree(dm,&refTree);
618:   DMPlexSetReferenceTree(dmNew,refTree);
619:   DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
620:   if (pSec) {
621:     PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
622:     PetscInt *childIDsShifted;
623:     PetscSection pSecShifted;

625:     PetscSectionGetChart(pSec,&pStart,&pEnd);
626:     DMPlexGetDepth(dm,&depth);
627:     pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
628:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
629:     PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
630:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
631:     PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
632:     for (p = pStartShifted; p < pEndShifted; p++) {
633:       /* start off assuming no children */
634:       PetscSectionSetDof(pSecShifted,p,0);
635:     }
636:     for (p = pStart; p < pEnd; p++) {
637:       PetscInt dof;
638:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

640:       PetscSectionGetDof(pSec,p,&dof);
641:       PetscSectionSetDof(pSecShifted,pNew,dof);
642:     }
643:     PetscSectionSetUp(pSecShifted);
644:     for (p = pStart; p < pEnd; p++) {
645:       PetscInt dof;
646:       PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);

648:       PetscSectionGetDof(pSec,p,&dof);
649:       if (dof) {
650:         PetscInt off, offNew;

652:         PetscSectionGetOffset(pSec,p,&off);
653:         PetscSectionGetOffset(pSecShifted,pNew,&offNew);
654:         parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
655:         childIDsShifted[offNew] = childIDs[off];
656:       }
657:     }
658:     DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
659:     PetscFree2(parentsShifted,childIDsShifted);
660:     PetscSectionDestroy(&pSecShifted);
661:   }
662:   return(0);
663: }

665: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
666: {
667:   PetscSF               sf;
668:   IS                    valueIS;
669:   const PetscInt       *values, *leaves;
670:   PetscInt             *depthShift;
671:   PetscInt              d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
672:   PetscBool             isper;
673:   const PetscReal      *maxCell, *L;
674:   const DMBoundaryType *bd;
675:   PetscErrorCode        ierr;

678:   DMGetPointSF(dm, &sf);
679:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
680:   nleaves = PetscMax(0, nleaves);
681:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
682:   /* Count ghost cells */
683:   DMLabelGetValueIS(label, &valueIS);
684:   ISGetLocalSize(valueIS, &numFS);
685:   ISGetIndices(valueIS, &values);
686:   Ng   = 0;
687:   for (fs = 0; fs < numFS; ++fs) {
688:     IS              faceIS;
689:     const PetscInt *faces;
690:     PetscInt        numFaces, f, numBdFaces = 0;

692:     DMLabelGetStratumIS(label, values[fs], &faceIS);
693:     ISGetLocalSize(faceIS, &numFaces);
694:     ISGetIndices(faceIS, &faces);
695:     for (f = 0; f < numFaces; ++f) {
696:       PetscInt numChildren;

698:       PetscFindInt(faces[f], nleaves, leaves, &loc);
699:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
700:       /* non-local and ancestors points don't get to register ghosts */
701:       if (loc >= 0 || numChildren) continue;
702:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
703:     }
704:     Ng += numBdFaces;
705:     ISDestroy(&faceIS);
706:   }
707:   DMPlexGetDepth(dm, &depth);
708:   PetscMalloc1(2*(depth+1), &depthShift);
709:   for (d = 0; d <= depth; d++) {
710:     PetscInt dEnd;

712:     DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
713:     depthShift[2*d]   = dEnd;
714:     depthShift[2*d+1] = 0;
715:   }
716:   if (depth >= 0) depthShift[2*depth+1] = Ng;
717:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
718:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
719:   /* Step 3: Set cone/support sizes for new points */
720:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
721:   DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
722:   for (c = cEnd; c < cEnd + Ng; ++c) {
723:     DMPlexSetConeSize(gdm, c, 1);
724:   }
725:   for (fs = 0; fs < numFS; ++fs) {
726:     IS              faceIS;
727:     const PetscInt *faces;
728:     PetscInt        numFaces, f;

730:     DMLabelGetStratumIS(label, values[fs], &faceIS);
731:     ISGetLocalSize(faceIS, &numFaces);
732:     ISGetIndices(faceIS, &faces);
733:     for (f = 0; f < numFaces; ++f) {
734:       PetscInt size, numChildren;

736:       PetscFindInt(faces[f], nleaves, leaves, &loc);
737:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
738:       if (loc >= 0 || numChildren) continue;
739:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
740:       DMPlexGetSupportSize(dm, faces[f], &size);
741:       if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
742:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
743:     }
744:     ISRestoreIndices(faceIS, &faces);
745:     ISDestroy(&faceIS);
746:   }
747:   /* Step 4: Setup ghosted DM */
748:   DMSetUp(gdm);
749:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
750:   /* Step 6: Set cones and supports for new points */
751:   ghostCell = cEnd;
752:   for (fs = 0; fs < numFS; ++fs) {
753:     IS              faceIS;
754:     const PetscInt *faces;
755:     PetscInt        numFaces, f;

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

763:       PetscFindInt(faces[f], nleaves, leaves, &loc);
764:       DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
765:       if (loc >= 0 || numChildren) continue;
766:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
767:       DMPlexSetCone(gdm, ghostCell, &newFace);
768:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
769:       ++ghostCell;
770:     }
771:     ISRestoreIndices(faceIS, &faces);
772:     ISDestroy(&faceIS);
773:   }
774:   ISRestoreIndices(valueIS, &values);
775:   ISDestroy(&valueIS);
776:   /* Step 7: Stratify */
777:   DMPlexStratify(gdm);
778:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
779:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
780:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
781:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
782:   PetscFree(depthShift);
783:   /* Step 7: Periodicity */
784:   DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
785:   DMSetPeriodicity(gdm, isper, maxCell,  L,  bd);
786:   if (numGhostCells) *numGhostCells = Ng;
787:   return(0);
788: }

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

793:   Collective on dm

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

799:   Output Parameters:
800: + numGhostCells - The number of ghost cells added to the DM
801: - dmGhosted - The new DM

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

805:   Level: developer

807: .seealso: DMCreate()
808: @*/
809: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
810: {
811:   DM             gdm;
812:   DMLabel        label;
813:   const char    *name = labelName ? labelName : "Face Sets";
814:   PetscInt       dim;
815:   PetscBool      flag;

822:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
823:   DMSetType(gdm, DMPLEX);
824:   DMGetDimension(dm, &dim);
825:   DMSetDimension(gdm, dim);
826:   DMPlexGetAdjacencyUseCone(dm, &flag);
827:   DMPlexSetAdjacencyUseCone(gdm, flag);
828:   DMPlexGetAdjacencyUseClosure(dm, &flag);
829:   DMPlexSetAdjacencyUseClosure(gdm, flag);
830:   DMGetLabel(dm, name, &label);
831:   if (!label) {
832:     /* Get label for boundary faces */
833:     DMCreateLabel(dm, name);
834:     DMGetLabel(dm, name, &label);
835:     DMPlexMarkBoundaryFaces(dm, label);
836:   }
837:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
838:   DMCopyBoundary(dm, gdm);
839:   *dmGhosted = gdm;
840:   return(0);
841: }

843: /*
844:   We are adding three kinds of points here:
845:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
846:     Non-replicated: Points which exist on the fault, but are not replicated
847:     Hybrid:         Entirely new points, such as cohesive cells

849:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
850:   each depth so that the new split/hybrid points can be inserted as a block.
851: */
852: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
853: {
854:   MPI_Comm         comm;
855:   IS               valueIS;
856:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
857:   const PetscInt  *values;          /* List of depths for which we have replicated points */
858:   IS              *splitIS;
859:   IS              *unsplitIS;
860:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
861:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
862:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
863:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
864:   const PetscInt **splitPoints;        /* Replicated points for each depth */
865:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
866:   PetscSection     coordSection;
867:   Vec              coordinates;
868:   PetscScalar     *coords;
869:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
870:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
871:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
872:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
873:   PetscInt        *coneNew, *coneONew, *supportNew;
874:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
875:   PetscErrorCode   ierr;

878:   PetscObjectGetComm((PetscObject)dm,&comm);
879:   DMGetDimension(dm, &dim);
880:   DMPlexGetDepth(dm, &depth);
881:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
882:   /* Count split points and add cohesive cells */
883:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
884:   PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
885:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
886:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
887:   for (d = 0; d <= depth; ++d) {
888:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
889:     depthEnd[d]           = pMaxNew[d];
890:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
891:     numSplitPoints[d]     = 0;
892:     numUnsplitPoints[d]   = 0;
893:     numHybridPoints[d]    = 0;
894:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
895:     splitPoints[d]        = NULL;
896:     unsplitPoints[d]      = NULL;
897:     splitIS[d]            = NULL;
898:     unsplitIS[d]          = NULL;
899:     /* we are shifting the existing hybrid points with the stratum behind them, so
900:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
901:     depthShift[2*d]       = depthMax[d];
902:     depthShift[2*d+1]     = 0;
903:   }
904:   if (label) {
905:     DMLabelGetValueIS(label, &valueIS);
906:     ISGetLocalSize(valueIS, &numSP);
907:     ISGetIndices(valueIS, &values);
908:   }
909:   for (sp = 0; sp < numSP; ++sp) {
910:     const PetscInt dep = values[sp];

912:     if ((dep < 0) || (dep > depth)) continue;
913:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
914:     if (splitIS[dep]) {
915:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
916:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
917:     }
918:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
919:     if (unsplitIS[dep]) {
920:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
921:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
922:     }
923:   }
924:   /* Calculate number of hybrid points */
925:   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   */
926:   for (d = 0; d <= depth; ++d) depthShift[2*d+1]      = numSplitPoints[d] + numHybridPoints[d];
927:   DMPlexShiftPointSetUp_Internal(depth,depthShift);
928:   /* the end of the points in this stratum that come before the new points:
929:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
930:    * added points */
931:   for (d = 0; d <= depth; ++d) pMaxNew[d]             = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
932:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
933:   /* Step 3: Set cone/support sizes for new points */
934:   for (dep = 0; dep <= depth; ++dep) {
935:     for (p = 0; p < numSplitPoints[dep]; ++p) {
936:       const PetscInt  oldp   = splitPoints[dep][p];
937:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
938:       const PetscInt  splitp = p    + pMaxNew[dep];
939:       const PetscInt *support;
940:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

942:       DMPlexGetConeSize(dm, oldp, &coneSize);
943:       DMPlexSetConeSize(sdm, splitp, coneSize);
944:       DMPlexGetSupportSize(dm, oldp, &supportSize);
945:       DMPlexSetSupportSize(sdm, splitp, supportSize);
946:       if (dep == depth-1) {
947:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

954:         DMPlexGetSupport(dm, oldp, &support);
955:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
956:           PetscInt val;

958:           DMLabelGetValue(label, support[e], &val);
959:           if (val == 1) ++qf;
960:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
961:           if ((val == 1) || (val == -(shift + 1))) ++qp;
962:         }
963:         /* Split old vertex: Edges into original vertex and new cohesive edge */
964:         DMPlexSetSupportSize(sdm, newp, qn+1);
965:         /* Split new vertex: Edges into split vertex and new cohesive edge */
966:         DMPlexSetSupportSize(sdm, splitp, qp+1);
967:         /* Add hybrid edge */
968:         DMPlexSetConeSize(sdm, hybedge, 2);
969:         DMPlexSetSupportSize(sdm, hybedge, qf);
970:       } else if (dep == dim-2) {
971:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

977:           DMLabelGetValue(label, support[e], &val);
978:           if (val == dim-1) ++qf;
979:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
980:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
981:         }
982:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
983:         DMPlexSetSupportSize(sdm, newp, qn+1);
984:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
985:         DMPlexSetSupportSize(sdm, splitp, qp+1);
986:         /* Add hybrid face */
987:         DMPlexSetConeSize(sdm, hybface, 4);
988:         DMPlexSetSupportSize(sdm, hybface, qf);
989:       }
990:     }
991:   }
992:   for (dep = 0; dep <= depth; ++dep) {
993:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
994:       const PetscInt  oldp   = unsplitPoints[dep][p];
995:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
996:       const PetscInt *support;
997:       PetscInt        coneSize, supportSize, qf, e, s;

999:       DMPlexGetConeSize(dm, oldp, &coneSize);
1000:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1001:       DMPlexGetSupport(dm, oldp, &support);
1002:       if (dep == 0) {
1003:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1005:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1006:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1007:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1008:           if (e >= 0) ++qf;
1009:         }
1010:         DMPlexSetSupportSize(sdm, newp, qf+2);
1011:         /* Add hybrid edge */
1012:         DMPlexSetConeSize(sdm, hybedge, 2);
1013:         for (e = 0, qf = 0; e < supportSize; ++e) {
1014:           PetscInt val;

1016:           DMLabelGetValue(label, support[e], &val);
1017:           /* Split and unsplit edges produce hybrid faces */
1018:           if (val == 1) ++qf;
1019:           if (val == (shift2 + 1)) ++qf;
1020:         }
1021:         DMPlexSetSupportSize(sdm, hybedge, qf);
1022:       } else if (dep == dim-2) {
1023:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1024:         PetscInt       val;

1026:         for (e = 0, qf = 0; e < supportSize; ++e) {
1027:           DMLabelGetValue(label, support[e], &val);
1028:           if (val == dim-1) qf += 2;
1029:           else              ++qf;
1030:         }
1031:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1032:         DMPlexSetSupportSize(sdm, newp, qf+2);
1033:         /* Add hybrid face */
1034:         for (e = 0, qf = 0; e < supportSize; ++e) {
1035:           DMLabelGetValue(label, support[e], &val);
1036:           if (val == dim-1) ++qf;
1037:         }
1038:         DMPlexSetConeSize(sdm, hybface, 4);
1039:         DMPlexSetSupportSize(sdm, hybface, qf);
1040:       }
1041:     }
1042:   }
1043:   /* Step 4: Setup split DM */
1044:   DMSetUp(sdm);
1045:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1046:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1047:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1048:   /* Step 6: Set cones and supports for new points */
1049:   for (dep = 0; dep <= depth; ++dep) {
1050:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1051:       const PetscInt  oldp   = splitPoints[dep][p];
1052:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1053:       const PetscInt  splitp = p    + pMaxNew[dep];
1054:       const PetscInt *cone, *support, *ornt;
1055:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1057:       DMPlexGetConeSize(dm, oldp, &coneSize);
1058:       DMPlexGetCone(dm, oldp, &cone);
1059:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1060:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1061:       DMPlexGetSupport(dm, oldp, &support);
1062:       if (dep == depth-1) {
1063:         PetscBool       hasUnsplit = PETSC_FALSE;
1064:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1065:         const PetscInt *supportF;

1067:         /* Split face:       copy in old face to new face to start */
1068:         DMPlexGetSupport(sdm, newp,  &supportF);
1069:         DMPlexSetSupport(sdm, splitp, supportF);
1070:         /* Split old face:   old vertices/edges in cone so no change */
1071:         /* Split new face:   new vertices/edges in cone */
1072:         for (q = 0; q < coneSize; ++q) {
1073:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1074:           if (v < 0) {
1075:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1076:             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);
1077:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1078:             hasUnsplit   = PETSC_TRUE;
1079:           } else {
1080:             coneNew[2+q] = v + pMaxNew[dep-1];
1081:             if (dep > 1) {
1082:               const PetscInt *econe;
1083:               PetscInt        econeSize, r, vs, vu;

1085:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1086:               DMPlexGetCone(dm, cone[q], &econe);
1087:               for (r = 0; r < econeSize; ++r) {
1088:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
1089:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1090:                 if (vs >= 0) continue;
1091:                 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);
1092:                 hasUnsplit   = PETSC_TRUE;
1093:               }
1094:             }
1095:           }
1096:         }
1097:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1098:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1099:         /* Face support */
1100:         for (s = 0; s < supportSize; ++s) {
1101:           PetscInt val;

1103:           DMLabelGetValue(label, support[s], &val);
1104:           if (val < 0) {
1105:             /* Split old face:   Replace negative side cell with cohesive cell */
1106:              DMPlexInsertSupport(sdm, newp, s, hybcell);
1107:           } else {
1108:             /* Split new face:   Replace positive side cell with cohesive cell */
1109:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
1110:             /* Get orientation for cohesive face */
1111:             {
1112:               const PetscInt *ncone, *nconeO;
1113:               PetscInt        nconeSize, nc;

1115:               DMPlexGetConeSize(dm, support[s], &nconeSize);
1116:               DMPlexGetCone(dm, support[s], &ncone);
1117:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
1118:               for (nc = 0; nc < nconeSize; ++nc) {
1119:                 if (ncone[nc] == oldp) {
1120:                   coneONew[0] = nconeO[nc];
1121:                   break;
1122:                 }
1123:               }
1124:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1125:             }
1126:           }
1127:         }
1128:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1129:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
1130:         coneNew[1]  = splitp;
1131:         coneONew[1] = coneONew[0];
1132:         for (q = 0; q < coneSize; ++q) {
1133:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1134:           if (v < 0) {
1135:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1136:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1137:             coneONew[2+q] = 0;
1138:           } else {
1139:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
1140:           }
1141:           coneONew[2+q] = 0;
1142:         }
1143:         DMPlexSetCone(sdm, hybcell, coneNew);
1144:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1145:         /* Label the hybrid cells on the boundary of the split */
1146:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1147:       } else if (dep == 0) {
1148:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

1154:           DMLabelGetValue(label, support[e], &val);
1155:           if ((val == 1) || (val == (shift + 1))) {
1156:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1157:           }
1158:         }
1159:         supportNew[qn] = hybedge;
1160:         DMPlexSetSupport(sdm, newp, supportNew);
1161:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1162:         for (e = 0, qp = 0; e < supportSize; ++e) {
1163:           PetscInt val, edge;

1165:           DMLabelGetValue(label, support[e], &val);
1166:           if (val == 1) {
1167:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1168:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1169:             supportNew[qp++] = edge + pMaxNew[dep+1];
1170:           } else if (val == -(shift + 1)) {
1171:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1172:           }
1173:         }
1174:         supportNew[qp] = hybedge;
1175:         DMPlexSetSupport(sdm, splitp, supportNew);
1176:         /* Hybrid edge:    Old and new split vertex */
1177:         coneNew[0] = newp;
1178:         coneNew[1] = splitp;
1179:         DMPlexSetCone(sdm, hybedge, coneNew);
1180:         for (e = 0, qf = 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[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1188:           }
1189:         }
1190:         DMPlexSetSupport(sdm, hybedge, supportNew);
1191:       } else if (dep == dim-2) {
1192:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

1194:         /* Split old edge:   old vertices in cone so no change */
1195:         /* Split new edge:   new vertices in cone */
1196:         for (q = 0; q < coneSize; ++q) {
1197:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1198:           if (v < 0) {
1199:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1200:             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);
1201:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1202:           } else {
1203:             coneNew[q] = v + pMaxNew[dep-1];
1204:           }
1205:         }
1206:         DMPlexSetCone(sdm, splitp, coneNew);
1207:         /* Split old edge: Faces in positive side cells and old split faces */
1208:         for (e = 0, q = 0; e < supportSize; ++e) {
1209:           PetscInt val;

1211:           DMLabelGetValue(label, support[e], &val);
1212:           if (val == dim-1) {
1213:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1214:           } else if (val == (shift + dim-1)) {
1215:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1216:           }
1217:         }
1218:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1219:         DMPlexSetSupport(sdm, newp, supportNew);
1220:         /* Split new edge: Faces in negative side cells and new split faces */
1221:         for (e = 0, q = 0; e < supportSize; ++e) {
1222:           PetscInt val, face;

1224:           DMLabelGetValue(label, support[e], &val);
1225:           if (val == dim-1) {
1226:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1227:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1228:             supportNew[q++] = face + pMaxNew[dep+1];
1229:           } else if (val == -(shift + dim-1)) {
1230:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1231:           }
1232:         }
1233:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1234:         DMPlexSetSupport(sdm, splitp, supportNew);
1235:         /* Hybrid face */
1236:         coneNew[0] = newp;
1237:         coneNew[1] = splitp;
1238:         for (v = 0; v < coneSize; ++v) {
1239:           PetscInt vertex;
1240:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1241:           if (vertex < 0) {
1242:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1243:             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);
1244:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1245:           } else {
1246:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1247:           }
1248:         }
1249:         DMPlexSetCone(sdm, hybface, coneNew);
1250:         for (e = 0, qf = 0; e < supportSize; ++e) {
1251:           PetscInt val, face;

1253:           DMLabelGetValue(label, support[e], &val);
1254:           if (val == dim-1) {
1255:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1256:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1257:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1258:           }
1259:         }
1260:         DMPlexSetSupport(sdm, hybface, supportNew);
1261:       }
1262:     }
1263:   }
1264:   for (dep = 0; dep <= depth; ++dep) {
1265:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1266:       const PetscInt  oldp   = unsplitPoints[dep][p];
1267:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1268:       const PetscInt *cone, *support, *ornt;
1269:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1271:       DMPlexGetConeSize(dm, oldp, &coneSize);
1272:       DMPlexGetCone(dm, oldp, &cone);
1273:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1274:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1275:       DMPlexGetSupport(dm, oldp, &support);
1276:       if (dep == 0) {
1277:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1279:         /* Unsplit vertex */
1280:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1281:         for (s = 0, q = 0; s < supportSize; ++s) {
1282:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1283:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1284:           if (e >= 0) {
1285:             supportNew[q++] = e + pMaxNew[dep+1];
1286:           }
1287:         }
1288:         supportNew[q++] = hybedge;
1289:         supportNew[q++] = hybedge;
1290:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1291:         DMPlexSetSupport(sdm, newp, supportNew);
1292:         /* Hybrid edge */
1293:         coneNew[0] = newp;
1294:         coneNew[1] = newp;
1295:         DMPlexSetCone(sdm, hybedge, coneNew);
1296:         for (e = 0, qf = 0; e < supportSize; ++e) {
1297:           PetscInt val, edge;

1299:           DMLabelGetValue(label, support[e], &val);
1300:           if (val == 1) {
1301:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1302:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1303:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1304:           } else if  (val ==  (shift2 + 1)) {
1305:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1306:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1307:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1308:           }
1309:         }
1310:         DMPlexSetSupport(sdm, hybedge, supportNew);
1311:       } else if (dep == dim-2) {
1312:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1318:           DMLabelGetValue(label, support[f], &val);
1319:           if (val == dim-1) {
1320:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1321:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1322:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1323:             supportNew[qf++] = face + pMaxNew[dep+1];
1324:           } else {
1325:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1326:           }
1327:         }
1328:         supportNew[qf++] = hybface;
1329:         supportNew[qf++] = hybface;
1330:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1331:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1332:         DMPlexSetSupport(sdm, newp, supportNew);
1333:         /* Add hybrid face */
1334:         coneNew[0] = newp;
1335:         coneNew[1] = newp;
1336:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1337:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1338:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1339:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1340:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1341:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1342:         DMPlexSetCone(sdm, hybface, coneNew);
1343:         for (f = 0, qf = 0; f < supportSize; ++f) {
1344:           PetscInt val, face;

1346:           DMLabelGetValue(label, support[f], &val);
1347:           if (val == dim-1) {
1348:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1349:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1350:           }
1351:         }
1352:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1353:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1354:         DMPlexSetSupport(sdm, hybface, supportNew);
1355:       }
1356:     }
1357:   }
1358:   /* Step 6b: Replace split points in negative side cones */
1359:   for (sp = 0; sp < numSP; ++sp) {
1360:     PetscInt        dep = values[sp];
1361:     IS              pIS;
1362:     PetscInt        numPoints;
1363:     const PetscInt *points;

1365:     if (dep >= 0) continue;
1366:     DMLabelGetStratumIS(label, dep, &pIS);
1367:     if (!pIS) continue;
1368:     dep  = -dep - shift;
1369:     ISGetLocalSize(pIS, &numPoints);
1370:     ISGetIndices(pIS, &points);
1371:     for (p = 0; p < numPoints; ++p) {
1372:       const PetscInt  oldp = points[p];
1373:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1374:       const PetscInt *cone;
1375:       PetscInt        coneSize, c;
1376:       /* PetscBool       replaced = PETSC_FALSE; */

1378:       /* Negative edge: replace split vertex */
1379:       /* Negative cell: replace split face */
1380:       DMPlexGetConeSize(sdm, newp, &coneSize);
1381:       DMPlexGetCone(sdm, newp, &cone);
1382:       for (c = 0; c < coneSize; ++c) {
1383:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1384:         PetscInt       csplitp, cp, val;

1386:         DMLabelGetValue(label, coldp, &val);
1387:         if (val == dep-1) {
1388:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1389:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1390:           csplitp  = pMaxNew[dep-1] + cp;
1391:           DMPlexInsertCone(sdm, newp, c, csplitp);
1392:           /* replaced = PETSC_TRUE; */
1393:         }
1394:       }
1395:       /* Cells with only a vertex or edge on the submesh have no replacement */
1396:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1397:     }
1398:     ISRestoreIndices(pIS, &points);
1399:     ISDestroy(&pIS);
1400:   }
1401:   /* Step 7: Stratify */
1402:   DMPlexStratify(sdm);
1403:   /* Step 8: Coordinates */
1404:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1405:   DMGetCoordinateSection(sdm, &coordSection);
1406:   DMGetCoordinatesLocal(sdm, &coordinates);
1407:   VecGetArray(coordinates, &coords);
1408:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1409:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1410:     const PetscInt splitp = pMaxNew[0] + v;
1411:     PetscInt       dof, off, soff, d;

1413:     PetscSectionGetDof(coordSection, newp, &dof);
1414:     PetscSectionGetOffset(coordSection, newp, &off);
1415:     PetscSectionGetOffset(coordSection, splitp, &soff);
1416:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1417:   }
1418:   VecRestoreArray(coordinates, &coords);
1419:   /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1420:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1421:   /* Step 10: Labels */
1422:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1423:   DMGetNumLabels(sdm, &numLabels);
1424:   for (dep = 0; dep <= depth; ++dep) {
1425:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1426:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1427:       const PetscInt splitp = pMaxNew[dep] + p;
1428:       PetscInt       l;

1430:       for (l = 0; l < numLabels; ++l) {
1431:         DMLabel     mlabel;
1432:         const char *lname;
1433:         PetscInt    val;
1434:         PetscBool   isDepth;

1436:         DMGetLabelName(sdm, l, &lname);
1437:         PetscStrcmp(lname, "depth", &isDepth);
1438:         if (isDepth) continue;
1439:         DMGetLabel(sdm, lname, &mlabel);
1440:         DMLabelGetValue(mlabel, newp, &val);
1441:         if (val >= 0) {
1442:           DMLabelSetValue(mlabel, splitp, val);
1443:         }
1444:       }
1445:     }
1446:   }
1447:   for (sp = 0; sp < numSP; ++sp) {
1448:     const PetscInt dep = values[sp];

1450:     if ((dep < 0) || (dep > depth)) continue;
1451:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1452:     ISDestroy(&splitIS[dep]);
1453:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1454:     ISDestroy(&unsplitIS[dep]);
1455:   }
1456:   if (label) {
1457:     ISRestoreIndices(valueIS, &values);
1458:     ISDestroy(&valueIS);
1459:   }
1460:   for (d = 0; d <= depth; ++d) {
1461:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1462:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1463:   }
1464:   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);
1465:   PetscFree3(coneNew, coneONew, supportNew);
1466:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1467:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1468:   return(0);
1469: }

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

1474:   Collective on dm

1476:   Input Parameters:
1477: + dm - The original DM
1478: - label - The label specifying the boundary faces (this could be auto-generated)

1480:   Output Parameters:
1481: - dmSplit - The new DM

1483:   Level: developer

1485: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1486: @*/
1487: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1488: {
1489:   DM             sdm;
1490:   PetscInt       dim;

1496:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1497:   DMSetType(sdm, DMPLEX);
1498:   DMGetDimension(dm, &dim);
1499:   DMSetDimension(sdm, dim);
1500:   switch (dim) {
1501:   case 2:
1502:   case 3:
1503:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1504:     break;
1505:   default:
1506:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1507:   }
1508:   *dmSplit = sdm;
1509:   return(0);
1510: }

1512: /* Returns the side of the surface for a given cell with a face on the surface */
1513: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1514: {
1515:   const PetscInt *cone, *ornt;
1516:   PetscInt        dim, coneSize, c;
1517:   PetscErrorCode  ierr;

1520:   *pos = PETSC_TRUE;
1521:   DMGetDimension(dm, &dim);
1522:   DMPlexGetConeSize(dm, cell, &coneSize);
1523:   DMPlexGetCone(dm, cell, &cone);
1524:   DMPlexGetConeOrientation(dm, cell, &ornt);
1525:   for (c = 0; c < coneSize; ++c) {
1526:     if (cone[c] == face) {
1527:       PetscInt o = ornt[c];

1529:       if (subdm) {
1530:         const PetscInt *subcone, *subornt;
1531:         PetscInt        subpoint, subface, subconeSize, sc;

1533:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1534:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1535:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1536:         DMPlexGetCone(subdm, subpoint, &subcone);
1537:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1538:         for (sc = 0; sc < subconeSize; ++sc) {
1539:           if (subcone[sc] == subface) {
1540:             o = subornt[0];
1541:             break;
1542:           }
1543:         }
1544:         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);
1545:       }
1546:       if (o >= 0) *pos = PETSC_TRUE;
1547:       else        *pos = PETSC_FALSE;
1548:       break;
1549:     }
1550:   }
1551:   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);
1552:   return(0);
1553: }

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

1559:   Input Parameters:
1560: + dm     - The DM
1561: . label  - A DMLabel marking the surface
1562: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1563: . flip   - Flag to flip the submesh normal and replace points on the other side
1564: - subdm  - The subDM associated with the label, or NULL

1566:   Output Parameter:
1567: . label - A DMLabel marking all surface points

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

1571:   Level: developer

1573: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1574: @*/
1575: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1576: {
1577:   DMLabel         depthLabel;
1578:   IS              dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1579:   const PetscInt *points, *subpoints;
1580:   const PetscInt  rev   = flip ? -1 : 1;
1581:   PetscInt       *pMax;
1582:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1583:   PetscErrorCode  ierr;

1586:   DMPlexGetDepth(dm, &depth);
1587:   DMGetDimension(dm, &dim);
1588:   pSize = PetscMax(depth, dim) + 1;
1589:   PetscMalloc1(pSize,&pMax);
1590:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1591:   DMPlexGetDepthLabel(dm, &depthLabel);
1592:   DMGetDimension(dm, &dim);
1593:   if (subdm) {
1594:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1595:     if (subpointIS) {
1596:       ISGetLocalSize(subpointIS, &numSubpoints);
1597:       ISGetIndices(subpointIS, &subpoints);
1598:     }
1599:   }
1600:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1601:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1602:   if (!dimIS) {
1603:     PetscFree(pMax);
1604:     ISDestroy(&subpointIS);
1605:     return(0);
1606:   }
1607:   ISGetLocalSize(dimIS, &numPoints);
1608:   ISGetIndices(dimIS, &points);
1609:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1610:     const PetscInt *support;
1611:     PetscInt        supportSize, s;

1613:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1614:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1615:     DMPlexGetSupport(dm, points[p], &support);
1616:     for (s = 0; s < supportSize; ++s) {
1617:       const PetscInt *cone;
1618:       PetscInt        coneSize, c;
1619:       PetscBool       pos;

1621:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1622:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1623:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1624:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1625:       /* Put faces touching the fault in the label */
1626:       DMPlexGetConeSize(dm, support[s], &coneSize);
1627:       DMPlexGetCone(dm, support[s], &cone);
1628:       for (c = 0; c < coneSize; ++c) {
1629:         const PetscInt point = cone[c];

1631:         DMLabelGetValue(label, point, &val);
1632:         if (val == -1) {
1633:           PetscInt *closure = NULL;
1634:           PetscInt  closureSize, cl;

1636:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1637:           for (cl = 0; cl < closureSize*2; cl += 2) {
1638:             const PetscInt clp  = closure[cl];
1639:             PetscInt       bval = -1;

1641:             DMLabelGetValue(label, clp, &val);
1642:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1643:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1644:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1645:               break;
1646:             }
1647:           }
1648:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1649:         }
1650:       }
1651:     }
1652:   }
1653:   ISRestoreIndices(dimIS, &points);
1654:   ISDestroy(&dimIS);
1655:   if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1656:   ISDestroy(&subpointIS);
1657:   /* Mark boundary points as unsplit */
1658:   if (blabel) {
1659:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1660:     ISGetLocalSize(dimIS, &numPoints);
1661:     ISGetIndices(dimIS, &points);
1662:     for (p = 0; p < numPoints; ++p) {
1663:       const PetscInt point = points[p];
1664:       PetscInt       val, bval;

1666:       DMLabelGetValue(blabel, point, &bval);
1667:       if (bval >= 0) {
1668:         DMLabelGetValue(label, point, &val);
1669:         if ((val < 0) || (val > dim)) {
1670:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1671:           DMLabelClearValue(blabel, point, bval);
1672:         }
1673:       }
1674:     }
1675:     for (p = 0; p < numPoints; ++p) {
1676:       const PetscInt point = points[p];
1677:       PetscInt       val, bval;

1679:       DMLabelGetValue(blabel, point, &bval);
1680:       if (bval >= 0) {
1681:         const PetscInt *cone,    *support;
1682:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1684:         /* Mark as unsplit */
1685:         DMLabelGetValue(label, point, &val);
1686:         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);
1687:         DMLabelClearValue(label, point, val);
1688:         DMLabelSetValue(label, point, shift2+val);
1689:         /* Check for cross-edge
1690:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1691:         if (val != 0) continue;
1692:         DMPlexGetSupport(dm, point, &support);
1693:         DMPlexGetSupportSize(dm, point, &supportSize);
1694:         for (s = 0; s < supportSize; ++s) {
1695:           DMPlexGetCone(dm, support[s], &cone);
1696:           DMPlexGetConeSize(dm, support[s], &coneSize);
1697:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1698:           DMLabelGetValue(blabel, cone[0], &valA);
1699:           DMLabelGetValue(blabel, cone[1], &valB);
1700:           DMLabelGetValue(blabel, support[s], &valE);
1701:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1702:         }
1703:       }
1704:     }
1705:     ISRestoreIndices(dimIS, &points);
1706:     ISDestroy(&dimIS);
1707:   }
1708:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1709:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1710:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1711:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1712:   cMax = cMax < 0 ? cEnd : cMax;
1713:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1714:   DMLabelGetStratumIS(label, 0, &dimIS);
1715:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1716:   if (dimIS && crossEdgeIS) {
1717:     IS vertIS = dimIS;

1719:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1720:     ISDestroy(&crossEdgeIS);
1721:     ISDestroy(&vertIS);
1722:   }
1723:   if (!dimIS) {
1724:     PetscFree(pMax);
1725:     return(0);
1726:   }
1727:   ISGetLocalSize(dimIS, &numPoints);
1728:   ISGetIndices(dimIS, &points);
1729:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1730:     PetscInt *star = NULL;
1731:     PetscInt  starSize, s;
1732:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1734:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1735:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1736:     while (again) {
1737:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1738:       again = 0;
1739:       for (s = 0; s < starSize*2; s += 2) {
1740:         const PetscInt  point = star[s];
1741:         const PetscInt *cone;
1742:         PetscInt        coneSize, c;

1744:         if ((point < cStart) || (point >= cMax)) continue;
1745:         DMLabelGetValue(label, point, &val);
1746:         if (val != -1) continue;
1747:         again = again == 1 ? 1 : 2;
1748:         DMPlexGetConeSize(dm, point, &coneSize);
1749:         DMPlexGetCone(dm, point, &cone);
1750:         for (c = 0; c < coneSize; ++c) {
1751:           DMLabelGetValue(label, cone[c], &val);
1752:           if (val != -1) {
1753:             const PetscInt *ccone;
1754:             PetscInt        cconeSize, cc, side;

1756:             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);
1757:             if (val > 0) side =  1;
1758:             else         side = -1;
1759:             DMLabelSetValue(label, point, side*(shift+dim));
1760:             /* Mark cell faces which touch the fault */
1761:             DMPlexGetConeSize(dm, point, &cconeSize);
1762:             DMPlexGetCone(dm, point, &ccone);
1763:             for (cc = 0; cc < cconeSize; ++cc) {
1764:               PetscInt *closure = NULL;
1765:               PetscInt  closureSize, cl;

1767:               DMLabelGetValue(label, ccone[cc], &val);
1768:               if (val != -1) continue;
1769:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1770:               for (cl = 0; cl < closureSize*2; cl += 2) {
1771:                 const PetscInt clp = closure[cl];

1773:                 DMLabelGetValue(label, clp, &val);
1774:                 if (val == -1) continue;
1775:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1776:                 break;
1777:               }
1778:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1779:             }
1780:             again = 1;
1781:             break;
1782:           }
1783:         }
1784:       }
1785:     }
1786:     /* Classify the rest by cell membership */
1787:     for (s = 0; s < starSize*2; s += 2) {
1788:       const PetscInt point = star[s];

1790:       DMLabelGetValue(label, point, &val);
1791:       if (val == -1) {
1792:         PetscInt  *sstar = NULL;
1793:         PetscInt   sstarSize, ss;
1794:         PetscBool  marked = PETSC_FALSE;

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

1800:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1801:           DMLabelGetValue(label, spoint, &val);
1802:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1803:           DMLabelGetValue(depthLabel, point, &dep);
1804:           if (val > 0) {
1805:             DMLabelSetValue(label, point,   shift+dep);
1806:           } else {
1807:             DMLabelSetValue(label, point, -(shift+dep));
1808:           }
1809:           marked = PETSC_TRUE;
1810:           break;
1811:         }
1812:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1813:         DMLabelGetValue(depthLabel, point, &dep);
1814:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1815:       }
1816:     }
1817:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1818:   }
1819:   ISRestoreIndices(dimIS, &points);
1820:   ISDestroy(&dimIS);
1821:   /* If any faces touching the fault divide cells on either side, split them */
1822:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1823:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1824:   ISExpand(facePosIS, faceNegIS, &dimIS);
1825:   ISDestroy(&facePosIS);
1826:   ISDestroy(&faceNegIS);
1827:   ISGetLocalSize(dimIS, &numPoints);
1828:   ISGetIndices(dimIS, &points);
1829:   for (p = 0; p < numPoints; ++p) {
1830:     const PetscInt  point = points[p];
1831:     const PetscInt *support;
1832:     PetscInt        supportSize, valA, valB;

1834:     DMPlexGetSupportSize(dm, point, &supportSize);
1835:     if (supportSize != 2) continue;
1836:     DMPlexGetSupport(dm, point, &support);
1837:     DMLabelGetValue(label, support[0], &valA);
1838:     DMLabelGetValue(label, support[1], &valB);
1839:     if ((valA == -1) || (valB == -1)) continue;
1840:     if (valA*valB > 0) continue;
1841:     /* Split the face */
1842:     DMLabelGetValue(label, point, &valA);
1843:     DMLabelClearValue(label, point, valA);
1844:     DMLabelSetValue(label, point, dim-1);
1845:     /* Label its closure:
1846:       unmarked: label as unsplit
1847:       incident: relabel as split
1848:       split:    do nothing
1849:     */
1850:     {
1851:       PetscInt *closure = NULL;
1852:       PetscInt  closureSize, cl;

1854:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1855:       for (cl = 0; cl < closureSize*2; cl += 2) {
1856:         DMLabelGetValue(label, closure[cl], &valA);
1857:         if (valA == -1) { /* Mark as unsplit */
1858:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1859:           DMLabelSetValue(label, closure[cl], shift2+dep);
1860:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1861:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1862:           DMLabelClearValue(label, closure[cl], valA);
1863:           DMLabelSetValue(label, closure[cl], dep);
1864:         }
1865:       }
1866:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1867:     }
1868:   }
1869:   ISRestoreIndices(dimIS, &points);
1870:   ISDestroy(&dimIS);
1871:   PetscFree(pMax);
1872:   return(0);
1873: }

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

1878:   Collective on dm

1880:   Input Parameters:
1881: + dm - The original DM
1882: - labelName - The label specifying the interface vertices

1884:   Output Parameters:
1885: + hybridLabel - The label fully marking the interface
1886: - dmHybrid - The new DM

1888:   Level: developer

1890: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1891: @*/
1892: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1893: {
1894:   DM             idm;
1895:   DMLabel        subpointMap, hlabel;
1896:   PetscInt       dim;

1903:   DMGetDimension(dm, &dim);
1904:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1905:   DMPlexOrient(idm);
1906:   DMPlexGetSubpointMap(idm, &subpointMap);
1907:   DMLabelDuplicate(subpointMap, &hlabel);
1908:   DMLabelClearStratum(hlabel, dim);
1909:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1910:   DMDestroy(&idm);
1911:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1912:   if (hybridLabel) *hybridLabel = hlabel;
1913:   else             {DMLabelDestroy(&hlabel);}
1914:   return(0);
1915: }

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

1919:      For any marked cell, the marked vertices constitute a single face
1920: */
1921: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1922: {
1923:   IS               subvertexIS = NULL;
1924:   const PetscInt  *subvertices;
1925:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1926:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1927:   PetscErrorCode   ierr;

1930:   *numFaces = 0;
1931:   *nFV      = 0;
1932:   DMPlexGetDepth(dm, &depth);
1933:   DMGetDimension(dm, &dim);
1934:   pSize = PetscMax(depth, dim) + 1;
1935:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1936:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1937:   for (d = 0; d <= depth; ++d) {
1938:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1939:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1940:   }
1941:   /* Loop over initial vertices and mark all faces in the collective star() */
1942:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1943:   if (subvertexIS) {
1944:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1945:     ISGetIndices(subvertexIS, &subvertices);
1946:   }
1947:   for (v = 0; v < numSubVerticesInitial; ++v) {
1948:     const PetscInt vertex = subvertices[v];
1949:     PetscInt      *star   = NULL;
1950:     PetscInt       starSize, s, numCells = 0, c;

1952:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1953:     for (s = 0; s < starSize*2; s += 2) {
1954:       const PetscInt point = star[s];
1955:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1956:     }
1957:     for (c = 0; c < numCells; ++c) {
1958:       const PetscInt cell    = star[c];
1959:       PetscInt      *closure = NULL;
1960:       PetscInt       closureSize, cl;
1961:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1963:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1964:       if (cellLoc == 2) continue;
1965:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1966:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1967:       for (cl = 0; cl < closureSize*2; cl += 2) {
1968:         const PetscInt point = closure[cl];
1969:         PetscInt       vertexLoc;

1971:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1972:           ++numCorners;
1973:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1974:           if (vertexLoc == value) closure[faceSize++] = point;
1975:         }
1976:       }
1977:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1978:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1979:       if (faceSize == *nFV) {
1980:         const PetscInt *cells = NULL;
1981:         PetscInt        numCells, nc;

1983:         ++(*numFaces);
1984:         for (cl = 0; cl < faceSize; ++cl) {
1985:           DMLabelSetValue(subpointMap, closure[cl], 0);
1986:         }
1987:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1988:         for (nc = 0; nc < numCells; ++nc) {
1989:           DMLabelSetValue(subpointMap, cells[nc], 2);
1990:         }
1991:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1992:       }
1993:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1994:     }
1995:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1996:   }
1997:   if (subvertexIS) {
1998:     ISRestoreIndices(subvertexIS, &subvertices);
1999:   }
2000:   ISDestroy(&subvertexIS);
2001:   PetscFree3(pStart,pEnd,pMax);
2002:   return(0);
2003: }

2005: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2006: {
2007:   IS               subvertexIS = NULL;
2008:   const PetscInt  *subvertices;
2009:   PetscInt        *pStart, *pEnd, *pMax;
2010:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
2011:   PetscErrorCode   ierr;

2014:   DMGetDimension(dm, &dim);
2015:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2016:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2017:   for (d = 0; d <= dim; ++d) {
2018:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2019:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2020:   }
2021:   /* Loop over initial vertices and mark all faces in the collective star() */
2022:   if (vertexLabel) {
2023:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2024:     if (subvertexIS) {
2025:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2026:       ISGetIndices(subvertexIS, &subvertices);
2027:     }
2028:   }
2029:   for (v = 0; v < numSubVerticesInitial; ++v) {
2030:     const PetscInt vertex = subvertices[v];
2031:     PetscInt      *star   = NULL;
2032:     PetscInt       starSize, s, numFaces = 0, f;

2034:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2035:     for (s = 0; s < starSize*2; s += 2) {
2036:       const PetscInt point = star[s];
2037:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2038:     }
2039:     for (f = 0; f < numFaces; ++f) {
2040:       const PetscInt face    = star[f];
2041:       PetscInt      *closure = NULL;
2042:       PetscInt       closureSize, c;
2043:       PetscInt       faceLoc;

2045:       DMLabelGetValue(subpointMap, face, &faceLoc);
2046:       if (faceLoc == dim-1) continue;
2047:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2048:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2049:       for (c = 0; c < closureSize*2; c += 2) {
2050:         const PetscInt point = closure[c];
2051:         PetscInt       vertexLoc;

2053:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2054:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2055:           if (vertexLoc != value) break;
2056:         }
2057:       }
2058:       if (c == closureSize*2) {
2059:         const PetscInt *support;
2060:         PetscInt        supportSize, s;

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

2065:           for (d = 0; d < dim; ++d) {
2066:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2067:               DMLabelSetValue(subpointMap, point, d);
2068:               break;
2069:             }
2070:           }
2071:         }
2072:         DMPlexGetSupportSize(dm, face, &supportSize);
2073:         DMPlexGetSupport(dm, face, &support);
2074:         for (s = 0; s < supportSize; ++s) {
2075:           DMLabelSetValue(subpointMap, support[s], dim);
2076:         }
2077:       }
2078:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2079:     }
2080:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2081:   }
2082:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2083:   ISDestroy(&subvertexIS);
2084:   PetscFree3(pStart,pEnd,pMax);
2085:   return(0);
2086: }

2088: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2089: {
2090:   DMLabel         label = NULL;
2091:   const PetscInt *cone;
2092:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2093:   PetscErrorCode  ierr;

2096:   *numFaces = 0;
2097:   *nFV = 0;
2098:   if (labelname) {DMGetLabel(dm, labelname, &label);}
2099:   *subCells = NULL;
2100:   DMGetDimension(dm, &dim);
2101:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2102:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2103:   if (cMax < 0) return(0);
2104:   if (label) {
2105:     for (c = cMax; c < cEnd; ++c) {
2106:       PetscInt val;

2108:       DMLabelGetValue(label, c, &val);
2109:       if (val == value) {
2110:         ++(*numFaces);
2111:         DMPlexGetConeSize(dm, c, &coneSize);
2112:       }
2113:     }
2114:   } else {
2115:     *numFaces = cEnd - cMax;
2116:     DMPlexGetConeSize(dm, cMax, &coneSize);
2117:   }
2118:   PetscMalloc1(*numFaces *2, subCells);
2119:   if (!(*numFaces)) return(0);
2120:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2121:   for (c = cMax; c < cEnd; ++c) {
2122:     const PetscInt *cells;
2123:     PetscInt        numCells;

2125:     if (label) {
2126:       PetscInt val;

2128:       DMLabelGetValue(label, c, &val);
2129:       if (val != value) continue;
2130:     }
2131:     DMPlexGetCone(dm, c, &cone);
2132:     for (p = 0; p < *nFV; ++p) {
2133:       DMLabelSetValue(subpointMap, cone[p], 0);
2134:     }
2135:     /* Negative face */
2136:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2137:     /* Not true in parallel
2138:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2139:     for (p = 0; p < numCells; ++p) {
2140:       DMLabelSetValue(subpointMap, cells[p], 2);
2141:       (*subCells)[subc++] = cells[p];
2142:     }
2143:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2144:     /* Positive face is not included */
2145:   }
2146:   return(0);
2147: }

2149: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2150: {
2151:   PetscInt      *pStart, *pEnd;
2152:   PetscInt       dim, cMax, cEnd, c, d;

2156:   DMGetDimension(dm, &dim);
2157:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2158:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2159:   if (cMax < 0) return(0);
2160:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2161:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2162:   for (c = cMax; c < cEnd; ++c) {
2163:     const PetscInt *cone;
2164:     PetscInt       *closure = NULL;
2165:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2167:     if (label) {
2168:       DMLabelGetValue(label, c, &val);
2169:       if (val != value) continue;
2170:     }
2171:     DMPlexGetConeSize(dm, c, &coneSize);
2172:     DMPlexGetCone(dm, c, &cone);
2173:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2174:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2175:     /* Negative face */
2176:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2177:     for (cl = 0; cl < closureSize*2; cl += 2) {
2178:       const PetscInt point = closure[cl];

2180:       for (d = 0; d <= dim; ++d) {
2181:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2182:           DMLabelSetValue(subpointMap, point, d);
2183:           break;
2184:         }
2185:       }
2186:     }
2187:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2188:     /* Cells -- positive face is not included */
2189:     for (cl = 0; cl < 1; ++cl) {
2190:       const PetscInt *support;
2191:       PetscInt        supportSize, s;

2193:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2194:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2195:       DMPlexGetSupport(dm, cone[cl], &support);
2196:       for (s = 0; s < supportSize; ++s) {
2197:         DMLabelSetValue(subpointMap, support[s], dim);
2198:       }
2199:     }
2200:   }
2201:   PetscFree2(pStart, pEnd);
2202:   return(0);
2203: }

2205: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2206: {
2207:   MPI_Comm       comm;
2208:   PetscBool      posOrient = PETSC_FALSE;
2209:   const PetscInt debug     = 0;
2210:   PetscInt       cellDim, faceSize, f;

2214:   PetscObjectGetComm((PetscObject)dm,&comm);
2215:   DMGetDimension(dm, &cellDim);
2216:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2218:   if (cellDim == 1 && numCorners == 2) {
2219:     /* Triangle */
2220:     faceSize  = numCorners-1;
2221:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2222:   } else if (cellDim == 2 && numCorners == 3) {
2223:     /* Triangle */
2224:     faceSize  = numCorners-1;
2225:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2226:   } else if (cellDim == 3 && numCorners == 4) {
2227:     /* Tetrahedron */
2228:     faceSize  = numCorners-1;
2229:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2230:   } else if (cellDim == 1 && numCorners == 3) {
2231:     /* Quadratic line */
2232:     faceSize  = 1;
2233:     posOrient = PETSC_TRUE;
2234:   } else if (cellDim == 2 && numCorners == 4) {
2235:     /* Quads */
2236:     faceSize = 2;
2237:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2238:       posOrient = PETSC_TRUE;
2239:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2240:       posOrient = PETSC_TRUE;
2241:     } else {
2242:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2243:         posOrient = PETSC_FALSE;
2244:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2245:     }
2246:   } else if (cellDim == 2 && numCorners == 6) {
2247:     /* Quadratic triangle (I hate this) */
2248:     /* Edges are determined by the first 2 vertices (corners of edges) */
2249:     const PetscInt faceSizeTri = 3;
2250:     PetscInt       sortedIndices[3], i, iFace;
2251:     PetscBool      found                    = PETSC_FALSE;
2252:     PetscInt       faceVerticesTriSorted[9] = {
2253:       0, 3,  4, /* bottom */
2254:       1, 4,  5, /* right */
2255:       2, 3,  5, /* left */
2256:     };
2257:     PetscInt       faceVerticesTri[9] = {
2258:       0, 3,  4, /* bottom */
2259:       1, 4,  5, /* right */
2260:       2, 5,  3, /* left */
2261:     };

2263:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2264:     PetscSortInt(faceSizeTri, sortedIndices);
2265:     for (iFace = 0; iFace < 3; ++iFace) {
2266:       const PetscInt ii = iFace*faceSizeTri;
2267:       PetscInt       fVertex, cVertex;

2269:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2270:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2271:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2272:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2273:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2274:               faceVertices[fVertex] = origVertices[cVertex];
2275:               break;
2276:             }
2277:           }
2278:         }
2279:         found = PETSC_TRUE;
2280:         break;
2281:       }
2282:     }
2283:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2284:     if (posOriented) *posOriented = PETSC_TRUE;
2285:     return(0);
2286:   } else if (cellDim == 2 && numCorners == 9) {
2287:     /* Quadratic quad (I hate this) */
2288:     /* Edges are determined by the first 2 vertices (corners of edges) */
2289:     const PetscInt faceSizeQuad = 3;
2290:     PetscInt       sortedIndices[3], i, iFace;
2291:     PetscBool      found                      = PETSC_FALSE;
2292:     PetscInt       faceVerticesQuadSorted[12] = {
2293:       0, 1,  4, /* bottom */
2294:       1, 2,  5, /* right */
2295:       2, 3,  6, /* top */
2296:       0, 3,  7, /* left */
2297:     };
2298:     PetscInt       faceVerticesQuad[12] = {
2299:       0, 1,  4, /* bottom */
2300:       1, 2,  5, /* right */
2301:       2, 3,  6, /* top */
2302:       3, 0,  7, /* left */
2303:     };

2305:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2306:     PetscSortInt(faceSizeQuad, sortedIndices);
2307:     for (iFace = 0; iFace < 4; ++iFace) {
2308:       const PetscInt ii = iFace*faceSizeQuad;
2309:       PetscInt       fVertex, cVertex;

2311:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2312:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2313:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2314:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2315:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2316:               faceVertices[fVertex] = origVertices[cVertex];
2317:               break;
2318:             }
2319:           }
2320:         }
2321:         found = PETSC_TRUE;
2322:         break;
2323:       }
2324:     }
2325:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2326:     if (posOriented) *posOriented = PETSC_TRUE;
2327:     return(0);
2328:   } else if (cellDim == 3 && numCorners == 8) {
2329:     /* Hexes
2330:        A hex is two oriented quads with the normal of the first
2331:        pointing up at the second.

2333:           7---6
2334:          /|  /|
2335:         4---5 |
2336:         | 1-|-2
2337:         |/  |/
2338:         0---3

2340:         Faces are determined by the first 4 vertices (corners of faces) */
2341:     const PetscInt faceSizeHex = 4;
2342:     PetscInt       sortedIndices[4], i, iFace;
2343:     PetscBool      found                     = PETSC_FALSE;
2344:     PetscInt       faceVerticesHexSorted[24] = {
2345:       0, 1, 2, 3,  /* bottom */
2346:       4, 5, 6, 7,  /* top */
2347:       0, 3, 4, 5,  /* front */
2348:       2, 3, 5, 6,  /* right */
2349:       1, 2, 6, 7,  /* back */
2350:       0, 1, 4, 7,  /* left */
2351:     };
2352:     PetscInt       faceVerticesHex[24] = {
2353:       1, 2, 3, 0,  /* bottom */
2354:       4, 5, 6, 7,  /* top */
2355:       0, 3, 5, 4,  /* front */
2356:       3, 2, 6, 5,  /* right */
2357:       2, 1, 7, 6,  /* back */
2358:       1, 0, 4, 7,  /* left */
2359:     };

2361:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2362:     PetscSortInt(faceSizeHex, sortedIndices);
2363:     for (iFace = 0; iFace < 6; ++iFace) {
2364:       const PetscInt ii = iFace*faceSizeHex;
2365:       PetscInt       fVertex, cVertex;

2367:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2368:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2369:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2370:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2371:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2372:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2373:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2374:               faceVertices[fVertex] = origVertices[cVertex];
2375:               break;
2376:             }
2377:           }
2378:         }
2379:         found = PETSC_TRUE;
2380:         break;
2381:       }
2382:     }
2383:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2384:     if (posOriented) *posOriented = PETSC_TRUE;
2385:     return(0);
2386:   } else if (cellDim == 3 && numCorners == 10) {
2387:     /* Quadratic tet */
2388:     /* Faces are determined by the first 3 vertices (corners of faces) */
2389:     const PetscInt faceSizeTet = 6;
2390:     PetscInt       sortedIndices[6], i, iFace;
2391:     PetscBool      found                     = PETSC_FALSE;
2392:     PetscInt       faceVerticesTetSorted[24] = {
2393:       0, 1, 2,  6, 7, 8, /* bottom */
2394:       0, 3, 4,  6, 7, 9,  /* front */
2395:       1, 4, 5,  7, 8, 9,  /* right */
2396:       2, 3, 5,  6, 8, 9,  /* left */
2397:     };
2398:     PetscInt       faceVerticesTet[24] = {
2399:       0, 1, 2,  6, 7, 8, /* bottom */
2400:       0, 4, 3,  6, 7, 9,  /* front */
2401:       1, 5, 4,  7, 8, 9,  /* right */
2402:       2, 3, 5,  8, 6, 9,  /* left */
2403:     };

2405:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2406:     PetscSortInt(faceSizeTet, sortedIndices);
2407:     for (iFace=0; iFace < 4; ++iFace) {
2408:       const PetscInt ii = iFace*faceSizeTet;
2409:       PetscInt       fVertex, cVertex;

2411:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2412:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2413:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2414:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2415:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2416:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2417:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2418:               faceVertices[fVertex] = origVertices[cVertex];
2419:               break;
2420:             }
2421:           }
2422:         }
2423:         found = PETSC_TRUE;
2424:         break;
2425:       }
2426:     }
2427:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2428:     if (posOriented) *posOriented = PETSC_TRUE;
2429:     return(0);
2430:   } else if (cellDim == 3 && numCorners == 27) {
2431:     /* Quadratic hexes (I hate this)
2432:        A hex is two oriented quads with the normal of the first
2433:        pointing up at the second.

2435:          7---6
2436:         /|  /|
2437:        4---5 |
2438:        | 3-|-2
2439:        |/  |/
2440:        0---1

2442:        Faces are determined by the first 4 vertices (corners of faces) */
2443:     const PetscInt faceSizeQuadHex = 9;
2444:     PetscInt       sortedIndices[9], i, iFace;
2445:     PetscBool      found                         = PETSC_FALSE;
2446:     PetscInt       faceVerticesQuadHexSorted[54] = {
2447:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2448:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2449:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2450:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2451:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2452:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2453:     };
2454:     PetscInt       faceVerticesQuadHex[54] = {
2455:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2456:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2457:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2458:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2459:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2460:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2461:     };

2463:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2464:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2465:     for (iFace = 0; iFace < 6; ++iFace) {
2466:       const PetscInt ii = iFace*faceSizeQuadHex;
2467:       PetscInt       fVertex, cVertex;

2469:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2470:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2471:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2472:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2473:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2474:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2475:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2476:               faceVertices[fVertex] = origVertices[cVertex];
2477:               break;
2478:             }
2479:           }
2480:         }
2481:         found = PETSC_TRUE;
2482:         break;
2483:       }
2484:     }
2485:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2486:     if (posOriented) *posOriented = PETSC_TRUE;
2487:     return(0);
2488:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2489:   if (!posOrient) {
2490:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2491:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2492:   } else {
2493:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2494:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2495:   }
2496:   if (posOriented) *posOriented = posOrient;
2497:   return(0);
2498: }

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

2504:   Not collective

2506:   Input Parameters:
2507: + dm           - The original mesh
2508: . cell         - The cell mesh point
2509: . faceSize     - The number of vertices on the face
2510: . face         - The face vertices
2511: . numCorners   - The number of vertices on the cell
2512: . indices      - Local numbering of face vertices in cell cone
2513: - origVertices - Original face vertices

2515:   Output Parameter:
2516: + faceVertices - The face vertices properly oriented
2517: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2519:   Level: developer

2521: .seealso: DMPlexGetCone()
2522: @*/
2523: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2524: {
2525:   const PetscInt *cone = NULL;
2526:   PetscInt        coneSize, v, f, v2;
2527:   PetscInt        oppositeVertex = -1;
2528:   PetscErrorCode  ierr;

2531:   DMPlexGetConeSize(dm, cell, &coneSize);
2532:   DMPlexGetCone(dm, cell, &cone);
2533:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2534:     PetscBool found = PETSC_FALSE;

2536:     for (f = 0; f < faceSize; ++f) {
2537:       if (face[f] == cone[v]) {
2538:         found = PETSC_TRUE; break;
2539:       }
2540:     }
2541:     if (found) {
2542:       indices[v2]      = v;
2543:       origVertices[v2] = cone[v];
2544:       ++v2;
2545:     } else {
2546:       oppositeVertex = v;
2547:     }
2548:   }
2549:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2550:   return(0);
2551: }

2553: /*
2554:   DMPlexInsertFace_Internal - Puts a face into the mesh

2556:   Not collective

2558:   Input Parameters:
2559:   + dm              - The DMPlex
2560:   . numFaceVertex   - The number of vertices in the face
2561:   . faceVertices    - The vertices in the face for dm
2562:   . subfaceVertices - The vertices in the face for subdm
2563:   . numCorners      - The number of vertices in the cell
2564:   . cell            - A cell in dm containing the face
2565:   . subcell         - A cell in subdm containing the face
2566:   . firstFace       - First face in the mesh
2567:   - newFacePoint    - Next face in the mesh

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

2572:   Level: developer
2573: */
2574: 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)
2575: {
2576:   MPI_Comm        comm;
2577:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2578:   const PetscInt *faces;
2579:   PetscInt        numFaces, coneSize;
2580:   PetscErrorCode  ierr;

2583:   PetscObjectGetComm((PetscObject)dm,&comm);
2584:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2585:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2586: #if 0
2587:   /* Cannot use this because support() has not been constructed yet */
2588:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2589: #else
2590:   {
2591:     PetscInt f;

2593:     numFaces = 0;
2594:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2595:     for (f = firstFace; f < *newFacePoint; ++f) {
2596:       PetscInt dof, off, d;

2598:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2599:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2600:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2601:       for (d = 0; d < dof; ++d) {
2602:         const PetscInt p = submesh->cones[off+d];
2603:         PetscInt       v;

2605:         for (v = 0; v < numFaceVertices; ++v) {
2606:           if (subfaceVertices[v] == p) break;
2607:         }
2608:         if (v == numFaceVertices) break;
2609:       }
2610:       if (d == dof) {
2611:         numFaces               = 1;
2612:         ((PetscInt*) faces)[0] = f;
2613:       }
2614:     }
2615:   }
2616: #endif
2617:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2618:   else if (numFaces == 1) {
2619:     /* Add the other cell neighbor for this face */
2620:     DMPlexSetCone(subdm, subcell, faces);
2621:   } else {
2622:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2623:     PetscBool posOriented;

2625:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2626:     origVertices        = &orientedVertices[numFaceVertices];
2627:     indices             = &orientedVertices[numFaceVertices*2];
2628:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2629:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2630:     /* TODO: I know that routine should return a permutation, not the indices */
2631:     for (v = 0; v < numFaceVertices; ++v) {
2632:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2633:       for (ov = 0; ov < numFaceVertices; ++ov) {
2634:         if (orientedVertices[ov] == vertex) {
2635:           orientedSubVertices[ov] = subvertex;
2636:           break;
2637:         }
2638:       }
2639:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2640:     }
2641:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2642:     DMPlexSetCone(subdm, subcell, newFacePoint);
2643:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2644:     ++(*newFacePoint);
2645:   }
2646: #if 0
2647:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2648: #else
2649:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2650: #endif
2651:   return(0);
2652: }

2654: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2655: {
2656:   MPI_Comm        comm;
2657:   DMLabel         subpointMap;
2658:   IS              subvertexIS,  subcellIS;
2659:   const PetscInt *subVertices, *subCells;
2660:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2661:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2662:   PetscInt        vStart, vEnd, c, f;
2663:   PetscErrorCode  ierr;

2666:   PetscObjectGetComm((PetscObject)dm,&comm);
2667:   /* Create subpointMap which marks the submesh */
2668:   DMLabelCreate("subpoint_map", &subpointMap);
2669:   DMPlexSetSubpointMap(subdm, subpointMap);
2670:   DMLabelDestroy(&subpointMap);
2671:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2672:   /* Setup chart */
2673:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2674:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2675:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2676:   DMPlexSetVTKCellHeight(subdm, 1);
2677:   /* Set cone sizes */
2678:   firstSubVertex = numSubCells;
2679:   firstSubFace   = numSubCells+numSubVertices;
2680:   newFacePoint   = firstSubFace;
2681:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2682:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2683:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2684:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2685:   for (c = 0; c < numSubCells; ++c) {
2686:     DMPlexSetConeSize(subdm, c, 1);
2687:   }
2688:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2689:     DMPlexSetConeSize(subdm, f, nFV);
2690:   }
2691:   DMSetUp(subdm);
2692:   /* Create face cones */
2693:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2694:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2695:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2696:   for (c = 0; c < numSubCells; ++c) {
2697:     const PetscInt cell    = subCells[c];
2698:     const PetscInt subcell = c;
2699:     PetscInt      *closure = NULL;
2700:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2702:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2703:     for (cl = 0; cl < closureSize*2; cl += 2) {
2704:       const PetscInt point = closure[cl];
2705:       PetscInt       subVertex;

2707:       if ((point >= vStart) && (point < vEnd)) {
2708:         ++numCorners;
2709:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2710:         if (subVertex >= 0) {
2711:           closure[faceSize] = point;
2712:           subface[faceSize] = firstSubVertex+subVertex;
2713:           ++faceSize;
2714:         }
2715:       }
2716:     }
2717:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2718:     if (faceSize == nFV) {
2719:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2720:     }
2721:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2722:   }
2723:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2724:   DMPlexSymmetrize(subdm);
2725:   DMPlexStratify(subdm);
2726:   /* Build coordinates */
2727:   {
2728:     PetscSection coordSection, subCoordSection;
2729:     Vec          coordinates, subCoordinates;
2730:     PetscScalar *coords, *subCoords;
2731:     PetscInt     numComp, coordSize, v;
2732:     const char  *name;

2734:     DMGetCoordinateSection(dm, &coordSection);
2735:     DMGetCoordinatesLocal(dm, &coordinates);
2736:     DMGetCoordinateSection(subdm, &subCoordSection);
2737:     PetscSectionSetNumFields(subCoordSection, 1);
2738:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2739:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2740:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2741:     for (v = 0; v < numSubVertices; ++v) {
2742:       const PetscInt vertex    = subVertices[v];
2743:       const PetscInt subvertex = firstSubVertex+v;
2744:       PetscInt       dof;

2746:       PetscSectionGetDof(coordSection, vertex, &dof);
2747:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2748:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2749:     }
2750:     PetscSectionSetUp(subCoordSection);
2751:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2752:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2753:     PetscObjectGetName((PetscObject)coordinates,&name);
2754:     PetscObjectSetName((PetscObject)subCoordinates,name);
2755:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2756:     VecSetType(subCoordinates,VECSTANDARD);
2757:     if (coordSize) {
2758:       VecGetArray(coordinates,    &coords);
2759:       VecGetArray(subCoordinates, &subCoords);
2760:       for (v = 0; v < numSubVertices; ++v) {
2761:         const PetscInt vertex    = subVertices[v];
2762:         const PetscInt subvertex = firstSubVertex+v;
2763:         PetscInt       dof, off, sdof, soff, d;

2765:         PetscSectionGetDof(coordSection, vertex, &dof);
2766:         PetscSectionGetOffset(coordSection, vertex, &off);
2767:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2768:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2769:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2770:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2771:       }
2772:       VecRestoreArray(coordinates,    &coords);
2773:       VecRestoreArray(subCoordinates, &subCoords);
2774:     }
2775:     DMSetCoordinatesLocal(subdm, subCoordinates);
2776:     VecDestroy(&subCoordinates);
2777:   }
2778:   /* Cleanup */
2779:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2780:   ISDestroy(&subvertexIS);
2781:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2782:   ISDestroy(&subcellIS);
2783:   return(0);
2784: }

2786: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2787: {
2788:   PetscInt       subPoint;

2791:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2792:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2793: }

2795: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2796: {
2797:   MPI_Comm         comm;
2798:   DMLabel          subpointMap;
2799:   IS              *subpointIS;
2800:   const PetscInt **subpoints;
2801:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2802:   PetscInt         totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2803:   PetscMPIInt      rank;
2804:   PetscErrorCode   ierr;

2807:   PetscObjectGetComm((PetscObject)dm,&comm);
2808:   MPI_Comm_rank(comm, &rank);
2809:   /* Create subpointMap which marks the submesh */
2810:   DMLabelCreate("subpoint_map", &subpointMap);
2811:   DMPlexSetSubpointMap(subdm, subpointMap);
2812:   if (cellHeight) {
2813:     if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2814:     else            {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2815:   } else {
2816:     DMLabel         depth;
2817:     IS              pointIS;
2818:     const PetscInt *points;
2819:     PetscInt        numPoints;

2821:     DMPlexGetDepthLabel(dm, &depth);
2822:     DMLabelGetStratumSize(label, value, &numPoints);
2823:     DMLabelGetStratumIS(label, value, &pointIS);
2824:     ISGetIndices(pointIS, &points);
2825:     for (p = 0; p < numPoints; ++p) {
2826:       PetscInt *closure = NULL;
2827:       PetscInt  closureSize, c, pdim;

2829:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2830:       for (c = 0; c < closureSize*2; c += 2) {
2831:         DMLabelGetValue(depth, closure[c], &pdim);
2832:         DMLabelSetValue(subpointMap, closure[c], pdim);
2833:       }
2834:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2835:     }
2836:     ISRestoreIndices(pointIS, &points);
2837:     ISDestroy(&pointIS);
2838:   }
2839:   DMLabelDestroy(&subpointMap);
2840:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2841:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2842:   cMax = (cMax < 0) ? cEnd : cMax;
2843:   /* Setup chart */
2844:   DMGetDimension(dm, &dim);
2845:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2846:   for (d = 0; d <= dim; ++d) {
2847:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2848:     totSubPoints += numSubPoints[d];
2849:   }
2850:   DMPlexSetChart(subdm, 0, totSubPoints);
2851:   DMPlexSetVTKCellHeight(subdm, cellHeight);
2852:   /* Set cone sizes */
2853:   firstSubPoint[dim] = 0;
2854:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2855:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2856:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2857:   for (d = 0; d <= dim; ++d) {
2858:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2859:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2860:   }
2861:   for (d = 0; d <= dim; ++d) {
2862:     for (p = 0; p < numSubPoints[d]; ++p) {
2863:       const PetscInt  point    = subpoints[d][p];
2864:       const PetscInt  subpoint = firstSubPoint[d] + p;
2865:       const PetscInt *cone;
2866:       PetscInt        coneSize, coneSizeNew, c, val;

2868:       DMPlexGetConeSize(dm, point, &coneSize);
2869:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2870:       if (cellHeight && (d == dim)) {
2871:         DMPlexGetCone(dm, point, &cone);
2872:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2873:           DMLabelGetValue(subpointMap, cone[c], &val);
2874:           if (val >= 0) coneSizeNew++;
2875:         }
2876:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2877:       }
2878:     }
2879:   }
2880:   DMSetUp(subdm);
2881:   /* Set cones */
2882:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2883:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2884:   for (d = 0; d <= dim; ++d) {
2885:     for (p = 0; p < numSubPoints[d]; ++p) {
2886:       const PetscInt  point    = subpoints[d][p];
2887:       const PetscInt  subpoint = firstSubPoint[d] + p;
2888:       const PetscInt *cone, *ornt;
2889:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

2891:       if (d == dim-1) {
2892:         const PetscInt *support, *cone, *ornt;
2893:         PetscInt        supportSize, coneSize, s, subc;

2895:         DMPlexGetSupport(dm, point, &support);
2896:         DMPlexGetSupportSize(dm, point, &supportSize);
2897:         for (s = 0; s < supportSize; ++s) {
2898:           if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2899:           PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2900:           if (subc >= 0) {
2901:             const PetscInt ccell = subpoints[d+1][subc];

2903:             DMPlexGetCone(dm, ccell, &cone);
2904:             DMPlexGetConeSize(dm, ccell, &coneSize);
2905:             DMPlexGetConeOrientation(dm, ccell, &ornt);
2906:             for (c = 0; c < coneSize; ++c) {
2907:               if (cone[c] == point) {
2908:                 fornt = ornt[c];
2909:                 break;
2910:               }
2911:             }
2912:             break;
2913:           }
2914:         }
2915:       }
2916:       DMPlexGetConeSize(dm, point, &coneSize);
2917:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2918:       DMPlexGetCone(dm, point, &cone);
2919:       DMPlexGetConeOrientation(dm, point, &ornt);
2920:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2921:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2922:         if (subc >= 0) {
2923:           coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2924:           orntNew[coneSizeNew] = ornt[c];
2925:           ++coneSizeNew;
2926:         }
2927:       }
2928:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2929:       if (fornt < 0) {
2930:         /* This should be replaced by a call to DMPlexReverseCell() */
2931: #if 0
2932:         DMPlexReverseCell(subdm, subpoint);
2933: #else
2934:         for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
2935:           PetscInt faceSize, tmp;

2937:           tmp        = coneNew[c];
2938:           coneNew[c] = coneNew[coneSizeNew-1-c];
2939:           coneNew[coneSizeNew-1-c] = tmp;
2940:           DMPlexGetConeSize(dm, cone[c], &faceSize);
2941:           tmp        = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2942:           orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2943:           orntNew[coneSizeNew-1-c] = tmp;
2944:         }
2945:       }
2946:       DMPlexSetCone(subdm, subpoint, coneNew);
2947:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2948: #endif
2949:     }
2950:   }
2951:   PetscFree2(coneNew,orntNew);
2952:   DMPlexSymmetrize(subdm);
2953:   DMPlexStratify(subdm);
2954:   /* Build coordinates */
2955:   {
2956:     PetscSection coordSection, subCoordSection;
2957:     Vec          coordinates, subCoordinates;
2958:     PetscScalar *coords, *subCoords;
2959:     PetscInt     cdim, numComp, coordSize;
2960:     const char  *name;

2962:     DMGetCoordinateDim(dm, &cdim);
2963:     DMGetCoordinateSection(dm, &coordSection);
2964:     DMGetCoordinatesLocal(dm, &coordinates);
2965:     DMGetCoordinateSection(subdm, &subCoordSection);
2966:     PetscSectionSetNumFields(subCoordSection, 1);
2967:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2968:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2969:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2970:     for (v = 0; v < numSubPoints[0]; ++v) {
2971:       const PetscInt vertex    = subpoints[0][v];
2972:       const PetscInt subvertex = firstSubPoint[0]+v;
2973:       PetscInt       dof;

2975:       PetscSectionGetDof(coordSection, vertex, &dof);
2976:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2977:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2978:     }
2979:     PetscSectionSetUp(subCoordSection);
2980:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2981:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
2982:     PetscObjectGetName((PetscObject)coordinates,&name);
2983:     PetscObjectSetName((PetscObject)subCoordinates,name);
2984:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2985:     VecSetBlockSize(subCoordinates, cdim);
2986:     VecSetType(subCoordinates,VECSTANDARD);
2987:     VecGetArray(coordinates,    &coords);
2988:     VecGetArray(subCoordinates, &subCoords);
2989:     for (v = 0; v < numSubPoints[0]; ++v) {
2990:       const PetscInt vertex    = subpoints[0][v];
2991:       const PetscInt subvertex = firstSubPoint[0]+v;
2992:       PetscInt dof, off, sdof, soff, d;

2994:       PetscSectionGetDof(coordSection, vertex, &dof);
2995:       PetscSectionGetOffset(coordSection, vertex, &off);
2996:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2997:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2998:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2999:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3000:     }
3001:     VecRestoreArray(coordinates,    &coords);
3002:     VecRestoreArray(subCoordinates, &subCoords);
3003:     DMSetCoordinatesLocal(subdm, subCoordinates);
3004:     VecDestroy(&subCoordinates);
3005:   }
3006:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3007:   {
3008:     PetscSF            sfPoint, sfPointSub;
3009:     IS                 subpIS;
3010:     const PetscSFNode *remotePoints;
3011:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3012:     const PetscInt    *localPoints, *subpoints;
3013:     PetscInt          *slocalPoints;
3014:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3015:     PetscMPIInt        rank;

3017:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3018:     DMGetPointSF(dm, &sfPoint);
3019:     DMGetPointSF(subdm, &sfPointSub);
3020:     DMPlexGetChart(dm, &pStart, &pEnd);
3021:     DMPlexGetChart(subdm, NULL, &numSubroots);
3022:     DMPlexCreateSubpointIS(subdm, &subpIS);
3023:     if (subpIS) {
3024:       ISGetIndices(subpIS, &subpoints);
3025:       ISGetLocalSize(subpIS, &numSubpoints);
3026:     }
3027:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3028:     if (numRoots >= 0) {
3029:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3030:       for (p = 0; p < pEnd-pStart; ++p) {
3031:         newLocalPoints[p].rank  = -2;
3032:         newLocalPoints[p].index = -2;
3033:       }
3034:       /* Set subleaves */
3035:       for (l = 0; l < numLeaves; ++l) {
3036:         const PetscInt point    = localPoints[l];
3037:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3039:         if (subpoint < 0) continue;
3040:         newLocalPoints[point-pStart].rank  = rank;
3041:         newLocalPoints[point-pStart].index = subpoint;
3042:         ++numSubleaves;
3043:       }
3044:       /* Must put in owned subpoints */
3045:       for (p = pStart; p < pEnd; ++p) {
3046:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3048:         if (subpoint < 0) {
3049:           newOwners[p-pStart].rank  = -3;
3050:           newOwners[p-pStart].index = -3;
3051:         } else {
3052:           newOwners[p-pStart].rank  = rank;
3053:           newOwners[p-pStart].index = subpoint;
3054:         }
3055:       }
3056:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3057:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3058:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3059:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3060:       PetscMalloc1(numSubleaves, &slocalPoints);
3061:       PetscMalloc1(numSubleaves, &sremotePoints);
3062:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3063:         const PetscInt point    = localPoints[l];
3064:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3066:         if (subpoint < 0) continue;
3067:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3068:         slocalPoints[sl]        = subpoint;
3069:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3070:         sremotePoints[sl].index = newLocalPoints[point].index;
3071:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3072:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3073:         ++sl;
3074:       }
3075:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3076:       PetscFree2(newLocalPoints,newOwners);
3077:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3078:     }
3079:     if (subpIS) {
3080:       ISRestoreIndices(subpIS, &subpoints);
3081:       ISDestroy(&subpIS);
3082:     }
3083:   }
3084:   /* Cleanup */
3085:   for (d = 0; d <= dim; ++d) {
3086:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3087:     ISDestroy(&subpointIS[d]);
3088:   }
3089:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3090:   return(0);
3091: }

3093: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3094: {

3098:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3099:   return(0);
3100: }

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

3105:   Input Parameters:
3106: + dm           - The original mesh
3107: . vertexLabel  - The DMLabel marking vertices contained in the surface
3108: - value        - The label value to use

3110:   Output Parameter:
3111: . subdm - The surface mesh

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

3115:   Level: developer

3117: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3118: @*/
3119: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3120: {
3121:   PetscInt       dim, cdim, depth;

3127:   DMGetDimension(dm, &dim);
3128:   DMPlexGetDepth(dm, &depth);
3129:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3130:   DMSetType(*subdm, DMPLEX);
3131:   DMSetDimension(*subdm, dim-1);
3132:   DMGetCoordinateDim(dm, &cdim);
3133:   DMSetCoordinateDim(*subdm, cdim);
3134:   if (depth == dim) {
3135:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3136:   } else {
3137:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3138:   }
3139:   return(0);
3140: }

3142: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3143: {
3144:   MPI_Comm        comm;
3145:   DMLabel         subpointMap;
3146:   IS              subvertexIS;
3147:   const PetscInt *subVertices;
3148:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3149:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3150:   PetscInt        cMax, c, f;
3151:   PetscErrorCode  ierr;

3154:   PetscObjectGetComm((PetscObject)dm, &comm);
3155:   /* Create subpointMap which marks the submesh */
3156:   DMLabelCreate("subpoint_map", &subpointMap);
3157:   DMPlexSetSubpointMap(subdm, subpointMap);
3158:   DMLabelDestroy(&subpointMap);
3159:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3160:   /* Setup chart */
3161:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3162:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3163:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3164:   DMPlexSetVTKCellHeight(subdm, 1);
3165:   /* Set cone sizes */
3166:   firstSubVertex = numSubCells;
3167:   firstSubFace   = numSubCells+numSubVertices;
3168:   newFacePoint   = firstSubFace;
3169:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3170:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3171:   for (c = 0; c < numSubCells; ++c) {
3172:     DMPlexSetConeSize(subdm, c, 1);
3173:   }
3174:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3175:     DMPlexSetConeSize(subdm, f, nFV);
3176:   }
3177:   DMSetUp(subdm);
3178:   /* Create face cones */
3179:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3180:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3181:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3182:   for (c = 0; c < numSubCells; ++c) {
3183:     const PetscInt  cell    = subCells[c];
3184:     const PetscInt  subcell = c;
3185:     const PetscInt *cone, *cells;
3186:     PetscInt        numCells, subVertex, p, v;

3188:     if (cell < cMax) continue;
3189:     DMPlexGetCone(dm, cell, &cone);
3190:     for (v = 0; v < nFV; ++v) {
3191:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3192:       subface[v] = firstSubVertex+subVertex;
3193:     }
3194:     DMPlexSetCone(subdm, newFacePoint, subface);
3195:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3196:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3197:     /* Not true in parallel
3198:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3199:     for (p = 0; p < numCells; ++p) {
3200:       PetscInt negsubcell;

3202:       if (cells[p] >= cMax) continue;
3203:       /* I know this is a crap search */
3204:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3205:         if (subCells[negsubcell] == cells[p]) break;
3206:       }
3207:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3208:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3209:     }
3210:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3211:     ++newFacePoint;
3212:   }
3213:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3214:   DMPlexSymmetrize(subdm);
3215:   DMPlexStratify(subdm);
3216:   /* Build coordinates */
3217:   {
3218:     PetscSection coordSection, subCoordSection;
3219:     Vec          coordinates, subCoordinates;
3220:     PetscScalar *coords, *subCoords;
3221:     PetscInt     cdim, numComp, coordSize, v;
3222:     const char  *name;

3224:     DMGetCoordinateDim(dm, &cdim);
3225:     DMGetCoordinateSection(dm, &coordSection);
3226:     DMGetCoordinatesLocal(dm, &coordinates);
3227:     DMGetCoordinateSection(subdm, &subCoordSection);
3228:     PetscSectionSetNumFields(subCoordSection, 1);
3229:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3230:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3231:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3232:     for (v = 0; v < numSubVertices; ++v) {
3233:       const PetscInt vertex    = subVertices[v];
3234:       const PetscInt subvertex = firstSubVertex+v;
3235:       PetscInt       dof;

3237:       PetscSectionGetDof(coordSection, vertex, &dof);
3238:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3239:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3240:     }
3241:     PetscSectionSetUp(subCoordSection);
3242:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3243:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3244:     PetscObjectGetName((PetscObject)coordinates,&name);
3245:     PetscObjectSetName((PetscObject)subCoordinates,name);
3246:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3247:     VecSetBlockSize(subCoordinates, cdim);
3248:     VecSetType(subCoordinates,VECSTANDARD);
3249:     VecGetArray(coordinates,    &coords);
3250:     VecGetArray(subCoordinates, &subCoords);
3251:     for (v = 0; v < numSubVertices; ++v) {
3252:       const PetscInt vertex    = subVertices[v];
3253:       const PetscInt subvertex = firstSubVertex+v;
3254:       PetscInt       dof, off, sdof, soff, d;

3256:       PetscSectionGetDof(coordSection, vertex, &dof);
3257:       PetscSectionGetOffset(coordSection, vertex, &off);
3258:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3259:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3260:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3261:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3262:     }
3263:     VecRestoreArray(coordinates,    &coords);
3264:     VecRestoreArray(subCoordinates, &subCoords);
3265:     DMSetCoordinatesLocal(subdm, subCoordinates);
3266:     VecDestroy(&subCoordinates);
3267:   }
3268:   /* Build SF */
3269:   CHKMEMQ;
3270:   {
3271:     PetscSF            sfPoint, sfPointSub;
3272:     const PetscSFNode *remotePoints;
3273:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3274:     const PetscInt    *localPoints;
3275:     PetscInt          *slocalPoints;
3276:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3277:     PetscMPIInt        rank;

3279:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3280:     DMGetPointSF(dm, &sfPoint);
3281:     DMGetPointSF(subdm, &sfPointSub);
3282:     DMPlexGetChart(dm, &pStart, &pEnd);
3283:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3284:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3285:     if (numRoots >= 0) {
3286:       /* Only vertices should be shared */
3287:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3288:       for (p = 0; p < pEnd-pStart; ++p) {
3289:         newLocalPoints[p].rank  = -2;
3290:         newLocalPoints[p].index = -2;
3291:       }
3292:       /* Set subleaves */
3293:       for (l = 0; l < numLeaves; ++l) {
3294:         const PetscInt point    = localPoints[l];
3295:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3297:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3298:         if (subPoint < 0) continue;
3299:         newLocalPoints[point-pStart].rank  = rank;
3300:         newLocalPoints[point-pStart].index = subPoint;
3301:         ++numSubLeaves;
3302:       }
3303:       /* Must put in owned subpoints */
3304:       for (p = pStart; p < pEnd; ++p) {
3305:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3307:         if (subPoint < 0) {
3308:           newOwners[p-pStart].rank  = -3;
3309:           newOwners[p-pStart].index = -3;
3310:         } else {
3311:           newOwners[p-pStart].rank  = rank;
3312:           newOwners[p-pStart].index = subPoint;
3313:         }
3314:       }
3315:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3316:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3317:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3318:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3319:       PetscMalloc1(numSubLeaves,    &slocalPoints);
3320:       PetscMalloc1(numSubLeaves, &sremotePoints);
3321:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3322:         const PetscInt point    = localPoints[l];
3323:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3325:         if (subPoint < 0) continue;
3326:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3327:         slocalPoints[sl]        = subPoint;
3328:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3329:         sremotePoints[sl].index = newLocalPoints[point].index;
3330:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3331:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3332:         ++sl;
3333:       }
3334:       PetscFree2(newLocalPoints,newOwners);
3335:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3336:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3337:     }
3338:   }
3339:   CHKMEMQ;
3340:   /* Cleanup */
3341:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3342:   ISDestroy(&subvertexIS);
3343:   PetscFree(subCells);
3344:   return(0);
3345: }

3347: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3348: {
3349:   DMLabel        label = NULL;

3353:   if (labelname) {DMGetLabel(dm, labelname, &label);}
3354:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3355:   return(0);
3356: }

3358: /*@C
3359:   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.

3361:   Input Parameters:
3362: + dm          - The original mesh
3363: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3364: . label       - A label name, or NULL
3365: - value  - A label value

3367:   Output Parameter:
3368: . subdm - The surface mesh

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

3372:   Level: developer

3374: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3375: @*/
3376: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3377: {
3378:   PetscInt       dim, cdim, depth;

3384:   DMGetDimension(dm, &dim);
3385:   DMPlexGetDepth(dm, &depth);
3386:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3387:   DMSetType(*subdm, DMPLEX);
3388:   DMSetDimension(*subdm, dim-1);
3389:   DMGetCoordinateDim(dm, &cdim);
3390:   DMSetCoordinateDim(*subdm, cdim);
3391:   if (depth == dim) {
3392:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3393:   } else {
3394:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3395:   }
3396:   return(0);
3397: }

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

3402:   Input Parameters:
3403: + dm        - The original mesh
3404: . cellLabel - The DMLabel marking cells contained in the new mesh
3405: - value     - The label value to use

3407:   Output Parameter:
3408: . subdm - The new mesh

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

3412:   Level: developer

3414: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3415: @*/
3416: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3417: {
3418:   PetscInt       dim;

3424:   DMGetDimension(dm, &dim);
3425:   DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3426:   DMSetType(*subdm, DMPLEX);
3427:   DMSetDimension(*subdm, dim);
3428:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3429:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3430:   return(0);
3431: }

3433: /*@
3434:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3436:   Input Parameter:
3437: . dm - The submesh DM

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

3442:   Level: developer

3444: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3445: @*/
3446: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3447: {
3451:   *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3452:   return(0);
3453: }

3455: /*@
3456:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3458:   Input Parameters:
3459: + dm - The submesh DM
3460: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3464:   Level: developer

3466: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3467: @*/
3468: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3469: {
3470:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3471:   DMLabel        tmp;

3476:   tmp  = mesh->subpointMap;
3477:   mesh->subpointMap = subpointMap;
3478:   ++mesh->subpointMap->refct;
3479:   DMLabelDestroy(&tmp);
3480:   return(0);
3481: }

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

3486:   Input Parameter:
3487: . dm - The submesh DM

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

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

3494:   Level: developer

3496: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3497: @*/
3498: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3499: {
3500:   MPI_Comm        comm;
3501:   DMLabel         subpointMap;
3502:   IS              is;
3503:   const PetscInt *opoints;
3504:   PetscInt       *points, *depths;
3505:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3506:   PetscErrorCode  ierr;

3511:   PetscObjectGetComm((PetscObject)dm,&comm);
3512:   *subpointIS = NULL;
3513:   DMPlexGetSubpointMap(dm, &subpointMap);
3514:   DMPlexGetDepth(dm, &depth);
3515:   if (subpointMap && depth >= 0) {
3516:     DMPlexGetChart(dm, &pStart, &pEnd);
3517:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3518:     DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3519:     depths[0] = depth;
3520:     depths[1] = 0;
3521:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3522:     PetscMalloc1(pEnd, &points);
3523:     for(d = 0, off = 0; d <= depth; ++d) {
3524:       const PetscInt dep = depths[d];

3526:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3527:       DMLabelGetStratumSize(subpointMap, dep, &n);
3528:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3529:         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);
3530:       } else {
3531:         if (!n) {
3532:           if (d == 0) {
3533:             /* Missing cells */
3534:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3535:           } else {
3536:             /* Missing faces */
3537:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3538:           }
3539:         }
3540:       }
3541:       if (n) {
3542:         DMLabelGetStratumIS(subpointMap, dep, &is);
3543:         ISGetIndices(is, &opoints);
3544:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3545:         ISRestoreIndices(is, &opoints);
3546:         ISDestroy(&is);
3547:       }
3548:     }
3549:     DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3550:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3551:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3552:   }
3553:   return(0);
3554: }