Actual source code: plexsubmesh.c

petsc-master 2015-01-29
Report Typos and Errors
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petscsf.h>

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

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

 17:     DMPlexGetSupportSize(dm, f, &supportSize);
 18:     if (supportSize == 1) {DMLabelSetValue(label, f, 1);}
 19:   }
 20:   return(0);
 21: }

 25: /*@
 26:   DMPlexMarkBoundaryFaces - Mark all faces on the boundary

 28:   Not Collective

 30:   Input Parameter:
 31: . dm - The original DM

 33:   Output Parameter:
 34: . label - The DMLabel marking boundary faces with value 1

 36:   Level: developer

 38: .seealso: DMLabelCreate(), DMPlexCreateLabel()
 39: @*/
 40: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
 41: {

 45:   DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
 46:   return(0);
 47: }

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

 54:   Input Parameters:
 55: + dm - The DM
 56: - label - A DMLabel marking the surface points

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

 61:   Level: developer

 63: .seealso: DMPlexLabelCohesiveComplete()
 64: @*/
 65: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
 66: {
 67:   IS              valueIS;
 68:   const PetscInt *values;
 69:   PetscInt        numValues, v;
 70:   PetscErrorCode  ierr;

 73:   DMLabelGetNumValues(label, &numValues);
 74:   DMLabelGetValueIS(label, &valueIS);
 75:   ISGetIndices(valueIS, &values);
 76:   for (v = 0; v < numValues; ++v) {
 77:     IS              pointIS;
 78:     const PetscInt *points;
 79:     PetscInt        numPoints, p;

 81:     DMLabelGetStratumSize(label, values[v], &numPoints);
 82:     DMLabelGetStratumIS(label, values[v], &pointIS);
 83:     ISGetIndices(pointIS, &points);
 84:     for (p = 0; p < numPoints; ++p) {
 85:       PetscInt *closure = NULL;
 86:       PetscInt  closureSize, c;

 88:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 89:       for (c = 0; c < closureSize*2; c += 2) {
 90:         DMLabelSetValue(label, closure[c], values[v]);
 91:       }
 92:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
 93:     }
 94:     ISRestoreIndices(pointIS, &points);
 95:     ISDestroy(&pointIS);
 96:   }
 97:   ISRestoreIndices(valueIS, &values);
 98:   ISDestroy(&valueIS);
 99:   return(0);
100: }

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

107:   Input Parameters:
108: + dm - The DM
109: - label - A DMLabel marking the surface points

111:   Output Parameter:
112: . label - A DMLabel incorporating cells

114:   Level: developer

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

118: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
119: @*/
120: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
121: {
122:   IS              valueIS;
123:   const PetscInt *values;
124:   PetscInt        numValues, v, cStart, cEnd;
125:   PetscErrorCode  ierr;

128:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
129:   DMLabelGetNumValues(label, &numValues);
130:   DMLabelGetValueIS(label, &valueIS);
131:   ISGetIndices(valueIS, &values);
132:   for (v = 0; v < numValues; ++v) {
133:     IS              pointIS;
134:     const PetscInt *points;
135:     PetscInt        numPoints, p;

137:     DMLabelGetStratumSize(label, values[v], &numPoints);
138:     DMLabelGetStratumIS(label, values[v], &pointIS);
139:     ISGetIndices(pointIS, &points);
140:     for (p = 0; p < numPoints; ++p) {
141:       PetscInt *closure = NULL;
142:       PetscInt  closureSize, point;

144:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
145:       point = closure[(closureSize-1)*2];
146:       if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]);}
147:       DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
148:     }
149:     ISRestoreIndices(pointIS, &points);
150:     ISDestroy(&pointIS);
151:   }
152:   ISRestoreIndices(valueIS, &values);
153:   ISDestroy(&valueIS);
154:   return(0);
155: }

159: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthMax[], PetscInt depthEnd[], PetscInt depthShift[])
160: {
161:   if (depth < 0) return p;
162:   /* Normal Cells                 */ if (p < depthMax[depth])                return p;
163:   /* Hybrid Cells+Normal Vertices */ if (p < depthMax[0])                    return p + depthShift[depth];
164:   /* Hybrid Vertices+Normal Faces */ if (depth < 2 || p < depthMax[depth-1]) return p + depthShift[depth] + depthShift[0];
165:   /* Hybrid Faces+Normal Edges    */ if (depth < 3 || p < depthMax[depth-2]) return p + depthShift[depth] + depthShift[0] + depthShift[depth-1];
166:   /* Hybrid Edges                 */                                         return p + depthShift[depth] + depthShift[0] + depthShift[depth-1] + depthShift[depth-2];
167: }

171: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
172: {
173:   PetscInt      *depthMax, *depthEnd;
174:   PetscInt       depth = 0, d, pStart, pEnd, p;

178:   DMPlexGetDepth(dm, &depth);
179:   if (depth < 0) return(0);
180:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
181:   /* Step 1: Expand chart */
182:   DMPlexGetChart(dm, &pStart, &pEnd);
183:   DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
184:   for (d = 0; d <= depth; ++d) {
185:     pEnd += depthShift[d];
186:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
187:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
188:   }
189:   DMPlexSetChart(dmNew, pStart, pEnd);
190:   /* Step 2: Set cone and support sizes */
191:   for (d = 0; d <= depth; ++d) {
192:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
193:     for (p = pStart; p < pEnd; ++p) {
194:       PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);
195:       PetscInt size;

197:       DMPlexGetConeSize(dm, p, &size);
198:       DMPlexSetConeSize(dmNew, newp, size);
199:       DMPlexGetSupportSize(dm, p, &size);
200:       DMPlexSetSupportSize(dmNew, newp, size);
201:     }
202:   }
203:   PetscFree2(depthMax,depthEnd);
204:   return(0);
205: }

209: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
210: {
211:   PetscInt      *depthEnd, *depthMax, *newpoints;
212:   PetscInt       depth = 0, d, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;

216:   DMPlexGetDepth(dm, &depth);
217:   if (depth < 0) return(0);
218:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
219:   DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
220:   PetscMalloc3(depth+1,&depthEnd,depth+1,&depthMax,PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
221:   DMPlexGetHybridBounds(dm, &depthMax[depth], depth > 0 ? &depthMax[depth-1] : NULL, &depthMax[1], &depthMax[0]);
222:   for (d = 0; d <= depth; ++d) {
223:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
224:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
225:   }
226:   /* Step 5: Set cones and supports */
227:   DMPlexGetChart(dm, &pStart, &pEnd);
228:   for (p = pStart; p < pEnd; ++p) {
229:     const PetscInt *points = NULL, *orientations = NULL;
230:     PetscInt        size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthMax, depthEnd, depthShift);

232:     DMPlexGetConeSize(dm, p, &size);
233:     DMPlexGetCone(dm, p, &points);
234:     DMPlexGetConeOrientation(dm, p, &orientations);
235:     for (i = 0; i < size; ++i) {
236:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
237:     }
238:     DMPlexSetCone(dmNew, newp, newpoints);
239:     DMPlexSetConeOrientation(dmNew, newp, orientations);
240:     DMPlexGetSupportSize(dm, p, &size);
241:     DMPlexGetSupportSize(dmNew, newp, &sizeNew);
242:     DMPlexGetSupport(dm, p, &points);
243:     for (i = 0; i < size; ++i) {
244:       newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthMax, depthEnd, depthShift);
245:     }
246:     for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
247:     DMPlexSetSupport(dmNew, newp, newpoints);
248:   }
249:   PetscFree3(depthEnd,depthMax,newpoints);
250:   return(0);
251: }

255: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
256: {
257:   PetscSection   coordSection, newCoordSection;
258:   Vec            coordinates, newCoordinates;
259:   PetscScalar   *coords, *newCoords;
260:   PetscInt      *depthEnd, coordSize;
261:   PetscInt       dim, depth = 0, d, vStart, vEnd, vStartNew, vEndNew, v;

265:   DMGetDimension(dm, &dim);
266:   DMPlexGetDepth(dm, &depth);
267:   PetscMalloc1(depth+1, &depthEnd);
268:   for (d = 0; d <= depth; ++d) {
269:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
270:   }
271:   /* Step 8: Convert coordinates */
272:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
273:   DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
274:   DMGetCoordinateSection(dm, &coordSection);
275:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
276:   PetscSectionSetNumFields(newCoordSection, 1);
277:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
278:   PetscSectionSetChart(newCoordSection, vStartNew, vEndNew);
279:   for (v = vStartNew; v < vEndNew; ++v) {
280:     PetscSectionSetDof(newCoordSection, v, dim);
281:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
282:   }
283:   PetscSectionSetUp(newCoordSection);
284:   DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
285:   PetscSectionGetStorageSize(newCoordSection, &coordSize);
286:   VecCreate(PetscObjectComm((PetscObject)dm), &newCoordinates);
287:   PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
288:   VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
289:   VecSetBlockSize(newCoordinates, dim);
290:   VecSetType(newCoordinates,VECSTANDARD);
291:   DMSetCoordinatesLocal(dmNew, newCoordinates);
292:   DMGetCoordinatesLocal(dm, &coordinates);
293:   VecGetArray(coordinates, &coords);
294:   VecGetArray(newCoordinates, &newCoords);
295:   for (v = vStart; v < vEnd; ++v) {
296:     PetscInt dof, off, noff, d;

298:     PetscSectionGetDof(coordSection, v, &dof);
299:     PetscSectionGetOffset(coordSection, v, &off);
300:     PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthEnd, depthEnd, depthShift), &noff);
301:     for (d = 0; d < dof; ++d) {
302:       newCoords[noff+d] = coords[off+d];
303:     }
304:   }
305:   VecRestoreArray(coordinates, &coords);
306:   VecRestoreArray(newCoordinates, &newCoords);
307:   VecDestroy(&newCoordinates);
308:   PetscSectionDestroy(&newCoordSection);
309:   PetscFree(depthEnd);
310:   return(0);
311: }

315: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
316: {
317:   PetscInt          *depthMax, *depthEnd;
318:   PetscInt           depth = 0, d;
319:   PetscSF            sfPoint, sfPointNew;
320:   const PetscSFNode *remotePoints;
321:   PetscSFNode       *gremotePoints;
322:   const PetscInt    *localPoints;
323:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
324:   PetscInt           numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
325:   PetscErrorCode     ierr;

328:   DMPlexGetDepth(dm, &depth);
329:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
330:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
331:   for (d = 0; d <= depth; ++d) {
332:     totShift += depthShift[d];
333:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
334:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
335:   }
336:   /* Step 9: Convert pointSF */
337:   DMGetPointSF(dm, &sfPoint);
338:   DMGetPointSF(dmNew, &sfPointNew);
339:   DMPlexGetChart(dm, &pStart, &pEnd);
340:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
341:   if (numRoots >= 0) {
342:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
343:     for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthMax, depthEnd, depthShift);
344:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
345:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
346:     PetscMalloc1(numLeaves,    &glocalPoints);
347:     PetscMalloc1(numLeaves, &gremotePoints);
348:     for (l = 0; l < numLeaves; ++l) {
349:       glocalPoints[l]        = DMPlexShiftPoint_Internal(localPoints[l], depth, depthMax, depthEnd, depthShift);
350:       gremotePoints[l].rank  = remotePoints[l].rank;
351:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
352:     }
353:     PetscFree2(newLocation,newRemoteLocation);
354:     PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
355:   }
356:   PetscFree2(depthMax,depthEnd);
357:   return(0);
358: }

362: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
363: {
364:   PetscSF            sfPoint;
365:   DMLabel            vtkLabel, ghostLabel;
366:   PetscInt          *depthMax, *depthEnd;
367:   const PetscSFNode *leafRemote;
368:   const PetscInt    *leafLocal;
369:   PetscInt           depth = 0, d, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
370:   PetscMPIInt        rank;
371:   PetscErrorCode     ierr;

374:   DMPlexGetDepth(dm, &depth);
375:   PetscMalloc2(depth+1,&depthMax,depth+1,&depthEnd);
376:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
377:   for (d = 0; d <= depth; ++d) {
378:     DMPlexGetDepthStratum(dm, d, NULL, &depthEnd[d]);
379:     depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
380:   }
381:   /* Step 10: Convert labels */
382:   DMPlexGetNumLabels(dm, &numLabels);
383:   for (l = 0; l < numLabels; ++l) {
384:     DMLabel         label, newlabel;
385:     const char     *lname;
386:     PetscBool       isDepth;
387:     IS              valueIS;
388:     const PetscInt *values;
389:     PetscInt        numValues, val;

391:     DMPlexGetLabelName(dm, l, &lname);
392:     PetscStrcmp(lname, "depth", &isDepth);
393:     if (isDepth) continue;
394:     DMPlexCreateLabel(dmNew, lname);
395:     DMPlexGetLabel(dm, lname, &label);
396:     DMPlexGetLabel(dmNew, lname, &newlabel);
397:     DMLabelGetValueIS(label, &valueIS);
398:     ISGetLocalSize(valueIS, &numValues);
399:     ISGetIndices(valueIS, &values);
400:     for (val = 0; val < numValues; ++val) {
401:       IS              pointIS;
402:       const PetscInt *points;
403:       PetscInt        numPoints, p;

405:       DMLabelGetStratumIS(label, values[val], &pointIS);
406:       ISGetLocalSize(pointIS, &numPoints);
407:       ISGetIndices(pointIS, &points);
408:       for (p = 0; p < numPoints; ++p) {
409:         const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthMax, depthEnd, depthShift);

411:         DMLabelSetValue(newlabel, newpoint, values[val]);
412:       }
413:       ISRestoreIndices(pointIS, &points);
414:       ISDestroy(&pointIS);
415:     }
416:     ISRestoreIndices(valueIS, &values);
417:     ISDestroy(&valueIS);
418:   }
419:   PetscFree2(depthMax,depthEnd);
420:   /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
421:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
422:   DMGetPointSF(dm, &sfPoint);
423:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
424:   PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
425:   DMPlexCreateLabel(dmNew, "vtk");
426:   DMPlexCreateLabel(dmNew, "ghost");
427:   DMPlexGetLabel(dmNew, "vtk", &vtkLabel);
428:   DMPlexGetLabel(dmNew, "ghost", &ghostLabel);
429:   for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
430:     for (; c < leafLocal[l] && c < cEnd; ++c) {
431:       DMLabelSetValue(vtkLabel, c, 1);
432:     }
433:     if (leafLocal[l] >= cEnd) break;
434:     if (leafRemote[l].rank == rank) {
435:       DMLabelSetValue(vtkLabel, c, 1);
436:     } else {
437:       DMLabelSetValue(ghostLabel, c, 2);
438:     }
439:   }
440:   for (; c < cEnd; ++c) {
441:     DMLabelSetValue(vtkLabel, c, 1);
442:   }
443:   if (0) {
444:     PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
445:     DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
446:     PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
447:   }
448:   DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
449:   for (f = fStart; f < fEnd; ++f) {
450:     PetscInt numCells;

452:     DMPlexGetSupportSize(dmNew, f, &numCells);
453:     if (numCells < 2) {
454:       DMLabelSetValue(ghostLabel, f, 1);
455:     } else {
456:       const PetscInt *cells = NULL;
457:       PetscInt        vA, vB;

459:       DMPlexGetSupport(dmNew, f, &cells);
460:       DMLabelGetValue(vtkLabel, cells[0], &vA);
461:       DMLabelGetValue(vtkLabel, cells[1], &vB);
462:       if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
463:     }
464:   }
465:   if (0) {
466:     PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
467:     DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
468:     PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
469:   }
470:   return(0);
471: }

