Actual source code: plexsubmesh.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/dmlabelimpl.h>
  3: #include <petscsf.h>

  5: static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid)
  6: {
  7:   DMPolytopeType ct;

  9:   DMPlexGetCellType(dm, p, &ct);
 10:   switch (ct) {
 11:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
 12:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
 13:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
 14:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
 15:     *isHybrid = PETSC_TRUE;
 16:     break;
 17:   default:
 18:     *isHybrid = PETSC_FALSE;
 19:     break;
 20:   }
 21:   return 0;
 22: }

 24: static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
 25: {
 26:   DMLabel ctLabel;

 28:   if (cStart) *cStart = -1;
 29:   if (cEnd) *cEnd = -1;
 30:   DMPlexGetCellTypeLabel(dm, &ctLabel);
 31:   switch (dim) {
 32:   case 1:
 33:     DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd);
 34:     break;
 35:   case 2:
 36:     DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd);
 37:     break;
 38:   case 3:
 39:     DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd);
 40:     if (*cStart < 0) DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd);
 41:     break;
 42:   default:
 43:     return 0;
 44:   }
 45:   return 0;
 46: }

 48: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
 49: {
 50:   PetscInt fStart, fEnd, f;

 52:   DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
 53:   for (f = fStart; f < fEnd; ++f) {
 54:     PetscInt supportSize;

 56:     DMPlexGetSupportSize(dm, f, &supportSize);
 57:     if (supportSize == 1) {
 58:       if (val < 0) {
 59:         PetscInt *closure = NULL;
 60:         PetscInt  clSize, cl, cval;

 62:         DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 63:         for (cl = 0; cl < clSize * 2; cl += 2) {
 64:           DMLabelGetValue(label, closure[cl], &cval);
 65:           if (cval < 0) continue;
 66:           DMLabelSetValue(label, f, cval);
 67:           break;
 68:         }
 69:         if (cl == clSize * 2) DMLabelSetValue(label, f, 1);
 70:         DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
 71:       } else {
 72:         DMLabelSetValue(label, f, val);
 73:       }
 74:     }
 75:   }
 76:   return 0;
 77: }

 79: /*@
 80:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 82:   Not Collective

 84:   Input Parameters:
 85: + dm - The original DM
 86: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)

 88:   Output Parameter:
 89: . label - The DMLabel marking boundary faces with the given value

 91:   Level: developer

 93: .seealso: `DMLabelCreate()`, `DMCreateLabel()`
 94: @*/
 95: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
 96: {
 97:   DMPlexInterpolatedFlag flg;

100:   DMPlexIsInterpolated(dm, &flg);
102:   DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
103:   return 0;
104: }

106: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
107: {
108:   IS              valueIS;
109:   PetscSF         sfPoint;
110:   const PetscInt *values;
111:   PetscInt        numValues, v, cStart, cEnd, nroots;

113:   DMLabelGetNumValues(label, &numValues);
114:   DMLabelGetValueIS(label, &valueIS);
115:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
116:   ISGetIndices(valueIS, &values);
117:   for (v = 0; v < numValues; ++v) {
118:     IS              pointIS;
119:     const PetscInt *points;
120:     PetscInt        numPoints, p;

122:     DMLabelGetStratumSize(label, values[v], &numPoints);
123:     DMLabelGetStratumIS(label, values[v], &pointIS);
124:     ISGetIndices(pointIS, &points);
125:     for (p = 0; p < numPoints; ++p) {
126:       PetscInt  q       = points[p];
127:       PetscInt *closure = NULL;
128:       PetscInt  closureSize, c;

130:       if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
131:         continue;
132:       }
133:       DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
134:       for (c = 0; c < closureSize * 2; c += 2) DMLabelSetValue(label, closure[c], values[v]);
135:       DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
136:     }
137:     ISRestoreIndices(pointIS, &points);
138:     ISDestroy(&pointIS);
139:   }
140:   ISRestoreIndices(valueIS, &values);
141:   ISDestroy(&valueIS);
142:   DMGetPointSF(dm, &sfPoint);
143:   PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
144:   if (nroots >= 0) {
145:     DMLabel         lblRoots, lblLeaves;
146:     IS              valueIS, pointIS;
147:     const PetscInt *values;
148:     PetscInt        numValues, v;

150:     /* Pull point contributions from remote leaves into local roots */
151:     DMLabelGather(label, sfPoint, &lblLeaves);
152:     DMLabelGetValueIS(lblLeaves, &valueIS);
153:     ISGetLocalSize(valueIS, &numValues);
154:     ISGetIndices(valueIS, &values);
155:     for (v = 0; v < numValues; ++v) {
156:       const PetscInt value = values[v];

158:       DMLabelGetStratumIS(lblLeaves, value, &pointIS);
159:       DMLabelInsertIS(label, pointIS, value);
160:       ISDestroy(&pointIS);
161:     }
162:     ISRestoreIndices(valueIS, &values);
163:     ISDestroy(&valueIS);
164:     DMLabelDestroy(&lblLeaves);
165:     /* Push point contributions from roots into remote leaves */
166:     DMLabelDistribute(label, sfPoint, &lblRoots);
167:     DMLabelGetValueIS(lblRoots, &valueIS);
168:     ISGetLocalSize(valueIS, &numValues);
169:     ISGetIndices(valueIS, &values);
170:     for (v = 0; v < numValues; ++v) {
171:       const PetscInt value = values[v];

173:       DMLabelGetStratumIS(lblRoots, value, &pointIS);
174:       DMLabelInsertIS(label, pointIS, value);
175:       ISDestroy(&pointIS);
176:     }
177:     ISRestoreIndices(valueIS, &values);
178:     ISDestroy(&valueIS);
179:     DMLabelDestroy(&lblRoots);
180:   }
181:   return 0;
182: }

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

187:   Input Parameters:
188: + dm - The DM
189: - label - A DMLabel marking the surface points

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

194:   Level: developer

196: .seealso: `DMPlexLabelCohesiveComplete()`
197: @*/
198: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
199: {
200:   DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
201:   return 0;
202: }

204: /*@
205:   DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point

207:   Input Parameters:
208: + dm - The DM
209: - label - A DMLabel marking the surface points

211:   Output Parameter:
212: . label - A DMLabel incorporating cells

214:   Level: developer

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

218: .seealso: `DMPlexLabelAddFaceCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
219: @*/
220: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
221: {
222:   IS              valueIS;
223:   const PetscInt *values;
224:   PetscInt        numValues, v, cStart, cEnd;

226:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
227:   DMLabelGetNumValues(label, &numValues);
228:   DMLabelGetValueIS(label, &valueIS);
229:   ISGetIndices(valueIS, &values);
230:   for (v = 0; v < numValues; ++v) {
231:     IS              pointIS;
232:     const PetscInt *points;
233:     PetscInt        numPoints, p;

235:     DMLabelGetStratumSize(label, values[v], &numPoints);
236:     DMLabelGetStratumIS(label, values[v], &pointIS);
237:     ISGetIndices(pointIS, &points);
238:     for (p = 0; p < numPoints; ++p) {
239:       PetscInt *closure = NULL;
240:       PetscInt  closureSize, cl;

242:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
243:       for (cl = closureSize - 1; cl > 0; --cl) {
244:         const PetscInt cell = closure[cl * 2];
245:         if ((cell >= cStart) && (cell < cEnd)) {
246:           DMLabelSetValue(label, cell, values[v]);
247:           break;
248:         }
249:       }
250:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
251:     }
252:     ISRestoreIndices(pointIS, &points);
253:     ISDestroy(&pointIS);
254:   }
255:   ISRestoreIndices(valueIS, &values);
256:   ISDestroy(&valueIS);
257:   return 0;
258: }

260: /*@
261:   DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face

263:   Input Parameters:
264: + dm - The DM
265: - label - A DMLabel marking the surface points

267:   Output Parameter:
268: . label - A DMLabel incorporating cells

270:   Level: developer

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

274: .seealso: `DMPlexLabelAddCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
275: @*/
276: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
277: {
278:   IS              valueIS;
279:   const PetscInt *values;
280:   PetscInt        numValues, v, cStart, cEnd, fStart, fEnd;

282:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
283:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
284:   DMLabelGetNumValues(label, &numValues);
285:   DMLabelGetValueIS(label, &valueIS);
286:   ISGetIndices(valueIS, &values);
287:   for (v = 0; v < numValues; ++v) {
288:     IS              pointIS;
289:     const PetscInt *points;
290:     PetscInt        numPoints, p;

292:     DMLabelGetStratumSize(label, values[v], &numPoints);
293:     DMLabelGetStratumIS(label, values[v], &pointIS);
294:     ISGetIndices(pointIS, &points);
295:     for (p = 0; p < numPoints; ++p) {
296:       const PetscInt face    = points[p];
297:       PetscInt      *closure = NULL;
298:       PetscInt       closureSize, cl;

300:       if ((face < fStart) || (face >= fEnd)) continue;
301:       DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
302:       for (cl = closureSize - 1; cl > 0; --cl) {
303:         const PetscInt cell = closure[cl * 2];
304:         if ((cell >= cStart) && (cell < cEnd)) {
305:           DMLabelSetValue(label, cell, values[v]);
306:           break;
307:         }
308:       }
309:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
310:     }
311:     ISRestoreIndices(pointIS, &points);
312:     ISDestroy(&pointIS);
313:   }
314:   ISRestoreIndices(valueIS, &values);
315:   ISDestroy(&valueIS);
316:   return 0;
317: }

319: /*@
320:   DMPlexLabelClearCells - Remove cells from a label

322:   Input Parameters:
323: + dm - The DM
324: - label - A DMLabel marking surface points and their adjacent cells

326:   Output Parameter:
327: . label - A DMLabel without cells

329:   Level: developer

331:   Note: This undoes DMPlexLabelAddCells() or DMPlexLabelAddFaceCells()

333: .seealso: `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`, `DMPlexLabelAddCells()`
334: @*/
335: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
336: {
337:   IS              valueIS;
338:   const PetscInt *values;
339:   PetscInt        numValues, v, cStart, cEnd;

341:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
342:   DMLabelGetNumValues(label, &numValues);
343:   DMLabelGetValueIS(label, &valueIS);
344:   ISGetIndices(valueIS, &values);
345:   for (v = 0; v < numValues; ++v) {
346:     IS              pointIS;
347:     const PetscInt *points;
348:     PetscInt        numPoints, p;

350:     DMLabelGetStratumSize(label, values[v], &numPoints);
351:     DMLabelGetStratumIS(label, values[v], &pointIS);
352:     ISGetIndices(pointIS, &points);
353:     for (p = 0; p < numPoints; ++p) {
354:       PetscInt point = points[p];

356:       if (point >= cStart && point < cEnd) DMLabelClearValue(label, point, values[v]);
357:     }
358:     ISRestoreIndices(pointIS, &points);
359:     ISDestroy(&pointIS);
360:   }
361:   ISRestoreIndices(valueIS, &values);
362:   ISDestroy(&valueIS);
363:   return 0;
364: }

366: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
367:  * index (skipping first, which is (0,0)) */
368: static inline PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
369: {
370:   PetscInt d, off = 0;

372:   /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
373:   for (d = 0; d < depth; d++) {
374:     PetscInt firstd     = d;
375:     PetscInt firstStart = depthShift[2 * d];
376:     PetscInt e;

378:     for (e = d + 1; e <= depth; e++) {
379:       if (depthShift[2 * e] < firstStart) {
380:         firstd     = e;
381:         firstStart = depthShift[2 * d];
382:       }
383:     }
384:     if (firstd != d) {
385:       PetscInt swap[2];

387:       e                     = firstd;
388:       swap[0]               = depthShift[2 * d];
389:       swap[1]               = depthShift[2 * d + 1];
390:       depthShift[2 * d]     = depthShift[2 * e];
391:       depthShift[2 * d + 1] = depthShift[2 * e + 1];
392:       depthShift[2 * e]     = swap[0];
393:       depthShift[2 * e + 1] = swap[1];
394:     }
395:   }
396:   /* convert (oldstart, added) to (oldstart, newstart) */
397:   for (d = 0; d <= depth; d++) {
398:     off += depthShift[2 * d + 1];
399:     depthShift[2 * d + 1] = depthShift[2 * d] + off;
400:   }
401:   return 0;
402: }

404: /* depthShift is a list of (old, new) pairs */
405: static inline PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
406: {
407:   PetscInt d;
408:   PetscInt newOff = 0;

410:   for (d = 0; d <= depth; d++) {
411:     if (p < depthShift[2 * d]) return p + newOff;
412:     else newOff = depthShift[2 * d + 1] - depthShift[2 * d];
413:   }
414:   return p + newOff;
415: }

417: /* depthShift is a list of (old, new) pairs */
418: static inline PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
419: {
420:   PetscInt d;
421:   PetscInt newOff = 0;

423:   for (d = 0; d <= depth; d++) {
424:     if (p < depthShift[2 * d + 1]) return p + newOff;
425:     else newOff = depthShift[2 * d] - depthShift[2 * d + 1];
426:   }
427:   return p + newOff;
428: }

430: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
431: {
432:   PetscInt depth = 0, d, pStart, pEnd, p;
433:   DMLabel  depthLabel;

435:   DMPlexGetDepth(dm, &depth);
436:   if (depth < 0) return 0;
437:   /* Step 1: Expand chart */
438:   DMPlexGetChart(dm, &pStart, &pEnd);
439:   pEnd = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
440:   DMPlexSetChart(dmNew, pStart, pEnd);
441:   DMCreateLabel(dmNew, "depth");
442:   DMPlexGetDepthLabel(dmNew, &depthLabel);
443:   DMCreateLabel(dmNew, "celltype");
444:   /* Step 2: Set cone and support sizes */
445:   for (d = 0; d <= depth; ++d) {
446:     PetscInt pStartNew, pEndNew;
447:     IS       pIS;

449:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
450:     pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
451:     pEndNew   = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
452:     ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
453:     DMLabelSetStratumIS(depthLabel, d, pIS);
454:     ISDestroy(&pIS);
455:     for (p = pStart; p < pEnd; ++p) {
456:       PetscInt       newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
457:       PetscInt       size;
458:       DMPolytopeType ct;

460:       DMPlexGetConeSize(dm, p, &size);
461:       DMPlexSetConeSize(dmNew, newp, size);
462:       DMPlexGetSupportSize(dm, p, &size);
463:       DMPlexSetSupportSize(dmNew, newp, size);
464:       DMPlexGetCellType(dm, p, &ct);
465:       DMPlexSetCellType(dmNew, newp, ct);
466:     }
467:   }
468:   return 0;
469: }

471: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
472: {
473:   PetscInt *newpoints;
474:   PetscInt  depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

476:   DMPlexGetDepth(dm, &depth);
477:   if (depth < 0) return 0;
478:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
479:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
480:   PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)), &newpoints);
481:   /* Step 5: Set cones and supports */
482:   DMPlexGetChart(dm, &pStart, &pEnd);
483:   for (p = pStart; p < pEnd; ++p) {
484:     const PetscInt *points = NULL, *orientations = NULL;
485:     PetscInt        size, sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);

487:     DMPlexGetConeSize(dm, p, &size);
488:     DMPlexGetCone(dm, p, &points);
489:     DMPlexGetConeOrientation(dm, p, &orientations);
490:     for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
491:     DMPlexSetCone(dmNew, newp, newpoints);
492:     DMPlexSetConeOrientation(dmNew, newp, orientations);
493:     DMPlexGetSupportSize(dm, p, &size);
494:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
495:     DMPlexGetSupport(dm, p, &points);
496:     for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
497:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
498:     DMPlexSetSupport(dmNew, newp, newpoints);
499:   }
500:   PetscFree(newpoints);
501:   return 0;
502: }

504: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
505: {
506:   PetscSection coordSection, newCoordSection;
507:   Vec          coordinates, newCoordinates;
508:   PetscScalar *coords, *newCoords;
509:   PetscInt     coordSize, sStart, sEnd;
510:   PetscInt     dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
511:   PetscBool    hasCells;

513:   DMGetCoordinateDim(dm, &dim);
514:   DMSetCoordinateDim(dmNew, dim);
515:   DMPlexGetDepth(dm, &depth);
516:   /* Step 8: Convert coordinates */
517:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
518:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
519:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
520:   DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
521:   DMGetCoordinateSection(dm, &coordSection);
522:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
523:   PetscSectionSetNumFields(newCoordSection, 1);
524:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
525:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
526:   hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
527:   PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
528:   if (hasCells) {
529:     for (c = cStart; c < cEnd; ++c) {
530:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;

532:       PetscSectionGetDof(coordSection, c, &dof);
533:       PetscSectionSetDof(newCoordSection, cNew, dof);
534:       PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
535:     }
536:   }
537:   for (v = vStartNew; v < vEndNew; ++v) {
538:     PetscSectionSetDof(newCoordSection, v, dim);
539:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
540:   }
541:   PetscSectionSetUp(newCoordSection);
542:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
543:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
544:   VecCreate(PETSC_COMM_SELF, &newCoordinates);
545:   PetscObjectSetName((PetscObject)newCoordinates, "coordinates");
546:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
547:   VecSetBlockSize(newCoordinates, dim);
548:   VecSetType(newCoordinates, VECSTANDARD);
549:   DMSetCoordinatesLocal(dmNew, newCoordinates);
550:   DMGetCoordinatesLocal(dm, &coordinates);
551:   VecGetArray(coordinates, &coords);
552:   VecGetArray(newCoordinates, &newCoords);
553:   if (hasCells) {
554:     for (c = cStart; c < cEnd; ++c) {
555:       PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;

557:       PetscSectionGetDof(coordSection, c, &dof);
558:       PetscSectionGetOffset(coordSection, c, &off);
559:       PetscSectionGetOffset(newCoordSection, cNew, &noff);
560:       for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
561:     }
562:   }
563:   for (v = vStart; v < vEnd; ++v) {
564:     PetscInt dof, off, noff, d;

566:     PetscSectionGetDof(coordSection, v, &dof);
567:     PetscSectionGetOffset(coordSection, v, &off);
568:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
569:     for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
570:   }
571:   VecRestoreArray(coordinates, &coords);
572:   VecRestoreArray(newCoordinates, &newCoords);
573:   VecDestroy(&newCoordinates);
574:   PetscSectionDestroy(&newCoordSection);
575:   return 0;
576: }