475: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
476: {
477:   PetscSF         sf;
478:   IS              valueIS;
479:   const PetscInt *values, *leaves;
480:   PetscInt       *depthShift;
481:   PetscInt        depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
482:   PetscErrorCode  ierr;

485:   DMGetPointSF(dm, &sf);
486:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
487:   nleaves = PetscMax(0, nleaves);
488:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
489:   /* Count ghost cells */
490:   DMLabelGetValueIS(label, &valueIS);
491:   ISGetLocalSize(valueIS, &numFS);
492:   ISGetIndices(valueIS, &values);
493:   Ng   = 0;
494:   for (fs = 0; fs < numFS; ++fs) {
495:     IS              faceIS;
496:     const PetscInt *faces;
497:     PetscInt        numFaces, f, numBdFaces = 0;

499:     DMLabelGetStratumIS(label, values[fs], &faceIS);
500:     ISGetLocalSize(faceIS, &numFaces);
501:     ISGetIndices(faceIS, &faces);
502:     for (f = 0; f < numFaces; ++f) {
503:       PetscFindInt(faces[f], nleaves, leaves, &loc);
504:       if (loc >= 0) continue;
505:       if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
506:     }
507:     Ng += numBdFaces;
508:     ISDestroy(&faceIS);
509:   }
510:   DMPlexGetDepth(dm, &depth);
511:   PetscMalloc1(depth+1, &depthShift);
512:   PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));
513:   if (depth >= 0) depthShift[depth] = Ng;
514:   DMPlexShiftSizes_Internal(dm, depthShift, gdm);
515:   /* Step 3: Set cone/support sizes for new points */
516:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
517:   DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
518:   for (c = cEnd; c < cEnd + Ng; ++c) {
519:     DMPlexSetConeSize(gdm, c, 1);
520:   }
521:   for (fs = 0; fs < numFS; ++fs) {
522:     IS              faceIS;
523:     const PetscInt *faces;
524:     PetscInt        numFaces, f;

526:     DMLabelGetStratumIS(label, values[fs], &faceIS);
527:     ISGetLocalSize(faceIS, &numFaces);
528:     ISGetIndices(faceIS, &faces);
529:     for (f = 0; f < numFaces; ++f) {
530:       PetscInt size;

532:       PetscFindInt(faces[f], nleaves, leaves, &loc);
533:       if (loc >= 0) continue;
534:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
535:       DMPlexGetSupportSize(dm, faces[f], &size);
536:       if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
537:       DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
538:     }
539:     ISRestoreIndices(faceIS, &faces);
540:     ISDestroy(&faceIS);
541:   }
542:   /* Step 4: Setup ghosted DM */
543:   DMSetUp(gdm);
544:   DMPlexShiftPoints_Internal(dm, depthShift, gdm);
545:   /* Step 6: Set cones and supports for new points */
546:   ghostCell = cEnd;
547:   for (fs = 0; fs < numFS; ++fs) {
548:     IS              faceIS;
549:     const PetscInt *faces;
550:     PetscInt        numFaces, f;

552:     DMLabelGetStratumIS(label, values[fs], &faceIS);
553:     ISGetLocalSize(faceIS, &numFaces);
554:     ISGetIndices(faceIS, &faces);
555:     for (f = 0; f < numFaces; ++f) {
556:       PetscInt newFace = faces[f] + Ng;

558:       PetscFindInt(faces[f], nleaves, leaves, &loc);
559:       if (loc >= 0) continue;
560:       if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
561:       DMPlexSetCone(gdm, ghostCell, &newFace);
562:       DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
563:       ++ghostCell;
564:     }
565:     ISRestoreIndices(faceIS, &faces);
566:     ISDestroy(&faceIS);
567:   }
568:   ISRestoreIndices(valueIS, &values);
569:   ISDestroy(&valueIS);
570:   /* Step 7: Stratify */
571:   DMPlexStratify(gdm);
572:   DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
573:   DMPlexShiftSF_Internal(dm, depthShift, gdm);
574:   DMPlexShiftLabels_Internal(dm, depthShift, gdm);
575:   PetscFree(depthShift);
576:   /* Step 7: Periodicity */
577:   if (dm->maxCell) {
578:     const PetscReal *maxCell, *L;
579:     DMGetPeriodicity(dm,  &maxCell, &L);
580:     DMSetPeriodicity(gdm,  maxCell,  L);
581:   }
582:   if (numGhostCells) *numGhostCells = Ng;
583:   return(0);
584: }

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

591:   Collective on dm

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

597:   Output Parameters:
598: + numGhostCells - The number of ghost cells added to the DM
599: - dmGhosted - The new DM

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

603:   Level: developer

605: .seealso: DMCreate()
606: @*/
607: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
608: {
609:   DM             gdm;
610:   DMLabel        label;
611:   const char    *name = labelName ? labelName : "Face Sets";
612:   PetscInt       dim;
613:   PetscBool      flag;

620:   DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
621:   DMSetType(gdm, DMPLEX);
622:   DMGetDimension(dm, &dim);
623:   DMSetDimension(gdm, dim);
624:   DMPlexGetAdjacencyUseCone(dm, &flag);
625:   DMPlexSetAdjacencyUseCone(gdm, flag);
626:   DMPlexGetAdjacencyUseClosure(dm, &flag);
627:   DMPlexSetAdjacencyUseClosure(gdm, flag);
628:   DMPlexGetLabel(dm, name, &label);
629:   if (!label) {
630:     /* Get label for boundary faces */
631:     DMPlexCreateLabel(dm, name);
632:     DMPlexGetLabel(dm, name, &label);
633:     DMPlexMarkBoundaryFaces(dm, label);
634:   }
635:   DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
636:   DMPlexCopyBoundary(dm, gdm);
637:   *dmGhosted = gdm;
638:   return(0);
639: }