578: static PetscErrorCode DMPlexShiftSF_Single(DM dm, PetscInt depthShift[], PetscSF sf, PetscSF sfNew)
579: {
580:   const PetscSFNode *remotePoints;
581:   PetscSFNode       *gremotePoints;
582:   const PetscInt    *localPoints;
583:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
584:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, depth = 0, totShift = 0;

586:   DMPlexGetDepth(dm, &depth);
587:   DMPlexGetChart(dm, &pStart, &pEnd);
588:   PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
589:   totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd;
590:   if (numRoots >= 0) {
591:     PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation);
592:     for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
593:     PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
594:     PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
595:     PetscMalloc1(numLeaves, &glocalPoints);
596:     PetscMalloc1(numLeaves, &gremotePoints);
597:     for (l = 0; l < numLeaves; ++l) {
598:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
599:       gremotePoints[l].rank  = remotePoints[l].rank;
600:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
601:     }
602:     PetscFree2(newLocation, newRemoteLocation);
603:     PetscSFSetGraph(sfNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
604:   }
605:   return 0;
606: }

608: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
609: {
610:   PetscSF   sfPoint, sfPointNew;
611:   PetscBool useNatural;

613:   /* Step 9: Convert pointSF */
614:   DMGetPointSF(dm, &sfPoint);
615:   DMGetPointSF(dmNew, &sfPointNew);
616:   DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew);
617:   /* Step 9b: Convert naturalSF */
618:   DMGetUseNatural(dm, &useNatural);
619:   if (useNatural) {
620:     PetscSF sfNat, sfNatNew;

622:     DMSetUseNatural(dmNew, useNatural);
623:     DMGetNaturalSF(dm, &sfNat);
624:     DMGetNaturalSF(dmNew, &sfNatNew);
625:     DMPlexShiftSF_Single(dm, depthShift, sfNat, sfNatNew);
626:   }
627:   return 0;
628: }

630: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
631: {
632:   PetscInt depth = 0, numLabels, l;

634:   DMPlexGetDepth(dm, &depth);
635:   /* Step 10: Convert labels */
636:   DMGetNumLabels(dm, &numLabels);
637:   for (l = 0; l < numLabels; ++l) {
638:     DMLabel         label, newlabel;
639:     const char     *lname;
640:     PetscBool       isDepth, isDim;
641:     IS              valueIS;
642:     const PetscInt *values;
643:     PetscInt        numValues, val;

645:     DMGetLabelName(dm, l, &lname);
646:     PetscStrcmp(lname, "depth", &isDepth);
647:     if (isDepth) continue;
648:     PetscStrcmp(lname, "dim", &isDim);
649:     if (isDim) continue;
650:     DMCreateLabel(dmNew, lname);
651:     DMGetLabel(dm, lname, &label);
652:     DMGetLabel(dmNew, lname, &newlabel);
653:     DMLabelGetDefaultValue(label, &val);
654:     DMLabelSetDefaultValue(newlabel, val);
655:     DMLabelGetValueIS(label, &valueIS);
656:     ISGetLocalSize(valueIS, &numValues);
657:     ISGetIndices(valueIS, &values);
658:     for (val = 0; val < numValues; ++val) {
659:       IS              pointIS;
660:       const PetscInt *points;
661:       PetscInt        numPoints, p;

663:       DMLabelGetStratumIS(label, values[val], &pointIS);
664:       ISGetLocalSize(pointIS, &numPoints);
665:       ISGetIndices(pointIS, &points);
666:       for (p = 0; p < numPoints; ++p) {
667:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);

669:         DMLabelSetValue(newlabel, newpoint, values[val]);
670:       }
671:       ISRestoreIndices(pointIS, &points);
672:       ISDestroy(&pointIS);
673:     }
674:     ISRestoreIndices(valueIS, &values);
675:     ISDestroy(&valueIS);
676:   }
677:   return 0;
678: }

680: static PetscErrorCode DMPlexCreateVTKLabel_Internal(DM dm, PetscBool createGhostLabel, DM dmNew)
681: {
682:   PetscSF            sfPoint;
683:   DMLabel            vtkLabel, ghostLabel = NULL;
684:   const PetscSFNode *leafRemote;
685:   const PetscInt    *leafLocal;
686:   PetscInt           cellHeight, cStart, cEnd, c, fStart, fEnd, f, numLeaves, l;
687:   PetscMPIInt        rank;

689:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
690:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
691:   DMGetPointSF(dm, &sfPoint);
692:   DMPlexGetVTKCellHeight(dmNew, &cellHeight);
693:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
694:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
695:   DMCreateLabel(dmNew, "vtk");
696:   DMGetLabel(dmNew, "vtk", &vtkLabel);
697:   if (createGhostLabel) {
698:     DMCreateLabel(dmNew, "ghost");
699:     DMGetLabel(dmNew, "ghost", &ghostLabel);
700:   }
701:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
702:     for (; c < leafLocal[l] && c < cEnd; ++c) DMLabelSetValue(vtkLabel, c, 1);
703:     if (leafLocal[l] >= cEnd) break;
704:     if (leafRemote[l].rank == rank) {
705:       DMLabelSetValue(vtkLabel, c, 1);
706:     } else if (ghostLabel) DMLabelSetValue(ghostLabel, c, 2);
707:   }
708:   for (; c < cEnd; ++c) DMLabelSetValue(vtkLabel, c, 1);
709:   if (ghostLabel) {
710:     DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
711:     for (f = fStart; f < fEnd; ++f) {
712:       PetscInt numCells;

714:       DMPlexGetSupportSize(dmNew, f, &numCells);
715:       if (numCells < 2) {
716:         DMLabelSetValue(ghostLabel, f, 1);
717:       } else {
718:         const PetscInt *cells = NULL;
719:         PetscInt        vA, vB;

721:         DMPlexGetSupport(dmNew, f, &cells);
722:         DMLabelGetValue(vtkLabel, cells[0], &vA);
723:         DMLabelGetValue(vtkLabel, cells[1], &vB);
724:         if (vA != 1 && vB != 1) DMLabelSetValue(ghostLabel, f, 1);
725:       }
726:     }
727:   }
728:   return 0;
729: }

731: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
732: {
733:   DM           refTree;
734:   PetscSection pSec;
735:   PetscInt    *parents, *childIDs;

737:   DMPlexGetReferenceTree(dm, &refTree);
738:   DMPlexSetReferenceTree(dmNew, refTree);
739:   DMPlexGetTree(dm, &pSec, &parents, &childIDs, NULL, NULL);
740:   if (pSec) {
741:     PetscInt     p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
742:     PetscInt    *childIDsShifted;
743:     PetscSection pSecShifted;

745:     PetscSectionGetChart(pSec, &pStart, &pEnd);
746:     DMPlexGetDepth(dm, &depth);
747:     pStartShifted = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
748:     pEndShifted   = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
749:     PetscMalloc2(pEndShifted - pStartShifted, &parentsShifted, pEndShifted - pStartShifted, &childIDsShifted);
750:     PetscSectionCreate(PetscObjectComm((PetscObject)dmNew), &pSecShifted);
751:     PetscSectionSetChart(pSecShifted, pStartShifted, pEndShifted);
752:     for (p = pStartShifted; p < pEndShifted; p++) {
753:       /* start off assuming no children */
754:       PetscSectionSetDof(pSecShifted, p, 0);
755:     }
756:     for (p = pStart; p < pEnd; p++) {
757:       PetscInt dof;
758:       PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);

760:       PetscSectionGetDof(pSec, p, &dof);
761:       PetscSectionSetDof(pSecShifted, pNew, dof);
762:     }
763:     PetscSectionSetUp(pSecShifted);
764:     for (p = pStart; p < pEnd; p++) {
765:       PetscInt dof;
766:       PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);

768:       PetscSectionGetDof(pSec, p, &dof);
769:       if (dof) {
770:         PetscInt off, offNew;

772:         PetscSectionGetOffset(pSec, p, &off);
773:         PetscSectionGetOffset(pSecShifted, pNew, &offNew);
774:         parentsShifted[offNew]  = DMPlexShiftPoint_Internal(parents[off], depth, depthShift);
775:         childIDsShifted[offNew] = childIDs[off];
776:       }
777:     }
778:     DMPlexSetTree(dmNew, pSecShifted, parentsShifted, childIDsShifted);
779:     PetscFree2(parentsShifted, childIDsShifted);
780:     PetscSectionDestroy(&pSecShifted);
781:   }
782:   return 0;
783: }

785: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
786: {
787:   PetscSF          sf;
788:   IS               valueIS;
789:   const PetscInt  *values, *leaves;
790:   PetscInt        *depthShift;
791:   PetscInt         d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
792:   const PetscReal *maxCell, *Lstart, *L;

794:   DMGetPointSF(dm, &sf);
795:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
796:   nleaves = PetscMax(0, nleaves);
797:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
798:   /* Count ghost cells */
799:   DMLabelGetValueIS(label, &valueIS);
800:   ISGetLocalSize(valueIS, &numFS);
801:   ISGetIndices(valueIS, &values);
802:   Ng = 0;
803:   for (fs = 0; fs < numFS; ++fs) {
804:     IS              faceIS;
805:     const PetscInt *faces;
806:     PetscInt        numFaces, f, numBdFaces = 0;

808:     DMLabelGetStratumIS(label, values[fs], &faceIS);
809:     ISGetLocalSize(faceIS, &numFaces);
810:     ISGetIndices(faceIS, &faces);
811:     for (f = 0; f < numFaces; ++f) {
812:       PetscInt numChildren;

814:       PetscFindInt(faces[f], nleaves, leaves, &loc);
815:       DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
816:       /* non-local and ancestors points don't get to register ghosts */
817:       if (loc >= 0 || numChildren) continue;
818:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
819:     }
820:     Ng += numBdFaces;
821:     ISRestoreIndices(faceIS, &faces);
822:     ISDestroy(&faceIS);
823:   }
824:   DMPlexGetDepth(dm, &depth);
825:   PetscMalloc1(2 * (depth + 1), &depthShift);
826:   for (d = 0; d <= depth; d++) {
827:     PetscInt dEnd;

829:     DMPlexGetDepthStratum(dm, d, NULL, &dEnd);
830:     depthShift[2 * d]     = dEnd;
831:     depthShift[2 * d + 1] = 0;
832:   }
833:   if (depth >= 0) depthShift[2 * depth + 1] = Ng;
834:   DMPlexShiftPointSetUp_Internal(depth, depthShift);
835:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
836:   /* Step 3: Set cone/support sizes for new points */
837:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
838:   for (c = cEnd; c < cEnd + Ng; ++c) DMPlexSetConeSize(gdm, c, 1);
839:   for (fs = 0; fs < numFS; ++fs) {
840:     IS              faceIS;
841:     const PetscInt *faces;
842:     PetscInt        numFaces, f;

844:     DMLabelGetStratumIS(label, values[fs], &faceIS);
845:     ISGetLocalSize(faceIS, &numFaces);
846:     ISGetIndices(faceIS, &faces);
847:     for (f = 0; f < numFaces; ++f) {
848:       PetscInt size, numChildren;

850:       PetscFindInt(faces[f], nleaves, leaves, &loc);
851:       DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
852:       if (loc >= 0 || numChildren) continue;
853:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
854:       DMPlexGetSupportSize(dm, faces[f], &size);
856:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
857:     }
858:     ISRestoreIndices(faceIS, &faces);
859:     ISDestroy(&faceIS);
860:   }
861:   /* Step 4: Setup ghosted DM */
862:   DMSetUp(gdm);
863:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
864:   /* Step 6: Set cones and supports for new points */
865:   ghostCell = cEnd;
866:   for (fs = 0; fs < numFS; ++fs) {
867:     IS              faceIS;
868:     const PetscInt *faces;
869:     PetscInt        numFaces, f;

871:     DMLabelGetStratumIS(label, values[fs], &faceIS);
872:     ISGetLocalSize(faceIS, &numFaces);
873:     ISGetIndices(faceIS, &faces);
874:     for (f = 0; f < numFaces; ++f) {
875:       PetscInt newFace = faces[f] + Ng, numChildren;

877:       PetscFindInt(faces[f], nleaves, leaves, &loc);
878:       DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
879:       if (loc >= 0 || numChildren) continue;
880:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
881:       DMPlexSetCone(gdm, ghostCell, &newFace);
882:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
883:       ++ghostCell;
884:     }
885:     ISRestoreIndices(faceIS, &faces);
886:     ISDestroy(&faceIS);
887:   }
888:   ISRestoreIndices(valueIS, &values);
889:   ISDestroy(&valueIS);
890:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
891:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
892:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
893:   DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm);
894:   DMPlexShiftTree_Internal(dm, depthShift, gdm);
895:   PetscFree(depthShift);
896:   for (c = cEnd; c < cEnd + Ng; ++c) DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST);
897:   /* Step 7: Periodicity */
898:   DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
899:   DMSetPeriodicity(gdm, maxCell, Lstart, L);
900:   if (numGhostCells) *numGhostCells = Ng;
901:   return 0;
902: }

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

907:   Collective on dm

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

913:   Output Parameters:
914: + numGhostCells - The number of ghost cells added to the DM
915: - dmGhosted - The new DM

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

919:   Level: developer

921: .seealso: `DMCreate()`
922: @*/
923: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
924: {
925:   DM          gdm;
926:   DMLabel     label;
927:   const char *name = labelName ? labelName : "Face Sets";
928:   PetscInt    dim, Ng = 0;
929:   PetscBool   useCone, useClosure;

934:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
935:   DMSetType(gdm, DMPLEX);
936:   DMGetDimension(dm, &dim);
937:   DMSetDimension(gdm, dim);
938:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
939:   DMSetBasicAdjacency(gdm, useCone, useClosure);
940:   DMGetLabel(dm, name, &label);
941:   if (!label) {
942:     /* Get label for boundary faces */
943:     DMCreateLabel(dm, name);
944:     DMGetLabel(dm, name, &label);
945:     DMPlexMarkBoundaryFaces(dm, 1, label);
946:   }
947:   DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
948:   DMCopyDisc(dm, gdm);
949:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, gdm);
950:   gdm->setfromoptionscalled = dm->setfromoptionscalled;
951:   if (numGhostCells) *numGhostCells = Ng;
952:   *dmGhosted = gdm;
953:   return 0;
954: }

956: static PetscErrorCode DivideCells_Private(DM dm, DMLabel label, DMPlexPointQueue queue)
957: {
958:   PetscInt dim, d, shift = 100, *pStart, *pEnd;

960:   DMGetDimension(dm, &dim);
961:   PetscMalloc2(dim, &pStart, dim, &pEnd);
962:   for (d = 0; d < dim; ++d) DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
963:   while (!DMPlexPointQueueEmpty(queue)) {
964:     PetscInt  cell    = -1;
965:     PetscInt *closure = NULL;
966:     PetscInt  closureSize, cl, cval;

968:     DMPlexPointQueueDequeue(queue, &cell);
969:     DMLabelGetValue(label, cell, &cval);
970:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
971:     /* Mark points in the cell closure that touch the fault */
972:     for (d = 0; d < dim; ++d) {
973:       for (cl = 0; cl < closureSize * 2; cl += 2) {
974:         const PetscInt clp = closure[cl];
975:         PetscInt       clval;

977:         if ((clp < pStart[d]) || (clp >= pEnd[d])) continue;
978:         DMLabelGetValue(label, clp, &clval);
979:         if (clval == -1) {
980:           const PetscInt *cone;
981:           PetscInt        coneSize, c;

983:           /* If a cone point touches the fault, then this point touches the fault */
984:           DMPlexGetCone(dm, clp, &cone);
985:           DMPlexGetConeSize(dm, clp, &coneSize);
986:           for (c = 0; c < coneSize; ++c) {
987:             PetscInt cpval;

989:             DMLabelGetValue(label, cone[c], &cpval);
990:             if (cpval != -1) {
991:               PetscInt dep;

993:               DMPlexGetPointDepth(dm, clp, &dep);
994:               clval = cval < 0 ? -(shift + dep) : shift + dep;
995:               DMLabelSetValue(label, clp, clval);
996:               break;
997:             }
998:           }
999:         }
1000:         /* Mark neighbor cells through marked faces (these cells must also touch the fault) */
1001:         if (d == dim - 1 && clval != -1) {
1002:           const PetscInt *support;
1003:           PetscInt        supportSize, s, nval;

1005:           DMPlexGetSupport(dm, clp, &support);
1006:           DMPlexGetSupportSize(dm, clp, &supportSize);
1007:           for (s = 0; s < supportSize; ++s) {
1008:             DMLabelGetValue(label, support[s], &nval);
1009:             if (nval == -1) {
1010:               DMLabelSetValue(label, support[s], clval < 0 ? clval - 1 : clval + 1);
1011:               DMPlexPointQueueEnqueue(queue, support[s]);
1012:             }
1013:           }
1014:         }
1015:       }
1016:     }
1017:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1018:   }
1019:   PetscFree2(pStart, pEnd);
1020:   return 0;
1021: }

1023: typedef struct {
1024:   DM               dm;
1025:   DMPlexPointQueue queue;
1026: } PointDivision;

1028: static PetscErrorCode divideCell(DMLabel label, PetscInt p, PetscInt val, void *ctx)
1029: {
1030:   PointDivision  *div  = (PointDivision *)ctx;
1031:   PetscInt        cval = val < 0 ? val - 1 : val + 1;
1032:   const PetscInt *support;
1033:   PetscInt        supportSize, s;

1035:   DMPlexGetSupport(div->dm, p, &support);
1036:   DMPlexGetSupportSize(div->dm, p, &supportSize);
1037:   for (s = 0; s < supportSize; ++s) {
1038:     DMLabelSetValue(label, support[s], cval);
1039:     DMPlexPointQueueEnqueue(div->queue, support[s]);
1040:   }
1041:   return 0;
1042: }

1044: /* Mark cells by label propagation */
1045: static PetscErrorCode DMPlexLabelFaultHalo(DM dm, DMLabel faultLabel)
1046: {
1047:   DMPlexPointQueue queue = NULL;
1048:   PointDivision    div;
1049:   PetscSF          pointSF;
1050:   IS               pointIS;
1051:   const PetscInt  *points;
1052:   PetscBool        empty;
1053:   PetscInt         dim, shift = 100, n, i;

1055:   DMGetDimension(dm, &dim);
1056:   DMPlexPointQueueCreate(1024, &queue);
1057:   div.dm    = dm;
1058:   div.queue = queue;
1059:   /* Enqueue cells on fault */
1060:   DMLabelGetStratumIS(faultLabel, shift + dim, &pointIS);
1061:   if (pointIS) {
1062:     ISGetLocalSize(pointIS, &n);
1063:     ISGetIndices(pointIS, &points);
1064:     for (i = 0; i < n; ++i) DMPlexPointQueueEnqueue(queue, points[i]);
1065:     ISRestoreIndices(pointIS, &points);
1066:     ISDestroy(&pointIS);
1067:   }
1068:   DMLabelGetStratumIS(faultLabel, -(shift + dim), &pointIS);
1069:   if (pointIS) {
1070:     ISGetLocalSize(pointIS, &n);
1071:     ISGetIndices(pointIS, &points);
1072:     for (i = 0; i < n; ++i) DMPlexPointQueueEnqueue(queue, points[i]);
1073:     ISRestoreIndices(pointIS, &points);
1074:     ISDestroy(&pointIS);
1075:   }

1077:   DMGetPointSF(dm, &pointSF);
1078:   DMLabelPropagateBegin(faultLabel, pointSF);
1079:   /* While edge queue is not empty: */
1080:   DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
1081:   while (!empty) {
1082:     DivideCells_Private(dm, faultLabel, queue);
1083:     DMLabelPropagatePush(faultLabel, pointSF, divideCell, &div);
1084:     DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
1085:   }
1086:   DMLabelPropagateEnd(faultLabel, pointSF);
1087:   DMPlexPointQueueDestroy(&queue);
1088:   return 0;
1089: }

1091: /*
1092:   We are adding three kinds of points here:
1093:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
1094:     Non-replicated: Points which exist on the fault, but are not replicated
1095:     Ghost:          These are shared fault faces which are not owned by this process. These do not produce hybrid cells and do not replicate
1096:     Hybrid:         Entirely new points, such as cohesive cells

1098:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
1099:   each depth so that the new split/hybrid points can be inserted as a block.
1100: */
1101: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
1102: {
1103:   MPI_Comm         comm;
1104:   IS               valueIS;
1105:   PetscInt         numSP = 0; /* The number of depths for which we have replicated points */
1106:   const PetscInt  *values;    /* List of depths for which we have replicated points */
1107:   IS              *splitIS;
1108:   IS              *unsplitIS;
1109:   IS               ghostIS;
1110:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
1111:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
1112:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
1113:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
1114:   PetscInt         numGhostPoints;     /* The number of unowned, shared fault faces */
1115:   const PetscInt **splitPoints;        /* Replicated points for each depth */
1116:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
1117:   const PetscInt  *ghostPoints;        /* Ghost fault faces */
1118:   PetscSection     coordSection;
1119:   Vec              coordinates;
1120:   PetscScalar     *coords;
1121:   PetscInt        *depthMax;   /* The first hybrid point at each depth in the original mesh */
1122:   PetscInt        *depthEnd;   /* The point limit at each depth in the original mesh */
1123:   PetscInt        *depthShift; /* Number of replicated+hybrid points at each depth */
1124:   PetscInt        *pMaxNew;    /* The first replicated point at each depth in the new mesh, hybrids come after this */
1125:   PetscInt        *coneNew, *coneONew, *supportNew;
1126:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;

1128:   PetscObjectGetComm((PetscObject)dm, &comm);
1129:   DMGetDimension(dm, &dim);
1130:   DMPlexGetDepth(dm, &depth);
1131:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1132:   /* We do not want this label automatically computed, instead we compute it here */
1133:   DMCreateLabel(sdm, "celltype");
1134:   /* Count split points and add cohesive cells */
1135:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
1136:   PetscMalloc5(depth + 1, &depthMax, depth + 1, &depthEnd, 2 * (depth + 1), &depthShift, depth + 1, &pMaxNew, depth + 1, &numHybridPointsOld);
1137:   PetscMalloc7(depth + 1, &splitIS, depth + 1, &unsplitIS, depth + 1, &numSplitPoints, depth + 1, &numUnsplitPoints, depth + 1, &numHybridPoints, depth + 1, &splitPoints, depth + 1, &unsplitPoints);
1138:   for (d = 0; d <= depth; ++d) {
1139:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
1140:     DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
1141:     depthEnd[d]           = pMaxNew[d];
1142:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1143:     numSplitPoints[d]     = 0;
1144:     numUnsplitPoints[d]   = 0;
1145:     numHybridPoints[d]    = 0;
1146:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1147:     splitPoints[d]        = NULL;
1148:     unsplitPoints[d]      = NULL;
1149:     splitIS[d]            = NULL;
1150:     unsplitIS[d]          = NULL;
1151:     /* we are shifting the existing hybrid points with the stratum behind them, so
1152:      * the split comes at the end of the normal points, i.e., at depthMax[d] */
1153:     depthShift[2 * d]     = depthMax[d];
1154:     depthShift[2 * d + 1] = 0;
1155:   }
1156:   numGhostPoints = 0;
1157:   ghostPoints    = NULL;
1158:   if (label) {
1159:     DMLabelGetValueIS(label, &valueIS);
1160:     ISGetLocalSize(valueIS, &numSP);
1161:     ISGetIndices(valueIS, &values);
1162:   }
1163:   for (sp = 0; sp < numSP; ++sp) {
1164:     const PetscInt dep = values[sp];

1166:     if ((dep < 0) || (dep > depth)) continue;
1167:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1168:     if (splitIS[dep]) {
1169:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1170:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
1171:     }
1172:     DMLabelGetStratumIS(label, shift2 + dep, &unsplitIS[dep]);
1173:     if (unsplitIS[dep]) {
1174:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1175:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1176:     }
1177:   }
1178:   DMLabelGetStratumIS(label, shift2 + dim - 1, &ghostIS);
1179:   if (ghostIS) {
1180:     ISGetLocalSize(ghostIS, &numGhostPoints);
1181:     ISGetIndices(ghostIS, &ghostPoints);
1182:   }
1183:   /* Calculate number of hybrid points */
1184:   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   */
1185:   for (d = 0; d <= depth; ++d) depthShift[2 * d + 1] = numSplitPoints[d] + numHybridPoints[d];
1186:   DMPlexShiftPointSetUp_Internal(depth, depthShift);
1187:   /* the end of the points in this stratum that come before the new points:
1188:    * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1189:    * added points */
1190:   for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d], depth, depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1191:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1192:   /* Step 3: Set cone/support sizes for new points */
1193:   for (dep = 0; dep <= depth; ++dep) {
1194:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1195:       const PetscInt  oldp   = splitPoints[dep][p];
1196:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1197:       const PetscInt  splitp = p + pMaxNew[dep];
1198:       const PetscInt *support;
1199:       DMPolytopeType  ct;
1200:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

1202:       DMPlexGetConeSize(dm, oldp, &coneSize);
1203:       DMPlexSetConeSize(sdm, splitp, coneSize);
1204:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1205:       DMPlexSetSupportSize(sdm, splitp, supportSize);
1206:       DMPlexGetCellType(dm, oldp, &ct);
1207:       DMPlexSetCellType(sdm, splitp, ct);
1208:       if (dep == depth - 1) {
1209:         const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];

1211:         /* Add cohesive cells, they are prisms */
1212:         DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
1213:         switch (coneSize) {
1214:         case 2:
1215:           DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);
1216:           break;
1217:         case 3:
1218:           DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);
1219:           break;
1220:         case 4:
1221:           DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);
1222:           break;
1223:         }
1224:         /* Shared fault faces with only one support cell now have two with the cohesive cell */
1225:         /*   TODO Check thaat oldp has rootdegree == 1 */
1226:         if (supportSize == 1) {
1227:           const PetscInt *support;
1228:           PetscInt        val;

1230:           DMPlexGetSupport(dm, oldp, &support);
1231:           DMLabelGetValue(label, support[0], &val);
1232:           if (val < 0) DMPlexSetSupportSize(sdm, splitp, 2);
1233:           else DMPlexSetSupportSize(sdm, newp, 2);
1234:         }
1235:       } else if (dep == 0) {
1236:         const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];

1238:         DMPlexGetSupport(dm, oldp, &support);
1239:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1240:           PetscInt val;

1242:           DMLabelGetValue(label, support[e], &val);
1243:           if (val == 1) ++qf;
1244:           if ((val == 1) || (val == (shift + 1))) ++qn;
1245:           if ((val == 1) || (val == -(shift + 1))) ++qp;
1246:         }
1247:         /* Split old vertex: Edges into original vertex and new cohesive edge */
1248:         DMPlexSetSupportSize(sdm, newp, qn + 1);
1249:         /* Split new vertex: Edges into split vertex and new cohesive edge */
1250:         DMPlexSetSupportSize(sdm, splitp, qp + 1);
1251:         /* Add hybrid edge */
1252:         DMPlexSetConeSize(sdm, hybedge, 2);
1253:         DMPlexSetSupportSize(sdm, hybedge, qf);
1254:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1255:       } else if (dep == dim - 2) {
1256:         const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];

1258:         DMPlexGetSupport(dm, oldp, &support);
1259:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1260:           PetscInt val;

1262:           DMLabelGetValue(label, support[e], &val);
1263:           if (val == dim - 1) ++qf;
1264:           if ((val == dim - 1) || (val == (shift + dim - 1))) ++qn;
1265:           if ((val == dim - 1) || (val == -(shift + dim - 1))) ++qp;
1266:         }
1267:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1268:         DMPlexSetSupportSize(sdm, newp, qn + 1);
1269:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1270:         DMPlexSetSupportSize(sdm, splitp, qp + 1);
1271:         /* Add hybrid face */
1272:         DMPlexSetConeSize(sdm, hybface, 4);
1273:         DMPlexSetSupportSize(sdm, hybface, qf);
1274:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1275:       }
1276:     }
1277:   }
1278:   for (dep = 0; dep <= depth; ++dep) {
1279:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1280:       const PetscInt  oldp = unsplitPoints[dep][p];
1281:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1282:       const PetscInt *support;
1283:       PetscInt        coneSize, supportSize, qf, e, s;

1285:       DMPlexGetConeSize(dm, oldp, &coneSize);
1286:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1287:       DMPlexGetSupport(dm, oldp, &support);
1288:       if (dep == 0) {
1289:         const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];

1291:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1292:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1293:           PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e);
1294:           if (e >= 0) ++qf;
1295:         }
1296:         DMPlexSetSupportSize(sdm, newp, qf + 2);
1297:         /* Add hybrid edge */
1298:         DMPlexSetConeSize(sdm, hybedge, 2);
1299:         for (e = 0, qf = 0; e < supportSize; ++e) {
1300:           PetscInt val;

1302:           DMLabelGetValue(label, support[e], &val);
1303:           /* Split and unsplit edges produce hybrid faces */
1304:           if (val == 1) ++qf;
1305:           if (val == (shift2 + 1)) ++qf;
1306:         }
1307:         DMPlexSetSupportSize(sdm, hybedge, qf);
1308:         DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1309:       } else if (dep == dim - 2) {
1310:         const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1311:         PetscInt       val;

1313:         for (e = 0, qf = 0; e < supportSize; ++e) {
1314:           DMLabelGetValue(label, support[e], &val);
1315:           if (val == dim - 1) qf += 2;
1316:           else ++qf;
1317:         }
1318:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1319:         DMPlexSetSupportSize(sdm, newp, qf + 2);
1320:         /* Add hybrid face */
1321:         for (e = 0, qf = 0; e < supportSize; ++e) {
1322:           DMLabelGetValue(label, support[e], &val);
1323:           if (val == dim - 1) ++qf;
1324:         }
1325:         DMPlexSetConeSize(sdm, hybface, 4);
1326:         DMPlexSetSupportSize(sdm, hybface, qf);
1327:         DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1328:       }
1329:     }
1330:   }
1331:   /* Step 4: Setup split DM */
1332:   DMSetUp(sdm);
1333:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1334:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1335:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneNew, PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneONew, PetscMax(maxSupportSize, maxSupportSizeNew), &supportNew);
1336:   /* Step 6: Set cones and supports for new points */
1337:   for (dep = 0; dep <= depth; ++dep) {
1338:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1339:       const PetscInt  oldp   = splitPoints[dep][p];
1340:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1341:       const PetscInt  splitp = p + pMaxNew[dep];
1342:       const PetscInt *cone, *support, *ornt;
1343:       DMPolytopeType  ct;
1344:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

1346:       DMPlexGetCellType(dm, oldp, &ct);
1347:       DMPlexGetConeSize(dm, oldp, &coneSize);
1348:       DMPlexGetCone(dm, oldp, &cone);
1349:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1350:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1351:       DMPlexGetSupport(dm, oldp, &support);
1352:       if (dep == depth - 1) {
1353:         PetscBool       hasUnsplit = PETSC_FALSE;
1354:         const PetscInt  hybcell    = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1355:         const PetscInt *supportF;

1357:         coneONew[0] = coneONew[1] = -1000;
1358:         /* Split face:       copy in old face to new face to start */
1359:         DMPlexGetSupport(sdm, newp, &supportF);
1360:         DMPlexSetSupport(sdm, splitp, supportF);
1361:         /* Split old face:   old vertices/edges in cone so no change */
1362:         /* Split new face:   new vertices/edges in cone */
1363:         for (q = 0; q < coneSize; ++q) {
1364:           PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1365:           if (v < 0) {
1366:             PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1368:             coneNew[2 + q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1369:             hasUnsplit     = PETSC_TRUE;
1370:           } else {
1371:             coneNew[2 + q] = v + pMaxNew[dep - 1];
1372:             if (dep > 1) {
1373:               const PetscInt *econe;
1374:               PetscInt        econeSize, r, vs, vu;

1376:               DMPlexGetConeSize(dm, cone[q], &econeSize);
1377:               DMPlexGetCone(dm, cone[q], &econe);
1378:               for (r = 0; r < econeSize; ++r) {
1379:                 PetscFindInt(econe[r], numSplitPoints[dep - 2], splitPoints[dep - 2], &vs);
1380:                 PetscFindInt(econe[r], numUnsplitPoints[dep - 2], unsplitPoints[dep - 2], &vu);
1381:                 if (vs >= 0) continue;
1383:                 hasUnsplit = PETSC_TRUE;
1384:               }
1385:             }
1386:           }
1387:         }
1388:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
1389:         DMPlexSetConeOrientation(sdm, splitp, ornt);
1390:         /* Face support */
1391:         PetscInt vals[2];

1393:         DMLabelGetValue(label, support[0], &vals[0]);
1394:         if (supportSize > 1) DMLabelGetValue(label, support[1], &vals[1]);
1395:         else vals[1] = -vals[0];

1398:         for (s = 0; s < 2; ++s) {
1399:           if (s >= supportSize) {
1400:             if (vals[s] < 0) {
1401:               /* Ghost old face:   Replace negative side cell with cohesive cell */
1402:               DMPlexInsertSupport(sdm, newp, s, hybcell);
1403:             } else {
1404:               /* Ghost new face:   Replace positive side cell with cohesive cell */
1405:               DMPlexInsertSupport(sdm, splitp, s, hybcell);
1406:             }
1407:           } else {
1408:             if (vals[s] < 0) {
1409:               /* Split old face:   Replace negative side cell with cohesive cell */
1410:               DMPlexInsertSupport(sdm, newp, s, hybcell);
1411:             } else {
1412:               /* Split new face:   Replace positive side cell with cohesive cell */
1413:               DMPlexInsertSupport(sdm, splitp, s, hybcell);
1414:             }
1415:           }
1416:         }
1417:         /* Get orientation for cohesive face using the positive side cell */
1418:         {
1419:           const PetscInt *ncone, *nconeO;
1420:           PetscInt        nconeSize, nc, ocell;
1421:           PetscBool       flip = PETSC_FALSE;

1423:           if (supportSize > 1) {
1424:             ocell = vals[0] < 0 ? support[1] : support[0];
1425:           } else {
1426:             ocell = support[0];
1427:             flip  = vals[0] < 0 ? PETSC_TRUE : PETSC_FALSE;
1428:           }
1429:           DMPlexGetConeSize(dm, ocell, &nconeSize);
1430:           DMPlexGetCone(dm, ocell, &ncone);
1431:           DMPlexGetConeOrientation(dm, ocell, &nconeO);
1432:           for (nc = 0; nc < nconeSize; ++nc) {
1433:             if (ncone[nc] == oldp) {
1434:               coneONew[0] = flip ? -(nconeO[nc] + 1) : nconeO[nc];
1435:               break;
1436:             }
1437:           }
1439:         }
1440:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
1441:         {
1442:           const PetscInt No = DMPolytopeTypeGetNumArrangments(ct) / 2;
1444:         }
1445:         const PetscInt *arr = DMPolytopeTypeGetArrangment(ct, coneONew[0]);

1447:         coneNew[0]  = newp; /* Extracted negative side orientation above */
1448:         coneNew[1]  = splitp;
1449:         coneONew[1] = coneONew[0];
1450:         for (q = 0; q < coneSize; ++q) {
1451:           /* Hybrid faces must follow order from oriented end face */
1452:           const PetscInt qa = arr[q * 2 + 0];
1453:           const PetscInt qo = arr[q * 2 + 1];
1454:           DMPolytopeType ft = dep == 2 ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;

1456:           PetscFindInt(cone[qa], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1457:           if (v < 0) {
1458:             PetscFindInt(cone[qa], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1459:             coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1460:           } else {
1461:             coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep];
1462:           }
1463:           coneONew[2 + q] = DMPolytopeTypeComposeOrientation(ft, qo, ornt[qa]);
1464:         }
1465:         DMPlexSetCone(sdm, hybcell, coneNew);
1466:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1467:         /* Label the hybrid cells on the boundary of the split */
1468:         if (hasUnsplit) DMLabelSetValue(label, -hybcell, dim);
1469:       } else if (dep == 0) {
1470:         const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];

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

1476:           DMLabelGetValue(label, support[e], &val);
1477:           if ((val == 1) || (val == (shift + 1))) supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1478:         }
1479:         supportNew[qn] = hybedge;
1480:         DMPlexSetSupport(sdm, newp, supportNew);
1481:         /* Split new vertex: Edges in new split faces and new cohesive edge */
1482:         for (e = 0, qp = 0; e < supportSize; ++e) {
1483:           PetscInt val, edge;

1485:           DMLabelGetValue(label, support[e], &val);
1486:           if (val == 1) {
1487:             PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1489:             supportNew[qp++] = edge + pMaxNew[dep + 1];
1490:           } else if (val == -(shift + 1)) {
1491:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1492:           }
1493:         }
1494:         supportNew[qp] = hybedge;
1495:         DMPlexSetSupport(sdm, splitp, supportNew);
1496:         /* Hybrid edge:    Old and new split vertex */
1497:         coneNew[0] = newp;
1498:         coneNew[1] = splitp;
1499:         DMPlexSetCone(sdm, hybedge, coneNew);
1500:         for (e = 0, qf = 0; e < supportSize; ++e) {
1501:           PetscInt val, edge;

1503:           DMLabelGetValue(label, support[e], &val);
1504:           if (val == 1) {
1505:             PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1507:             supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1508:           }
1509:         }
1510:         DMPlexSetSupport(sdm, hybedge, supportNew);
1511:       } else if (dep == dim - 2) {
1512:         const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];

1514:         /* Split old edge:   old vertices in cone so no change */
1515:         /* Split new edge:   new vertices in cone */
1516:         for (q = 0; q < coneSize; ++q) {
1517:           PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1518:           if (v < 0) {
1519:             PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1521:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1522:           } else {
1523:             coneNew[q] = v + pMaxNew[dep - 1];
1524:           }
1525:         }
1526:         DMPlexSetCone(sdm, splitp, coneNew);
1527:         /* Split old edge: Faces in positive side cells and old split faces */
1528:         for (e = 0, q = 0; e < supportSize; ++e) {
1529:           PetscInt val;

1531:           DMLabelGetValue(label, support[e], &val);
1532:           if (val == dim - 1) {
1533:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1534:           } else if (val == (shift + dim - 1)) {
1535:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1536:           }
1537:         }
1538:         supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1539:         DMPlexSetSupport(sdm, newp, supportNew);
1540:         /* Split new edge: Faces in negative side cells and new split faces */
1541:         for (e = 0, q = 0; e < supportSize; ++e) {
1542:           PetscInt val, face;

1544:           DMLabelGetValue(label, support[e], &val);
1545:           if (val == dim - 1) {
1546:             PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1548:             supportNew[q++] = face + pMaxNew[dep + 1];
1549:           } else if (val == -(shift + dim - 1)) {
1550:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1551:           }
1552:         }
1553:         supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1554:         DMPlexSetSupport(sdm, splitp, supportNew);
1555:         /* Hybrid face */
1556:         coneNew[0] = newp;
1557:         coneNew[1] = splitp;
1558:         for (v = 0; v < coneSize; ++v) {
1559:           PetscInt vertex;
1560:           PetscFindInt(cone[v], numSplitPoints[dep - 1], splitPoints[dep - 1], &vertex);
1561:           if (vertex < 0) {
1562:             PetscFindInt(cone[v], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &vertex);
1564:             coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1565:           } else {
1566:             coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1567:           }
1568:         }
1569:         DMPlexSetCone(sdm, hybface, coneNew);
1570:         for (e = 0, qf = 0; e < supportSize; ++e) {
1571:           PetscInt val, face;

1573:           DMLabelGetValue(label, support[e], &val);
1574:           if (val == dim - 1) {
1575:             PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1577:             supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1578:           }
1579:         }
1580:         DMPlexSetSupport(sdm, hybface, supportNew);
1581:       }
1582:     }
1583:   }
1584:   for (dep = 0; dep <= depth; ++dep) {
1585:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1586:       const PetscInt  oldp = unsplitPoints[dep][p];
1587:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1588:       const PetscInt *cone, *support;
1589:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1591:       DMPlexGetConeSize(dm, oldp, &coneSize);
1592:       DMPlexGetCone(dm, oldp, &cone);
1593:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1594:       DMPlexGetSupport(dm, oldp, &support);
1595:       if (dep == 0) {
1596:         const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];

1598:         /* Unsplit vertex */
1599:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1600:         for (s = 0, q = 0; s < supportSize; ++s) {
1601:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1602:           PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e);
1603:           if (e >= 0) supportNew[q++] = e + pMaxNew[dep + 1];
1604:         }
1605:         supportNew[q++] = hybedge;
1606:         supportNew[q++] = hybedge;
1608:         DMPlexSetSupport(sdm, newp, supportNew);
1609:         /* Hybrid edge */
1610:         coneNew[0] = newp;
1611:         coneNew[1] = newp;
1612:         DMPlexSetCone(sdm, hybedge, coneNew);
1613:         for (e = 0, qf = 0; e < supportSize; ++e) {
1614:           PetscInt val, edge;

1616:           DMLabelGetValue(label, support[e], &val);
1617:           if (val == 1) {
1618:             PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1620:             supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1621:           } else if (val == (shift2 + 1)) {
1622:             PetscFindInt(support[e], numUnsplitPoints[dep + 1], unsplitPoints[dep + 1], &edge);
1624:             supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2] + numSplitPoints[dep + 1];
1625:           }
1626:         }
1627:         DMPlexSetSupport(sdm, hybedge, supportNew);
1628:       } else if (dep == dim - 2) {
1629:         const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];

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