643: /*
644:   We are adding three kinds of points here:
645:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
646:     Non-replicated: Points which exist on the fault, but are not replicated
647:     Hybrid:         Entirely new points, such as cohesive cells

649:   When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
650:   each depth so that the new split/hybrid points can be inserted as a block.
651: */
652: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
653: {
654:   MPI_Comm         comm;
655:   IS               valueIS;
656:   PetscInt         numSP = 0;       /* The number of depths for which we have replicated points */
657:   const PetscInt  *values;          /* List of depths for which we have replicated points */
658:   IS              *splitIS;
659:   IS              *unsplitIS;
660:   PetscInt        *numSplitPoints;     /* The number of replicated points at each depth */
661:   PetscInt        *numUnsplitPoints;   /* The number of non-replicated points at each depth which still give rise to hybrid points */
662:   PetscInt        *numHybridPoints;    /* The number of new hybrid points at each depth */
663:   PetscInt        *numHybridPointsOld; /* The number of existing hybrid points at each depth */
664:   const PetscInt **splitPoints;        /* Replicated points for each depth */
665:   const PetscInt **unsplitPoints;      /* Non-replicated points for each depth */
666:   PetscSection     coordSection;
667:   Vec              coordinates;
668:   PetscScalar     *coords;
669:   PetscInt         depths[4];          /* Depths in the order that plex points are numbered */
670:   PetscInt        *depthMax;           /* The first hybrid point at each depth in the original mesh */
671:   PetscInt        *depthEnd;           /* The point limit at each depth in the original mesh */
672:   PetscInt        *depthShift;         /* Number of replicated+hybrid points at each depth */
673:   PetscInt        *depthOffset;        /* Prefix sums of depthShift */
674:   PetscInt        *pMaxNew;            /* The first replicated point at each depth in the new mesh, hybrids come after this */
675:   PetscInt        *coneNew, *coneONew, *supportNew;
676:   PetscInt         shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
677:   PetscErrorCode   ierr;

680:   PetscObjectGetComm((PetscObject)dm,&comm);
681:   DMGetDimension(dm, &dim);
682:   DMPlexGetDepth(dm, &depth);
683:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
684:   depths[0] = depth;
685:   depths[1] = 0;
686:   depths[2] = depth-1;
687:   depths[3] = 1;
688:   /* Count split points and add cohesive cells */
689:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
690:   PetscMalloc6(depth+1,&depthMax,depth+1,&depthEnd,depth+1,&depthShift,depth+1,&depthOffset,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
691:   PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
692:   PetscMemzero(depthShift,  (depth+1) * sizeof(PetscInt));
693:   PetscMemzero(depthOffset, (depth+1) * sizeof(PetscInt));
694:   DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
695:   for (d = 0; d <= depth; ++d) {
696:     DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
697:     depthEnd[d]           = pMaxNew[d];
698:     depthMax[d]           = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
699:     numSplitPoints[d]     = 0;
700:     numUnsplitPoints[d]   = 0;
701:     numHybridPoints[d]    = 0;
702:     numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
703:     splitPoints[d]        = NULL;
704:     unsplitPoints[d]      = NULL;
705:     splitIS[d]            = NULL;
706:     unsplitIS[d]          = NULL;
707:   }
708:   if (label) {
709:     DMLabelGetValueIS(label, &valueIS);
710:     ISGetLocalSize(valueIS, &numSP);
711:     ISGetIndices(valueIS, &values);
712:   }
713:   for (sp = 0; sp < numSP; ++sp) {
714:     const PetscInt dep = values[sp];

716:     if ((dep < 0) || (dep > depth)) continue;
717:     DMLabelGetStratumIS(label, dep, &splitIS[dep]);
718:     if (splitIS[dep]) {
719:       ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
720:       ISGetIndices(splitIS[dep], &splitPoints[dep]);
721:     }
722:     DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
723:     if (unsplitIS[dep]) {
724:       ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
725:       ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
726:     }
727:   }
728:   /* Calculate number of hybrid points */
729:   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   */
730:   for (d = 0; d <= depth; ++d) depthShift[d]          = numSplitPoints[d] + numHybridPoints[d];
731:   for (d = 1; d <= depth; ++d) depthOffset[depths[d]] = depthOffset[depths[d-1]] + depthShift[depths[d-1]];
732:   for (d = 0; d <= depth; ++d) pMaxNew[d]            += depthOffset[d] - numHybridPointsOld[d];
733:   DMPlexShiftSizes_Internal(dm, depthShift, sdm);
734:   /* Step 3: Set cone/support sizes for new points */
735:   for (dep = 0; dep <= depth; ++dep) {
736:     for (p = 0; p < numSplitPoints[dep]; ++p) {
737:       const PetscInt  oldp   = splitPoints[dep][p];
738:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
739:       const PetscInt  splitp = p    + pMaxNew[dep];
740:       const PetscInt *support;
741:       PetscInt        coneSize, supportSize, qf, qn, qp, e;

743:       DMPlexGetConeSize(dm, oldp, &coneSize);
744:       DMPlexSetConeSize(sdm, splitp, coneSize);
745:       DMPlexGetSupportSize(dm, oldp, &supportSize);
746:       DMPlexSetSupportSize(sdm, splitp, supportSize);
747:       if (dep == depth-1) {
748:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

755:         DMPlexGetSupport(dm, oldp, &support);
756:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
757:           PetscInt val;

759:           DMLabelGetValue(label, support[e], &val);
760:           if (val == 1) ++qf;
761:           if ((val == 1) || (val ==  (shift + 1))) ++qn;
762:           if ((val == 1) || (val == -(shift + 1))) ++qp;
763:         }
764:         /* Split old vertex: Edges into original vertex and new cohesive edge */
765:         DMPlexSetSupportSize(sdm, newp, qn+1);
766:         /* Split new vertex: Edges into split vertex and new cohesive edge */
767:         DMPlexSetSupportSize(sdm, splitp, qp+1);
768:         /* Add hybrid edge */
769:         DMPlexSetConeSize(sdm, hybedge, 2);
770:         DMPlexSetSupportSize(sdm, hybedge, qf);
771:       } else if (dep == dim-2) {
772:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

774:         DMPlexGetSupport(dm, oldp, &support);
775:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
776:           PetscInt val;

778:           DMLabelGetValue(label, support[e], &val);
779:           if (val == dim-1) ++qf;
780:           if ((val == dim-1) || (val ==  (shift + dim-1))) ++qn;
781:           if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
782:         }
783:         /* Split old edge: Faces into original edge and cohesive face (positive side?) */
784:         DMPlexSetSupportSize(sdm, newp, qn+1);
785:         /* Split new edge: Faces into split edge and cohesive face (negative side?) */
786:         DMPlexSetSupportSize(sdm, splitp, qp+1);
787:         /* Add hybrid face */
788:         DMPlexSetConeSize(sdm, hybface, 4);
789:         DMPlexSetSupportSize(sdm, hybface, qf);
790:       }
791:     }
792:   }
793:   for (dep = 0; dep <= depth; ++dep) {
794:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
795:       const PetscInt  oldp   = unsplitPoints[dep][p];
796:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
797:       const PetscInt *support;
798:       PetscInt        coneSize, supportSize, qf, e, s;

800:       DMPlexGetConeSize(dm, oldp, &coneSize);
801:       DMPlexGetSupportSize(dm, oldp, &supportSize);
802:       DMPlexGetSupport(dm, oldp, &support);
803:       if (dep == 0) {
804:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

806:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
807:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
808:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
809:           if (e >= 0) ++qf;
810:         }
811:         DMPlexSetSupportSize(sdm, newp, qf+2);
812:         /* Add hybrid edge */
813:         DMPlexSetConeSize(sdm, hybedge, 2);
814:         for (e = 0, qf = 0; e < supportSize; ++e) {
815:           PetscInt val;

817:           DMLabelGetValue(label, support[e], &val);
818:           /* Split and unsplit edges produce hybrid faces */
819:           if (val == 1) ++qf;
820:           if (val == (shift2 + 1)) ++qf;
821:         }
822:         DMPlexSetSupportSize(sdm, hybedge, qf);
823:       } else if (dep == dim-2) {
824:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
825:         PetscInt       val;

827:         for (e = 0, qf = 0; e < supportSize; ++e) {
828:           DMLabelGetValue(label, support[e], &val);
829:           if (val == dim-1) qf += 2;
830:           else              ++qf;
831:         }
832:         /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
833:         DMPlexSetSupportSize(sdm, newp, qf+2);
834:         /* Add hybrid face */
835:         for (e = 0, qf = 0; e < supportSize; ++e) {
836:           DMLabelGetValue(label, support[e], &val);
837:           if (val == dim-1) ++qf;
838:         }
839:         DMPlexSetConeSize(sdm, hybface, 4);
840:         DMPlexSetSupportSize(sdm, hybface, qf);
841:       }
842:     }
843:   }
844:   /* Step 4: Setup split DM */
845:   DMSetUp(sdm);
846:   DMPlexShiftPoints_Internal(dm, depthShift, sdm);
847:   DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
848:   PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
849:   /* Step 6: Set cones and supports for new points */
850:   for (dep = 0; dep <= depth; ++dep) {
851:     for (p = 0; p < numSplitPoints[dep]; ++p) {
852:       const PetscInt  oldp   = splitPoints[dep][p];
853:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
854:       const PetscInt  splitp = p    + pMaxNew[dep];
855:       const PetscInt *cone, *support, *ornt;
856:       PetscInt        coneSize, supportSize, q, qf, qn, qp, v, e, s;

858:       DMPlexGetConeSize(dm, oldp, &coneSize);
859:       DMPlexGetCone(dm, oldp, &cone);
860:       DMPlexGetConeOrientation(dm, oldp, &ornt);
861:       DMPlexGetSupportSize(dm, oldp, &supportSize);
862:       DMPlexGetSupport(dm, oldp, &support);
863:       if (dep == depth-1) {
864:         PetscBool       hasUnsplit = PETSC_FALSE;
865:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
866:         const PetscInt *supportF;

868:         /* Split face:       copy in old face to new face to start */
869:         DMPlexGetSupport(sdm, newp,  &supportF);
870:         DMPlexSetSupport(sdm, splitp, supportF);
871:         /* Split old face:   old vertices/edges in cone so no change */
872:         /* Split new face:   new vertices/edges in cone */
873:         for (q = 0; q < coneSize; ++q) {
874:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
875:           if (v < 0) {
876:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
877:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
878:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
879:             hasUnsplit   = PETSC_TRUE;
880:           } else {
881:             coneNew[2+q] = v + pMaxNew[dep-1];
882:             if (dep > 1) {
883:               const PetscInt *econe;
884:               PetscInt        econeSize, r, vs, vu;

886:               DMPlexGetConeSize(dm, cone[q], &econeSize);
887:               DMPlexGetCone(dm, cone[q], &econe);
888:               for (r = 0; r < econeSize; ++r) {
889:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
890:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
891:                 if (vs >= 0) continue;
892:                 if (vu < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", econe[r], dep-2);
893:                 hasUnsplit   = PETSC_TRUE;
894:               }
895:             }
896:           }
897:         }
898:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
899:         DMPlexSetConeOrientation(sdm, splitp, ornt);
900:         /* Face support */
901:         for (s = 0; s < supportSize; ++s) {
902:           PetscInt val;

904:           DMLabelGetValue(label, support[s], &val);
905:           if (val < 0) {
906:             /* Split old face:   Replace negative side cell with cohesive cell */
907:              DMPlexInsertSupport(sdm, newp, s, hybcell);
908:           } else {
909:             /* Split new face:   Replace positive side cell with cohesive cell */
910:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
911:             /* Get orientation for cohesive face */
912:             {
913:               const PetscInt *ncone, *nconeO;
914:               PetscInt        nconeSize, nc;

916:               DMPlexGetConeSize(dm, support[s], &nconeSize);
917:               DMPlexGetCone(dm, support[s], &ncone);
918:               DMPlexGetConeOrientation(dm, support[s], &nconeO);
919:               for (nc = 0; nc < nconeSize; ++nc) {
920:                 if (ncone[nc] == oldp) {
921:                   coneONew[0] = nconeO[nc];
922:                   break;
923:                 }
924:               }
925:               if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
926:             }
927:           }
928:         }
929:         /* Cohesive cell:    Old and new split face, then new cohesive faces */
930:         coneNew[0]  = newp;   /* Extracted negative side orientation above */
931:         coneNew[1]  = splitp;
932:         coneONew[1] = coneONew[0];
933:         for (q = 0; q < coneSize; ++q) {
934:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
935:           if (v < 0) {
936:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
937:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
938:             coneONew[2+q] = 0;
939:           } else {
940:             coneNew[2+q]  = v + pMaxNew[dep] + numSplitPoints[dep];
941:           }
942:           coneONew[2+q] = 0;
943:         }
944:         DMPlexSetCone(sdm, hybcell, coneNew);
945:         DMPlexSetConeOrientation(sdm, hybcell, coneONew);
946:         /* Label the hybrid cells on the boundary of the split */
947:         if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
948:       } else if (dep == 0) {
949:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

955:           DMLabelGetValue(label, support[e], &val);
956:           if ((val == 1) || (val == (shift + 1))) {
957:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
958:           }
959:         }
960:         supportNew[qn] = hybedge;
961:         DMPlexSetSupport(sdm, newp, supportNew);
962:         /* Split new vertex: Edges in new split faces and new cohesive edge */
963:         for (e = 0, qp = 0; e < supportSize; ++e) {
964:           PetscInt val, edge;

966:           DMLabelGetValue(label, support[e], &val);
967:           if (val == 1) {
968:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
969:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
970:             supportNew[qp++] = edge + pMaxNew[dep+1];
971:           } else if (val == -(shift + 1)) {
972:             supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
973:           }
974:         }
975:         supportNew[qp] = hybedge;
976:         DMPlexSetSupport(sdm, splitp, supportNew);
977:         /* Hybrid edge:    Old and new split vertex */
978:         coneNew[0] = newp;
979:         coneNew[1] = splitp;
980:         DMPlexSetCone(sdm, hybedge, coneNew);
981:         for (e = 0, qf = 0; e < supportSize; ++e) {
982:           PetscInt val, edge;

984:           DMLabelGetValue(label, support[e], &val);
985:           if (val == 1) {
986:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
987:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
988:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
989:           }
990:         }
991:         DMPlexSetSupport(sdm, hybedge, supportNew);
992:       } else if (dep == dim-2) {
993:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

995:         /* Split old edge:   old vertices in cone so no change */
996:         /* Split new edge:   new vertices in cone */
997:         for (q = 0; q < coneSize; ++q) {
998:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
999:           if (v < 0) {
1000:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1001:             if (v < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[q], dep-1);
1002:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1003:           } else {
1004:             coneNew[q] = v + pMaxNew[dep-1];
1005:           }
1006:         }
1007:         DMPlexSetCone(sdm, splitp, coneNew);
1008:         /* Split old edge: Faces in positive side cells and old split faces */
1009:         for (e = 0, q = 0; e < supportSize; ++e) {
1010:           PetscInt val;

1012:           DMLabelGetValue(label, support[e], &val);
1013:           if (val == dim-1) {
1014:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1015:           } else if (val == (shift + dim-1)) {
1016:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1017:           }
1018:         }
1019:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1020:         DMPlexSetSupport(sdm, newp, supportNew);
1021:         /* Split new edge: Faces in negative side cells and new split faces */
1022:         for (e = 0, q = 0; e < supportSize; ++e) {
1023:           PetscInt val, face;

1025:           DMLabelGetValue(label, support[e], &val);
1026:           if (val == dim-1) {
1027:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1028:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1029:             supportNew[q++] = face + pMaxNew[dep+1];
1030:           } else if (val == -(shift + dim-1)) {
1031:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1032:           }
1033:         }
1034:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1035:         DMPlexSetSupport(sdm, splitp, supportNew);
1036:         /* Hybrid face */
1037:         coneNew[0] = newp;
1038:         coneNew[1] = splitp;
1039:         for (v = 0; v < coneSize; ++v) {
1040:           PetscInt vertex;
1041:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1042:           if (vertex < 0) {
1043:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1044:             if (vertex < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %d in split or unsplit points of depth %d", cone[v], dep-1);
1045:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1046:           } else {
1047:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1048:           }
1049:         }
1050:         DMPlexSetCone(sdm, hybface, coneNew);
1051:         for (e = 0, qf = 0; e < supportSize; ++e) {
1052:           PetscInt val, face;

1054:           DMLabelGetValue(label, support[e], &val);
1055:           if (val == dim-1) {
1056:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1057:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1058:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1059:           }
1060:         }
1061:         DMPlexSetSupport(sdm, hybface, supportNew);
1062:       }
1063:     }
1064:   }
1065:   for (dep = 0; dep <= depth; ++dep) {
1066:     for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1067:       const PetscInt  oldp   = unsplitPoints[dep][p];
1068:       const PetscInt  newp   = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*oldp + depthOffset[dep]*/;
1069:       const PetscInt *cone, *support, *ornt;
1070:       PetscInt        coneSize, supportSize, supportSizeNew, q, qf, e, f, s;

1072:       DMPlexGetConeSize(dm, oldp, &coneSize);
1073:       DMPlexGetCone(dm, oldp, &cone);
1074:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1075:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1076:       DMPlexGetSupport(dm, oldp, &support);
1077:       if (dep == 0) {
1078:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

1080:         /* Unsplit vertex */
1081:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1082:         for (s = 0, q = 0; s < supportSize; ++s) {
1083:           supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthMax, depthEnd, depthShift) /*support[s] + depthOffset[dep+1]*/;
1084:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1085:           if (e >= 0) {
1086:             supportNew[q++] = e + pMaxNew[dep+1];
1087:           }
1088:         }
1089:         supportNew[q++] = hybedge;
1090:         supportNew[q++] = hybedge;
1091:         if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1092:         DMPlexSetSupport(sdm, newp, supportNew);
1093:         /* Hybrid edge */
1094:         coneNew[0] = newp;
1095:         coneNew[1] = newp;
1096:         DMPlexSetCone(sdm, hybedge, coneNew);
1097:         for (e = 0, qf = 0; e < supportSize; ++e) {
1098:           PetscInt val, edge;

1100:           DMLabelGetValue(label, support[e], &val);
1101:           if (val == 1) {
1102:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1103:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1104:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1105:           } else if  (val ==  (shift2 + 1)) {
1106:             PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1107:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1108:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1109:           }
1110:         }
1111:         DMPlexSetSupport(sdm, hybedge, supportNew);
1112:       } else if (dep == dim-2) {
1113:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

1119:           DMLabelGetValue(label, support[f], &val);
1120:           if (val == dim-1) {
1121:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1122:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1123:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1124:             supportNew[qf++] = face + pMaxNew[dep+1];
1125:           } else {
1126:             supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthMax, depthEnd, depthShift) /*support[f] + depthOffset[dep+1]*/;
1127:           }
1128:         }
1129:         supportNew[qf++] = hybface;
1130:         supportNew[qf++] = hybface;
1131:         DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1132:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1133:         DMPlexSetSupport(sdm, newp, supportNew);
1134:         /* Add hybrid face */
1135:         coneNew[0] = newp;
1136:         coneNew[1] = newp;
1137:         PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1138:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1139:         coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1140:         PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1141:         if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1142:         coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1143:         DMPlexSetCone(sdm, hybface, coneNew);
1144:         for (f = 0, qf = 0; f < supportSize; ++f) {
1145:           PetscInt val, face;

1147:           DMLabelGetValue(label, support[f], &val);
1148:           if (val == dim-1) {
1149:             PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1150:             supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1151:           }
1152:         }
1153:         DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1154:         if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1155:         DMPlexSetSupport(sdm, hybface, supportNew);
1156:       }
1157:     }
1158:   }
1159:   /* Step 6b: Replace split points in negative side cones */
1160:   for (sp = 0; sp < numSP; ++sp) {
1161:     PetscInt        dep = values[sp];
1162:     IS              pIS;
1163:     PetscInt        numPoints;
1164:     const PetscInt *points;

1166:     if (dep >= 0) continue;
1167:     DMLabelGetStratumIS(label, dep, &pIS);
1168:     if (!pIS) continue;
1169:     dep  = -dep - shift;
1170:     ISGetLocalSize(pIS, &numPoints);
1171:     ISGetIndices(pIS, &points);
1172:     for (p = 0; p < numPoints; ++p) {
1173:       const PetscInt  oldp = points[p];
1174:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + oldp*/;
1175:       const PetscInt *cone;
1176:       PetscInt        coneSize, c;
1177:       /* PetscBool       replaced = PETSC_FALSE; */

1179:       /* Negative edge: replace split vertex */
1180:       /* Negative cell: replace split face */
1181:       DMPlexGetConeSize(sdm, newp, &coneSize);
1182:       DMPlexGetCone(sdm, newp, &cone);
1183:       for (c = 0; c < coneSize; ++c) {
1184:         const PetscInt coldp = cone[c] - depthOffset[dep-1];
1185:         PetscInt       csplitp, cp, val;

1187:         DMLabelGetValue(label, coldp, &val);
1188:         if (val == dep-1) {
1189:           PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1190:           if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1191:           csplitp  = pMaxNew[dep-1] + cp;
1192:           DMPlexInsertCone(sdm, newp, c, csplitp);
1193:           /* replaced = PETSC_TRUE; */
1194:         }
1195:       }
1196:       /* Cells with only a vertex or edge on the submesh have no replacement */
1197:       /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1198:     }
1199:     ISRestoreIndices(pIS, &points);
1200:     ISDestroy(&pIS);
1201:   }
1202:   /* Step 7: Stratify */
1203:   DMPlexStratify(sdm);
1204:   /* Step 8: Coordinates */
1205:   DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1206:   DMGetCoordinateSection(sdm, &coordSection);
1207:   DMGetCoordinatesLocal(sdm, &coordinates);
1208:   VecGetArray(coordinates, &coords);
1209:   for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1210:     const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthMax, depthEnd, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1211:     const PetscInt splitp = pMaxNew[0] + v;
1212:     PetscInt       dof, off, soff, d;