1635:           DMLabelGetValue(label, support[f], &val);
1636:           if (val == dim - 1) {
1637:             PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1639:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1640:             supportNew[qf++] = face + pMaxNew[dep + 1];
1641:           } else {
1642:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1643:           }
1644:         }
1645:         supportNew[qf++] = hybface;
1646:         supportNew[qf++] = hybface;
1647:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1649:         DMPlexSetSupport(sdm, newp, supportNew);
1650:         /* Add hybrid face */
1651:         coneNew[0] = newp;
1652:         coneNew[1] = newp;
1653:         PetscFindInt(cone[0], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1655:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1656:         PetscFindInt(cone[1], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1658:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1659:         DMPlexSetCone(sdm, hybface, coneNew);
1660:         for (f = 0, qf = 0; f < supportSize; ++f) {
1661:           PetscInt val, face;

1663:           DMLabelGetValue(label, support[f], &val);
1664:           if (val == dim - 1) {
1665:             PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1666:             supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1667:           }
1668:         }
1669:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1671:         DMPlexSetSupport(sdm, hybface, supportNew);
1672:       }
1673:     }
1674:   }
1675:   /* Step 6b: Replace split points in negative side cones */
1676:   for (sp = 0; sp < numSP; ++sp) {
1677:     PetscInt        dep = values[sp];
1678:     IS              pIS;
1679:     PetscInt        numPoints;
1680:     const PetscInt *points;

1682:     if (dep >= 0) continue;
1683:     DMLabelGetStratumIS(label, dep, &pIS);
1684:     if (!pIS) continue;
1685:     dep = -dep - shift;
1686:     ISGetLocalSize(pIS, &numPoints);
1687:     ISGetIndices(pIS, &points);
1688:     for (p = 0; p < numPoints; ++p) {
1689:       const PetscInt  oldp = points[p];
1690:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1691:       const PetscInt *cone;
1692:       PetscInt        coneSize, c;
1693:       /* PetscBool       replaced = PETSC_FALSE; */

1695:       /* Negative edge: replace split vertex */
1696:       /* Negative cell: replace split face */
1697:       DMPlexGetConeSize(sdm, newp, &coneSize);
1698:       DMPlexGetCone(sdm, newp, &cone);
1699:       for (c = 0; c < coneSize; ++c) {
1700:         const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c], depth, depthShift);
1701:         PetscInt       csplitp, cp, val;

1703:         DMLabelGetValue(label, coldp, &val);
1704:         if (val == dep - 1) {
1705:           PetscFindInt(coldp, numSplitPoints[dep - 1], splitPoints[dep - 1], &cp);
1707:           csplitp = pMaxNew[dep - 1] + cp;
1708:           DMPlexInsertCone(sdm, newp, c, csplitp);
1709:           /* replaced = PETSC_TRUE; */
1710:         }
1711:       }
1712:       /* Cells with only a vertex or edge on the submesh have no replacement */
1714:     }
1715:     ISRestoreIndices(pIS, &points);
1716:     ISDestroy(&pIS);
1717:   }
1718:   /* Step 7: Coordinates */
1719:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1720:   DMGetCoordinateSection(sdm, &coordSection);
1721:   DMGetCoordinatesLocal(sdm, &coordinates);
1722:   VecGetArray(coordinates, &coords);
1723:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1724:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1725:     const PetscInt splitp = pMaxNew[0] + v;
1726:     PetscInt       dof, off, soff, d;

1728:     PetscSectionGetDof(coordSection, newp, &dof);
1729:     PetscSectionGetOffset(coordSection, newp, &off);
1730:     PetscSectionGetOffset(coordSection, splitp, &soff);
1731:     for (d = 0; d < dof; ++d) coords[soff + d] = coords[off + d];
1732:   }
1733:   VecRestoreArray(coordinates, &coords);
1734:   /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1735:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1736:   /*   TODO We need to associate the ghost points with the correct replica */
1737:   /* Step 9: Labels */
1738:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1739:   DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm);
1740:   DMGetNumLabels(sdm, &numLabels);
1741:   for (dep = 0; dep <= depth; ++dep) {
1742:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1743:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1744:       const PetscInt splitp = pMaxNew[dep] + p;
1745:       PetscInt       l;

1747:       if (splitLabel) {
1748:         const PetscInt val = 100 + dep;

1750:         DMLabelSetValue(splitLabel, newp, val);
1751:         DMLabelSetValue(splitLabel, splitp, -val);
1752:       }
1753:       for (l = 0; l < numLabels; ++l) {
1754:         DMLabel     mlabel;
1755:         const char *lname;
1756:         PetscInt    val;
1757:         PetscBool   isDepth;

1759:         DMGetLabelName(sdm, l, &lname);
1760:         PetscStrcmp(lname, "depth", &isDepth);
1761:         if (isDepth) continue;
1762:         DMGetLabel(sdm, lname, &mlabel);
1763:         DMLabelGetValue(mlabel, newp, &val);
1764:         if (val >= 0) DMLabelSetValue(mlabel, splitp, val);
1765:       }
1766:     }
1767:   }
1768:   for (sp = 0; sp < numSP; ++sp) {
1769:     const PetscInt dep = values[sp];

1771:     if ((dep < 0) || (dep > depth)) continue;
1772:     if (splitIS[dep]) ISRestoreIndices(splitIS[dep], &splitPoints[dep]);
1773:     ISDestroy(&splitIS[dep]);
1774:     if (unsplitIS[dep]) ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);
1775:     ISDestroy(&unsplitIS[dep]);
1776:   }
1777:   if (ghostIS) ISRestoreIndices(ghostIS, &ghostPoints);
1778:   ISDestroy(&ghostIS);
1779:   if (label) {
1780:     ISRestoreIndices(valueIS, &values);
1781:     ISDestroy(&valueIS);
1782:   }
1783:   for (d = 0; d <= depth; ++d) {
1784:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1785:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1786:   }
1787:   PetscFree3(coneNew, coneONew, supportNew);
1788:   PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1789:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1790:   return 0;
1791: }

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

1796:   Collective on dm

1798:   Input Parameters:
1799: + dm - The original DM
1800: - label - The label specifying the boundary faces (this could be auto-generated)

1802:   Output Parameters:
1803: + splitLabel - The label containing the split points, or NULL if no output is desired
1804: - dmSplit - The new DM

1806:   Level: developer

1808: .seealso: `DMCreate()`, `DMPlexLabelCohesiveComplete()`
1809: @*/
1810: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1811: {
1812:   DM       sdm;
1813:   PetscInt dim;

1817:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1818:   DMSetType(sdm, DMPLEX);
1819:   DMGetDimension(dm, &dim);
1820:   DMSetDimension(sdm, dim);
1821:   switch (dim) {
1822:   case 2:
1823:   case 3:
1824:     DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1825:     break;
1826:   default:
1827:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %" PetscInt_FMT, dim);
1828:   }
1829:   *dmSplit = sdm;
1830:   return 0;
1831: }

1833: /* Returns the side of the surface for a given cell with a face on the surface */
1834: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1835: {
1836:   const PetscInt *cone, *ornt;
1837:   PetscInt        dim, coneSize, c;

1839:   *pos = PETSC_TRUE;
1840:   DMGetDimension(dm, &dim);
1841:   DMPlexGetConeSize(dm, cell, &coneSize);
1842:   DMPlexGetCone(dm, cell, &cone);
1843:   DMPlexGetConeOrientation(dm, cell, &ornt);
1844:   for (c = 0; c < coneSize; ++c) {
1845:     if (cone[c] == face) {
1846:       PetscInt o = ornt[c];

1848:       if (subdm) {
1849:         const PetscInt *subcone, *subornt;
1850:         PetscInt        subpoint, subface, subconeSize, sc;

1852:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1853:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1854:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1855:         DMPlexGetCone(subdm, subpoint, &subcone);
1856:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1857:         for (sc = 0; sc < subconeSize; ++sc) {
1858:           if (subcone[sc] == subface) {
1859:             o = subornt[0];
1860:             break;
1861:           }
1862:         }
1864:       }
1865:       if (o >= 0) *pos = PETSC_TRUE;
1866:       else *pos = PETSC_FALSE;
1867:       break;
1868:     }
1869:   }
1871:   return 0;
1872: }

1874: static PetscErrorCode CheckFaultEdge_Private(DM dm, DMLabel label)
1875: {
1876:   IS              facePosIS, faceNegIS, dimIS;
1877:   const PetscInt *points;
1878:   PetscInt        dim, numPoints, p, shift = 100, shift2 = 200;

1880:   DMGetDimension(dm, &dim);
1881:   /* If any faces touching the fault divide cells on either side, split them */
1882:   DMLabelGetStratumIS(label, shift + dim - 1, &facePosIS);
1883:   DMLabelGetStratumIS(label, -(shift + dim - 1), &faceNegIS);
1884:   if (!facePosIS || !faceNegIS) {
1885:     ISDestroy(&facePosIS);
1886:     ISDestroy(&faceNegIS);
1887:     return 0;
1888:   }
1889:   ISExpand(facePosIS, faceNegIS, &dimIS);
1890:   ISDestroy(&facePosIS);
1891:   ISDestroy(&faceNegIS);
1892:   ISGetLocalSize(dimIS, &numPoints);
1893:   ISGetIndices(dimIS, &points);
1894:   for (p = 0; p < numPoints; ++p) {
1895:     const PetscInt  point = points[p];
1896:     const PetscInt *support;
1897:     PetscInt        supportSize, valA, valB;

1899:     DMPlexGetSupportSize(dm, point, &supportSize);
1900:     if (supportSize != 2) continue;
1901:     DMPlexGetSupport(dm, point, &support);
1902:     DMLabelGetValue(label, support[0], &valA);
1903:     DMLabelGetValue(label, support[1], &valB);
1904:     if ((valA == -1) || (valB == -1)) continue;
1905:     if (valA * valB > 0) continue;
1906:     /* Check that this face is not incident on only unsplit faces, meaning has at least one split face */
1907:     {
1908:       PetscInt *closure = NULL;
1909:       PetscBool split   = PETSC_FALSE;
1910:       PetscInt  closureSize, cl;

1912:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1913:       for (cl = 0; cl < closureSize * 2; cl += 2) {
1914:         DMLabelGetValue(label, closure[cl], &valA);
1915:         if ((valA >= 0) && (valA <= dim)) {
1916:           split = PETSC_TRUE;
1917:           break;
1918:         }
1919:       }
1920:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1921:       if (!split) continue;
1922:     }
1923:     /* Split the face */
1924:     DMLabelGetValue(label, point, &valA);
1925:     DMLabelClearValue(label, point, valA);
1926:     DMLabelSetValue(label, point, dim - 1);
1927:     /* Label its closure:
1928:       unmarked: label as unsplit
1929:       incident: relabel as split
1930:       split:    do nothing
1931:     */
1932:     {
1933:       PetscInt *closure = NULL;
1934:       PetscInt  closureSize, cl, dep;

1936:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1937:       for (cl = 0; cl < closureSize * 2; cl += 2) {
1938:         DMLabelGetValue(label, closure[cl], &valA);
1939:         if (valA == -1) { /* Mark as unsplit */
1940:           DMPlexGetPointDepth(dm, closure[cl], &dep);
1941:           DMLabelSetValue(label, closure[cl], shift2 + dep);
1942:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1943:           DMPlexGetPointDepth(dm, closure[cl], &dep);
1944:           DMLabelClearValue(label, closure[cl], valA);
1945:           DMLabelSetValue(label, closure[cl], dep);
1946:         }
1947:       }
1948:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1949:     }
1950:   }
1951:   ISRestoreIndices(dimIS, &points);
1952:   ISDestroy(&dimIS);
1953:   return 0;
1954: }

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

1960:   Input Parameters:
1961: + dm     - The DM
1962: . label  - A DMLabel marking the surface
1963: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1964: . bvalue - Value of DMLabel marking the vertices on the boundary
1965: . flip   - Flag to flip the submesh normal and replace points on the other side
1966: - subdm  - The subDM associated with the label, or NULL

1968:   Output Parameter:
1969: . label - A DMLabel marking all surface points

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

1973:   Level: developer