1214:     PetscSectionGetDof(coordSection, newp, &dof);
1215:     PetscSectionGetOffset(coordSection, newp, &off);
1216:     PetscSectionGetOffset(coordSection, splitp, &soff);
1217:     for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1218:   }
1219:   VecRestoreArray(coordinates, &coords);
1220:   /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1221:   DMPlexShiftSF_Internal(dm, depthShift, sdm);
1222:   /* Step 10: Labels */
1223:   DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1224:   DMPlexGetNumLabels(sdm, &numLabels);
1225:   for (dep = 0; dep <= depth; ++dep) {
1226:     for (p = 0; p < numSplitPoints[dep]; ++p) {
1227:       const PetscInt newp   = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1228:       const PetscInt splitp = pMaxNew[dep] + p;
1229:       PetscInt       l;

1231:       for (l = 0; l < numLabels; ++l) {
1232:         DMLabel     mlabel;
1233:         const char *lname;
1234:         PetscInt    val;
1235:         PetscBool   isDepth;

1237:         DMPlexGetLabelName(sdm, l, &lname);
1238:         PetscStrcmp(lname, "depth", &isDepth);
1239:         if (isDepth) continue;
1240:         DMPlexGetLabel(sdm, lname, &mlabel);
1241:         DMLabelGetValue(mlabel, newp, &val);
1242:         if (val >= 0) {
1243:           DMLabelSetValue(mlabel, splitp, val);
1244: #if 0
1245:           /* Do not put cohesive edges into the label */
1246:           if (dep == 0) {
1247:             const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1248:             DMLabelSetValue(mlabel, cedge, val);
1249:           } else if (dep == dim-2) {
1250:             const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1251:             DMLabelSetValue(mlabel, cface, val);
1252:           }
1253:           /* Do not put cohesive faces into the label */
1254: #endif
1255:         }
1256:       }
1257:     }
1258:   }
1259:   for (sp = 0; sp < numSP; ++sp) {
1260:     const PetscInt dep = values[sp];

1262:     if ((dep < 0) || (dep > depth)) continue;
1263:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1264:     ISDestroy(&splitIS[dep]);
1265:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1266:     ISDestroy(&unsplitIS[dep]);
1267:   }
1268:   if (label) {
1269:     ISRestoreIndices(valueIS, &values);
1270:     ISDestroy(&valueIS);
1271:   }
1272:   for (d = 0; d <= depth; ++d) {
1273:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1274:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1275:   }
1276:   DMPlexSetHybridBounds(sdm, depth >= 0 ? pMaxNew[depth] : PETSC_DETERMINE, depth>1 ? pMaxNew[depth-1] : PETSC_DETERMINE, depth>2 ? pMaxNew[1] : PETSC_DETERMINE, depth >= 0 ? pMaxNew[0] : PETSC_DETERMINE);
1277:   PetscFree3(coneNew, coneONew, supportNew);
1278:   PetscFree6(depthMax, depthEnd, depthShift, depthOffset, pMaxNew, numHybridPointsOld);
1279:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1280:   return(0);
1281: }

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

1288:   Collective on dm

1290:   Input Parameters:
1291: + dm - The original DM
1292: - label - The label specifying the boundary faces (this could be auto-generated)

1294:   Output Parameters:
1295: - dmSplit - The new DM

1297:   Level: developer

1299: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1300: @*/
1301: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1302: {
1303:   DM             sdm;
1304:   PetscInt       dim;

1310:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1311:   DMSetType(sdm, DMPLEX);
1312:   DMGetDimension(dm, &dim);
1313:   DMSetDimension(sdm, dim);
1314:   switch (dim) {
1315:   case 2:
1316:   case 3:
1317:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1318:     break;
1319:   default:
1320:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1321:   }
1322:   *dmSplit = sdm;
1323:   return(0);
1324: }

1328: /* Returns the side of the surface for a given cell with a face on the surface */
1329: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1330: {
1331:   const PetscInt *cone, *ornt;
1332:   PetscInt        dim, coneSize, c;
1333:   PetscErrorCode  ierr;

1336:   *pos = PETSC_TRUE;
1337:   DMGetDimension(dm, &dim);
1338:   DMPlexGetConeSize(dm, cell, &coneSize);
1339:   DMPlexGetCone(dm, cell, &cone);
1340:   DMPlexGetConeOrientation(dm, cell, &ornt);
1341:   for (c = 0; c < coneSize; ++c) {
1342:     if (cone[c] == face) {
1343:       PetscInt o = ornt[c];

1345:       if (subdm) {
1346:         const PetscInt *subcone, *subornt;
1347:         PetscInt        subpoint, subface, subconeSize, sc;

1349:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1350:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1351:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1352:         DMPlexGetCone(subdm, subpoint, &subcone);
1353:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1354:         for (sc = 0; sc < subconeSize; ++sc) {
1355:           if (subcone[sc] == subface) {
1356:             o = subornt[0];
1357:             break;
1358:           }
1359:         }
1360:         if (sc >= subconeSize) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find subpoint %d (%d) in cone for subpoint %d (%d)", subface, face, subpoint, cell);
1361:       }
1362:       if (o >= 0) *pos = PETSC_TRUE;
1363:       else        *pos = PETSC_FALSE;
1364:       break;
1365:     }
1366:   }
1367:   if (c == coneSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell %d in split face %d support does not have it in the cone", cell, face);
1368:   return(0);
1369: }

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

1377:   Input Parameters:
1378: + dm     - The DM
1379: . label  - A DMLabel marking the surface
1380: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1381: . flip   - Flag to flip the submesh normal and replace points on the other side
1382: - subdm  - The subDM associated with the label, or NULL

1384:   Output Parameter:
1385: . label - A DMLabel marking all surface points

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

1389:   Level: developer

1391: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1392: @*/
1393: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1394: {
1395:   DMLabel         depthLabel;
1396:   IS              dimIS, subpointIS, facePosIS, faceNegIS, crossEdgeIS = NULL;
1397:   const PetscInt *points, *subpoints;
1398:   const PetscInt  rev   = flip ? -1 : 1;
1399:   PetscInt       *pMax;
1400:   PetscInt        shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1401:   PetscErrorCode  ierr;

1404:   DMPlexGetDepth(dm, &depth);
1405:   DMGetDimension(dm, &dim);
1406:   pSize = PetscMax(depth, dim) + 1;
1407:   PetscMalloc1(pSize,&pMax);
1408:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1409:   DMPlexGetDepthLabel(dm, &depthLabel);
1410:   DMGetDimension(dm, &dim);
1411:   if (subdm) {
1412:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1413:     if (subpointIS) {
1414:       ISGetLocalSize(subpointIS, &numSubpoints);
1415:       ISGetIndices(subpointIS, &subpoints);
1416:     }
1417:   }
1418:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1419:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1420:   if (!dimIS) {
1421:     PetscFree(pMax);
1422:     return(0);
1423:   }
1424:   ISGetLocalSize(dimIS, &numPoints);
1425:   ISGetIndices(dimIS, &points);
1426:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1427:     const PetscInt *support;
1428:     PetscInt        supportSize, s;

1430:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1431:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1432:     DMPlexGetSupport(dm, points[p], &support);
1433:     for (s = 0; s < supportSize; ++s) {
1434:       const PetscInt *cone;
1435:       PetscInt        coneSize, c;
1436:       PetscBool       pos;

1438:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1439:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1440:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1441:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1442:       /* Put faces touching the fault in the label */
1443:       DMPlexGetConeSize(dm, support[s], &coneSize);
1444:       DMPlexGetCone(dm, support[s], &cone);
1445:       for (c = 0; c < coneSize; ++c) {
1446:         const PetscInt point = cone[c];

1448:         DMLabelGetValue(label, point, &val);
1449:         if (val == -1) {
1450:           PetscInt *closure = NULL;
1451:           PetscInt  closureSize, cl;

1453:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1454:           for (cl = 0; cl < closureSize*2; cl += 2) {
1455:             const PetscInt clp  = closure[cl];
1456:             PetscInt       bval = -1;

1458:             DMLabelGetValue(label, clp, &val);
1459:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1460:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1461:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1462:               break;
1463:             }
1464:           }
1465:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1466:         }
1467:       }
1468:     }
1469:   }
1470:   ISRestoreIndices(dimIS, &points);
1471:   ISDestroy(&dimIS);
1472:   if (subdm) {
1473:     if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1474:     ISDestroy(&subpointIS);
1475:   }
1476:   /* Mark boundary points as unsplit */
1477:   if (blabel) {
1478:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1479:     ISGetLocalSize(dimIS, &numPoints);
1480:     ISGetIndices(dimIS, &points);
1481:     for (p = 0; p < numPoints; ++p) {
1482:       const PetscInt point = points[p];
1483:       PetscInt       val, bval;

1485:       DMLabelGetValue(blabel, point, &bval);
1486:       if (bval >= 0) {
1487:         DMLabelGetValue(label, point, &val);
1488:         if ((val < 0) || (val > dim)) {
1489:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1490:           DMLabelClearValue(blabel, point, bval);
1491:         }
1492:       }
1493:     }
1494:     for (p = 0; p < numPoints; ++p) {
1495:       const PetscInt point = points[p];
1496:       PetscInt       val, bval;

1498:       DMLabelGetValue(blabel, point, &bval);
1499:       if (bval >= 0) {
1500:         const PetscInt *cone,    *support;
1501:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1503:         /* Mark as unsplit */
1504:         DMLabelGetValue(label, point, &val);
1505:         if ((val < 0) || (val > dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d has label value %d, should be part of the fault", point, val);
1506:         DMLabelClearValue(label, point, val);
1507:         DMLabelSetValue(label, point, shift2+val);
1508:         /* Check for cross-edge
1509:              A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1510:         if (val != 0) continue;
1511:         DMPlexGetSupport(dm, point, &support);
1512:         DMPlexGetSupportSize(dm, point, &supportSize);
1513:         for (s = 0; s < supportSize; ++s) {
1514:           DMPlexGetCone(dm, support[s], &cone);
1515:           DMPlexGetConeSize(dm, support[s], &coneSize);
1516:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1517:           DMLabelGetValue(blabel, cone[0], &valA);
1518:           DMLabelGetValue(blabel, cone[1], &valB);
1519:           DMLabelGetValue(blabel, support[s], &valE);
1520:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1521:         }
1522:       }
1523:     }
1524:     ISRestoreIndices(dimIS, &points);
1525:     ISDestroy(&dimIS);
1526:   }
1527:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1528:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1529:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1530:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1531:   cMax = cMax < 0 ? cEnd : cMax;
1532:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1533:   DMLabelGetStratumIS(label, 0, &dimIS);
1534:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1535:   if (dimIS && crossEdgeIS) {
1536:     IS vertIS = dimIS;

1538:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1539:     ISDestroy(&crossEdgeIS);
1540:     ISDestroy(&vertIS);
1541:   }
1542:   if (!dimIS) {
1543:     PetscFree(pMax);
1544:     return(0);
1545:   }
1546:   ISGetLocalSize(dimIS, &numPoints);
1547:   ISGetIndices(dimIS, &points);
1548:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1549:     PetscInt *star = NULL;
1550:     PetscInt  starSize, s;
1551:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1553:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1554:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1555:     while (again) {
1556:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1557:       again = 0;
1558:       for (s = 0; s < starSize*2; s += 2) {
1559:         const PetscInt  point = star[s];
1560:         const PetscInt *cone;
1561:         PetscInt        coneSize, c;

1563:         if ((point < cStart) || (point >= cMax)) continue;
1564:         DMLabelGetValue(label, point, &val);
1565:         if (val != -1) continue;
1566:         again = again == 1 ? 1 : 2;
1567:         DMPlexGetConeSize(dm, point, &coneSize);
1568:         DMPlexGetCone(dm, point, &cone);
1569:         for (c = 0; c < coneSize; ++c) {
1570:           DMLabelGetValue(label, cone[c], &val);
1571:           if (val != -1) {
1572:             const PetscInt *ccone;
1573:             PetscInt        cconeSize, cc, side;

1575:             if (abs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
1576:             if (val > 0) side =  1;
1577:             else         side = -1;
1578:             DMLabelSetValue(label, point, side*(shift+dim));
1579:             /* Mark cell faces which touch the fault */
1580:             DMPlexGetConeSize(dm, point, &cconeSize);
1581:             DMPlexGetCone(dm, point, &ccone);
1582:             for (cc = 0; cc < cconeSize; ++cc) {
1583:               PetscInt *closure = NULL;
1584:               PetscInt  closureSize, cl;

1586:               DMLabelGetValue(label, ccone[cc], &val);
1587:               if (val != -1) continue;
1588:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1589:               for (cl = 0; cl < closureSize*2; cl += 2) {
1590:                 const PetscInt clp = closure[cl];

1592:                 DMLabelGetValue(label, clp, &val);
1593:                 if (val == -1) continue;
1594:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1595:                 break;
1596:               }
1597:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1598:             }
1599:             again = 1;
1600:             break;
1601:           }
1602:         }
1603:       }
1604:     }
1605:     /* Classify the rest by cell membership */
1606:     for (s = 0; s < starSize*2; s += 2) {
1607:       const PetscInt point = star[s];

1609:       DMLabelGetValue(label, point, &val);
1610:       if (val == -1) {
1611:         PetscInt  *sstar = NULL;
1612:         PetscInt   sstarSize, ss;
1613:         PetscBool  marked = PETSC_FALSE;

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

1619:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1620:           DMLabelGetValue(label, spoint, &val);
1621:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1622:           DMLabelGetValue(depthLabel, point, &dep);
1623:           if (val > 0) {
1624:             DMLabelSetValue(label, point,   shift+dep);
1625:           } else {
1626:             DMLabelSetValue(label, point, -(shift+dep));
1627:           }
1628:           marked = PETSC_TRUE;
1629:           break;
1630:         }
1631:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1632:         DMLabelGetValue(depthLabel, point, &dep);
1633:         if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1634:       }
1635:     }
1636:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1637:   }
1638:   ISRestoreIndices(dimIS, &points);
1639:   ISDestroy(&dimIS);
1640:   /* If any faces touching the fault divide cells on either side, split them */
1641:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1642:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1643:   ISExpand(facePosIS, faceNegIS, &dimIS);
1644:   ISDestroy(&facePosIS);
1645:   ISDestroy(&faceNegIS);
1646:   ISGetLocalSize(dimIS, &numPoints);
1647:   ISGetIndices(dimIS, &points);
1648:   for (p = 0; p < numPoints; ++p) {
1649:     const PetscInt  point = points[p];
1650:     const PetscInt *support;
1651:     PetscInt        supportSize, valA, valB;

1653:     DMPlexGetSupportSize(dm, point, &supportSize);
1654:     if (supportSize != 2) continue;
1655:     DMPlexGetSupport(dm, point, &support);
1656:     DMLabelGetValue(label, support[0], &valA);
1657:     DMLabelGetValue(label, support[1], &valB);
1658:     if ((valA == -1) || (valB == -1)) continue;
1659:     if (valA*valB > 0) continue;
1660:     /* Split the face */
1661:     DMLabelGetValue(label, point, &valA);
1662:     DMLabelClearValue(label, point, valA);
1663:     DMLabelSetValue(label, point, dim-1);
1664:     /* Label its closure:
1665:       unmarked: label as unsplit
1666:       incident: relabel as split
1667:       split:    do nothing
1668:     */
1669:     {
1670:       PetscInt *closure = NULL;
1671:       PetscInt  closureSize, cl;

1673:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1674:       for (cl = 0; cl < closureSize*2; cl += 2) {
1675:         DMLabelGetValue(label, closure[cl], &valA);
1676:         if (valA == -1) { /* Mark as unsplit */
1677:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1678:           DMLabelSetValue(label, closure[cl], shift2+dep);
1679:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1680:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1681:           DMLabelClearValue(label, closure[cl], valA);
1682:           DMLabelSetValue(label, closure[cl], dep);
1683:         }
1684:       }
1685:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1686:     }
1687:   }
1688:   ISRestoreIndices(dimIS, &points);
1689:   ISDestroy(&dimIS);
1690:   PetscFree(pMax);
1691:   return(0);
1692: }

1696: /*@C
1697:   DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface

1699:   Collective on dm

1701:   Input Parameters:
1702: + dm - The original DM
1703: - labelName - The label specifying the interface vertices

1705:   Output Parameters:
1706: + hybridLabel - The label fully marking the interface
1707: - dmHybrid - The new DM

1709:   Level: developer

1711: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1712: @*/
1713: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1714: {
1715:   DM             idm;
1716:   DMLabel        subpointMap, hlabel;
1717:   PetscInt       dim;

1724:   DMGetDimension(dm, &dim);
1725:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1726:   DMPlexOrient(idm);
1727:   DMPlexGetSubpointMap(idm, &subpointMap);
1728:   DMLabelDuplicate(subpointMap, &hlabel);
1729:   DMLabelClearStratum(hlabel, dim);
1730:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1731:   DMDestroy(&idm);
1732:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1733:   if (hybridLabel) *hybridLabel = hlabel;
1734:   else             {DMLabelDestroy(&hlabel);}
1735:   return(0);
1736: }

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

1742:      For any marked cell, the marked vertices constitute a single face
1743: */
1744: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1745: {
1746:   IS               subvertexIS = NULL;
1747:   const PetscInt  *subvertices;
1748:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1749:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1750:   PetscErrorCode   ierr;

1753:   *numFaces = 0;
1754:   *nFV      = 0;
1755:   DMPlexGetDepth(dm, &depth);
1756:   DMGetDimension(dm, &dim);
1757:   pSize = PetscMax(depth, dim) + 1;
1758:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1759:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1760:   for (d = 0; d <= depth; ++d) {
1761:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1762:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1763:   }
1764:   /* Loop over initial vertices and mark all faces in the collective star() */
1765:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1766:   if (subvertexIS) {
1767:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1768:     ISGetIndices(subvertexIS, &subvertices);
1769:   }
1770:   for (v = 0; v < numSubVerticesInitial; ++v) {
1771:     const PetscInt vertex = subvertices[v];
1772:     PetscInt      *star   = NULL;
1773:     PetscInt       starSize, s, numCells = 0, c;

1775:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1776:     for (s = 0; s < starSize*2; s += 2) {
1777:       const PetscInt point = star[s];
1778:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1779:     }
1780:     for (c = 0; c < numCells; ++c) {
1781:       const PetscInt cell    = star[c];
1782:       PetscInt      *closure = NULL;
1783:       PetscInt       closureSize, cl;
1784:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1786:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1787:       if (cellLoc == 2) continue;
1788:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1789:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1790:       for (cl = 0; cl < closureSize*2; cl += 2) {
1791:         const PetscInt point = closure[cl];
1792:         PetscInt       vertexLoc;

1794:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1795:           ++numCorners;
1796:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1797:           if (vertexLoc == value) closure[faceSize++] = point;
1798:         }
1799:       }
1800:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1801:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1802:       if (faceSize == *nFV) {
1803:         const PetscInt *cells = NULL;
1804:         PetscInt        numCells, nc;

1806:         ++(*numFaces);
1807:         for (cl = 0; cl < faceSize; ++cl) {
1808:           DMLabelSetValue(subpointMap, closure[cl], 0);
1809:         }
1810:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1811:         for (nc = 0; nc < numCells; ++nc) {
1812:           DMLabelSetValue(subpointMap, cells[nc], 2);
1813:         }
1814:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1815:       }
1816:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1817:     }
1818:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1819:   }
1820:   if (subvertexIS) {
1821:     ISRestoreIndices(subvertexIS, &subvertices);
1822:   }
1823:   ISDestroy(&subvertexIS);
1824:   PetscFree3(pStart,pEnd,pMax);
1825:   return(0);
1826: }

1830: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1831: {
1832:   IS               subvertexIS;
1833:   const PetscInt  *subvertices;
1834:   PetscInt        *pStart, *pEnd, *pMax;
1835:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
1836:   PetscErrorCode   ierr;

1839:   DMGetDimension(dm, &dim);
1840:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1841:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1842:   for (d = 0; d <= dim; ++d) {
1843:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1844:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1845:   }
1846:   /* Loop over initial vertices and mark all faces in the collective star() */
1847:   DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1848:   if (subvertexIS) {
1849:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1850:     ISGetIndices(subvertexIS, &subvertices);
1851:   }
1852:   for (v = 0; v < numSubVerticesInitial; ++v) {
1853:     const PetscInt vertex = subvertices[v];
1854:     PetscInt      *star   = NULL;
1855:     PetscInt       starSize, s, numFaces = 0, f;

1857:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1858:     for (s = 0; s < starSize*2; s += 2) {
1859:       const PetscInt point = star[s];
1860:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1861:     }
1862:     for (f = 0; f < numFaces; ++f) {
1863:       const PetscInt face    = star[f];
1864:       PetscInt      *closure = NULL;
1865:       PetscInt       closureSize, c;
1866:       PetscInt       faceLoc;

1868:       DMLabelGetValue(subpointMap, face, &faceLoc);
1869:       if (faceLoc == dim-1) continue;
1870:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1871:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1872:       for (c = 0; c < closureSize*2; c += 2) {
1873:         const PetscInt point = closure[c];
1874:         PetscInt       vertexLoc;

1876:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1877:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1878:           if (vertexLoc != value) break;
1879:         }
1880:       }
1881:       if (c == closureSize*2) {
1882:         const PetscInt *support;
1883:         PetscInt        supportSize, s;

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

1888:           for (d = 0; d < dim; ++d) {
1889:             if ((point >= pStart[d]) && (point < pEnd[d])) {
1890:               DMLabelSetValue(subpointMap, point, d);
1891:               break;
1892:             }
1893:           }
1894:         }
1895:         DMPlexGetSupportSize(dm, face, &supportSize);
1896:         DMPlexGetSupport(dm, face, &support);
1897:         for (s = 0; s < supportSize; ++s) {
1898:           DMLabelSetValue(subpointMap, support[s], dim);
1899:         }
1900:       }
1901:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1902:     }
1903:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1904:   }
1905:   if (subvertexIS) {
1906:     ISRestoreIndices(subvertexIS, &subvertices);
1907:   }
1908:   ISDestroy(&subvertexIS);
1909:   PetscFree3(pStart,pEnd,pMax);
1910:   return(0);
1911: }

1915: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
1916: {
1917:   DMLabel         label = NULL;
1918:   const PetscInt *cone;
1919:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize;
1920:   PetscErrorCode  ierr;

1923:   *numFaces = 0;
1924:   *nFV = 0;
1925:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1926:   *subCells = NULL;
1927:   DMGetDimension(dm, &dim);
1928:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1929:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1930:   if (cMax < 0) return(0);
1931:   if (label) {
1932:     for (c = cMax; c < cEnd; ++c) {
1933:       PetscInt val;

1935:       DMLabelGetValue(label, c, &val);
1936:       if (val == value) {
1937:         ++(*numFaces);
1938:         DMPlexGetConeSize(dm, c, &coneSize);
1939:       }
1940:     }
1941:   } else {
1942:     *numFaces = cEnd - cMax;
1943:     DMPlexGetConeSize(dm, cMax, &coneSize);
1944:   }
1945:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
1946:   PetscMalloc1(*numFaces *2, subCells);
1947:   for (c = cMax; c < cEnd; ++c) {
1948:     const PetscInt *cells;
1949:     PetscInt        numCells;

1951:     if (label) {
1952:       PetscInt val;

1954:       DMLabelGetValue(label, c, &val);
1955:       if (val != value) continue;
1956:     }
1957:     DMPlexGetCone(dm, c, &cone);
1958:     for (p = 0; p < *nFV; ++p) {
1959:       DMLabelSetValue(subpointMap, cone[p], 0);
1960:     }
1961:     /* Negative face */
1962:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
1963:     /* Not true in parallel
1964:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
1965:     for (p = 0; p < numCells; ++p) {
1966:       DMLabelSetValue(subpointMap, cells[p], 2);
1967:       (*subCells)[subc++] = cells[p];
1968:     }
1969:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
1970:     /* Positive face is not included */
1971:   }
1972:   return(0);
1973: }

1977: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DMLabel subpointMap, DM subdm)
1978: {
1979:   DMLabel        label = NULL;
1980:   PetscInt      *pStart, *pEnd;
1981:   PetscInt       dim, cMax, cEnd, c, d;

1985:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1986:   DMGetDimension(dm, &dim);
1987:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1988:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1989:   if (cMax < 0) return(0);
1990:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
1991:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1992:   for (c = cMax; c < cEnd; ++c) {
1993:     const PetscInt *cone;
1994:     PetscInt       *closure = NULL;
1995:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

1997:     if (label) {
1998:       DMLabelGetValue(label, c, &val);
1999:       if (val != value) continue;
2000:     }
2001:     DMPlexGetConeSize(dm, c, &coneSize);
2002:     DMPlexGetCone(dm, c, &cone);
2003:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
2004:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2005:     /* Negative face */
2006:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2007:     for (cl = 0; cl < closureSize*2; cl += 2) {
2008:       const PetscInt point = closure[cl];

2010:       for (d = 0; d <= dim; ++d) {
2011:         if ((point >= pStart[d]) && (point < pEnd[d])) {
2012:           DMLabelSetValue(subpointMap, point, d);
2013:           break;
2014:         }
2015:       }
2016:     }
2017:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2018:     /* Cells -- positive face is not included */
2019:     for (cl = 0; cl < 1; ++cl) {
2020:       const PetscInt *support;
2021:       PetscInt        supportSize, s;

2023:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2024:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2025:       DMPlexGetSupport(dm, cone[cl], &support);
2026:       for (s = 0; s < supportSize; ++s) {
2027:         DMLabelSetValue(subpointMap, support[s], dim);
2028:       }
2029:     }
2030:   }
2031:   PetscFree2(pStart, pEnd);
2032:   return(0);
2033: }