1975: .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexLabelComplete()`
1976: @*/
1977: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscInt bvalue, PetscBool flip, DM subdm)
1978: {
1979:   DMLabel         depthLabel;
1980:   IS              dimIS, subpointIS = NULL;
1981:   const PetscInt *points, *subpoints;
1982:   const PetscInt  rev   = flip ? -1 : 1;
1983:   PetscInt        shift = 100, shift2 = 200, shift3 = 300, dim, depth, numPoints, numSubpoints, p, val;

1985:   DMPlexGetDepth(dm, &depth);
1986:   DMGetDimension(dm, &dim);
1987:   DMPlexGetDepthLabel(dm, &depthLabel);
1988:   if (subdm) {
1989:     DMPlexGetSubpointIS(subdm, &subpointIS);
1990:     if (subpointIS) {
1991:       ISGetLocalSize(subpointIS, &numSubpoints);
1992:       ISGetIndices(subpointIS, &subpoints);
1993:     }
1994:   }
1995:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1996:   DMLabelGetStratumIS(label, dim - 1, &dimIS);
1997:   if (!dimIS) goto divide;
1998:   ISGetLocalSize(dimIS, &numPoints);
1999:   ISGetIndices(dimIS, &points);
2000:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
2001:     const PetscInt *support;
2002:     PetscInt        supportSize, s;

2004:     DMPlexGetSupportSize(dm, points[p], &supportSize);
2005: #if 0
2006:     if (supportSize != 2) {
2007:       const PetscInt *lp;
2008:       PetscInt        Nlp, pind;

2010:       /* Check that for a cell with a single support face, that face is in the SF */
2011:       /*   THis check only works for the remote side. We would need root side information */
2012:       PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
2013:       PetscFindInt(points[p], Nlp, lp, &pind);
2015:     }
2016: #endif
2017:     DMPlexGetSupport(dm, points[p], &support);
2018:     for (s = 0; s < supportSize; ++s) {
2019:       const PetscInt *cone;
2020:       PetscInt        coneSize, c;
2021:       PetscBool       pos;

2023:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
2024:       if (pos) DMLabelSetValue(label, support[s], rev * (shift + dim));
2025:       else DMLabelSetValue(label, support[s], -rev * (shift + dim));
2026:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
2027:       /* Put faces touching the fault in the label */
2028:       DMPlexGetConeSize(dm, support[s], &coneSize);
2029:       DMPlexGetCone(dm, support[s], &cone);
2030:       for (c = 0; c < coneSize; ++c) {
2031:         const PetscInt point = cone[c];

2033:         DMLabelGetValue(label, point, &val);
2034:         if (val == -1) {
2035:           PetscInt *closure = NULL;
2036:           PetscInt  closureSize, cl;

2038:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2039:           for (cl = 0; cl < closureSize * 2; cl += 2) {
2040:             const PetscInt clp  = closure[cl];
2041:             PetscInt       bval = -1;

2043:             DMLabelGetValue(label, clp, &val);
2044:             if (blabel) DMLabelGetValue(blabel, clp, &bval);
2045:             if ((val >= 0) && (val < dim - 1) && (bval < 0)) {
2046:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift + dim - 1 : -(shift + dim - 1));
2047:               break;
2048:             }
2049:           }
2050:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2051:         }
2052:       }
2053:     }
2054:   }
2055:   ISRestoreIndices(dimIS, &points);
2056:   ISDestroy(&dimIS);
2057:   /* Mark boundary points as unsplit */
2058:   if (blabel) {
2059:     IS bdIS;

2061:     DMLabelGetStratumIS(blabel, bvalue, &bdIS);
2062:     ISGetLocalSize(bdIS, &numPoints);
2063:     ISGetIndices(bdIS, &points);
2064:     for (p = 0; p < numPoints; ++p) {
2065:       const PetscInt point = points[p];
2066:       PetscInt       val, bval;

2068:       DMLabelGetValue(blabel, point, &bval);
2069:       if (bval >= 0) {
2070:         DMLabelGetValue(label, point, &val);
2071:         if ((val < 0) || (val > dim)) {
2072:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
2073:           DMLabelClearValue(blabel, point, bval);
2074:         }
2075:       }
2076:     }
2077:     for (p = 0; p < numPoints; ++p) {
2078:       const PetscInt point = points[p];
2079:       PetscInt       val, bval;

2081:       DMLabelGetValue(blabel, point, &bval);
2082:       if (bval >= 0) {
2083:         const PetscInt *cone, *support;
2084:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

2086:         /* Mark as unsplit */
2087:         DMLabelGetValue(label, point, &val);
2089:         DMLabelClearValue(label, point, val);
2090:         DMLabelSetValue(label, point, shift2 + val);
2091:         /* Check for cross-edge
2092:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
2093:         if (val != 0) continue;
2094:         DMPlexGetSupport(dm, point, &support);
2095:         DMPlexGetSupportSize(dm, point, &supportSize);
2096:         for (s = 0; s < supportSize; ++s) {
2097:           DMPlexGetCone(dm, support[s], &cone);
2098:           DMPlexGetConeSize(dm, support[s], &coneSize);
2100:           DMLabelGetValue(blabel, cone[0], &valA);
2101:           DMLabelGetValue(blabel, cone[1], &valB);
2102:           DMLabelGetValue(blabel, support[s], &valE);
2103:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) DMLabelSetValue(blabel, support[s], 2);
2104:         }
2105:       }
2106:     }
2107:     ISRestoreIndices(bdIS, &points);
2108:     ISDestroy(&bdIS);
2109:   }
2110:   /* Mark ghost fault cells */
2111:   {
2112:     PetscSF         sf;
2113:     const PetscInt *leaves;
2114:     PetscInt        Nl, l;

2116:     DMGetPointSF(dm, &sf);
2117:     PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL);
2118:     DMLabelGetStratumIS(label, dim - 1, &dimIS);
2119:     if (!dimIS) goto divide;
2120:     ISGetLocalSize(dimIS, &numPoints);
2121:     ISGetIndices(dimIS, &points);
2122:     if (Nl > 0) {
2123:       for (p = 0; p < numPoints; ++p) {
2124:         const PetscInt point = points[p];
2125:         PetscInt       val;

2127:         PetscFindInt(point, Nl, leaves, &l);
2128:         if (l >= 0) {
2129:           PetscInt *closure = NULL;
2130:           PetscInt  closureSize, cl;

2132:           DMLabelGetValue(label, point, &val);
2134:           DMLabelClearValue(label, point, val);
2135:           DMLabelSetValue(label, point, shift3 + val);
2136:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2137:           for (cl = 2; cl < closureSize * 2; cl += 2) {
2138:             const PetscInt clp = closure[cl];

2140:             DMLabelGetValue(label, clp, &val);
2142:             DMLabelClearValue(label, clp, val);
2143:             DMLabelSetValue(label, clp, shift3 + val);
2144:           }
2145:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2146:         }
2147:       }
2148:     }
2149:     ISRestoreIndices(dimIS, &points);
2150:     ISDestroy(&dimIS);
2151:   }
2152: divide:
2153:   if (subpointIS) ISRestoreIndices(subpointIS, &subpoints);
2154:   DMPlexLabelFaultHalo(dm, label);
2155:   CheckFaultEdge_Private(dm, label);
2156:   return 0;
2157: }

2159: /* Check that no cell have all vertices on the fault */
2160: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2161: {
2162:   IS              subpointIS;
2163:   const PetscInt *dmpoints;
2164:   PetscInt        defaultValue, cStart, cEnd, c, vStart, vEnd;

2166:   if (!label) return 0;
2167:   DMLabelGetDefaultValue(label, &defaultValue);
2168:   DMPlexGetSubpointIS(subdm, &subpointIS);
2169:   if (!subpointIS) return 0;
2170:   DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2171:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2172:   ISGetIndices(subpointIS, &dmpoints);
2173:   for (c = cStart; c < cEnd; ++c) {
2174:     PetscBool invalidCell = PETSC_TRUE;
2175:     PetscInt *closure     = NULL;
2176:     PetscInt  closureSize, cl;

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

2182:       if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2183:       DMLabelGetValue(label, closure[cl], &value);
2184:       if (value == defaultValue) {
2185:         invalidCell = PETSC_FALSE;
2186:         break;
2187:       }
2188:     }
2189:     DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2190:     if (invalidCell) {
2191:       ISRestoreIndices(subpointIS, &dmpoints);
2192:       ISDestroy(&subpointIS);
2193:       DMDestroy(&subdm);
2194:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %" PetscInt_FMT " has all of its vertices on the submesh.", dmpoints[c]);
2195:     }
2196:   }
2197:   ISRestoreIndices(subpointIS, &dmpoints);
2198:   return 0;
2199: }

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

2204:   Collective on dm

2206:   Input Parameters:
2207: + dm - The original DM
2208: . label - The label specifying the interface vertices
2209: . bdlabel - The optional label specifying the interface boundary vertices
2210: - bdvalue - Value of optional label specifying the interface boundary vertices

2212:   Output Parameters:
2213: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2214: . splitLabel - The label containing the split points, or NULL if no output is desired
2215: . dmInterface - The new interface DM, or NULL
2216: - dmHybrid - The new DM with cohesive cells

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

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

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

2231:   Level: developer

2233: .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexLabelCohesiveComplete()`, `DMPlexGetSubpointMap()`, `DMCreate()`
2234: @*/
2235: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, PetscInt bdvalue, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2236: {
2237:   DM       idm;
2238:   DMLabel  subpointMap, hlabel, slabel = NULL;
2239:   PetscInt dim;

2248:   DMGetDimension(dm, &dim);
2249:   DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2250:   DMPlexCheckValidSubmesh_Private(dm, label, idm);
2251:   DMPlexOrient(idm);
2252:   DMPlexGetSubpointMap(idm, &subpointMap);
2253:   DMLabelDuplicate(subpointMap, &hlabel);
2254:   DMLabelClearStratum(hlabel, dim);
2255:   if (splitLabel) {
2256:     const char *name;
2257:     char        sname[PETSC_MAX_PATH_LEN];

2259:     PetscObjectGetName((PetscObject)hlabel, &name);
2260:     PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2261:     PetscStrcat(sname, " split");
2262:     DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2263:   }
2264:   DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, bdvalue, PETSC_FALSE, idm);
2265:   if (dmInterface) {
2266:     *dmInterface = idm;
2267:   } else DMDestroy(&idm);
2268:   DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2269:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmHybrid);
2270:   if (hybridLabel) *hybridLabel = hlabel;
2271:   else DMLabelDestroy(&hlabel);
2272:   if (splitLabel) *splitLabel = slabel;
2273:   {
2274:     DM      cdm;
2275:     DMLabel ctLabel;

2277:     /* We need to somehow share the celltype label with the coordinate dm */
2278:     DMGetCoordinateDM(*dmHybrid, &cdm);
2279:     DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel);
2280:     DMSetLabel(cdm, ctLabel);
2281:   }
2282:   return 0;
2283: }

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

2287:      For any marked cell, the marked vertices constitute a single face
2288: */
2289: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2290: {
2291:   IS              subvertexIS = NULL;
2292:   const PetscInt *subvertices;
2293:   PetscInt       *pStart, *pEnd, pSize;
2294:   PetscInt        depth, dim, d, numSubVerticesInitial = 0, v;

2296:   *numFaces = 0;
2297:   *nFV      = 0;
2298:   DMPlexGetDepth(dm, &depth);
2299:   DMGetDimension(dm, &dim);
2300:   pSize = PetscMax(depth, dim) + 1;
2301:   PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2302:   for (d = 0; d <= depth; ++d) DMPlexGetSimplexOrBoxCells(dm, depth - d, &pStart[d], &pEnd[d]);
2303:   /* Loop over initial vertices and mark all faces in the collective star() */
2304:   if (vertexLabel) DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2305:   if (subvertexIS) {
2306:     ISGetSize(subvertexIS, &numSubVerticesInitial);
2307:     ISGetIndices(subvertexIS, &subvertices);
2308:   }
2309:   for (v = 0; v < numSubVerticesInitial; ++v) {
2310:     const PetscInt vertex = subvertices[v];
2311:     PetscInt      *star   = NULL;
2312:     PetscInt       starSize, s, numCells = 0, c;

2314:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2315:     for (s = 0; s < starSize * 2; s += 2) {
2316:       const PetscInt point = star[s];
2317:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2318:     }
2319:     for (c = 0; c < numCells; ++c) {
2320:       const PetscInt cell    = star[c];
2321:       PetscInt      *closure = NULL;
2322:       PetscInt       closureSize, cl;
2323:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

2325:       DMLabelGetValue(subpointMap, cell, &cellLoc);
2326:       if (cellLoc == 2) continue;
2328:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2329:       for (cl = 0; cl < closureSize * 2; cl += 2) {
2330:         const PetscInt point = closure[cl];
2331:         PetscInt       vertexLoc;

2333:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2334:           ++numCorners;
2335:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2336:           if (vertexLoc == value) closure[faceSize++] = point;
2337:         }
2338:       }
2339:       if (!(*nFV)) DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);
2341:       if (faceSize == *nFV) {
2342:         const PetscInt *cells = NULL;
2343:         PetscInt        numCells, nc;

2345:         ++(*numFaces);
2346:         for (cl = 0; cl < faceSize; ++cl) DMLabelSetValue(subpointMap, closure[cl], 0);
2347:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2348:         for (nc = 0; nc < numCells; ++nc) DMLabelSetValue(subpointMap, cells[nc], 2);
2349:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2350:       }
2351:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2352:     }
2353:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2354:   }
2355:   if (subvertexIS) ISRestoreIndices(subvertexIS, &subvertices);
2356:   ISDestroy(&subvertexIS);
2357:   PetscFree2(pStart, pEnd);
2358:   return 0;
2359: }

2361: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2362: {
2363:   IS              subvertexIS = NULL;
2364:   const PetscInt *subvertices;
2365:   PetscInt       *pStart, *pEnd;
2366:   PetscInt        dim, d, numSubVerticesInitial = 0, v;

2368:   DMGetDimension(dm, &dim);
2369:   PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd);
2370:   for (d = 0; d <= dim; ++d) DMPlexGetSimplexOrBoxCells(dm, dim - d, &pStart[d], &pEnd[d]);
2371:   /* Loop over initial vertices and mark all faces in the collective star() */
2372:   if (vertexLabel) {
2373:     DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2374:     if (subvertexIS) {
2375:       ISGetSize(subvertexIS, &numSubVerticesInitial);
2376:       ISGetIndices(subvertexIS, &subvertices);
2377:     }
2378:   }
2379:   for (v = 0; v < numSubVerticesInitial; ++v) {
2380:     const PetscInt vertex = subvertices[v];
2381:     PetscInt      *star   = NULL;
2382:     PetscInt       starSize, s, numFaces = 0, f;

2384:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2385:     for (s = 0; s < starSize * 2; s += 2) {
2386:       const PetscInt point = star[s];
2387:       PetscInt       faceLoc;

2389:       if ((point >= pStart[dim - 1]) && (point < pEnd[dim - 1])) {
2390:         if (markedFaces) {
2391:           DMLabelGetValue(vertexLabel, point, &faceLoc);
2392:           if (faceLoc < 0) continue;
2393:         }
2394:         star[numFaces++] = point;
2395:       }
2396:     }
2397:     for (f = 0; f < numFaces; ++f) {
2398:       const PetscInt face    = star[f];
2399:       PetscInt      *closure = NULL;
2400:       PetscInt       closureSize, c;
2401:       PetscInt       faceLoc;

2403:       DMLabelGetValue(subpointMap, face, &faceLoc);
2404:       if (faceLoc == dim - 1) continue;
2406:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2407:       for (c = 0; c < closureSize * 2; c += 2) {
2408:         const PetscInt point = closure[c];
2409:         PetscInt       vertexLoc;

2411:         if ((point >= pStart[0]) && (point < pEnd[0])) {
2412:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
2413:           if (vertexLoc != value) break;
2414:         }
2415:       }
2416:       if (c == closureSize * 2) {
2417:         const PetscInt *support;
2418:         PetscInt        supportSize, s;

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

2423:           for (d = 0; d < dim; ++d) {
2424:             if ((point >= pStart[d]) && (point < pEnd[d])) {
2425:               DMLabelSetValue(subpointMap, point, d);
2426:               break;
2427:             }
2428:           }
2429:         }
2430:         DMPlexGetSupportSize(dm, face, &supportSize);
2431:         DMPlexGetSupport(dm, face, &support);
2432:         for (s = 0; s < supportSize; ++s) DMLabelSetValue(subpointMap, support[s], dim);
2433:       }
2434:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2435:     }
2436:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2437:   }
2438:   if (subvertexIS) ISRestoreIndices(subvertexIS, &subvertices);
2439:   ISDestroy(&subvertexIS);
2440:   PetscFree2(pStart, pEnd);
2441:   return 0;
2442: }

2444: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2445: {
2446:   DMLabel         label = NULL;
2447:   const PetscInt *cone;
2448:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;

2450:   *numFaces = 0;
2451:   *nFV      = 0;
2452:   if (labelname) DMGetLabel(dm, labelname, &label);
2453:   *subCells = NULL;
2454:   DMGetDimension(dm, &dim);
2455:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2456:   if (cMax < 0) return 0;
2457:   if (label) {
2458:     for (c = cMax; c < cEnd; ++c) {
2459:       PetscInt val;

2461:       DMLabelGetValue(label, c, &val);
2462:       if (val == value) {
2463:         ++(*numFaces);
2464:         DMPlexGetConeSize(dm, c, &coneSize);
2465:       }
2466:     }
2467:   } else {
2468:     *numFaces = cEnd - cMax;
2469:     DMPlexGetConeSize(dm, cMax, &coneSize);
2470:   }
2471:   PetscMalloc1(*numFaces * 2, subCells);
2472:   if (!(*numFaces)) return 0;
2473:   *nFV = hasLagrange ? coneSize / 3 : coneSize / 2;
2474:   for (c = cMax; c < cEnd; ++c) {
2475:     const PetscInt *cells;
2476:     PetscInt        numCells;

2478:     if (label) {
2479:       PetscInt val;

2481:       DMLabelGetValue(label, c, &val);
2482:       if (val != value) continue;
2483:     }
2484:     DMPlexGetCone(dm, c, &cone);
2485:     for (p = 0; p < *nFV; ++p) DMLabelSetValue(subpointMap, cone[p], 0);
2486:     /* Negative face */
2487:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2488:     /* Not true in parallel
2490:     for (p = 0; p < numCells; ++p) {
2491:       DMLabelSetValue(subpointMap, cells[p], 2);
2492:       (*subCells)[subc++] = cells[p];
2493:     }
2494:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2495:     /* Positive face is not included */
2496:   }
2497:   return 0;
2498: }

2500: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2501: {
2502:   PetscInt *pStart, *pEnd;
2503:   PetscInt  dim, cMax, cEnd, c, d;

2505:   DMGetDimension(dm, &dim);
2506:   DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2507:   if (cMax < 0) return 0;
2508:   PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd);
2509:   for (d = 0; d <= dim; ++d) DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2510:   for (c = cMax; c < cEnd; ++c) {
2511:     const PetscInt *cone;
2512:     PetscInt       *closure = NULL;
2513:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

2515:     if (label) {
2516:       DMLabelGetValue(label, c, &val);
2517:       if (val != value) continue;
2518:     }
2519:     DMPlexGetConeSize(dm, c, &coneSize);
2520:     DMPlexGetCone(dm, c, &cone);
2521:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2523:     /* Negative face */
2524:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2525:     for (cl = 0; cl < closureSize * 2; cl += 2) {
2526:       const PetscInt point = closure[cl];

2528:       for (d = 0; d <= dim; ++d) {
2529:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2530:           DMLabelSetValue(subpointMap, point, d);
2531:           break;
2532:         }
2533:       }
2534:     }
2535:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2536:     /* Cells -- positive face is not included */
2537:     for (cl = 0; cl < 1; ++cl) {
2538:       const PetscInt *support;
2539:       PetscInt        supportSize, s;

2541:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2543:       DMPlexGetSupport(dm, cone[cl], &support);
2544:       for (s = 0; s < supportSize; ++s) DMLabelSetValue(subpointMap, support[s], dim);
2545:     }
2546:   }
2547:   PetscFree2(pStart, pEnd);
2548:   return 0;
2549: }