2037: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2038: {
2039:   MPI_Comm       comm;
2040:   PetscBool      posOrient = PETSC_FALSE;
2041:   const PetscInt debug     = 0;
2042:   PetscInt       cellDim, faceSize, f;

2046:   PetscObjectGetComm((PetscObject)dm,&comm);
2047:   DMGetDimension(dm, &cellDim);
2048:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2050:   if (cellDim == 1 && numCorners == 2) {
2051:     /* Triangle */
2052:     faceSize  = numCorners-1;
2053:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2054:   } else if (cellDim == 2 && numCorners == 3) {
2055:     /* Triangle */
2056:     faceSize  = numCorners-1;
2057:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2058:   } else if (cellDim == 3 && numCorners == 4) {
2059:     /* Tetrahedron */
2060:     faceSize  = numCorners-1;
2061:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2062:   } else if (cellDim == 1 && numCorners == 3) {
2063:     /* Quadratic line */
2064:     faceSize  = 1;
2065:     posOrient = PETSC_TRUE;
2066:   } else if (cellDim == 2 && numCorners == 4) {
2067:     /* Quads */
2068:     faceSize = 2;
2069:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2070:       posOrient = PETSC_TRUE;
2071:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2072:       posOrient = PETSC_TRUE;
2073:     } else {
2074:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2075:         posOrient = PETSC_FALSE;
2076:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2077:     }
2078:   } else if (cellDim == 2 && numCorners == 6) {
2079:     /* Quadratic triangle (I hate this) */
2080:     /* Edges are determined by the first 2 vertices (corners of edges) */
2081:     const PetscInt faceSizeTri = 3;
2082:     PetscInt       sortedIndices[3], i, iFace;
2083:     PetscBool      found                    = PETSC_FALSE;
2084:     PetscInt       faceVerticesTriSorted[9] = {
2085:       0, 3,  4, /* bottom */
2086:       1, 4,  5, /* right */
2087:       2, 3,  5, /* left */
2088:     };
2089:     PetscInt       faceVerticesTri[9] = {
2090:       0, 3,  4, /* bottom */
2091:       1, 4,  5, /* right */
2092:       2, 5,  3, /* left */
2093:     };

2095:     faceSize = faceSizeTri;
2096:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2097:     PetscSortInt(faceSizeTri, sortedIndices);
2098:     for (iFace = 0; iFace < 3; ++iFace) {
2099:       const PetscInt ii = iFace*faceSizeTri;
2100:       PetscInt       fVertex, cVertex;

2102:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2103:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2104:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2105:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2106:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2107:               faceVertices[fVertex] = origVertices[cVertex];
2108:               break;
2109:             }
2110:           }
2111:         }
2112:         found = PETSC_TRUE;
2113:         break;
2114:       }
2115:     }
2116:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2117:     if (posOriented) *posOriented = PETSC_TRUE;
2118:     return(0);
2119:   } else if (cellDim == 2 && numCorners == 9) {
2120:     /* Quadratic quad (I hate this) */
2121:     /* Edges are determined by the first 2 vertices (corners of edges) */
2122:     const PetscInt faceSizeQuad = 3;
2123:     PetscInt       sortedIndices[3], i, iFace;
2124:     PetscBool      found                      = PETSC_FALSE;
2125:     PetscInt       faceVerticesQuadSorted[12] = {
2126:       0, 1,  4, /* bottom */
2127:       1, 2,  5, /* right */
2128:       2, 3,  6, /* top */
2129:       0, 3,  7, /* left */
2130:     };
2131:     PetscInt       faceVerticesQuad[12] = {
2132:       0, 1,  4, /* bottom */
2133:       1, 2,  5, /* right */
2134:       2, 3,  6, /* top */
2135:       3, 0,  7, /* left */
2136:     };

2138:     faceSize = faceSizeQuad;
2139:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2140:     PetscSortInt(faceSizeQuad, sortedIndices);
2141:     for (iFace = 0; iFace < 4; ++iFace) {
2142:       const PetscInt ii = iFace*faceSizeQuad;
2143:       PetscInt       fVertex, cVertex;

2145:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2146:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2147:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2148:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2149:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2150:               faceVertices[fVertex] = origVertices[cVertex];
2151:               break;
2152:             }
2153:           }
2154:         }
2155:         found = PETSC_TRUE;
2156:         break;
2157:       }
2158:     }
2159:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2160:     if (posOriented) *posOriented = PETSC_TRUE;
2161:     return(0);
2162:   } else if (cellDim == 3 && numCorners == 8) {
2163:     /* Hexes
2164:        A hex is two oriented quads with the normal of the first
2165:        pointing up at the second.

2167:           7---6
2168:          /|  /|
2169:         4---5 |
2170:         | 1-|-2
2171:         |/  |/
2172:         0---3

2174:         Faces are determined by the first 4 vertices (corners of faces) */
2175:     const PetscInt faceSizeHex = 4;
2176:     PetscInt       sortedIndices[4], i, iFace;
2177:     PetscBool      found                     = PETSC_FALSE;
2178:     PetscInt       faceVerticesHexSorted[24] = {
2179:       0, 1, 2, 3,  /* bottom */
2180:       4, 5, 6, 7,  /* top */
2181:       0, 3, 4, 5,  /* front */
2182:       2, 3, 5, 6,  /* right */
2183:       1, 2, 6, 7,  /* back */
2184:       0, 1, 4, 7,  /* left */
2185:     };
2186:     PetscInt       faceVerticesHex[24] = {
2187:       1, 2, 3, 0,  /* bottom */
2188:       4, 5, 6, 7,  /* top */
2189:       0, 3, 5, 4,  /* front */
2190:       3, 2, 6, 5,  /* right */
2191:       2, 1, 7, 6,  /* back */
2192:       1, 0, 4, 7,  /* left */
2193:     };

2195:     faceSize = faceSizeHex;
2196:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2197:     PetscSortInt(faceSizeHex, sortedIndices);
2198:     for (iFace = 0; iFace < 6; ++iFace) {
2199:       const PetscInt ii = iFace*faceSizeHex;
2200:       PetscInt       fVertex, cVertex;

2202:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2203:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2204:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2205:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2206:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2207:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2208:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2209:               faceVertices[fVertex] = origVertices[cVertex];
2210:               break;
2211:             }
2212:           }
2213:         }
2214:         found = PETSC_TRUE;
2215:         break;
2216:       }
2217:     }
2218:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2219:     if (posOriented) *posOriented = PETSC_TRUE;
2220:     return(0);
2221:   } else if (cellDim == 3 && numCorners == 10) {
2222:     /* Quadratic tet */
2223:     /* Faces are determined by the first 3 vertices (corners of faces) */
2224:     const PetscInt faceSizeTet = 6;
2225:     PetscInt       sortedIndices[6], i, iFace;
2226:     PetscBool      found                     = PETSC_FALSE;
2227:     PetscInt       faceVerticesTetSorted[24] = {
2228:       0, 1, 2,  6, 7, 8, /* bottom */
2229:       0, 3, 4,  6, 7, 9,  /* front */
2230:       1, 4, 5,  7, 8, 9,  /* right */
2231:       2, 3, 5,  6, 8, 9,  /* left */
2232:     };
2233:     PetscInt       faceVerticesTet[24] = {
2234:       0, 1, 2,  6, 7, 8, /* bottom */
2235:       0, 4, 3,  6, 7, 9,  /* front */
2236:       1, 5, 4,  7, 8, 9,  /* right */
2237:       2, 3, 5,  8, 6, 9,  /* left */
2238:     };

2240:     faceSize = faceSizeTet;
2241:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2242:     PetscSortInt(faceSizeTet, sortedIndices);
2243:     for (iFace=0; iFace < 4; ++iFace) {
2244:       const PetscInt ii = iFace*faceSizeTet;
2245:       PetscInt       fVertex, cVertex;

2247:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2248:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2249:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2250:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2251:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2252:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2253:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2254:               faceVertices[fVertex] = origVertices[cVertex];
2255:               break;
2256:             }
2257:           }
2258:         }
2259:         found = PETSC_TRUE;
2260:         break;
2261:       }
2262:     }
2263:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2264:     if (posOriented) *posOriented = PETSC_TRUE;
2265:     return(0);
2266:   } else if (cellDim == 3 && numCorners == 27) {
2267:     /* Quadratic hexes (I hate this)
2268:        A hex is two oriented quads with the normal of the first
2269:        pointing up at the second.

2271:          7---6
2272:         /|  /|
2273:        4---5 |
2274:        | 3-|-2
2275:        |/  |/
2276:        0---1

2278:        Faces are determined by the first 4 vertices (corners of faces) */
2279:     const PetscInt faceSizeQuadHex = 9;
2280:     PetscInt       sortedIndices[9], i, iFace;
2281:     PetscBool      found                         = PETSC_FALSE;
2282:     PetscInt       faceVerticesQuadHexSorted[54] = {
2283:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2284:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2285:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2286:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2287:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2288:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2289:     };
2290:     PetscInt       faceVerticesQuadHex[54] = {
2291:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2292:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2293:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2294:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2295:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2296:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2297:     };

2299:     faceSize = faceSizeQuadHex;
2300:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2301:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2302:     for (iFace = 0; iFace < 6; ++iFace) {
2303:       const PetscInt ii = iFace*faceSizeQuadHex;
2304:       PetscInt       fVertex, cVertex;

2306:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2307:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2308:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2309:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2310:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2311:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2312:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2313:               faceVertices[fVertex] = origVertices[cVertex];
2314:               break;
2315:             }
2316:           }
2317:         }
2318:         found = PETSC_TRUE;
2319:         break;
2320:       }
2321:     }
2322:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2323:     if (posOriented) *posOriented = PETSC_TRUE;
2324:     return(0);
2325:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2326:   if (!posOrient) {
2327:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2328:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2329:   } else {
2330:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2331:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2332:   }
2333:   if (posOriented) *posOriented = posOrient;
2334:   return(0);
2335: }

2339: /*
2340:     Given a cell and a face, as a set of vertices,
2341:       return the oriented face, as a set of vertices, in faceVertices
2342:     The orientation is such that the face normal points out of the cell
2343: */
2344: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2345: {
2346:   const PetscInt *cone = NULL;
2347:   PetscInt        coneSize, v, f, v2;
2348:   PetscInt        oppositeVertex = -1;
2349:   PetscErrorCode  ierr;

2352:   DMPlexGetConeSize(dm, cell, &coneSize);
2353:   DMPlexGetCone(dm, cell, &cone);
2354:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2355:     PetscBool found = PETSC_FALSE;

2357:     for (f = 0; f < faceSize; ++f) {
2358:       if (face[f] == cone[v]) {
2359:         found = PETSC_TRUE; break;
2360:       }
2361:     }
2362:     if (found) {
2363:       indices[v2]      = v;
2364:       origVertices[v2] = cone[v];
2365:       ++v2;
2366:     } else {
2367:       oppositeVertex = v;
2368:     }
2369:   }
2370:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2371:   return(0);
2372: }

2376: /*
2377:   DMPlexInsertFace_Internal - Puts a face into the mesh

2379:   Not collective

2381:   Input Parameters:
2382:   + dm              - The DMPlex
2383:   . numFaceVertex   - The number of vertices in the face
2384:   . faceVertices    - The vertices in the face for dm
2385:   . subfaceVertices - The vertices in the face for subdm
2386:   . numCorners      - The number of vertices in the cell
2387:   . cell            - A cell in dm containing the face
2388:   . subcell         - A cell in subdm containing the face
2389:   . firstFace       - First face in the mesh
2390:   - newFacePoint    - Next face in the mesh

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

2395:   Level: developer
2396: */
2397: 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)
2398: {
2399:   MPI_Comm        comm;
2400:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2401:   const PetscInt *faces;
2402:   PetscInt        numFaces, coneSize;
2403:   PetscErrorCode  ierr;

2406:   PetscObjectGetComm((PetscObject)dm,&comm);
2407:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2408:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2409: #if 0
2410:   /* Cannot use this because support() has not been constructed yet */
2411:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2412: #else
2413:   {
2414:     PetscInt f;

2416:     numFaces = 0;
2417:     DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2418:     for (f = firstFace; f < *newFacePoint; ++f) {
2419:       PetscInt dof, off, d;

2421:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2422:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2423:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2424:       for (d = 0; d < dof; ++d) {
2425:         const PetscInt p = submesh->cones[off+d];
2426:         PetscInt       v;

2428:         for (v = 0; v < numFaceVertices; ++v) {
2429:           if (subfaceVertices[v] == p) break;
2430:         }
2431:         if (v == numFaceVertices) break;
2432:       }
2433:       if (d == dof) {
2434:         numFaces               = 1;
2435:         ((PetscInt*) faces)[0] = f;
2436:       }
2437:     }
2438:   }
2439: #endif
2440:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2441:   else if (numFaces == 1) {
2442:     /* Add the other cell neighbor for this face */
2443:     DMPlexSetCone(subdm, subcell, faces);
2444:   } else {
2445:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2446:     PetscBool posOriented;

2448:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2449:     origVertices        = &orientedVertices[numFaceVertices];
2450:     indices             = &orientedVertices[numFaceVertices*2];
2451:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2452:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2453:     /* TODO: I know that routine should return a permutation, not the indices */
2454:     for (v = 0; v < numFaceVertices; ++v) {
2455:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2456:       for (ov = 0; ov < numFaceVertices; ++ov) {
2457:         if (orientedVertices[ov] == vertex) {
2458:           orientedSubVertices[ov] = subvertex;
2459:           break;
2460:         }
2461:       }
2462:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2463:     }
2464:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2465:     DMPlexSetCone(subdm, subcell, newFacePoint);
2466:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2467:     ++(*newFacePoint);
2468:   }
2469: #if 0
2470:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2471: #else
2472:   DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2473: #endif
2474:   return(0);
2475: }

2479: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2480: {
2481:   MPI_Comm        comm;
2482:   DMLabel         subpointMap;
2483:   IS              subvertexIS,  subcellIS;
2484:   const PetscInt *subVertices, *subCells;
2485:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2486:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2487:   PetscInt        vStart, vEnd, c, f;
2488:   PetscErrorCode  ierr;

2491:   PetscObjectGetComm((PetscObject)dm,&comm);
2492:   /* Create subpointMap which marks the submesh */
2493:   DMLabelCreate("subpoint_map", &subpointMap);
2494:   DMPlexSetSubpointMap(subdm, subpointMap);
2495:   DMLabelDestroy(&subpointMap);
2496:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2497:   /* Setup chart */
2498:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2499:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2500:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2501:   DMPlexSetVTKCellHeight(subdm, 1);
2502:   /* Set cone sizes */
2503:   firstSubVertex = numSubCells;
2504:   firstSubFace   = numSubCells+numSubVertices;
2505:   newFacePoint   = firstSubFace;
2506:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2507:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2508:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2509:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2510:   for (c = 0; c < numSubCells; ++c) {
2511:     DMPlexSetConeSize(subdm, c, 1);
2512:   }
2513:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2514:     DMPlexSetConeSize(subdm, f, nFV);
2515:   }
2516:   DMSetUp(subdm);
2517:   /* Create face cones */
2518:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2519:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2520:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2521:   for (c = 0; c < numSubCells; ++c) {
2522:     const PetscInt cell    = subCells[c];
2523:     const PetscInt subcell = c;
2524:     PetscInt      *closure = NULL;
2525:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2527:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2528:     for (cl = 0; cl < closureSize*2; cl += 2) {
2529:       const PetscInt point = closure[cl];
2530:       PetscInt       subVertex;

2532:       if ((point >= vStart) && (point < vEnd)) {
2533:         ++numCorners;
2534:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2535:         if (subVertex >= 0) {
2536:           closure[faceSize] = point;
2537:           subface[faceSize] = firstSubVertex+subVertex;
2538:           ++faceSize;
2539:         }
2540:       }
2541:     }
2542:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2543:     if (faceSize == nFV) {
2544:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2545:     }
2546:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2547:   }
2548:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2549:   DMPlexSymmetrize(subdm);
2550:   DMPlexStratify(subdm);
2551:   /* Build coordinates */
2552:   {
2553:     PetscSection coordSection, subCoordSection;
2554:     Vec          coordinates, subCoordinates;
2555:     PetscScalar *coords, *subCoords;
2556:     PetscInt     numComp, coordSize, v;

2558:     DMGetCoordinateSection(dm, &coordSection);
2559:     DMGetCoordinatesLocal(dm, &coordinates);
2560:     DMGetCoordinateSection(subdm, &subCoordSection);
2561:     PetscSectionSetNumFields(subCoordSection, 1);
2562:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2563:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2564:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2565:     for (v = 0; v < numSubVertices; ++v) {
2566:       const PetscInt vertex    = subVertices[v];
2567:       const PetscInt subvertex = firstSubVertex+v;
2568:       PetscInt       dof;

2570:       PetscSectionGetDof(coordSection, vertex, &dof);
2571:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2572:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2573:     }
2574:     PetscSectionSetUp(subCoordSection);
2575:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2576:     VecCreate(comm, &subCoordinates);
2577:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2578:     VecSetType(subCoordinates,VECSTANDARD);
2579:     if (coordSize) {
2580:       VecGetArray(coordinates,    &coords);
2581:       VecGetArray(subCoordinates, &subCoords);
2582:       for (v = 0; v < numSubVertices; ++v) {
2583:         const PetscInt vertex    = subVertices[v];
2584:         const PetscInt subvertex = firstSubVertex+v;
2585:         PetscInt       dof, off, sdof, soff, d;

2587:         PetscSectionGetDof(coordSection, vertex, &dof);
2588:         PetscSectionGetOffset(coordSection, vertex, &off);
2589:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2590:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2591:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2592:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2593:       }
2594:       VecRestoreArray(coordinates,    &coords);
2595:       VecRestoreArray(subCoordinates, &subCoords);
2596:     }
2597:     DMSetCoordinatesLocal(subdm, subCoordinates);
2598:     VecDestroy(&subCoordinates);
2599:   }
2600:   /* Cleanup */
2601:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2602:   ISDestroy(&subvertexIS);
2603:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2604:   ISDestroy(&subcellIS);
2605:   return(0);
2606: }

2610: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2611: {
2612:   MPI_Comm         comm;
2613:   DMLabel          subpointMap;
2614:   IS              *subpointIS;
2615:   const PetscInt **subpoints;
2616:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *coneONew;
2617:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2618:   PetscErrorCode   ierr;

2621:   PetscObjectGetComm((PetscObject)dm,&comm);
2622:   /* Create subpointMap which marks the submesh */
2623:   DMLabelCreate("subpoint_map", &subpointMap);
2624:   DMPlexSetSubpointMap(subdm, subpointMap);
2625:   DMLabelDestroy(&subpointMap);
2626:   if (vertexLabel) {DMPlexMarkSubmesh_Interpolated(dm, vertexLabel, value, subpointMap, subdm);}
2627:   /* Setup chart */
2628:   DMGetDimension(dm, &dim);
2629:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2630:   for (d = 0; d <= dim; ++d) {
2631:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2632:     totSubPoints += numSubPoints[d];
2633:   }
2634:   DMPlexSetChart(subdm, 0, totSubPoints);
2635:   DMPlexSetVTKCellHeight(subdm, 1);
2636:   /* Set cone sizes */
2637:   firstSubPoint[dim] = 0;
2638:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2639:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2640:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2641:   for (d = 0; d <= dim; ++d) {
2642:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2643:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2644:   }
2645:   for (d = 0; d <= dim; ++d) {
2646:     for (p = 0; p < numSubPoints[d]; ++p) {
2647:       const PetscInt  point    = subpoints[d][p];
2648:       const PetscInt  subpoint = firstSubPoint[d] + p;
2649:       const PetscInt *cone;
2650:       PetscInt        coneSize, coneSizeNew, c, val;

2652:       DMPlexGetConeSize(dm, point, &coneSize);
2653:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2654:       if (d == dim) {
2655:         DMPlexGetCone(dm, point, &cone);
2656:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2657:           DMLabelGetValue(subpointMap, cone[c], &val);
2658:           if (val >= 0) coneSizeNew++;
2659:         }
2660:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2661:       }
2662:     }
2663:   }
2664:   DMSetUp(subdm);
2665:   /* Set cones */
2666:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2667:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&coneONew);
2668:   for (d = 0; d <= dim; ++d) {
2669:     for (p = 0; p < numSubPoints[d]; ++p) {
2670:       const PetscInt  point    = subpoints[d][p];
2671:       const PetscInt  subpoint = firstSubPoint[d] + p;
2672:       const PetscInt *cone, *ornt;
2673:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

2675:       DMPlexGetConeSize(dm, point, &coneSize);
2676:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2677:       DMPlexGetCone(dm, point, &cone);
2678:       DMPlexGetConeOrientation(dm, point, &ornt);
2679:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2680:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2681:         if (subc >= 0) {
2682:           coneNew[coneSizeNew]  = firstSubPoint[d-1] + subc;
2683:           coneONew[coneSizeNew] = ornt[c];
2684:           ++coneSizeNew;
2685:         }
2686:       }
2687:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2688:       DMPlexSetCone(subdm, subpoint, coneNew);
2689:       DMPlexSetConeOrientation(subdm, subpoint, coneONew);
2690:     }
2691:   }
2692:   PetscFree2(coneNew,coneONew);
2693:   DMPlexSymmetrize(subdm);
2694:   DMPlexStratify(subdm);
2695:   /* Build coordinates */
2696:   {
2697:     PetscSection coordSection, subCoordSection;
2698:     Vec          coordinates, subCoordinates;
2699:     PetscScalar *coords, *subCoords;
2700:     PetscInt     numComp, coordSize;

2702:     DMGetCoordinateSection(dm, &coordSection);
2703:     DMGetCoordinatesLocal(dm, &coordinates);
2704:     DMGetCoordinateSection(subdm, &subCoordSection);
2705:     PetscSectionSetNumFields(subCoordSection, 1);
2706:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2707:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2708:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2709:     for (v = 0; v < numSubPoints[0]; ++v) {
2710:       const PetscInt vertex    = subpoints[0][v];
2711:       const PetscInt subvertex = firstSubPoint[0]+v;
2712:       PetscInt       dof;

2714:       PetscSectionGetDof(coordSection, vertex, &dof);
2715:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2716:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2717:     }
2718:     PetscSectionSetUp(subCoordSection);
2719:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2720:     VecCreate(comm, &subCoordinates);
2721:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2722:     VecSetType(subCoordinates,VECSTANDARD);
2723:     VecGetArray(coordinates,    &coords);
2724:     VecGetArray(subCoordinates, &subCoords);
2725:     for (v = 0; v < numSubPoints[0]; ++v) {
2726:       const PetscInt vertex    = subpoints[0][v];
2727:       const PetscInt subvertex = firstSubPoint[0]+v;
2728:       PetscInt dof, off, sdof, soff, d;

2730:       PetscSectionGetDof(coordSection, vertex, &dof);
2731:       PetscSectionGetOffset(coordSection, vertex, &off);
2732:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2733:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2734:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2735:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2736:     }
2737:     VecRestoreArray(coordinates,    &coords);
2738:     VecRestoreArray(subCoordinates, &subCoords);
2739:     DMSetCoordinatesLocal(subdm, subCoordinates);
2740:     VecDestroy(&subCoordinates);
2741:   }
2742:   /* Cleanup */
2743:   for (d = 0; d <= dim; ++d) {
2744:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
2745:     ISDestroy(&subpointIS[d]);
2746:   }
2747:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2748:   return(0);
2749: }

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

2756:   Input Parameters:
2757: + dm           - The original mesh
2758: . vertexLabel  - The DMLabel marking vertices contained in the surface
2759: - value        - The label value to use

2761:   Output Parameter:
2762: . subdm - The surface mesh

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

2766:   Level: developer

2768: .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue()
2769: @*/
2770: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
2771: {
2772:   PetscInt       dim, depth;

2778:   DMGetDimension(dm, &dim);
2779:   DMPlexGetDepth(dm, &depth);
2780:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2781:   DMSetType(*subdm, DMPLEX);
2782:   DMSetDimension(*subdm, dim-1);
2783:   if (depth == dim) {
2784:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
2785:   } else {
2786:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
2787:   }
2788:   return(0);
2789: }

2793: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2794: {
2795:   PetscInt       subPoint;

2798:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2799:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2800: }