2551: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2552: {
2553:   MPI_Comm       comm;
2554:   PetscBool      posOrient = PETSC_FALSE;
2555:   const PetscInt debug     = 0;
2556:   PetscInt       cellDim, faceSize, f;

2558:   PetscObjectGetComm((PetscObject)dm, &comm);
2559:   DMGetDimension(dm, &cellDim);
2560:   if (debug) PetscPrintf(comm, "cellDim: %" PetscInt_FMT " numCorners: %" PetscInt_FMT "\n", cellDim, numCorners);

2562:   if (cellDim == 1 && numCorners == 2) {
2563:     /* Triangle */
2564:     faceSize  = numCorners - 1;
2565:     posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2566:   } else if (cellDim == 2 && numCorners == 3) {
2567:     /* Triangle */
2568:     faceSize  = numCorners - 1;
2569:     posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2570:   } else if (cellDim == 3 && numCorners == 4) {
2571:     /* Tetrahedron */
2572:     faceSize  = numCorners - 1;
2573:     posOrient = (oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2574:   } else if (cellDim == 1 && numCorners == 3) {
2575:     /* Quadratic line */
2576:     faceSize  = 1;
2577:     posOrient = PETSC_TRUE;
2578:   } else if (cellDim == 2 && numCorners == 4) {
2579:     /* Quads */
2580:     faceSize = 2;
2581:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2582:       posOrient = PETSC_TRUE;
2583:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2584:       posOrient = PETSC_TRUE;
2585:     } else {
2586:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2587:         posOrient = PETSC_FALSE;
2588:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2589:     }
2590:   } else if (cellDim == 2 && numCorners == 6) {
2591:     /* Quadratic triangle (I hate this) */
2592:     /* Edges are determined by the first 2 vertices (corners of edges) */
2593:     const PetscInt faceSizeTri = 3;
2594:     PetscInt       sortedIndices[3], i, iFace;
2595:     PetscBool      found                    = PETSC_FALSE;
2596:     PetscInt       faceVerticesTriSorted[9] = {
2597:       0, 3, 4, /* bottom */
2598:       1, 4, 5, /* right */
2599:       2, 3, 5, /* left */
2600:     };
2601:     PetscInt faceVerticesTri[9] = {
2602:       0, 3, 4, /* bottom */
2603:       1, 4, 5, /* right */
2604:       2, 5, 3, /* left */
2605:     };

2607:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2608:     PetscSortInt(faceSizeTri, sortedIndices);
2609:     for (iFace = 0; iFace < 3; ++iFace) {
2610:       const PetscInt ii = iFace * faceSizeTri;
2611:       PetscInt       fVertex, cVertex;

2613:       if ((sortedIndices[0] == faceVerticesTriSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTriSorted[ii + 1])) {
2614:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2615:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2616:             if (indices[cVertex] == faceVerticesTri[ii + fVertex]) {
2617:               faceVertices[fVertex] = origVertices[cVertex];
2618:               break;
2619:             }
2620:           }
2621:         }
2622:         found = PETSC_TRUE;
2623:         break;
2624:       }
2625:     }
2627:     if (posOriented) *posOriented = PETSC_TRUE;
2628:     return 0;
2629:   } else if (cellDim == 2 && numCorners == 9) {
2630:     /* Quadratic quad (I hate this) */
2631:     /* Edges are determined by the first 2 vertices (corners of edges) */
2632:     const PetscInt faceSizeQuad = 3;
2633:     PetscInt       sortedIndices[3], i, iFace;
2634:     PetscBool      found                      = PETSC_FALSE;
2635:     PetscInt       faceVerticesQuadSorted[12] = {
2636:       0, 1, 4, /* bottom */
2637:       1, 2, 5, /* right */
2638:       2, 3, 6, /* top */
2639:       0, 3, 7, /* left */
2640:     };
2641:     PetscInt faceVerticesQuad[12] = {
2642:       0, 1, 4, /* bottom */
2643:       1, 2, 5, /* right */
2644:       2, 3, 6, /* top */
2645:       3, 0, 7, /* left */
2646:     };

2648:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2649:     PetscSortInt(faceSizeQuad, sortedIndices);
2650:     for (iFace = 0; iFace < 4; ++iFace) {
2651:       const PetscInt ii = iFace * faceSizeQuad;
2652:       PetscInt       fVertex, cVertex;

2654:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadSorted[ii + 1])) {
2655:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2656:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2657:             if (indices[cVertex] == faceVerticesQuad[ii + fVertex]) {
2658:               faceVertices[fVertex] = origVertices[cVertex];
2659:               break;
2660:             }
2661:           }
2662:         }
2663:         found = PETSC_TRUE;
2664:         break;
2665:       }
2666:     }
2668:     if (posOriented) *posOriented = PETSC_TRUE;
2669:     return 0;
2670:   } else if (cellDim == 3 && numCorners == 8) {
2671:     /* Hexes
2672:        A hex is two oriented quads with the normal of the first
2673:        pointing up at the second.

2675:           7---6
2676:          /|  /|
2677:         4---5 |
2678:         | 1-|-2
2679:         |/  |/
2680:         0---3

2682:         Faces are determined by the first 4 vertices (corners of faces) */
2683:     const PetscInt faceSizeHex = 4;
2684:     PetscInt       sortedIndices[4], i, iFace;
2685:     PetscBool      found                     = PETSC_FALSE;
2686:     PetscInt       faceVerticesHexSorted[24] = {
2687:       0, 1, 2, 3, /* bottom */
2688:       4, 5, 6, 7, /* top */
2689:       0, 3, 4, 5, /* front */
2690:       2, 3, 5, 6, /* right */
2691:       1, 2, 6, 7, /* back */
2692:       0, 1, 4, 7, /* left */
2693:     };
2694:     PetscInt faceVerticesHex[24] = {
2695:       1, 2, 3, 0, /* bottom */
2696:       4, 5, 6, 7, /* top */
2697:       0, 3, 5, 4, /* front */
2698:       3, 2, 6, 5, /* right */
2699:       2, 1, 7, 6, /* back */
2700:       1, 0, 4, 7, /* left */
2701:     };

2703:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2704:     PetscSortInt(faceSizeHex, sortedIndices);
2705:     for (iFace = 0; iFace < 6; ++iFace) {
2706:       const PetscInt ii = iFace * faceSizeHex;
2707:       PetscInt       fVertex, cVertex;

2709:       if ((sortedIndices[0] == faceVerticesHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesHexSorted[ii + 3])) {
2710:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2711:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2712:             if (indices[cVertex] == faceVerticesHex[ii + fVertex]) {
2713:               faceVertices[fVertex] = origVertices[cVertex];
2714:               break;
2715:             }
2716:           }
2717:         }
2718:         found = PETSC_TRUE;
2719:         break;
2720:       }
2721:     }
2723:     if (posOriented) *posOriented = PETSC_TRUE;
2724:     return 0;
2725:   } else if (cellDim == 3 && numCorners == 10) {
2726:     /* Quadratic tet */
2727:     /* Faces are determined by the first 3 vertices (corners of faces) */
2728:     const PetscInt faceSizeTet = 6;
2729:     PetscInt       sortedIndices[6], i, iFace;
2730:     PetscBool      found                     = PETSC_FALSE;
2731:     PetscInt       faceVerticesTetSorted[24] = {
2732:       0, 1, 2, 6, 7, 8, /* bottom */
2733:       0, 3, 4, 6, 7, 9, /* front */
2734:       1, 4, 5, 7, 8, 9, /* right */
2735:       2, 3, 5, 6, 8, 9, /* left */
2736:     };
2737:     PetscInt faceVerticesTet[24] = {
2738:       0, 1, 2, 6, 7, 8, /* bottom */
2739:       0, 4, 3, 6, 7, 9, /* front */
2740:       1, 5, 4, 7, 8, 9, /* right */
2741:       2, 3, 5, 8, 6, 9, /* left */
2742:     };

2744:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2745:     PetscSortInt(faceSizeTet, sortedIndices);
2746:     for (iFace = 0; iFace < 4; ++iFace) {
2747:       const PetscInt ii = iFace * faceSizeTet;
2748:       PetscInt       fVertex, cVertex;

2750:       if ((sortedIndices[0] == faceVerticesTetSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTetSorted[ii + 1]) && (sortedIndices[2] == faceVerticesTetSorted[ii + 2]) && (sortedIndices[3] == faceVerticesTetSorted[ii + 3])) {
2751:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2752:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2753:             if (indices[cVertex] == faceVerticesTet[ii + fVertex]) {
2754:               faceVertices[fVertex] = origVertices[cVertex];
2755:               break;
2756:             }
2757:           }
2758:         }
2759:         found = PETSC_TRUE;
2760:         break;
2761:       }
2762:     }
2764:     if (posOriented) *posOriented = PETSC_TRUE;
2765:     return 0;
2766:   } else if (cellDim == 3 && numCorners == 27) {
2767:     /* Quadratic hexes (I hate this)
2768:        A hex is two oriented quads with the normal of the first
2769:        pointing up at the second.

2771:          7---6
2772:         /|  /|
2773:        4---5 |
2774:        | 3-|-2
2775:        |/  |/
2776:        0---1

2778:        Faces are determined by the first 4 vertices (corners of faces) */
2779:     const PetscInt faceSizeQuadHex = 9;
2780:     PetscInt       sortedIndices[9], i, iFace;
2781:     PetscBool      found                         = PETSC_FALSE;
2782:     PetscInt       faceVerticesQuadHexSorted[54] = {
2783:       0, 1, 2, 3, 8,  9,  10, 11, 24, /* bottom */
2784:       4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2785:       0, 1, 4, 5, 8,  12, 16, 17, 22, /* front */
2786:       1, 2, 5, 6, 9,  13, 17, 18, 21, /* right */
2787:       2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2788:       0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2789:     };
2790:     PetscInt faceVerticesQuadHex[54] = {
2791:       3, 2, 1, 0, 10, 9,  8,  11, 24, /* bottom */
2792:       4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2793:       0, 1, 5, 4, 8,  17, 12, 16, 22, /* front */
2794:       1, 2, 6, 5, 9,  18, 13, 17, 21, /* right */
2795:       2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2796:       3, 0, 4, 7, 11, 16, 15, 19, 20  /* left */
2797:     };

2799:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2800:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2801:     for (iFace = 0; iFace < 6; ++iFace) {
2802:       const PetscInt ii = iFace * faceSizeQuadHex;
2803:       PetscInt       fVertex, cVertex;

2805:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesQuadHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesQuadHexSorted[ii + 3])) {
2806:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2807:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2808:             if (indices[cVertex] == faceVerticesQuadHex[ii + fVertex]) {
2809:               faceVertices[fVertex] = origVertices[cVertex];
2810:               break;
2811:             }
2812:           }
2813:         }
2814:         found = PETSC_TRUE;
2815:         break;
2816:       }
2817:     }
2819:     if (posOriented) *posOriented = PETSC_TRUE;
2820:     return 0;
2821:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2822:   if (!posOrient) {
2823:     if (debug) PetscPrintf(comm, "  Reversing initial face orientation\n");
2824:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize - 1 - f];
2825:   } else {
2826:     if (debug) PetscPrintf(comm, "  Keeping initial face orientation\n");
2827:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2828:   }
2829:   if (posOriented) *posOriented = posOrient;
2830:   return 0;
2831: }

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

2837:   Not collective

2839:   Input Parameters:
2840: + dm           - The original mesh
2841: . cell         - The cell mesh point
2842: . faceSize     - The number of vertices on the face
2843: . face         - The face vertices
2844: . numCorners   - The number of vertices on the cell
2845: . indices      - Local numbering of face vertices in cell cone
2846: - origVertices - Original face vertices

2848:   Output Parameters:
2849: + faceVertices - The face vertices properly oriented
2850: - posOriented  - PETSC_TRUE if the face was oriented with outward normal

2852:   Level: developer

2854: .seealso: `DMPlexGetCone()`
2855: @*/
2856: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2857: {
2858:   const PetscInt *cone = NULL;
2859:   PetscInt        coneSize, v, f, v2;
2860:   PetscInt        oppositeVertex = -1;

2862:   DMPlexGetConeSize(dm, cell, &coneSize);
2863:   DMPlexGetCone(dm, cell, &cone);
2864:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2865:     PetscBool found = PETSC_FALSE;

2867:     for (f = 0; f < faceSize; ++f) {
2868:       if (face[f] == cone[v]) {
2869:         found = PETSC_TRUE;
2870:         break;
2871:       }
2872:     }
2873:     if (found) {
2874:       indices[v2]      = v;
2875:       origVertices[v2] = cone[v];
2876:       ++v2;
2877:     } else {
2878:       oppositeVertex = v;
2879:     }
2880:   }
2881:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2882:   return 0;
2883: }

2885: /*
2886:   DMPlexInsertFace_Internal - Puts a face into the mesh

2888:   Not collective

2890:   Input Parameters:
2891:   + dm              - The DMPlex
2892:   . numFaceVertex   - The number of vertices in the face
2893:   . faceVertices    - The vertices in the face for dm
2894:   . subfaceVertices - The vertices in the face for subdm
2895:   . numCorners      - The number of vertices in the cell
2896:   . cell            - A cell in dm containing the face
2897:   . subcell         - A cell in subdm containing the face
2898:   . firstFace       - First face in the mesh
2899:   - newFacePoint    - Next face in the mesh

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

2904:   Level: developer
2905: */
2906: 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)
2907: {
2908:   MPI_Comm        comm;
2909:   DM_Plex        *submesh = (DM_Plex *)subdm->data;
2910:   const PetscInt *faces;
2911:   PetscInt        numFaces, coneSize;

2913:   PetscObjectGetComm((PetscObject)dm, &comm);
2914:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2916: #if 0
2917:   /* Cannot use this because support() has not been constructed yet */
2918:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2919: #else
2920:   {
2921:     PetscInt f;

2923:     numFaces = 0;
2924:     DMGetWorkArray(subdm, 1, MPIU_INT, (void **)&faces);
2925:     for (f = firstFace; f < *newFacePoint; ++f) {
2926:       PetscInt dof, off, d;

2928:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2929:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2930:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2931:       for (d = 0; d < dof; ++d) {
2932:         const PetscInt p = submesh->cones[off + d];
2933:         PetscInt       v;

2935:         for (v = 0; v < numFaceVertices; ++v) {
2936:           if (subfaceVertices[v] == p) break;
2937:         }
2938:         if (v == numFaceVertices) break;
2939:       }
2940:       if (d == dof) {
2941:         numFaces               = 1;
2942:         ((PetscInt *)faces)[0] = f;
2943:       }
2944:     }
2945:   }
2946: #endif
2948:   if (numFaces == 1) {
2949:     /* Add the other cell neighbor for this face */
2950:     DMPlexSetCone(subdm, subcell, faces);
2951:   } else {
2952:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2953:     PetscBool posOriented;

2955:     DMGetWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2956:     origVertices        = &orientedVertices[numFaceVertices];
2957:     indices             = &orientedVertices[numFaceVertices * 2];
2958:     orientedSubVertices = &orientedVertices[numFaceVertices * 3];
2959:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2960:     /* TODO: I know that routine should return a permutation, not the indices */
2961:     for (v = 0; v < numFaceVertices; ++v) {
2962:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2963:       for (ov = 0; ov < numFaceVertices; ++ov) {
2964:         if (orientedVertices[ov] == vertex) {
2965:           orientedSubVertices[ov] = subvertex;
2966:           break;
2967:         }
2968:       }
2970:     }
2971:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2972:     DMPlexSetCone(subdm, subcell, newFacePoint);
2973:     DMRestoreWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2974:     ++(*newFacePoint);
2975:   }
2976: #if 0
2977:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2978: #else
2979:   DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **)&faces);
2980: #endif
2981:   return 0;
2982: }

2984: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2985: {
2986:   MPI_Comm        comm;
2987:   DMLabel         subpointMap;
2988:   IS              subvertexIS, subcellIS;
2989:   const PetscInt *subVertices, *subCells;
2990:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2991:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2992:   PetscInt        vStart, vEnd, c, f;

2994:   PetscObjectGetComm((PetscObject)dm, &comm);
2995:   /* Create subpointMap which marks the submesh */
2996:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2997:   DMPlexSetSubpointMap(subdm, subpointMap);
2998:   DMLabelDestroy(&subpointMap);
2999:   if (vertexLabel) DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);
3000:   /* Setup chart */
3001:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3002:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3003:   DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices);
3004:   DMPlexSetVTKCellHeight(subdm, 1);
3005:   /* Set cone sizes */
3006:   firstSubVertex = numSubCells;
3007:   firstSubFace   = numSubCells + numSubVertices;
3008:   newFacePoint   = firstSubFace;
3009:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3010:   if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
3011:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
3012:   if (subcellIS) ISGetIndices(subcellIS, &subCells);
3013:   for (c = 0; c < numSubCells; ++c) DMPlexSetConeSize(subdm, c, 1);
3014:   for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) DMPlexSetConeSize(subdm, f, nFV);
3015:   DMSetUp(subdm);
3016:   /* Create face cones */
3017:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3018:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3019:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3020:   for (c = 0; c < numSubCells; ++c) {
3021:     const PetscInt cell    = subCells[c];
3022:     const PetscInt subcell = c;
3023:     PetscInt      *closure = NULL;
3024:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

3026:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
3027:     for (cl = 0; cl < closureSize * 2; cl += 2) {
3028:       const PetscInt point = closure[cl];
3029:       PetscInt       subVertex;

3031:       if ((point >= vStart) && (point < vEnd)) {
3032:         ++numCorners;
3033:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
3034:         if (subVertex >= 0) {
3035:           closure[faceSize] = point;
3036:           subface[faceSize] = firstSubVertex + subVertex;
3037:           ++faceSize;
3038:         }
3039:       }
3040:     }
3042:     if (faceSize == nFV) DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
3043:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
3044:   }
3045:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3046:   DMPlexSymmetrize(subdm);
3047:   DMPlexStratify(subdm);
3048:   /* Build coordinates */
3049:   {
3050:     PetscSection coordSection, subCoordSection;
3051:     Vec          coordinates, subCoordinates;
3052:     PetscScalar *coords, *subCoords;
3053:     PetscInt     numComp, coordSize, v;
3054:     const char  *name;

3056:     DMGetCoordinateSection(dm, &coordSection);
3057:     DMGetCoordinatesLocal(dm, &coordinates);
3058:     DMGetCoordinateSection(subdm, &subCoordSection);
3059:     PetscSectionSetNumFields(subCoordSection, 1);
3060:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3061:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3062:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices);
3063:     for (v = 0; v < numSubVertices; ++v) {
3064:       const PetscInt vertex    = subVertices[v];
3065:       const PetscInt subvertex = firstSubVertex + v;
3066:       PetscInt       dof;

3068:       PetscSectionGetDof(coordSection, vertex, &dof);
3069:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3070:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3071:     }
3072:     PetscSectionSetUp(subCoordSection);
3073:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3074:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3075:     PetscObjectGetName((PetscObject)coordinates, &name);
3076:     PetscObjectSetName((PetscObject)subCoordinates, name);
3077:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3078:     VecSetType(subCoordinates, VECSTANDARD);
3079:     if (coordSize) {
3080:       VecGetArray(coordinates, &coords);
3081:       VecGetArray(subCoordinates, &subCoords);
3082:       for (v = 0; v < numSubVertices; ++v) {
3083:         const PetscInt vertex    = subVertices[v];
3084:         const PetscInt subvertex = firstSubVertex + v;
3085:         PetscInt       dof, off, sdof, soff, d;

3087:         PetscSectionGetDof(coordSection, vertex, &dof);
3088:         PetscSectionGetOffset(coordSection, vertex, &off);
3089:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3090:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3092:         for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3093:       }
3094:       VecRestoreArray(coordinates, &coords);
3095:       VecRestoreArray(subCoordinates, &subCoords);
3096:     }
3097:     DMSetCoordinatesLocal(subdm, subCoordinates);
3098:     VecDestroy(&subCoordinates);
3099:   }
3100:   /* Cleanup */
3101:   if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
3102:   ISDestroy(&subvertexIS);
3103:   if (subcellIS) ISRestoreIndices(subcellIS, &subCells);
3104:   ISDestroy(&subcellIS);
3105:   return 0;
3106: }

3108: /* TODO: Fix this to properly propagate up error conditions it may find */
3109: static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3110: {
3111:   PetscInt       subPoint;

3114:   PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3115:   if (ierr) return -1;
3116:   return subPoint < 0 ? subPoint : firstSubPoint + subPoint;
3117: }

3119: /* TODO: Fix this to properly propagate up error conditions it may find */
3120: static inline PetscInt DMPlexFilterPointPerm_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[], const PetscInt subIndices[])
3121: {
3122:   PetscInt       subPoint;

3125:   PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3126:   if (ierr) return -1;
3127:   return subPoint < 0 ? subPoint : firstSubPoint + (subIndices ? subIndices[subPoint] : subPoint);
3128: }

3130: static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm)
3131: {
3132:   DMLabel  depthLabel;
3133:   PetscInt Nl, l, d;

3135:   // Reset depth label for fast lookup
3136:   DMPlexGetDepthLabel(dm, &depthLabel);
3137:   DMLabelMakeAllInvalid_Internal(depthLabel);
3138:   DMGetNumLabels(dm, &Nl);
3139:   for (l = 0; l < Nl; ++l) {
3140:     DMLabel         label, newlabel;
3141:     const char     *lname;
3142:     PetscBool       isDepth, isDim, isCelltype, isVTK;
3143:     IS              valueIS;
3144:     const PetscInt *values;
3145:     PetscInt        Nv, v;

3147:     DMGetLabelName(dm, l, &lname);
3148:     PetscStrcmp(lname, "depth", &isDepth);
3149:     PetscStrcmp(lname, "dim", &isDim);
3150:     PetscStrcmp(lname, "celltype", &isCelltype);
3151:     PetscStrcmp(lname, "vtk", &isVTK);
3152:     if (isDepth || isDim || isCelltype || isVTK) continue;
3153:     DMCreateLabel(subdm, lname);
3154:     DMGetLabel(dm, lname, &label);
3155:     DMGetLabel(subdm, lname, &newlabel);
3156:     DMLabelGetDefaultValue(label, &v);
3157:     DMLabelSetDefaultValue(newlabel, v);
3158:     DMLabelGetValueIS(label, &valueIS);
3159:     ISGetLocalSize(valueIS, &Nv);
3160:     ISGetIndices(valueIS, &values);
3161:     for (v = 0; v < Nv; ++v) {
3162:       IS              pointIS;
3163:       const PetscInt *points;
3164:       PetscInt        Np, p;

3166:       DMLabelGetStratumIS(label, values[v], &pointIS);
3167:       ISGetLocalSize(pointIS, &Np);
3168:       ISGetIndices(pointIS, &points);
3169:       for (p = 0; p < Np; ++p) {
3170:         const PetscInt point = points[p];
3171:         PetscInt       subp;

3173:         DMPlexGetPointDepth(dm, point, &d);
3174:         subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]);
3175:         if (subp >= 0) DMLabelSetValue(newlabel, subp, values[v]);
3176:       }
3177:       ISRestoreIndices(pointIS, &points);
3178:       ISDestroy(&pointIS);
3179:     }
3180:     ISRestoreIndices(valueIS, &values);
3181:     ISDestroy(&valueIS);
3182:   }
3183:   return 0;
3184: }

3186: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3187: {
3188:   MPI_Comm         comm;
3189:   DMLabel          subpointMap;
3190:   IS              *subpointIS;
3191:   const PetscInt **subpoints;
3192:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3193:   PetscInt         totSubPoints = 0, maxConeSize, dim, sdim, cdim, p, d, v;
3194:   PetscMPIInt      rank;

3196:   PetscObjectGetComm((PetscObject)dm, &comm);
3197:   MPI_Comm_rank(comm, &rank);
3198:   /* Create subpointMap which marks the submesh */
3199:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3200:   DMPlexSetSubpointMap(subdm, subpointMap);
3201:   if (cellHeight) {
3202:     if (isCohesive) DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
3203:     else DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);
3204:   } else {
3205:     DMLabel         depth;
3206:     IS              pointIS;
3207:     const PetscInt *points;
3208:     PetscInt        numPoints = 0;

3210:     DMPlexGetDepthLabel(dm, &depth);
3211:     DMLabelGetStratumIS(label, value, &pointIS);
3212:     if (pointIS) {
3213:       ISGetIndices(pointIS, &points);
3214:       ISGetLocalSize(pointIS, &numPoints);
3215:     }
3216:     for (p = 0; p < numPoints; ++p) {
3217:       PetscInt *closure = NULL;
3218:       PetscInt  closureSize, c, pdim;

3220:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3221:       for (c = 0; c < closureSize * 2; c += 2) {
3222:         DMLabelGetValue(depth, closure[c], &pdim);
3223:         DMLabelSetValue(subpointMap, closure[c], pdim);
3224:       }
3225:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3226:     }
3227:     if (pointIS) ISRestoreIndices(pointIS, &points);
3228:     ISDestroy(&pointIS);
3229:   }
3230:   /* Setup chart */
3231:   DMGetDimension(dm, &dim);
3232:   DMGetCoordinateDim(dm, &cdim);
3233:   PetscMalloc4(dim + 1, &numSubPoints, dim + 1, &firstSubPoint, dim + 1, &subpointIS, dim + 1, &subpoints);
3234:   for (d = 0; d <= dim; ++d) {
3235:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3236:     totSubPoints += numSubPoints[d];
3237:   }
3238:   // Determine submesh dimension
3239:   DMGetDimension(subdm, &sdim);
3240:   if (sdim > 0) {
3241:     // Calling function knows what dimension to use, and we include neighboring cells as well
3242:     sdim = dim;
3243:   } else {
3244:     // We reset the subdimension based on what is being selected
3245:     PetscInt lsdim;
3246:     for (lsdim = dim; lsdim >= 0; --lsdim)
3247:       if (numSubPoints[lsdim]) break;
3248:     MPI_Allreduce(&lsdim, &sdim, 1, MPIU_INT, MPIU_MAX, comm);
3249:     DMSetDimension(subdm, sdim);
3250:     DMSetCoordinateDim(subdm, cdim);
3251:   }
3252:   DMPlexSetChart(subdm, 0, totSubPoints);
3253:   DMPlexSetVTKCellHeight(subdm, cellHeight);
3254:   /* Set cone sizes */
3255:   firstSubPoint[sdim] = 0;
3256:   firstSubPoint[0]    = firstSubPoint[sdim] + numSubPoints[sdim];
3257:   if (sdim > 1) firstSubPoint[sdim - 1] = firstSubPoint[0] + numSubPoints[0];
3258:   if (sdim > 2) firstSubPoint[sdim - 2] = firstSubPoint[sdim - 1] + numSubPoints[sdim - 1];
3259:   for (d = 0; d <= sdim; ++d) {
3260:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3261:     if (subpointIS[d]) ISGetIndices(subpointIS[d], &subpoints[d]);
3262:   }
3263:   /* We do not want this label automatically computed, instead we compute it here */
3264:   DMCreateLabel(subdm, "celltype");
3265:   for (d = 0; d <= sdim; ++d) {
3266:     for (p = 0; p < numSubPoints[d]; ++p) {
3267:       const PetscInt  point    = subpoints[d][p];
3268:       const PetscInt  subpoint = firstSubPoint[d] + p;
3269:       const PetscInt *cone;
3270:       PetscInt        coneSize;

3272:       DMPlexGetConeSize(dm, point, &coneSize);
3273:       if (cellHeight && (d == sdim)) {
3274:         PetscInt coneSizeNew, c, val;

3276:         DMPlexGetCone(dm, point, &cone);
3277:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3278:           DMLabelGetValue(subpointMap, cone[c], &val);
3279:           if (val >= 0) coneSizeNew++;
3280:         }
3281:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3282:         DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3283:       } else {
3284:         DMPolytopeType ct;

3286:         DMPlexSetConeSize(subdm, subpoint, coneSize);
3287:         DMPlexGetCellType(dm, point, &ct);
3288:         DMPlexSetCellType(subdm, subpoint, ct);
3289:       }
3290:     }
3291:   }
3292:   DMLabelDestroy(&subpointMap);
3293:   DMSetUp(subdm);
3294:   /* Set cones */
3295:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3296:   PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
3297:   for (d = 0; d <= sdim; ++d) {
3298:     for (p = 0; p < numSubPoints[d]; ++p) {
3299:       const PetscInt  point    = subpoints[d][p];
3300:       const PetscInt  subpoint = firstSubPoint[d] + p;
3301:       const PetscInt *cone, *ornt;
3302:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;

3304:       if (d == sdim - 1) {
3305:         const PetscInt *support, *cone, *ornt;
3306:         PetscInt        supportSize, coneSize, s, subc;

3308:         DMPlexGetSupport(dm, point, &support);
3309:         DMPlexGetSupportSize(dm, point, &supportSize);
3310:         for (s = 0; s < supportSize; ++s) {
3311:           PetscBool isHybrid = PETSC_FALSE;

3313:           DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3314:           if (!isHybrid) continue;
3315:           PetscFindInt(support[s], numSubPoints[d + 1], subpoints[d + 1], &subc);
3316:           if (subc >= 0) {
3317:             const PetscInt ccell = subpoints[d + 1][subc];

3319:             DMPlexGetCone(dm, ccell, &cone);
3320:             DMPlexGetConeSize(dm, ccell, &coneSize);
3321:             DMPlexGetConeOrientation(dm, ccell, &ornt);
3322:             for (c = 0; c < coneSize; ++c) {
3323:               if (cone[c] == point) {
3324:                 fornt = ornt[c];
3325:                 break;
3326:               }
3327:             }
3328:             break;
3329:           }
3330:         }
3331:       }
3332:       DMPlexGetConeSize(dm, point, &coneSize);
3333:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3334:       DMPlexGetCone(dm, point, &cone);
3335:       DMPlexGetConeOrientation(dm, point, &ornt);
3336:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3337:         PetscFindInt(cone[c], numSubPoints[d - 1], subpoints[d - 1], &subc);
3338:         if (subc >= 0) {
3339:           coneNew[coneSizeNew] = firstSubPoint[d - 1] + subc;
3340:           orntNew[coneSizeNew] = ornt[c];
3341:           ++coneSizeNew;
3342:         }
3343:       }
3345:       DMPlexSetCone(subdm, subpoint, coneNew);
3346:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3347:       if (fornt < 0) DMPlexOrientPoint(subdm, subpoint, fornt);
3348:     }
3349:   }
3350:   PetscFree2(coneNew, orntNew);
3351:   DMPlexSymmetrize(subdm);
3352:   DMPlexStratify(subdm);
3353:   /* Build coordinates */
3354:   {
3355:     PetscSection coordSection, subCoordSection;
3356:     Vec          coordinates, subCoordinates;
3357:     PetscScalar *coords, *subCoords;
3358:     PetscInt     cdim, numComp, coordSize;
3359:     const char  *name;

3361:     DMGetCoordinateDim(dm, &cdim);
3362:     DMGetCoordinateSection(dm, &coordSection);
3363:     DMGetCoordinatesLocal(dm, &coordinates);
3364:     DMGetCoordinateSection(subdm, &subCoordSection);
3365:     PetscSectionSetNumFields(subCoordSection, 1);
3366:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3367:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3368:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0] + numSubPoints[0]);
3369:     for (v = 0; v < numSubPoints[0]; ++v) {
3370:       const PetscInt vertex    = subpoints[0][v];
3371:       const PetscInt subvertex = firstSubPoint[0] + v;
3372:       PetscInt       dof;

3374:       PetscSectionGetDof(coordSection, vertex, &dof);
3375:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3376:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3377:     }
3378:     PetscSectionSetUp(subCoordSection);
3379:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3380:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3381:     PetscObjectGetName((PetscObject)coordinates, &name);
3382:     PetscObjectSetName((PetscObject)subCoordinates, name);
3383:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3384:     VecSetBlockSize(subCoordinates, cdim);
3385:     VecSetType(subCoordinates, VECSTANDARD);
3386:     VecGetArray(coordinates, &coords);
3387:     VecGetArray(subCoordinates, &subCoords);
3388:     for (v = 0; v < numSubPoints[0]; ++v) {
3389:       const PetscInt vertex    = subpoints[0][v];
3390:       const PetscInt subvertex = firstSubPoint[0] + v;
3391:       PetscInt       dof, off, sdof, soff, d;

3393:       PetscSectionGetDof(coordSection, vertex, &dof);
3394:       PetscSectionGetOffset(coordSection, vertex, &off);
3395:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3396:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3398:       for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3399:     }
3400:     VecRestoreArray(coordinates, &coords);
3401:     VecRestoreArray(subCoordinates, &subCoords);
3402:     DMSetCoordinatesLocal(subdm, subCoordinates);
3403:     VecDestroy(&subCoordinates);
3404:   }
3405:   /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3406:   {
3407:     PetscSF            sfPoint, sfPointSub;
3408:     IS                 subpIS;
3409:     const PetscSFNode *remotePoints;
3410:     PetscSFNode       *sremotePoints = NULL, *newLocalPoints = NULL, *newOwners = NULL;
3411:     const PetscInt    *localPoints, *subpoints, *rootdegree;
3412:     PetscInt          *slocalPoints = NULL, *sortedPoints = NULL, *sortedIndices = NULL;
3413:     PetscInt           numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl = 0, ll = 0, pStart, pEnd, p;
3414:     PetscMPIInt        rank, size;

3416:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
3417:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
3418:     DMGetPointSF(dm, &sfPoint);
3419:     DMGetPointSF(subdm, &sfPointSub);
3420:     DMPlexGetChart(dm, &pStart, &pEnd);
3421:     DMPlexGetChart(subdm, NULL, &numSubroots);
3422:     DMPlexGetSubpointIS(subdm, &subpIS);
3423:     if (subpIS) {
3424:       PetscBool sorted = PETSC_TRUE;

3426:       ISGetIndices(subpIS, &subpoints);
3427:       ISGetLocalSize(subpIS, &numSubpoints);
3428:       for (p = 1; p < numSubpoints; ++p) sorted = sorted && (subpoints[p] >= subpoints[p - 1]) ? PETSC_TRUE : PETSC_FALSE;
3429:       if (!sorted) {
3430:         PetscMalloc2(numSubpoints, &sortedPoints, numSubpoints, &sortedIndices);
3431:         for (p = 0; p < numSubpoints; ++p) sortedIndices[p] = p;
3432:         PetscArraycpy(sortedPoints, subpoints, numSubpoints);
3433:         PetscSortIntWithArray(numSubpoints, sortedPoints, sortedIndices);
3434:       }
3435:     }
3436:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3437:     if (numRoots >= 0) {
3438:       PetscSFComputeDegreeBegin(sfPoint, &rootdegree);
3439:       PetscSFComputeDegreeEnd(sfPoint, &rootdegree);
3440:       PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners);
3441:       for (p = 0; p < pEnd - pStart; ++p) {
3442:         newLocalPoints[p].rank  = -2;
3443:         newLocalPoints[p].index = -2;
3444:       }
3445:       /* Set subleaves */
3446:       for (l = 0; l < numLeaves; ++l) {
3447:         const PetscInt point    = localPoints[l];
3448:         const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);

3450:         if (subpoint < 0) continue;
3451:         newLocalPoints[point - pStart].rank  = rank;
3452:         newLocalPoints[point - pStart].index = subpoint;
3453:         ++numSubleaves;
3454:       }
3455:       /* Must put in owned subpoints */
3456:       for (p = pStart; p < pEnd; ++p) {
3457:         newOwners[p - pStart].rank  = -3;
3458:         newOwners[p - pStart].index = -3;
3459:       }
3460:       for (p = 0; p < numSubpoints; ++p) {
3461:         /* Hold on to currently owned points */
3462:         if (rootdegree[subpoints[p] - pStart]) newOwners[subpoints[p] - pStart].rank = rank + size;
3463:         else newOwners[subpoints[p] - pStart].rank = rank;
3464:         newOwners[subpoints[p] - pStart].index = p;
3465:       }
3466:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3467:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3468:       for (p = pStart; p < pEnd; ++p)
3469:         if (newOwners[p - pStart].rank >= size) newOwners[p - pStart].rank -= size;
3470:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3471:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3472:       PetscMalloc1(numSubleaves, &slocalPoints);
3473:       PetscMalloc1(numSubleaves, &sremotePoints);
3474:       for (l = 0; l < numLeaves; ++l) {
3475:         const PetscInt point    = localPoints[l];
3476:         const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);

3478:         if (subpoint < 0) continue;
3479:         if (newLocalPoints[point].rank == rank) {
3480:           ++ll;
3481:           continue;
3482:         }
3483:         slocalPoints[sl]        = subpoint;
3484:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3485:         sremotePoints[sl].index = newLocalPoints[point].index;
3488:         ++sl;
3489:       }
3491:       PetscFree2(newLocalPoints, newOwners);
3492:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3493:     }
3494:     if (subpIS) ISRestoreIndices(subpIS, &subpoints);
3495:     PetscFree2(sortedPoints, sortedIndices);
3496:   }
3497:   /* Filter labels */
3498:   DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm);
3499:   /* Cleanup */
3500:   for (d = 0; d <= sdim; ++d) {
3501:     if (subpointIS[d]) ISRestoreIndices(subpointIS[d], &subpoints[d]);
3502:     ISDestroy(&subpointIS[d]);
3503:   }
3504:   PetscFree4(numSubPoints, firstSubPoint, subpointIS, subpoints);
3505:   return 0;
3506: }