2804: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
2805: {
2806:   MPI_Comm        comm;
2807:   DMLabel         subpointMap;
2808:   IS              subvertexIS;
2809:   const PetscInt *subVertices;
2810:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
2811:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2812:   PetscInt        cMax, c, f;
2813:   PetscErrorCode  ierr;

2816:   PetscObjectGetComm((PetscObject)dm, &comm);
2817:   /* Create subpointMap which marks the submesh */
2818:   DMLabelCreate("subpoint_map", &subpointMap);
2819:   DMPlexSetSubpointMap(subdm, subpointMap);
2820:   DMLabelDestroy(&subpointMap);
2821:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
2822:   /* Setup chart */
2823:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2824:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2825:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2826:   DMPlexSetVTKCellHeight(subdm, 1);
2827:   /* Set cone sizes */
2828:   firstSubVertex = numSubCells;
2829:   firstSubFace   = numSubCells+numSubVertices;
2830:   newFacePoint   = firstSubFace;
2831:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2832:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2833:   for (c = 0; c < numSubCells; ++c) {
2834:     DMPlexSetConeSize(subdm, c, 1);
2835:   }
2836:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2837:     DMPlexSetConeSize(subdm, f, nFV);
2838:   }
2839:   DMSetUp(subdm);
2840:   /* Create face cones */
2841:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2842:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2843:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2844:   for (c = 0; c < numSubCells; ++c) {
2845:     const PetscInt  cell    = subCells[c];
2846:     const PetscInt  subcell = c;
2847:     const PetscInt *cone, *cells;
2848:     PetscInt        numCells, subVertex, p, v;

2850:     if (cell < cMax) continue;
2851:     DMPlexGetCone(dm, cell, &cone);
2852:     for (v = 0; v < nFV; ++v) {
2853:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
2854:       subface[v] = firstSubVertex+subVertex;
2855:     }
2856:     DMPlexSetCone(subdm, newFacePoint, subface);
2857:     DMPlexSetCone(subdm, subcell, &newFacePoint);
2858:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
2859:     /* Not true in parallel
2860:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2861:     for (p = 0; p < numCells; ++p) {
2862:       PetscInt negsubcell;

2864:       if (cells[p] >= cMax) continue;
2865:       /* I know this is a crap search */
2866:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2867:         if (subCells[negsubcell] == cells[p]) break;
2868:       }
2869:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2870:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
2871:     }
2872:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
2873:     ++newFacePoint;
2874:   }
2875:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2876:   DMPlexSymmetrize(subdm);
2877:   DMPlexStratify(subdm);
2878:   /* Build coordinates */
2879:   {
2880:     PetscSection coordSection, subCoordSection;
2881:     Vec          coordinates, subCoordinates;
2882:     PetscScalar *coords, *subCoords;
2883:     PetscInt     numComp, coordSize, v;

2885:     DMGetCoordinateSection(dm, &coordSection);
2886:     DMGetCoordinatesLocal(dm, &coordinates);
2887:     DMGetCoordinateSection(subdm, &subCoordSection);
2888:     PetscSectionSetNumFields(subCoordSection, 1);
2889:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2890:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2891:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2892:     for (v = 0; v < numSubVertices; ++v) {
2893:       const PetscInt vertex    = subVertices[v];
2894:       const PetscInt subvertex = firstSubVertex+v;
2895:       PetscInt       dof;

2897:       PetscSectionGetDof(coordSection, vertex, &dof);
2898:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2899:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2900:     }
2901:     PetscSectionSetUp(subCoordSection);
2902:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2903:     VecCreate(comm, &subCoordinates);
2904:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2905:     VecSetType(subCoordinates,VECSTANDARD);
2906:     VecGetArray(coordinates,    &coords);
2907:     VecGetArray(subCoordinates, &subCoords);
2908:     for (v = 0; v < numSubVertices; ++v) {
2909:       const PetscInt vertex    = subVertices[v];
2910:       const PetscInt subvertex = firstSubVertex+v;
2911:       PetscInt       dof, off, sdof, soff, d;

2913:       PetscSectionGetDof(coordSection, vertex, &dof);
2914:       PetscSectionGetOffset(coordSection, vertex, &off);
2915:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2916:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2917:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2918:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2919:     }
2920:     VecRestoreArray(coordinates,    &coords);
2921:     VecRestoreArray(subCoordinates, &subCoords);
2922:     DMSetCoordinatesLocal(subdm, subCoordinates);
2923:     VecDestroy(&subCoordinates);
2924:   }
2925:   /* Build SF */
2926:   CHKMEMQ;
2927:   {
2928:     PetscSF            sfPoint, sfPointSub;
2929:     const PetscSFNode *remotePoints;
2930:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
2931:     const PetscInt    *localPoints;
2932:     PetscInt          *slocalPoints;
2933:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
2934:     PetscMPIInt        rank;

2936:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2937:     DMGetPointSF(dm, &sfPoint);
2938:     DMGetPointSF(subdm, &sfPointSub);
2939:     DMPlexGetChart(dm, &pStart, &pEnd);
2940:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2941:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2942:     if (numRoots >= 0) {
2943:       /* Only vertices should be shared */
2944:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2945:       for (p = 0; p < pEnd-pStart; ++p) {
2946:         newLocalPoints[p].rank  = -2;
2947:         newLocalPoints[p].index = -2;
2948:       }
2949:       /* Set subleaves */
2950:       for (l = 0; l < numLeaves; ++l) {
2951:         const PetscInt point    = localPoints[l];
2952:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2954:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
2955:         if (subPoint < 0) continue;
2956:         newLocalPoints[point-pStart].rank  = rank;
2957:         newLocalPoints[point-pStart].index = subPoint;
2958:         ++numSubLeaves;
2959:       }
2960:       /* Must put in owned subpoints */
2961:       for (p = pStart; p < pEnd; ++p) {
2962:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

2964:         if (subPoint < 0) {
2965:           newOwners[p-pStart].rank  = -3;
2966:           newOwners[p-pStart].index = -3;
2967:         } else {
2968:           newOwners[p-pStart].rank  = rank;
2969:           newOwners[p-pStart].index = subPoint;
2970:         }
2971:       }
2972:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2973:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2974:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2975:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2976:       PetscMalloc1(numSubLeaves,    &slocalPoints);
2977:       PetscMalloc1(numSubLeaves, &sremotePoints);
2978:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
2979:         const PetscInt point    = localPoints[l];
2980:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2982:         if (subPoint < 0) continue;
2983:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
2984:         slocalPoints[sl]        = subPoint;
2985:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
2986:         sremotePoints[sl].index = newLocalPoints[point].index;
2987:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
2988:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
2989:         ++sl;
2990:       }
2991:       PetscFree2(newLocalPoints,newOwners);
2992:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
2993:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
2994:     }
2995:   }
2996:   CHKMEMQ;
2997:   /* Cleanup */
2998:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2999:   ISDestroy(&subvertexIS);
3000:   PetscFree(subCells);
3001:   return(0);
3002: }

3006: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char label[], PetscInt value, DM subdm)
3007: {
3008:   MPI_Comm         comm;
3009:   DMLabel          subpointMap;
3010:   IS              *subpointIS;
3011:   const PetscInt **subpoints;
3012:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3013:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
3014:   PetscErrorCode   ierr;

3017:   PetscObjectGetComm((PetscObject)dm,&comm);
3018:   /* Create subpointMap which marks the submesh */
3019:   DMLabelCreate("subpoint_map", &subpointMap);
3020:   DMPlexSetSubpointMap(subdm, subpointMap);
3021:   DMLabelDestroy(&subpointMap);
3022:   DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
3023:   /* Setup chart */
3024:   DMGetDimension(dm, &dim);
3025:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3026:   for (d = 0; d <= dim; ++d) {
3027:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3028:     totSubPoints += numSubPoints[d];
3029:   }
3030:   DMPlexSetChart(subdm, 0, totSubPoints);
3031:   DMPlexSetVTKCellHeight(subdm, 1);
3032:   /* Set cone sizes */
3033:   firstSubPoint[dim] = 0;
3034:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
3035:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
3036:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3037:   for (d = 0; d <= dim; ++d) {
3038:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3039:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3040:   }
3041:   for (d = 0; d <= dim; ++d) {
3042:     for (p = 0; p < numSubPoints[d]; ++p) {
3043:       const PetscInt  point    = subpoints[d][p];
3044:       const PetscInt  subpoint = firstSubPoint[d] + p;
3045:       const PetscInt *cone;
3046:       PetscInt        coneSize, coneSizeNew, c, val;

3048:       DMPlexGetConeSize(dm, point, &coneSize);
3049:       DMPlexSetConeSize(subdm, subpoint, coneSize);
3050:       if (d == dim) {
3051:         DMPlexGetCone(dm, point, &cone);
3052:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3053:           DMLabelGetValue(subpointMap, cone[c], &val);
3054:           if (val >= 0) coneSizeNew++;
3055:         }
3056:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3057:       }
3058:     }
3059:   }
3060:   DMSetUp(subdm);
3061:   /* Set cones */
3062:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3063:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3064:   for (d = 0; d <= dim; ++d) {
3065:     for (p = 0; p < numSubPoints[d]; ++p) {
3066:       const PetscInt  point    = subpoints[d][p];
3067:       const PetscInt  subpoint = firstSubPoint[d] + p;
3068:       const PetscInt *cone, *ornt;
3069:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

3071:       DMPlexGetConeSize(dm, point, &coneSize);
3072:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3073:       DMPlexGetCone(dm, point, &cone);
3074:       DMPlexGetConeOrientation(dm, point, &ornt);
3075:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3076:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3077:         if (subc >= 0) {
3078:           coneNew[coneSizeNew]   = firstSubPoint[d-1] + subc;
3079:           orntNew[coneSizeNew++] = ornt[c];
3080:         }
3081:       }
3082:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3083:       DMPlexSetCone(subdm, subpoint, coneNew);
3084:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3085:     }
3086:   }
3087:   PetscFree2(coneNew,orntNew);
3088:   DMPlexSymmetrize(subdm);
3089:   DMPlexStratify(subdm);
3090:   /* Build coordinates */
3091:   {
3092:     PetscSection coordSection, subCoordSection;
3093:     Vec          coordinates, subCoordinates;
3094:     PetscScalar *coords, *subCoords;
3095:     PetscInt     numComp, coordSize;

3097:     DMGetCoordinateSection(dm, &coordSection);
3098:     DMGetCoordinatesLocal(dm, &coordinates);
3099:     DMGetCoordinateSection(subdm, &subCoordSection);
3100:     PetscSectionSetNumFields(subCoordSection, 1);
3101:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3102:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3103:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3104:     for (v = 0; v < numSubPoints[0]; ++v) {
3105:       const PetscInt vertex    = subpoints[0][v];
3106:       const PetscInt subvertex = firstSubPoint[0]+v;
3107:       PetscInt       dof;

3109:       PetscSectionGetDof(coordSection, vertex, &dof);
3110:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3111:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3112:     }
3113:     PetscSectionSetUp(subCoordSection);
3114:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3115:     VecCreate(comm, &subCoordinates);
3116:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3117:     VecSetType(subCoordinates,VECSTANDARD);
3118:     VecGetArray(coordinates,    &coords);
3119:     VecGetArray(subCoordinates, &subCoords);
3120:     for (v = 0; v < numSubPoints[0]; ++v) {
3121:       const PetscInt vertex    = subpoints[0][v];
3122:       const PetscInt subvertex = firstSubPoint[0]+v;
3123:       PetscInt dof, off, sdof, soff, d;

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

3148:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3149:     DMGetPointSF(dm, &sfPoint);
3150:     DMGetPointSF(subdm, &sfPointSub);
3151:     DMPlexGetChart(dm, &pStart, &pEnd);
3152:     DMPlexGetChart(subdm, NULL, &numSubroots);
3153:     DMPlexCreateSubpointIS(subdm, &subpIS);
3154:     if (subpIS) {
3155:       ISGetIndices(subpIS, &subpoints);
3156:       ISGetLocalSize(subpIS, &numSubpoints);
3157:     }
3158:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3159:     if (numRoots >= 0) {
3160:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3161:       for (p = 0; p < pEnd-pStart; ++p) {
3162:         newLocalPoints[p].rank  = -2;
3163:         newLocalPoints[p].index = -2;
3164:       }
3165:       /* Set subleaves */
3166:       for (l = 0; l < numLeaves; ++l) {
3167:         const PetscInt point    = localPoints[l];
3168:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3170:         if (subpoint < 0) continue;
3171:         newLocalPoints[point-pStart].rank  = rank;
3172:         newLocalPoints[point-pStart].index = subpoint;
3173:         ++numSubleaves;
3174:       }
3175:       /* Must put in owned subpoints */
3176:       for (p = pStart; p < pEnd; ++p) {
3177:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3179:         if (subpoint < 0) {
3180:           newOwners[p-pStart].rank  = -3;
3181:           newOwners[p-pStart].index = -3;
3182:         } else {
3183:           newOwners[p-pStart].rank  = rank;
3184:           newOwners[p-pStart].index = subpoint;
3185:         }
3186:       }
3187:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3188:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3189:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3190:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3191:       PetscMalloc1(numSubleaves, &slocalPoints);
3192:       PetscMalloc1(numSubleaves, &sremotePoints);
3193:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3194:         const PetscInt point    = localPoints[l];
3195:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3197:         if (subpoint < 0) continue;
3198:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3199:         slocalPoints[sl]        = subpoint;
3200:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3201:         sremotePoints[sl].index = newLocalPoints[point].index;
3202:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3203:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3204:         ++sl;
3205:       }
3206:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3207:       PetscFree2(newLocalPoints,newOwners);
3208:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3209:     }
3210:     if (subpIS) {
3211:       ISRestoreIndices(subpIS, &subpoints);
3212:       ISDestroy(&subpIS);
3213:     }
3214:   }
3215:   /* Cleanup */
3216:   for (d = 0; d <= dim; ++d) {
3217:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3218:     ISDestroy(&subpointIS[d]);
3219:   }
3220:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3221:   return(0);
3222: }

3226: /*
3227:   DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label an be given to restrict the cells.

3229:   Input Parameters:
3230: + dm          - The original mesh
3231: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3232: . label       - A label name, or NULL
3233: - value  - A label value

3235:   Output Parameter:
3236: . subdm - The surface mesh

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

3240:   Level: developer

3242: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3243: */
3244: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3245: {
3246:   PetscInt       dim, depth;

3252:   DMGetDimension(dm, &dim);
3253:   DMPlexGetDepth(dm, &depth);
3254:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3255:   DMSetType(*subdm, DMPLEX);
3256:   DMSetDimension(*subdm, dim-1);
3257:   if (depth == dim) {
3258:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3259:   } else {
3260:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3261:   }
3262:   return(0);
3263: }

3267: /*@
3268:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3270:   Input Parameter:
3271: . dm - The submesh DM

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

3276:   Level: developer

3278: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3279: @*/
3280: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3281: {
3282:   DM_Plex *mesh = (DM_Plex*) dm->data;

3287:   *subpointMap = mesh->subpointMap;
3288:   return(0);
3289: }

3293: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3294: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3295: {
3296:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3297:   DMLabel        tmp;

3302:   tmp  = mesh->subpointMap;
3303:   mesh->subpointMap = subpointMap;
3304:   ++mesh->subpointMap->refct;
3305:   DMLabelDestroy(&tmp);
3306:   return(0);
3307: }

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

3314:   Input Parameter:
3315: . dm - The submesh DM

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

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

3322:   Level: developer

3324: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3325: @*/
3326: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3327: {
3328:   MPI_Comm        comm;
3329:   DMLabel         subpointMap;
3330:   IS              is;
3331:   const PetscInt *opoints;
3332:   PetscInt       *points, *depths;
3333:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3334:   PetscErrorCode  ierr;

3339:   PetscObjectGetComm((PetscObject)dm,&comm);
3340:   *subpointIS = NULL;
3341:   DMPlexGetSubpointMap(dm, &subpointMap);
3342:   DMPlexGetDepth(dm, &depth);
3343:   if (subpointMap && depth >= 0) {
3344:     DMPlexGetChart(dm, &pStart, &pEnd);
3345:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3346:     DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3347:     depths[0] = depth;
3348:     depths[1] = 0;
3349:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3350:     PetscMalloc1(pEnd, &points);
3351:     for(d = 0, off = 0; d <= depth; ++d) {
3352:       const PetscInt dep = depths[d];

3354:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3355:       DMLabelGetStratumSize(subpointMap, dep, &n);
3356:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3357:         if (n != depEnd-depStart) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d at depth %d should be %d", n, dep, depEnd-depStart);
3358:       } else {
3359:         if (!n) {
3360:           if (d == 0) {
3361:             /* Missing cells */
3362:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3363:           } else {
3364:             /* Missing faces */
3365:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3366:           }
3367:         }
3368:       }
3369:       if (n) {
3370:         DMLabelGetStratumIS(subpointMap, dep, &is);
3371:         ISGetIndices(is, &opoints);
3372:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3373:         ISRestoreIndices(is, &opoints);
3374:         ISDestroy(&is);
3375:       }
3376:     }
3377:     DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3378:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3379:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3380:   }
3381:   return(0);
3382: }