3508: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3509: {
3510:   DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3511:   return 0;
3512: }

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

3517:   Input Parameters:
3518: + dm           - The original mesh
3519: . vertexLabel  - The DMLabel marking points contained in the surface
3520: . value        - The label value to use
3521: - markedFaces  - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked

3523:   Output Parameter:
3524: . subdm - The surface mesh

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

3528:   Level: developer

3530: .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`
3531: @*/
3532: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3533: {
3534:   DMPlexInterpolatedFlag interpolated;
3535:   PetscInt               dim, cdim;

3539:   DMGetDimension(dm, &dim);
3540:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3541:   DMSetType(*subdm, DMPLEX);
3542:   DMSetDimension(*subdm, dim - 1);
3543:   DMGetCoordinateDim(dm, &cdim);
3544:   DMSetCoordinateDim(*subdm, cdim);
3545:   DMPlexIsInterpolated(dm, &interpolated);
3547:   if (interpolated) {
3548:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3549:   } else {
3550:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3551:   }
3552:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3553:   return 0;
3554: }

3556: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3557: {
3558:   MPI_Comm        comm;
3559:   DMLabel         subpointMap;
3560:   IS              subvertexIS;
3561:   const PetscInt *subVertices;
3562:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3563:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3564:   PetscInt        c, f;

3566:   PetscObjectGetComm((PetscObject)dm, &comm);
3567:   /* Create subpointMap which marks the submesh */
3568:   DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3569:   DMPlexSetSubpointMap(subdm, subpointMap);
3570:   DMLabelDestroy(&subpointMap);
3571:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3572:   /* Setup chart */
3573:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3574:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3575:   DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices);
3576:   DMPlexSetVTKCellHeight(subdm, 1);
3577:   /* Set cone sizes */
3578:   firstSubVertex = numSubCells;
3579:   firstSubFace   = numSubCells + numSubVertices;
3580:   newFacePoint   = firstSubFace;
3581:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3582:   if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
3583:   for (c = 0; c < numSubCells; ++c) DMPlexSetConeSize(subdm, c, 1);
3584:   for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) DMPlexSetConeSize(subdm, f, nFV);
3585:   DMSetUp(subdm);
3586:   /* Create face cones */
3587:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3588:   DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3589:   for (c = 0; c < numSubCells; ++c) {
3590:     const PetscInt  cell    = subCells[c];
3591:     const PetscInt  subcell = c;
3592:     const PetscInt *cone, *cells;
3593:     PetscBool       isHybrid = PETSC_FALSE;
3594:     PetscInt        numCells, subVertex, p, v;

3596:     DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3597:     if (!isHybrid) continue;
3598:     DMPlexGetCone(dm, cell, &cone);
3599:     for (v = 0; v < nFV; ++v) {
3600:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3601:       subface[v] = firstSubVertex + subVertex;
3602:     }
3603:     DMPlexSetCone(subdm, newFacePoint, subface);
3604:     DMPlexSetCone(subdm, subcell, &newFacePoint);
3605:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3606:     /* Not true in parallel
3608:     for (p = 0; p < numCells; ++p) {
3609:       PetscInt  negsubcell;
3610:       PetscBool isHybrid = PETSC_FALSE;

3612:       DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3613:       if (isHybrid) continue;
3614:       /* I know this is a crap search */
3615:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3616:         if (subCells[negsubcell] == cells[p]) break;
3617:       }
3619:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3620:     }
3621:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3622:     ++newFacePoint;
3623:   }
3624:   DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3625:   DMPlexSymmetrize(subdm);
3626:   DMPlexStratify(subdm);
3627:   /* Build coordinates */
3628:   {
3629:     PetscSection coordSection, subCoordSection;
3630:     Vec          coordinates, subCoordinates;
3631:     PetscScalar *coords, *subCoords;
3632:     PetscInt     cdim, numComp, coordSize, v;
3633:     const char  *name;

3635:     DMGetCoordinateDim(dm, &cdim);
3636:     DMGetCoordinateSection(dm, &coordSection);
3637:     DMGetCoordinatesLocal(dm, &coordinates);
3638:     DMGetCoordinateSection(subdm, &subCoordSection);
3639:     PetscSectionSetNumFields(subCoordSection, 1);
3640:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3641:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3642:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices);
3643:     for (v = 0; v < numSubVertices; ++v) {
3644:       const PetscInt vertex    = subVertices[v];
3645:       const PetscInt subvertex = firstSubVertex + v;
3646:       PetscInt       dof;

3648:       PetscSectionGetDof(coordSection, vertex, &dof);
3649:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3650:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3651:     }
3652:     PetscSectionSetUp(subCoordSection);
3653:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3654:     VecCreate(PETSC_COMM_SELF, &subCoordinates);
3655:     PetscObjectGetName((PetscObject)coordinates, &name);
3656:     PetscObjectSetName((PetscObject)subCoordinates, name);
3657:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3658:     VecSetBlockSize(subCoordinates, cdim);
3659:     VecSetType(subCoordinates, VECSTANDARD);
3660:     VecGetArray(coordinates, &coords);
3661:     VecGetArray(subCoordinates, &subCoords);
3662:     for (v = 0; v < numSubVertices; ++v) {
3663:       const PetscInt vertex    = subVertices[v];
3664:       const PetscInt subvertex = firstSubVertex + v;
3665:       PetscInt       dof, off, sdof, soff, d;

3667:       PetscSectionGetDof(coordSection, vertex, &dof);
3668:       PetscSectionGetOffset(coordSection, vertex, &off);
3669:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3670:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3672:       for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3673:     }
3674:     VecRestoreArray(coordinates, &coords);
3675:     VecRestoreArray(subCoordinates, &subCoords);
3676:     DMSetCoordinatesLocal(subdm, subCoordinates);
3677:     VecDestroy(&subCoordinates);
3678:   }
3679:   /* Build SF */
3680:   CHKMEMQ;
3681:   {
3682:     PetscSF            sfPoint, sfPointSub;
3683:     const PetscSFNode *remotePoints;
3684:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
3685:     const PetscInt    *localPoints;
3686:     PetscInt          *slocalPoints;
3687:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells + numSubFaces + numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3688:     PetscMPIInt        rank;

3690:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
3691:     DMGetPointSF(dm, &sfPoint);
3692:     DMGetPointSF(subdm, &sfPointSub);
3693:     DMPlexGetChart(dm, &pStart, &pEnd);
3694:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3695:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3696:     if (numRoots >= 0) {
3697:       /* Only vertices should be shared */
3698:       PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners);
3699:       for (p = 0; p < pEnd - pStart; ++p) {
3700:         newLocalPoints[p].rank  = -2;
3701:         newLocalPoints[p].index = -2;
3702:       }
3703:       /* Set subleaves */
3704:       for (l = 0; l < numLeaves; ++l) {
3705:         const PetscInt point    = localPoints[l];
3706:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3709:         if (subPoint < 0) continue;
3710:         newLocalPoints[point - pStart].rank  = rank;
3711:         newLocalPoints[point - pStart].index = subPoint;
3712:         ++numSubLeaves;
3713:       }
3714:       /* Must put in owned subpoints */
3715:       for (p = pStart; p < pEnd; ++p) {
3716:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

3718:         if (subPoint < 0) {
3719:           newOwners[p - pStart].rank  = -3;
3720:           newOwners[p - pStart].index = -3;
3721:         } else {
3722:           newOwners[p - pStart].rank  = rank;
3723:           newOwners[p - pStart].index = subPoint;
3724:         }
3725:       }
3726:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3727:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3728:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3729:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3730:       PetscMalloc1(numSubLeaves, &slocalPoints);
3731:       PetscMalloc1(numSubLeaves, &sremotePoints);
3732:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3733:         const PetscInt point    = localPoints[l];
3734:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

3736:         if (subPoint < 0) continue;
3737:         if (newLocalPoints[point].rank == rank) {
3738:           ++ll;
3739:           continue;
3740:         }
3741:         slocalPoints[sl]        = subPoint;
3742:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3743:         sremotePoints[sl].index = newLocalPoints[point].index;
3746:         ++sl;
3747:       }
3748:       PetscFree2(newLocalPoints, newOwners);
3750:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3751:     }
3752:   }
3753:   CHKMEMQ;
3754:   /* Cleanup */
3755:   if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
3756:   ISDestroy(&subvertexIS);
3757:   PetscFree(subCells);
3758:   return 0;
3759: }

3761: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3762: {
3763:   DMLabel label = NULL;

3765:   if (labelname) DMGetLabel(dm, labelname, &label);
3766:   DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3767:   return 0;
3768: }

3770: /*@C
3771:   DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label can be given to restrict the cells.

3773:   Input Parameters:
3774: + dm          - The original mesh
3775: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3776: . label       - A label name, or NULL
3777: - value  - A label value

3779:   Output Parameter:
3780: . subdm - The surface mesh

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

3784:   Level: developer

3786: .seealso: `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()`
3787: @*/
3788: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3789: {
3790:   PetscInt dim, cdim, depth;

3794:   DMGetDimension(dm, &dim);
3795:   DMPlexGetDepth(dm, &depth);
3796:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3797:   DMSetType(*subdm, DMPLEX);
3798:   DMSetDimension(*subdm, dim - 1);
3799:   DMGetCoordinateDim(dm, &cdim);
3800:   DMSetCoordinateDim(*subdm, cdim);
3801:   if (depth == dim) {
3802:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3803:   } else {
3804:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3805:   }
3806:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3807:   return 0;
3808: }

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

3813:   Input Parameters:
3814: + dm        - The original mesh
3815: . cellLabel - The DMLabel marking cells contained in the new mesh
3816: - value     - The label value to use

3818:   Output Parameter:
3819: . subdm - The new mesh

3821:   Notes:
3822:   This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`.

3824:   Level: developer

3826: .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`, `DMPlexCreateSubmesh()`
3827: @*/
3828: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3829: {
3830:   PetscInt dim, overlap;

3834:   DMGetDimension(dm, &dim);
3835:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3836:   DMSetType(*subdm, DMPLEX);
3837:   /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3838:   DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3839:   DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3840:   // It is possible to obtain a surface mesh where some faces are in SF
3841:   //   We should either mark the mesh as having an overlap, or delete these from the SF
3842:   DMPlexGetOverlap(dm, &overlap);
3843:   if (!overlap) {
3844:     PetscSF         sf;
3845:     const PetscInt *leaves;
3846:     PetscInt        cStart, cEnd, Nl;
3847:     PetscBool       hasSubcell = PETSC_FALSE, ghasSubcell;

3849:     DMPlexGetHeightStratum(*subdm, 0, &cStart, &cEnd);
3850:     DMGetPointSF(*subdm, &sf);
3851:     PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL);
3852:     for (PetscInt l = 0; l < Nl; ++l) {
3853:       const PetscInt point = leaves ? leaves[l] : l;

3855:       if (point >= cStart && point < cEnd) {
3856:         hasSubcell = PETSC_TRUE;
3857:         break;
3858:       }
3859:     }
3860:     MPIU_Allreduce(&hasSubcell, &ghasSubcell, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
3861:     if (ghasSubcell) DMPlexSetOverlap(*subdm, NULL, 1);
3862:   }
3863:   return 0;
3864: }

3866: /*@
3867:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3869:   Input Parameter:
3870: . dm - The submesh DM

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

3875:   Level: developer

3877: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
3878: @*/
3879: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3880: {
3883:   *subpointMap = ((DM_Plex *)dm->data)->subpointMap;
3884:   return 0;
3885: }

3887: /*@
3888:   DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values

3890:   Input Parameters:
3891: + dm - The submesh DM
3892: - subpointMap - The DMLabel of all the points from the original mesh in this submesh

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

3896:   Level: developer

3898: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
3899: @*/
3900: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3901: {
3902:   DM_Plex *mesh = (DM_Plex *)dm->data;
3903:   DMLabel  tmp;

3906:   tmp               = mesh->subpointMap;
3907:   mesh->subpointMap = subpointMap;
3908:   PetscObjectReference((PetscObject)mesh->subpointMap);
3909:   DMLabelDestroy(&tmp);
3910:   return 0;
3911: }

3913: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
3914: {
3915:   DMLabel  spmap;
3916:   PetscInt depth, d;

3918:   DMPlexGetSubpointMap(dm, &spmap);
3919:   DMPlexGetDepth(dm, &depth);
3920:   if (spmap && depth >= 0) {
3921:     DM_Plex  *mesh = (DM_Plex *)dm->data;
3922:     PetscInt *points, *depths;
3923:     PetscInt  pStart, pEnd, p, off;

3925:     DMPlexGetChart(dm, &pStart, &pEnd);
3927:     PetscMalloc1(pEnd, &points);
3928:     DMGetWorkArray(dm, depth + 1, MPIU_INT, &depths);
3929:     depths[0] = depth;
3930:     depths[1] = 0;
3931:     for (d = 2; d <= depth; ++d) depths[d] = depth + 1 - d;
3932:     for (d = 0, off = 0; d <= depth; ++d) {
3933:       const PetscInt dep = depths[d];
3934:       PetscInt       depStart, depEnd, n;

3936:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3937:       DMLabelGetStratumSize(spmap, dep, &n);
3938:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3940:       } else {
3941:         if (!n) {
3942:           if (d == 0) {
3943:             /* Missing cells */
3944:             for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = -1;
3945:           } else {
3946:             /* Missing faces */
3947:             for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3948:           }
3949:         }
3950:       }
3951:       if (n) {
3952:         IS              is;
3953:         const PetscInt *opoints;

3955:         DMLabelGetStratumIS(spmap, dep, &is);
3956:         ISGetIndices(is, &opoints);
3957:         for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3958:         ISRestoreIndices(is, &opoints);
3959:         ISDestroy(&is);
3960:       }
3961:     }
3962:     DMRestoreWorkArray(dm, depth + 1, MPIU_INT, &depths);
3964:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3965:     PetscObjectStateGet((PetscObject)spmap, &mesh->subpointState);
3966:   }
3967:   return 0;
3968: }

3970: /*@
3971:   DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data

3973:   Input Parameter:
3974: . dm - The submesh DM

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

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

3981:   Level: developer

3983: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()`
3984: @*/
3985: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
3986: {
3987:   DM_Plex         *mesh = (DM_Plex *)dm->data;
3988:   DMLabel          spmap;
3989:   PetscObjectState state;

3993:   DMPlexGetSubpointMap(dm, &spmap);
3994:   PetscObjectStateGet((PetscObject)spmap, &state);
3995:   if (state != mesh->subpointState || !mesh->subpointIS) DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);
3996:   *subpointIS = mesh->subpointIS;
3997:   return 0;
3998: }

4000: /*@
4001:   DMGetEnclosureRelation - Get the relationship between dmA and dmB

4003:   Input Parameters:
4004: + dmA - The first DM
4005: - dmB - The second DM

4007:   Output Parameter:
4008: . rel - The relation of dmA to dmB

4010:   Level: intermediate

4012: .seealso: `DMGetEnclosurePoint()`
4013: @*/
4014: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
4015: {
4016:   DM       plexA, plexB, sdm;
4017:   DMLabel  spmap;
4018:   PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;

4021:   *rel = DM_ENC_NONE;
4022:   if (!dmA || !dmB) return 0;
4025:   if (dmA == dmB) {
4026:     *rel = DM_ENC_EQUALITY;
4027:     return 0;
4028:   }
4029:   DMConvert(dmA, DMPLEX, &plexA);
4030:   DMConvert(dmB, DMPLEX, &plexB);
4031:   DMPlexGetChart(plexA, &pStartA, &pEndA);
4032:   DMPlexGetChart(plexB, &pStartB, &pEndB);
4033:   /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
4034:     - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
4035:   if ((pStartA == pStartB) && (pEndA == pEndB)) {
4036:     *rel = DM_ENC_EQUALITY;
4037:     goto end;
4038:   }
4039:   NpA = pEndA - pStartA;
4040:   NpB = pEndB - pStartB;
4041:   if (NpA == NpB) goto end;
4042:   sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
4043:   DMPlexGetSubpointMap(sdm, &spmap);
4044:   if (!spmap) goto end;
4045:   /* TODO Check the space mapped to by subpointMap is same size as dm */
4046:   if (NpA > NpB) {
4047:     *rel = DM_ENC_SUPERMESH;
4048:   } else {
4049:     *rel = DM_ENC_SUBMESH;
4050:   }
4051: end:
4052:   DMDestroy(&plexA);
4053:   DMDestroy(&plexB);
4054:   return 0;
4055: }

4057: /*@
4058:   DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB

4060:   Input Parameters:
4061: + dmA   - The first DM
4062: . dmB   - The second DM
4063: . etype - The type of enclosure relation that dmA has to dmB
4064: - pB    - A point of dmB

4066:   Output Parameter:
4067: . pA    - The corresponding point of dmA

4069:   Level: intermediate

4071: .seealso: `DMGetEnclosureRelation()`
4072: @*/
4073: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
4074: {
4075:   DM              sdm;
4076:   IS              subpointIS;
4077:   const PetscInt *subpoints;
4078:   PetscInt        numSubpoints;

4080:   /* TODO Cache the IS, making it look like an index */
4081:   switch (etype) {
4082:   case DM_ENC_SUPERMESH:
4083:     sdm = dmB;
4084:     DMPlexGetSubpointIS(sdm, &subpointIS);
4085:     ISGetIndices(subpointIS, &subpoints);
4086:     *pA = subpoints[pB];
4087:     ISRestoreIndices(subpointIS, &subpoints);
4088:     break;
4089:   case DM_ENC_SUBMESH:
4090:     sdm = dmA;
4091:     DMPlexGetSubpointIS(sdm, &subpointIS);
4092:     ISGetLocalSize(subpointIS, &numSubpoints);
4093:     ISGetIndices(subpointIS, &subpoints);
4094:     PetscFindInt(pB, numSubpoints, subpoints, pA);
4095:     if (*pA < 0) {
4096:       DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
4097:       DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
4098:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB);
4099:     }
4100:     ISRestoreIndices(subpointIS, &subpoints);
4101:     break;
4102:   case DM_ENC_EQUALITY:
4103:   case DM_ENC_NONE:
4104:     *pA = pB;
4105:     break;
4106:   case DM_ENC_UNKNOWN: {
4107:     DMEnclosureType enc;

4109:     DMGetEnclosureRelation(dmA, dmB, &enc);
4110:     DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
4111:   } break;
4112:   default:
4113:     SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int)etype);
4114:   }
4115:   return 0;
4116: }