Actual source code: plexsubmesh.c

petsc-master 2014-10-23
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:   VecSetType(newCoordinates,VECSTANDARD);
290:   DMSetCoordinatesLocal(dmNew, newCoordinates);
291:   DMGetCoordinatesLocal(dm, &coordinates);
292:   VecGetArray(coordinates, &coords);
293:   VecGetArray(newCoordinates, &newCoords);
294:   for (v = vStart; v < vEnd; ++v) {
295:     PetscInt dof, off, noff, d;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

584:   Collective on dm

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

590:   Output Parameters:
591: + numGhostCells - The number of ghost cells added to the DM
592: - dmGhosted - The new DM

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

596:   Level: developer

598: .seealso: DMCreate()
599: @*/
600: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
601: {
602:   DM             gdm;
603:   DMLabel        label;
604:   const char    *name = labelName ? labelName : "Face Sets";
605:   PetscInt       dim;
606:   PetscBool      flag;

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

636: /*
637:   We are adding three kinds of points here:
638:     Replicated:     Copies of points which exist in the mesh, such as vertices identified across a fault
639:     Non-replicated: Points which exist on the fault, but are not replicated
640:     Hybrid:         Entirely new points, such as cohesive cells

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

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

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

736:       DMPlexGetConeSize(dm, oldp, &coneSize);
737:       DMPlexSetConeSize(sdm, splitp, coneSize);
738:       DMPlexGetSupportSize(dm, oldp, &supportSize);
739:       DMPlexSetSupportSize(sdm, splitp, supportSize);
740:       if (dep == depth-1) {
741:         const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

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

748:         DMPlexGetSupport(dm, oldp, &support);
749:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
750:           PetscInt val;

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

767:         DMPlexGetSupport(dm, oldp, &support);
768:         for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
769:           PetscInt val;

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

793:       DMPlexGetConeSize(dm, oldp, &coneSize);
794:       DMPlexGetSupportSize(dm, oldp, &supportSize);
795:       DMPlexGetSupport(dm, oldp, &support);
796:       if (dep == 0) {
797:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

799:         /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
800:         for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
801:           PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
802:           if (e >= 0) ++qf;
803:         }
804:         DMPlexSetSupportSize(sdm, newp, qf+2);
805:         /* Add hybrid edge */
806:         DMPlexSetConeSize(sdm, hybedge, 2);
807:         for (e = 0, qf = 0; e < supportSize; ++e) {
808:           PetscInt val;

810:           DMLabelGetValue(label, support[e], &val);
811:           /* Split and unsplit edges produce hybrid faces */
812:           if (val == 1) ++qf;
813:           if (val == (shift2 + 1)) ++qf;
814:         }
815:         DMPlexSetSupportSize(sdm, hybedge, qf);
816:       } else if (dep == dim-2) {
817:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
818:         PetscInt       val;

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

851:       DMPlexGetConeSize(dm, oldp, &coneSize);
852:       DMPlexGetCone(dm, oldp, &cone);
853:       DMPlexGetConeOrientation(dm, oldp, &ornt);
854:       DMPlexGetSupportSize(dm, oldp, &supportSize);
855:       DMPlexGetSupport(dm, oldp, &support);
856:       if (dep == depth-1) {
857:         PetscBool       hasUnsplit = PETSC_FALSE;
858:         const PetscInt  hybcell    = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
859:         const PetscInt *supportF;

861:         /* Split face:       copy in old face to new face to start */
862:         DMPlexGetSupport(sdm, newp,  &supportF);
863:         DMPlexSetSupport(sdm, splitp, supportF);
864:         /* Split old face:   old vertices/edges in cone so no change */
865:         /* Split new face:   new vertices/edges in cone */
866:         for (q = 0; q < coneSize; ++q) {
867:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
868:           if (v < 0) {
869:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
870:             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);
871:             coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
872:             hasUnsplit   = PETSC_TRUE;
873:           } else {
874:             coneNew[2+q] = v + pMaxNew[dep-1];
875:             if (dep > 1) {
876:               const PetscInt *econe;
877:               PetscInt        econeSize, r, vs, vu;

879:               DMPlexGetConeSize(dm, cone[q], &econeSize);
880:               DMPlexGetCone(dm, cone[q], &econe);
881:               for (r = 0; r < econeSize; ++r) {
882:                 PetscFindInt(econe[r], numSplitPoints[dep-2],   splitPoints[dep-2],   &vs);
883:                 PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
884:                 if (vs >= 0) continue;
885:                 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);
886:                 hasUnsplit   = PETSC_TRUE;
887:               }
888:             }
889:           }
890:         }
891:         DMPlexSetCone(sdm, splitp, &coneNew[2]);
892:         DMPlexSetConeOrientation(sdm, splitp, ornt);
893:         /* Face support */
894:         for (s = 0; s < supportSize; ++s) {
895:           PetscInt val;

897:           DMLabelGetValue(label, support[s], &val);
898:           if (val < 0) {
899:             /* Split old face:   Replace negative side cell with cohesive cell */
900:              DMPlexInsertSupport(sdm, newp, s, hybcell);
901:           } else {
902:             /* Split new face:   Replace positive side cell with cohesive cell */
903:             DMPlexInsertSupport(sdm, splitp, s, hybcell);
904:             /* Get orientation for cohesive face */
905:             {
906:               const PetscInt *ncone, *nconeO;
907:               PetscInt        nconeSize, nc;

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

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

948:           DMLabelGetValue(label, support[e], &val);
949:           if ((val == 1) || (val == (shift + 1))) {
950:             supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
951:           }
952:         }
953:         supportNew[qn] = hybedge;
954:         DMPlexSetSupport(sdm, newp, supportNew);
955:         /* Split new vertex: Edges in new split faces and new cohesive edge */
956:         for (e = 0, qp = 0; e < supportSize; ++e) {
957:           PetscInt val, edge;

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

977:           DMLabelGetValue(label, support[e], &val);
978:           if (val == 1) {
979:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
980:             if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
981:             supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
982:           }
983:         }
984:         DMPlexSetSupport(sdm, hybedge, supportNew);
985:       } else if (dep == dim-2) {
986:         const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];

988:         /* Split old edge:   old vertices in cone so no change */
989:         /* Split new edge:   new vertices in cone */
990:         for (q = 0; q < coneSize; ++q) {
991:           PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
992:           if (v < 0) {
993:             PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
994:             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);
995:             coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthMax, depthEnd, depthShift) /*cone[q] + depthOffset[dep-1]*/;
996:           } else {
997:             coneNew[q] = v + pMaxNew[dep-1];
998:           }
999:         }
1000:         DMPlexSetCone(sdm, splitp, coneNew);
1001:         /* Split old edge: Faces in positive side cells and old split faces */
1002:         for (e = 0, q = 0; e < supportSize; ++e) {
1003:           PetscInt val;

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

1018:           DMLabelGetValue(label, support[e], &val);
1019:           if (val == dim-1) {
1020:             PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1021:             if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1022:             supportNew[q++] = face + pMaxNew[dep+1];
1023:           } else if (val == -(shift + dim-1)) {
1024:             supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthMax, depthEnd, depthShift) /*support[e] + depthOffset[dep+1]*/;
1025:           }
1026:         }
1027:         supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1028:         DMPlexSetSupport(sdm, splitp, supportNew);
1029:         /* Hybrid face */
1030:         coneNew[0] = newp;
1031:         coneNew[1] = splitp;
1032:         for (v = 0; v < coneSize; ++v) {
1033:           PetscInt vertex;
1034:           PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1035:           if (vertex < 0) {
1036:             PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1037:             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);
1038:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1039:           } else {
1040:             coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1041:           }
1042:         }
1043:         DMPlexSetCone(sdm, hybface, coneNew);
1044:         for (e = 0, qf = 0; e < supportSize; ++e) {
1045:           PetscInt val, face;

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

1065:       DMPlexGetConeSize(dm, oldp, &coneSize);
1066:       DMPlexGetCone(dm, oldp, &cone);
1067:       DMPlexGetConeOrientation(dm, oldp, &ornt);
1068:       DMPlexGetSupportSize(dm, oldp, &supportSize);
1069:       DMPlexGetSupport(dm, oldp, &support);
1070:       if (dep == 0) {
1071:         const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];

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

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

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

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

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

1159:     if (dep >= 0) continue;
1160:     DMLabelGetStratumIS(label, dep, &pIS);
1161:     if (!pIS) continue;
1162:     dep  = -dep - shift;
1163:     ISGetLocalSize(pIS, &numPoints);
1164:     ISGetIndices(pIS, &points);
1165:     for (p = 0; p < numPoints; ++p) {
1166:       const PetscInt  oldp = points[p];
1167:       const PetscInt  newp = DMPlexShiftPoint_Internal(oldp, depth, depthMax, depthEnd, depthShift) /*depthOffset[dep] + oldp*/;
1168:       const PetscInt *cone;
1169:       PetscInt        coneSize, c;
1170:       /* PetscBool       replaced = PETSC_FALSE; */

1172:       /* Negative edge: replace split vertex */
1173:       /* Negative cell: replace split face */
1174:       DMPlexGetConeSize(sdm, newp, &coneSize);
1175:       DMPlexGetCone(sdm, newp, &cone);
1176:       for (c = 0; c < coneSize; ++c) {
1177:         const PetscInt coldp = cone[c] - depthOffset[dep-1];
1178:         PetscInt       csplitp, cp, val;

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

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

1224:       for (l = 0; l < numLabels; ++l) {
1225:         DMLabel     mlabel;
1226:         const char *lname;
1227:         PetscInt    val;
1228:         PetscBool   isDepth;

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

1255:     if ((dep < 0) || (dep > depth)) continue;
1256:     if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1257:     ISDestroy(&splitIS[dep]);
1258:     if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1259:     ISDestroy(&unsplitIS[dep]);
1260:   }
1261:   if (label) {
1262:     ISRestoreIndices(valueIS, &values);
1263:     ISDestroy(&valueIS);
1264:   }
1265:   for (d = 0; d <= depth; ++d) {
1266:     DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1267:     pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1268:   }
1269:   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);
1270:   PetscFree3(coneNew, coneONew, supportNew);
1271:   PetscFree6(depthMax, depthEnd, depthShift, depthOffset, pMaxNew, numHybridPointsOld);
1272:   PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1273:   return(0);
1274: }

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

1281:   Collective on dm

1283:   Input Parameters:
1284: + dm - The original DM
1285: - label - The label specifying the boundary faces (this could be auto-generated)

1287:   Output Parameters:
1288: - dmSplit - The new DM

1290:   Level: developer

1292: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1293: @*/
1294: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1295: {
1296:   DM             sdm;
1297:   PetscInt       dim;

1303:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1304:   DMSetType(sdm, DMPLEX);
1305:   DMGetDimension(dm, &dim);
1306:   DMSetDimension(sdm, dim);
1307:   switch (dim) {
1308:   case 2:
1309:   case 3:
1310:     DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1311:     break;
1312:   default:
1313:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1314:   }
1315:   *dmSplit = sdm;
1316:   return(0);
1317: }

1321: /* Returns the side of the surface for a given cell with a face on the surface */
1322: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1323: {
1324:   const PetscInt *cone, *ornt;
1325:   PetscInt        dim, coneSize, c;
1326:   PetscErrorCode  ierr;

1329:   *pos = PETSC_TRUE;
1330:   DMGetDimension(dm, &dim);
1331:   DMPlexGetConeSize(dm, cell, &coneSize);
1332:   DMPlexGetCone(dm, cell, &cone);
1333:   DMPlexGetConeOrientation(dm, cell, &ornt);
1334:   for (c = 0; c < coneSize; ++c) {
1335:     if (cone[c] == face) {
1336:       PetscInt o = ornt[c];

1338:       if (subdm) {
1339:         const PetscInt *subcone, *subornt;
1340:         PetscInt        subpoint, subface, subconeSize, sc;

1342:         PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1343:         PetscFindInt(face, numSubpoints, subpoints, &subface);
1344:         DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1345:         DMPlexGetCone(subdm, subpoint, &subcone);
1346:         DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1347:         for (sc = 0; sc < subconeSize; ++sc) {
1348:           if (subcone[sc] == subface) {
1349:             o = subornt[0];
1350:             break;
1351:           }
1352:         }
1353:         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);
1354:       }
1355:       if (o >= 0) *pos = PETSC_TRUE;
1356:       else        *pos = PETSC_FALSE;
1357:       break;
1358:     }
1359:   }
1360:   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);
1361:   return(0);
1362: }

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

1370:   Input Parameters:
1371: + dm     - The DM
1372: . label  - A DMLabel marking the surface
1373: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1374: . flip   - Flag to flip the submesh normal and replace points on the other side
1375: - subdm  - The subDM associated with the label, or NULL

1377:   Output Parameter:
1378: . label - A DMLabel marking all surface points

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

1382:   Level: developer

1384: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1385: @*/
1386: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1387: {
1388:   DMLabel         depthLabel;
1389:   IS              dimIS, subpointIS, facePosIS, faceNegIS, crossEdgeIS = NULL;
1390:   const PetscInt *points, *subpoints;
1391:   const PetscInt  rev   = flip ? -1 : 1;
1392:   PetscInt        shift = 100, shift2 = 200, dim, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1393:   PetscErrorCode  ierr;

1396:   DMPlexGetDepthLabel(dm, &depthLabel);
1397:   DMGetDimension(dm, &dim);
1398:   if (subdm) {
1399:     DMPlexCreateSubpointIS(subdm, &subpointIS);
1400:     if (subpointIS) {
1401:       ISGetLocalSize(subpointIS, &numSubpoints);
1402:       ISGetIndices(subpointIS, &subpoints);
1403:     }
1404:   }
1405:   /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1406:   DMLabelGetStratumIS(label, dim-1, &dimIS);
1407:   if (!dimIS) return(0);
1408:   ISGetLocalSize(dimIS, &numPoints);
1409:   ISGetIndices(dimIS, &points);
1410:   for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1411:     const PetscInt *support;
1412:     PetscInt        supportSize, s;

1414:     DMPlexGetSupportSize(dm, points[p], &supportSize);
1415:     if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1416:     DMPlexGetSupport(dm, points[p], &support);
1417:     for (s = 0; s < supportSize; ++s) {
1418:       const PetscInt *cone;
1419:       PetscInt        coneSize, c;
1420:       PetscBool       pos;

1422:       GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1423:       if (pos) {DMLabelSetValue(label, support[s],  rev*(shift+dim));}
1424:       else     {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1425:       if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1426:       /* Put faces touching the fault in the label */
1427:       DMPlexGetConeSize(dm, support[s], &coneSize);
1428:       DMPlexGetCone(dm, support[s], &cone);
1429:       for (c = 0; c < coneSize; ++c) {
1430:         const PetscInt point = cone[c];

1432:         DMLabelGetValue(label, point, &val);
1433:         if (val == -1) {
1434:           PetscInt *closure = NULL;
1435:           PetscInt  closureSize, cl;

1437:           DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1438:           for (cl = 0; cl < closureSize*2; cl += 2) {
1439:             const PetscInt clp  = closure[cl];
1440:             PetscInt       bval = -1;

1442:             DMLabelGetValue(label, clp, &val);
1443:             if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1444:             if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1445:               DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1446:               break;
1447:             }
1448:           }
1449:           DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1450:         }
1451:       }
1452:     }
1453:   }
1454:   ISRestoreIndices(dimIS, &points);
1455:   ISDestroy(&dimIS);
1456:   if (subdm) {
1457:     if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1458:     ISDestroy(&subpointIS);
1459:   }
1460:   /* Mark boundary points as unsplit */
1461:   if (blabel) {
1462:     DMLabelGetStratumIS(blabel, 1, &dimIS);
1463:     ISGetLocalSize(dimIS, &numPoints);
1464:     ISGetIndices(dimIS, &points);
1465:     for (p = 0; p < numPoints; ++p) {
1466:       const PetscInt point = points[p];
1467:       PetscInt       val, bval;

1469:       DMLabelGetValue(blabel, point, &bval);
1470:       if (bval >= 0) {
1471:         DMLabelGetValue(label, point, &val);
1472:         if ((val < 0) || (val > dim)) {
1473:           /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1474:           DMLabelClearValue(blabel, point, bval);
1475:         }
1476:       }
1477:     }
1478:     for (p = 0; p < numPoints; ++p) {
1479:       const PetscInt point = points[p];
1480:       PetscInt       val, bval;

1482:       DMLabelGetValue(blabel, point, &bval);
1483:       if (bval >= 0) {
1484:         const PetscInt *cone,    *support;
1485:         PetscInt        coneSize, supportSize, s, valA, valB, valE;

1487:         /* Mark as unsplit */
1488:         DMLabelGetValue(label, point, &val);
1489:         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);
1490:         DMLabelClearValue(label, point, val);
1491:         DMLabelSetValue(label, point, shift2+val);
1492:         /* Check for cross-edge */
1493:         if (val != 0) continue;
1494:         DMPlexGetSupport(dm, point, &support);
1495:         DMPlexGetSupportSize(dm, point, &supportSize);
1496:         for (s = 0; s < supportSize; ++s) {
1497:           DMPlexGetCone(dm, support[s], &cone);
1498:           DMPlexGetConeSize(dm, support[s], &coneSize);
1499:           if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1500:           DMLabelGetValue(blabel, cone[0], &valA);
1501:           DMLabelGetValue(blabel, cone[1], &valB);
1502:           DMLabelGetValue(blabel, support[s], &valE);
1503:           if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1504:         }
1505:       }
1506:     }
1507:     ISRestoreIndices(dimIS, &points);
1508:     ISDestroy(&dimIS);
1509:   }
1510:   /* Search for other cells/faces/edges connected to the fault by a vertex */
1511:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1512:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1513:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1514:   cMax = cMax < 0 ? cEnd : cMax;
1515:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1516:   DMLabelGetStratumIS(label, 0, &dimIS);
1517:   if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1518:   if (dimIS && crossEdgeIS) {
1519:     IS vertIS = dimIS;

1521:     ISExpand(vertIS, crossEdgeIS, &dimIS);
1522:     ISDestroy(&crossEdgeIS);
1523:     ISDestroy(&vertIS);
1524:   }
1525:   if (!dimIS) return(0);
1526:   ISGetLocalSize(dimIS, &numPoints);
1527:   ISGetIndices(dimIS, &points);
1528:   for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1529:     PetscInt *star = NULL;
1530:     PetscInt  starSize, s;
1531:     PetscInt  again = 1;  /* 0: Finished 1: Keep iterating after a change 2: No change */

1533:     /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1534:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1535:     while (again) {
1536:       if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1537:       again = 0;
1538:       for (s = 0; s < starSize*2; s += 2) {
1539:         const PetscInt  point = star[s];
1540:         const PetscInt *cone;
1541:         PetscInt        coneSize, c;

1543:         if ((point < cStart) || (point >= cMax)) continue;
1544:         DMLabelGetValue(label, point, &val);
1545:         if (val != -1) continue;
1546:         again = again == 1 ? 1 : 2;
1547:         DMPlexGetConeSize(dm, point, &coneSize);
1548:         DMPlexGetCone(dm, point, &cone);
1549:         for (c = 0; c < coneSize; ++c) {
1550:           DMLabelGetValue(label, cone[c], &val);
1551:           if (val != -1) {
1552:             const PetscInt *ccone;
1553:             PetscInt        cconeSize, cc, side;

1555:             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);
1556:             if (val > 0) side =  1;
1557:             else         side = -1;
1558:             DMLabelSetValue(label, point, side*(shift+dim));
1559:             /* Mark cell faces which touch the fault */
1560:             DMPlexGetConeSize(dm, point, &cconeSize);
1561:             DMPlexGetCone(dm, point, &ccone);
1562:             for (cc = 0; cc < cconeSize; ++cc) {
1563:               PetscInt *closure = NULL;
1564:               PetscInt  closureSize, cl;

1566:               DMLabelGetValue(label, ccone[cc], &val);
1567:               if (val != -1) continue;
1568:               DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1569:               for (cl = 0; cl < closureSize*2; cl += 2) {
1570:                 const PetscInt clp = closure[cl];

1572:                 DMLabelGetValue(label, clp, &val);
1573:                 if (val == -1) continue;
1574:                 DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1575:                 break;
1576:               }
1577:               DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1578:             }
1579:             again = 1;
1580:             break;
1581:           }
1582:         }
1583:       }
1584:     }
1585:     /* Classify the rest by cell membership */
1586:     for (s = 0; s < starSize*2; s += 2) {
1587:       const PetscInt point = star[s];

1589:       DMLabelGetValue(label, point, &val);
1590:       if (val == -1) {
1591:         PetscInt  *sstar = NULL;
1592:         PetscInt   sstarSize, ss;
1593:         PetscBool  marked = PETSC_FALSE;

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

1599:           if ((spoint < cStart) || (spoint >= cMax)) continue;
1600:           DMLabelGetValue(label, spoint, &val);
1601:           if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1602:           DMLabelGetValue(depthLabel, point, &dep);
1603:           if (val > 0) {
1604:             DMLabelSetValue(label, point,   shift+dep);
1605:           } else {
1606:             DMLabelSetValue(label, point, -(shift+dep));
1607:           }
1608:           marked = PETSC_TRUE;
1609:           break;
1610:         }
1611:         DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1612:         if (!marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1613:       }
1614:     }
1615:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1616:   }
1617:   ISRestoreIndices(dimIS, &points);
1618:   ISDestroy(&dimIS);
1619:   /* If any faces touching the fault divide cells on either side, split them */
1620:   DMLabelGetStratumIS(label,   shift+dim-1,  &facePosIS);
1621:   DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1622:   ISExpand(facePosIS, faceNegIS, &dimIS);
1623:   ISDestroy(&facePosIS);
1624:   ISDestroy(&faceNegIS);
1625:   ISGetLocalSize(dimIS, &numPoints);
1626:   ISGetIndices(dimIS, &points);
1627:   for (p = 0; p < numPoints; ++p) {
1628:     const PetscInt  point = points[p];
1629:     const PetscInt *support;
1630:     PetscInt        supportSize, valA, valB;

1632:     DMPlexGetSupportSize(dm, point, &supportSize);
1633:     if (supportSize != 2) continue;
1634:     DMPlexGetSupport(dm, point, &support);
1635:     DMLabelGetValue(label, support[0], &valA);
1636:     DMLabelGetValue(label, support[1], &valB);
1637:     if ((valA == -1) || (valB == -1)) continue;
1638:     if (valA*valB > 0) continue;
1639:     /* Split the face */
1640:     DMLabelGetValue(label, point, &valA);
1641:     DMLabelClearValue(label, point, valA);
1642:     DMLabelSetValue(label, point, dim-1);
1643:     /* Label its closure:
1644:       unmarked: label as unsplit
1645:       incident: relabel as split
1646:       split:    do nothing
1647:     */
1648:     {
1649:       PetscInt *closure = NULL;
1650:       PetscInt  closureSize, cl;

1652:       DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1653:       for (cl = 0; cl < closureSize*2; cl += 2) {
1654:         DMLabelGetValue(label, closure[cl], &valA);
1655:         if (valA == -1) { /* Mark as unsplit */
1656:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1657:           DMLabelSetValue(label, closure[cl], shift2+dep);
1658:         } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1659:           DMLabelGetValue(depthLabel, closure[cl], &dep);
1660:           DMLabelClearValue(label, closure[cl], valA);
1661:           DMLabelSetValue(label, closure[cl], dep);
1662:         }
1663:       }
1664:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1665:     }
1666:   }
1667:   ISRestoreIndices(dimIS, &points);
1668:   ISDestroy(&dimIS);
1669:   return(0);
1670: }

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

1677:   Collective on dm

1679:   Input Parameters:
1680: + dm - The original DM
1681: - labelName - The label specifying the interface vertices

1683:   Output Parameters:
1684: + hybridLabel - The label fully marking the interface
1685: - dmHybrid - The new DM

1687:   Level: developer

1689: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1690: @*/
1691: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1692: {
1693:   DM             idm;
1694:   DMLabel        subpointMap, hlabel;
1695:   PetscInt       dim;

1702:   DMGetDimension(dm, &dim);
1703:   DMPlexCreateSubmesh(dm, label, 1, &idm);
1704:   DMPlexOrient(idm);
1705:   DMPlexGetSubpointMap(idm, &subpointMap);
1706:   DMLabelDuplicate(subpointMap, &hlabel);
1707:   DMLabelClearStratum(hlabel, dim);
1708:   DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1709:   DMDestroy(&idm);
1710:   DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1711:   if (hybridLabel) *hybridLabel = hlabel;
1712:   else             {DMLabelDestroy(&hlabel);}
1713:   return(0);
1714: }

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

1720:      For any marked cell, the marked vertices constitute a single face
1721: */
1722: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1723: {
1724:   IS               subvertexIS = NULL;
1725:   const PetscInt  *subvertices;
1726:   PetscInt        *pStart, *pEnd, *pMax, pSize;
1727:   PetscInt         depth, dim, d, numSubVerticesInitial = 0, v;
1728:   PetscErrorCode   ierr;

1731:   *numFaces = 0;
1732:   *nFV      = 0;
1733:   DMPlexGetDepth(dm, &depth);
1734:   DMGetDimension(dm, &dim);
1735:   pSize = PetscMax(depth, dim) + 1;
1736:   PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1737:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1738:   for (d = 0; d <= depth; ++d) {
1739:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1740:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1741:   }
1742:   /* Loop over initial vertices and mark all faces in the collective star() */
1743:   if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1744:   if (subvertexIS) {
1745:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1746:     ISGetIndices(subvertexIS, &subvertices);
1747:   }
1748:   for (v = 0; v < numSubVerticesInitial; ++v) {
1749:     const PetscInt vertex = subvertices[v];
1750:     PetscInt      *star   = NULL;
1751:     PetscInt       starSize, s, numCells = 0, c;

1753:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1754:     for (s = 0; s < starSize*2; s += 2) {
1755:       const PetscInt point = star[s];
1756:       if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1757:     }
1758:     for (c = 0; c < numCells; ++c) {
1759:       const PetscInt cell    = star[c];
1760:       PetscInt      *closure = NULL;
1761:       PetscInt       closureSize, cl;
1762:       PetscInt       cellLoc, numCorners = 0, faceSize = 0;

1764:       DMLabelGetValue(subpointMap, cell, &cellLoc);
1765:       if (cellLoc == 2) continue;
1766:       if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1767:       DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1768:       for (cl = 0; cl < closureSize*2; cl += 2) {
1769:         const PetscInt point = closure[cl];
1770:         PetscInt       vertexLoc;

1772:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1773:           ++numCorners;
1774:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1775:           if (vertexLoc == value) closure[faceSize++] = point;
1776:         }
1777:       }
1778:       if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1779:       if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1780:       if (faceSize == *nFV) {
1781:         const PetscInt *cells = NULL;
1782:         PetscInt        numCells, nc;

1784:         ++(*numFaces);
1785:         for (cl = 0; cl < faceSize; ++cl) {
1786:           DMLabelSetValue(subpointMap, closure[cl], 0);
1787:         }
1788:         DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
1789:         for (nc = 0; nc < numCells; ++nc) {
1790:           DMLabelSetValue(subpointMap, cells[nc], 2);
1791:         }
1792:         DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
1793:       }
1794:       DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1795:     }
1796:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1797:   }
1798:   if (subvertexIS) {
1799:     ISRestoreIndices(subvertexIS, &subvertices);
1800:   }
1801:   ISDestroy(&subvertexIS);
1802:   PetscFree3(pStart,pEnd,pMax);
1803:   return(0);
1804: }

1808: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
1809: {
1810:   IS               subvertexIS;
1811:   const PetscInt  *subvertices;
1812:   PetscInt        *pStart, *pEnd, *pMax;
1813:   PetscInt         dim, d, numSubVerticesInitial = 0, v;
1814:   PetscErrorCode   ierr;

1817:   DMGetDimension(dm, &dim);
1818:   PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
1819:   DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
1820:   for (d = 0; d <= dim; ++d) {
1821:     DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1822:     if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1823:   }
1824:   /* Loop over initial vertices and mark all faces in the collective star() */
1825:   DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
1826:   if (subvertexIS) {
1827:     ISGetSize(subvertexIS, &numSubVerticesInitial);
1828:     ISGetIndices(subvertexIS, &subvertices);
1829:   }
1830:   for (v = 0; v < numSubVerticesInitial; ++v) {
1831:     const PetscInt vertex = subvertices[v];
1832:     PetscInt      *star   = NULL;
1833:     PetscInt       starSize, s, numFaces = 0, f;

1835:     DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1836:     for (s = 0; s < starSize*2; s += 2) {
1837:       const PetscInt point = star[s];
1838:       if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
1839:     }
1840:     for (f = 0; f < numFaces; ++f) {
1841:       const PetscInt face    = star[f];
1842:       PetscInt      *closure = NULL;
1843:       PetscInt       closureSize, c;
1844:       PetscInt       faceLoc;

1846:       DMLabelGetValue(subpointMap, face, &faceLoc);
1847:       if (faceLoc == dim-1) continue;
1848:       if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
1849:       DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1850:       for (c = 0; c < closureSize*2; c += 2) {
1851:         const PetscInt point = closure[c];
1852:         PetscInt       vertexLoc;

1854:         if ((point >= pStart[0]) && (point < pEnd[0])) {
1855:           DMLabelGetValue(vertexLabel, point, &vertexLoc);
1856:           if (vertexLoc != value) break;
1857:         }
1858:       }
1859:       if (c == closureSize*2) {
1860:         const PetscInt *support;
1861:         PetscInt        supportSize, s;

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

1866:           for (d = 0; d < dim; ++d) {
1867:             if ((point >= pStart[d]) && (point < pEnd[d])) {
1868:               DMLabelSetValue(subpointMap, point, d);
1869:               break;
1870:             }
1871:           }
1872:         }
1873:         DMPlexGetSupportSize(dm, face, &supportSize);
1874:         DMPlexGetSupport(dm, face, &support);
1875:         for (s = 0; s < supportSize; ++s) {
1876:           DMLabelSetValue(subpointMap, support[s], dim);
1877:         }
1878:       }
1879:       DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
1880:     }
1881:     DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1882:   }
1883:   if (subvertexIS) {
1884:     ISRestoreIndices(subvertexIS, &subvertices);
1885:   }
1886:   ISDestroy(&subvertexIS);
1887:   PetscFree3(pStart,pEnd,pMax);
1888:   return(0);
1889: }

1893: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
1894: {
1895:   DMLabel         label = NULL;
1896:   const PetscInt *cone;
1897:   PetscInt        dim, cMax, cEnd, c, subc = 0, p, coneSize;
1898:   PetscErrorCode  ierr;

1901:   *numFaces = 0;
1902:   *nFV = 0;
1903:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1904:   *subCells = NULL;
1905:   DMGetDimension(dm, &dim);
1906:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1907:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1908:   if (cMax < 0) return(0);
1909:   if (label) {
1910:     for (c = cMax; c < cEnd; ++c) {
1911:       PetscInt val;

1913:       DMLabelGetValue(label, c, &val);
1914:       if (val == value) {
1915:         ++(*numFaces);
1916:         DMPlexGetConeSize(dm, c, &coneSize);
1917:       }
1918:     }
1919:   } else {
1920:     *numFaces = cEnd - cMax;
1921:     DMPlexGetConeSize(dm, cMax, &coneSize);
1922:   }
1923:   *nFV = hasLagrange ? coneSize/3 : coneSize/2;
1924:   PetscMalloc1(*numFaces *2, subCells);
1925:   for (c = cMax; c < cEnd; ++c) {
1926:     const PetscInt *cells;
1927:     PetscInt        numCells;

1929:     if (label) {
1930:       PetscInt val;

1932:       DMLabelGetValue(label, c, &val);
1933:       if (val != value) continue;
1934:     }
1935:     DMPlexGetCone(dm, c, &cone);
1936:     for (p = 0; p < *nFV; ++p) {
1937:       DMLabelSetValue(subpointMap, cone[p], 0);
1938:     }
1939:     /* Negative face */
1940:     DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
1941:     /* Not true in parallel
1942:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
1943:     for (p = 0; p < numCells; ++p) {
1944:       DMLabelSetValue(subpointMap, cells[p], 2);
1945:       (*subCells)[subc++] = cells[p];
1946:     }
1947:     DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
1948:     /* Positive face is not included */
1949:   }
1950:   return(0);
1951: }

1955: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DMLabel subpointMap, DM subdm)
1956: {
1957:   DMLabel        label = NULL;
1958:   PetscInt      *pStart, *pEnd;
1959:   PetscInt       dim, cMax, cEnd, c, d;

1963:   if (labelname) {DMPlexGetLabel(dm, labelname, &label);}
1964:   DMGetDimension(dm, &dim);
1965:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
1966:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1967:   if (cMax < 0) return(0);
1968:   PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
1969:   for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1970:   for (c = cMax; c < cEnd; ++c) {
1971:     const PetscInt *cone;
1972:     PetscInt       *closure = NULL;
1973:     PetscInt        fconeSize, coneSize, closureSize, cl, val;

1975:     if (label) {
1976:       DMLabelGetValue(label, c, &val);
1977:       if (val != value) continue;
1978:     }
1979:     DMPlexGetConeSize(dm, c, &coneSize);
1980:     DMPlexGetCone(dm, c, &cone);
1981:     DMPlexGetConeSize(dm, cone[0], &fconeSize);
1982:     if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
1983:     /* Negative face */
1984:     DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1985:     for (cl = 0; cl < closureSize*2; cl += 2) {
1986:       const PetscInt point = closure[cl];

1988:       for (d = 0; d <= dim; ++d) {
1989:         if ((point >= pStart[d]) && (point < pEnd[d])) {
1990:           DMLabelSetValue(subpointMap, point, d);
1991:           break;
1992:         }
1993:       }
1994:     }
1995:     DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
1996:     /* Cells -- positive face is not included */
1997:     for (cl = 0; cl < 1; ++cl) {
1998:       const PetscInt *support;
1999:       PetscInt        supportSize, s;

2001:       DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2002:       /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2003:       DMPlexGetSupport(dm, cone[cl], &support);
2004:       for (s = 0; s < supportSize; ++s) {
2005:         DMLabelSetValue(subpointMap, support[s], dim);
2006:       }
2007:     }
2008:   }
2009:   PetscFree2(pStart, pEnd);
2010:   return(0);
2011: }

2015: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2016: {
2017:   MPI_Comm       comm;
2018:   PetscBool      posOrient = PETSC_FALSE;
2019:   const PetscInt debug     = 0;
2020:   PetscInt       cellDim, faceSize, f;

2024:   PetscObjectGetComm((PetscObject)dm,&comm);
2025:   DMGetDimension(dm, &cellDim);
2026:   if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}

2028:   if (cellDim == 1 && numCorners == 2) {
2029:     /* Triangle */
2030:     faceSize  = numCorners-1;
2031:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2032:   } else if (cellDim == 2 && numCorners == 3) {
2033:     /* Triangle */
2034:     faceSize  = numCorners-1;
2035:     posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2036:   } else if (cellDim == 3 && numCorners == 4) {
2037:     /* Tetrahedron */
2038:     faceSize  = numCorners-1;
2039:     posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2040:   } else if (cellDim == 1 && numCorners == 3) {
2041:     /* Quadratic line */
2042:     faceSize  = 1;
2043:     posOrient = PETSC_TRUE;
2044:   } else if (cellDim == 2 && numCorners == 4) {
2045:     /* Quads */
2046:     faceSize = 2;
2047:     if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2048:       posOrient = PETSC_TRUE;
2049:     } else if ((indices[0] == 3) && (indices[1] == 0)) {
2050:       posOrient = PETSC_TRUE;
2051:     } else {
2052:       if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2053:         posOrient = PETSC_FALSE;
2054:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2055:     }
2056:   } else if (cellDim == 2 && numCorners == 6) {
2057:     /* Quadratic triangle (I hate this) */
2058:     /* Edges are determined by the first 2 vertices (corners of edges) */
2059:     const PetscInt faceSizeTri = 3;
2060:     PetscInt       sortedIndices[3], i, iFace;
2061:     PetscBool      found                    = PETSC_FALSE;
2062:     PetscInt       faceVerticesTriSorted[9] = {
2063:       0, 3,  4, /* bottom */
2064:       1, 4,  5, /* right */
2065:       2, 3,  5, /* left */
2066:     };
2067:     PetscInt       faceVerticesTri[9] = {
2068:       0, 3,  4, /* bottom */
2069:       1, 4,  5, /* right */
2070:       2, 5,  3, /* left */
2071:     };

2073:     faceSize = faceSizeTri;
2074:     for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2075:     PetscSortInt(faceSizeTri, sortedIndices);
2076:     for (iFace = 0; iFace < 3; ++iFace) {
2077:       const PetscInt ii = iFace*faceSizeTri;
2078:       PetscInt       fVertex, cVertex;

2080:       if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2081:           (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2082:         for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2083:           for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2084:             if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2085:               faceVertices[fVertex] = origVertices[cVertex];
2086:               break;
2087:             }
2088:           }
2089:         }
2090:         found = PETSC_TRUE;
2091:         break;
2092:       }
2093:     }
2094:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2095:     if (posOriented) *posOriented = PETSC_TRUE;
2096:     return(0);
2097:   } else if (cellDim == 2 && numCorners == 9) {
2098:     /* Quadratic quad (I hate this) */
2099:     /* Edges are determined by the first 2 vertices (corners of edges) */
2100:     const PetscInt faceSizeQuad = 3;
2101:     PetscInt       sortedIndices[3], i, iFace;
2102:     PetscBool      found                      = PETSC_FALSE;
2103:     PetscInt       faceVerticesQuadSorted[12] = {
2104:       0, 1,  4, /* bottom */
2105:       1, 2,  5, /* right */
2106:       2, 3,  6, /* top */
2107:       0, 3,  7, /* left */
2108:     };
2109:     PetscInt       faceVerticesQuad[12] = {
2110:       0, 1,  4, /* bottom */
2111:       1, 2,  5, /* right */
2112:       2, 3,  6, /* top */
2113:       3, 0,  7, /* left */
2114:     };

2116:     faceSize = faceSizeQuad;
2117:     for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2118:     PetscSortInt(faceSizeQuad, sortedIndices);
2119:     for (iFace = 0; iFace < 4; ++iFace) {
2120:       const PetscInt ii = iFace*faceSizeQuad;
2121:       PetscInt       fVertex, cVertex;

2123:       if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2124:           (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2125:         for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2126:           for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2127:             if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2128:               faceVertices[fVertex] = origVertices[cVertex];
2129:               break;
2130:             }
2131:           }
2132:         }
2133:         found = PETSC_TRUE;
2134:         break;
2135:       }
2136:     }
2137:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2138:     if (posOriented) *posOriented = PETSC_TRUE;
2139:     return(0);
2140:   } else if (cellDim == 3 && numCorners == 8) {
2141:     /* Hexes
2142:        A hex is two oriented quads with the normal of the first
2143:        pointing up at the second.

2145:           7---6
2146:          /|  /|
2147:         4---5 |
2148:         | 1-|-2
2149:         |/  |/
2150:         0---3

2152:         Faces are determined by the first 4 vertices (corners of faces) */
2153:     const PetscInt faceSizeHex = 4;
2154:     PetscInt       sortedIndices[4], i, iFace;
2155:     PetscBool      found                     = PETSC_FALSE;
2156:     PetscInt       faceVerticesHexSorted[24] = {
2157:       0, 1, 2, 3,  /* bottom */
2158:       4, 5, 6, 7,  /* top */
2159:       0, 3, 4, 5,  /* front */
2160:       2, 3, 5, 6,  /* right */
2161:       1, 2, 6, 7,  /* back */
2162:       0, 1, 4, 7,  /* left */
2163:     };
2164:     PetscInt       faceVerticesHex[24] = {
2165:       1, 2, 3, 0,  /* bottom */
2166:       4, 5, 6, 7,  /* top */
2167:       0, 3, 5, 4,  /* front */
2168:       3, 2, 6, 5,  /* right */
2169:       2, 1, 7, 6,  /* back */
2170:       1, 0, 4, 7,  /* left */
2171:     };

2173:     faceSize = faceSizeHex;
2174:     for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2175:     PetscSortInt(faceSizeHex, sortedIndices);
2176:     for (iFace = 0; iFace < 6; ++iFace) {
2177:       const PetscInt ii = iFace*faceSizeHex;
2178:       PetscInt       fVertex, cVertex;

2180:       if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2181:           (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2182:           (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2183:           (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2184:         for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2185:           for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2186:             if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2187:               faceVertices[fVertex] = origVertices[cVertex];
2188:               break;
2189:             }
2190:           }
2191:         }
2192:         found = PETSC_TRUE;
2193:         break;
2194:       }
2195:     }
2196:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2197:     if (posOriented) *posOriented = PETSC_TRUE;
2198:     return(0);
2199:   } else if (cellDim == 3 && numCorners == 10) {
2200:     /* Quadratic tet */
2201:     /* Faces are determined by the first 3 vertices (corners of faces) */
2202:     const PetscInt faceSizeTet = 6;
2203:     PetscInt       sortedIndices[6], i, iFace;
2204:     PetscBool      found                     = PETSC_FALSE;
2205:     PetscInt       faceVerticesTetSorted[24] = {
2206:       0, 1, 2,  6, 7, 8, /* bottom */
2207:       0, 3, 4,  6, 7, 9,  /* front */
2208:       1, 4, 5,  7, 8, 9,  /* right */
2209:       2, 3, 5,  6, 8, 9,  /* left */
2210:     };
2211:     PetscInt       faceVerticesTet[24] = {
2212:       0, 1, 2,  6, 7, 8, /* bottom */
2213:       0, 4, 3,  6, 7, 9,  /* front */
2214:       1, 5, 4,  7, 8, 9,  /* right */
2215:       2, 3, 5,  8, 6, 9,  /* left */
2216:     };

2218:     faceSize = faceSizeTet;
2219:     for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2220:     PetscSortInt(faceSizeTet, sortedIndices);
2221:     for (iFace=0; iFace < 4; ++iFace) {
2222:       const PetscInt ii = iFace*faceSizeTet;
2223:       PetscInt       fVertex, cVertex;

2225:       if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2226:           (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2227:           (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2228:           (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2229:         for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2230:           for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2231:             if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2232:               faceVertices[fVertex] = origVertices[cVertex];
2233:               break;
2234:             }
2235:           }
2236:         }
2237:         found = PETSC_TRUE;
2238:         break;
2239:       }
2240:     }
2241:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2242:     if (posOriented) *posOriented = PETSC_TRUE;
2243:     return(0);
2244:   } else if (cellDim == 3 && numCorners == 27) {
2245:     /* Quadratic hexes (I hate this)
2246:        A hex is two oriented quads with the normal of the first
2247:        pointing up at the second.

2249:          7---6
2250:         /|  /|
2251:        4---5 |
2252:        | 3-|-2
2253:        |/  |/
2254:        0---1

2256:        Faces are determined by the first 4 vertices (corners of faces) */
2257:     const PetscInt faceSizeQuadHex = 9;
2258:     PetscInt       sortedIndices[9], i, iFace;
2259:     PetscBool      found                         = PETSC_FALSE;
2260:     PetscInt       faceVerticesQuadHexSorted[54] = {
2261:       0, 1, 2, 3,  8, 9, 10, 11,  24, /* bottom */
2262:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2263:       0, 1, 4, 5,  8, 12, 16, 17,  22, /* front */
2264:       1, 2, 5, 6,  9, 13, 17, 18,  21, /* right */
2265:       2, 3, 6, 7,  10, 14, 18, 19,  23, /* back */
2266:       0, 3, 4, 7,  11, 15, 16, 19,  20, /* left */
2267:     };
2268:     PetscInt       faceVerticesQuadHex[54] = {
2269:       3, 2, 1, 0,  10, 9, 8, 11,  24, /* bottom */
2270:       4, 5, 6, 7,  12, 13, 14, 15,  25, /* top */
2271:       0, 1, 5, 4,  8, 17, 12, 16,  22, /* front */
2272:       1, 2, 6, 5,  9, 18, 13, 17,  21, /* right */
2273:       2, 3, 7, 6,  10, 19, 14, 18,  23, /* back */
2274:       3, 0, 4, 7,  11, 16, 15, 19,  20 /* left */
2275:     };

2277:     faceSize = faceSizeQuadHex;
2278:     for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2279:     PetscSortInt(faceSizeQuadHex, sortedIndices);
2280:     for (iFace = 0; iFace < 6; ++iFace) {
2281:       const PetscInt ii = iFace*faceSizeQuadHex;
2282:       PetscInt       fVertex, cVertex;

2284:       if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2285:           (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2286:           (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2287:           (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2288:         for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2289:           for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2290:             if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2291:               faceVertices[fVertex] = origVertices[cVertex];
2292:               break;
2293:             }
2294:           }
2295:         }
2296:         found = PETSC_TRUE;
2297:         break;
2298:       }
2299:     }
2300:     if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2301:     if (posOriented) *posOriented = PETSC_TRUE;
2302:     return(0);
2303:   } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2304:   if (!posOrient) {
2305:     if (debug) {PetscPrintf(comm, "  Reversing initial face orientation\n");}
2306:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2307:   } else {
2308:     if (debug) {PetscPrintf(comm, "  Keeping initial face orientation\n");}
2309:     for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2310:   }
2311:   if (posOriented) *posOriented = posOrient;
2312:   return(0);
2313: }

2317: /*
2318:     Given a cell and a face, as a set of vertices,
2319:       return the oriented face, as a set of vertices, in faceVertices
2320:     The orientation is such that the face normal points out of the cell
2321: */
2322: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2323: {
2324:   const PetscInt *cone = NULL;
2325:   PetscInt        coneSize, v, f, v2;
2326:   PetscInt        oppositeVertex = -1;
2327:   PetscErrorCode  ierr;

2330:   DMPlexGetConeSize(dm, cell, &coneSize);
2331:   DMPlexGetCone(dm, cell, &cone);
2332:   for (v = 0, v2 = 0; v < coneSize; ++v) {
2333:     PetscBool found = PETSC_FALSE;

2335:     for (f = 0; f < faceSize; ++f) {
2336:       if (face[f] == cone[v]) {
2337:         found = PETSC_TRUE; break;
2338:       }
2339:     }
2340:     if (found) {
2341:       indices[v2]      = v;
2342:       origVertices[v2] = cone[v];
2343:       ++v2;
2344:     } else {
2345:       oppositeVertex = v;
2346:     }
2347:   }
2348:   DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2349:   return(0);
2350: }

2354: /*
2355:   DMPlexInsertFace_Internal - Puts a face into the mesh

2357:   Not collective

2359:   Input Parameters:
2360:   + dm              - The DMPlex
2361:   . numFaceVertex   - The number of vertices in the face
2362:   . faceVertices    - The vertices in the face for dm
2363:   . subfaceVertices - The vertices in the face for subdm
2364:   . numCorners      - The number of vertices in the cell
2365:   . cell            - A cell in dm containing the face
2366:   . subcell         - A cell in subdm containing the face
2367:   . firstFace       - First face in the mesh
2368:   - newFacePoint    - Next face in the mesh

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

2373:   Level: developer
2374: */
2375: 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)
2376: {
2377:   MPI_Comm        comm;
2378:   DM_Plex        *submesh = (DM_Plex*) subdm->data;
2379:   const PetscInt *faces;
2380:   PetscInt        numFaces, coneSize;
2381:   PetscErrorCode  ierr;

2384:   PetscObjectGetComm((PetscObject)dm,&comm);
2385:   DMPlexGetConeSize(subdm, subcell, &coneSize);
2386:   if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2387: #if 0
2388:   /* Cannot use this because support() has not been constructed yet */
2389:   DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2390: #else
2391:   {
2392:     PetscInt f;

2394:     numFaces = 0;
2395:     DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2396:     for (f = firstFace; f < *newFacePoint; ++f) {
2397:       PetscInt dof, off, d;

2399:       PetscSectionGetDof(submesh->coneSection, f, &dof);
2400:       PetscSectionGetOffset(submesh->coneSection, f, &off);
2401:       /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2402:       for (d = 0; d < dof; ++d) {
2403:         const PetscInt p = submesh->cones[off+d];
2404:         PetscInt       v;

2406:         for (v = 0; v < numFaceVertices; ++v) {
2407:           if (subfaceVertices[v] == p) break;
2408:         }
2409:         if (v == numFaceVertices) break;
2410:       }
2411:       if (d == dof) {
2412:         numFaces               = 1;
2413:         ((PetscInt*) faces)[0] = f;
2414:       }
2415:     }
2416:   }
2417: #endif
2418:   if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2419:   else if (numFaces == 1) {
2420:     /* Add the other cell neighbor for this face */
2421:     DMPlexSetCone(subdm, subcell, faces);
2422:   } else {
2423:     PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2424:     PetscBool posOriented;

2426:     DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2427:     origVertices        = &orientedVertices[numFaceVertices];
2428:     indices             = &orientedVertices[numFaceVertices*2];
2429:     orientedSubVertices = &orientedVertices[numFaceVertices*3];
2430:     DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2431:     /* TODO: I know that routine should return a permutation, not the indices */
2432:     for (v = 0; v < numFaceVertices; ++v) {
2433:       const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2434:       for (ov = 0; ov < numFaceVertices; ++ov) {
2435:         if (orientedVertices[ov] == vertex) {
2436:           orientedSubVertices[ov] = subvertex;
2437:           break;
2438:         }
2439:       }
2440:       if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2441:     }
2442:     DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2443:     DMPlexSetCone(subdm, subcell, newFacePoint);
2444:     DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2445:     ++(*newFacePoint);
2446:   }
2447: #if 0
2448:   DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2449: #else
2450:   DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2451: #endif
2452:   return(0);
2453: }

2457: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2458: {
2459:   MPI_Comm        comm;
2460:   DMLabel         subpointMap;
2461:   IS              subvertexIS,  subcellIS;
2462:   const PetscInt *subVertices, *subCells;
2463:   PetscInt        numSubVertices, firstSubVertex, numSubCells;
2464:   PetscInt       *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2465:   PetscInt        vStart, vEnd, c, f;
2466:   PetscErrorCode  ierr;

2469:   PetscObjectGetComm((PetscObject)dm,&comm);
2470:   /* Create subpointMap which marks the submesh */
2471:   DMLabelCreate("subpoint_map", &subpointMap);
2472:   DMPlexSetSubpointMap(subdm, subpointMap);
2473:   DMLabelDestroy(&subpointMap);
2474:   if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2475:   /* Setup chart */
2476:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2477:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2478:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2479:   DMPlexSetVTKCellHeight(subdm, 1);
2480:   /* Set cone sizes */
2481:   firstSubVertex = numSubCells;
2482:   firstSubFace   = numSubCells+numSubVertices;
2483:   newFacePoint   = firstSubFace;
2484:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2485:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2486:   DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2487:   if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2488:   for (c = 0; c < numSubCells; ++c) {
2489:     DMPlexSetConeSize(subdm, c, 1);
2490:   }
2491:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2492:     DMPlexSetConeSize(subdm, f, nFV);
2493:   }
2494:   DMSetUp(subdm);
2495:   /* Create face cones */
2496:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2497:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2498:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2499:   for (c = 0; c < numSubCells; ++c) {
2500:     const PetscInt cell    = subCells[c];
2501:     const PetscInt subcell = c;
2502:     PetscInt      *closure = NULL;
2503:     PetscInt       closureSize, cl, numCorners = 0, faceSize = 0;

2505:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2506:     for (cl = 0; cl < closureSize*2; cl += 2) {
2507:       const PetscInt point = closure[cl];
2508:       PetscInt       subVertex;

2510:       if ((point >= vStart) && (point < vEnd)) {
2511:         ++numCorners;
2512:         PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2513:         if (subVertex >= 0) {
2514:           closure[faceSize] = point;
2515:           subface[faceSize] = firstSubVertex+subVertex;
2516:           ++faceSize;
2517:         }
2518:       }
2519:     }
2520:     if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2521:     if (faceSize == nFV) {
2522:       DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2523:     }
2524:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2525:   }
2526:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2527:   DMPlexSymmetrize(subdm);
2528:   DMPlexStratify(subdm);
2529:   /* Build coordinates */
2530:   {
2531:     PetscSection coordSection, subCoordSection;
2532:     Vec          coordinates, subCoordinates;
2533:     PetscScalar *coords, *subCoords;
2534:     PetscInt     numComp, coordSize, v;

2536:     DMGetCoordinateSection(dm, &coordSection);
2537:     DMGetCoordinatesLocal(dm, &coordinates);
2538:     DMGetCoordinateSection(subdm, &subCoordSection);
2539:     PetscSectionSetNumFields(subCoordSection, 1);
2540:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2541:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2542:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2543:     for (v = 0; v < numSubVertices; ++v) {
2544:       const PetscInt vertex    = subVertices[v];
2545:       const PetscInt subvertex = firstSubVertex+v;
2546:       PetscInt       dof;

2548:       PetscSectionGetDof(coordSection, vertex, &dof);
2549:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2550:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2551:     }
2552:     PetscSectionSetUp(subCoordSection);
2553:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2554:     VecCreate(comm, &subCoordinates);
2555:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2556:     VecSetType(subCoordinates,VECSTANDARD);
2557:     if (coordSize) {
2558:       VecGetArray(coordinates,    &coords);
2559:       VecGetArray(subCoordinates, &subCoords);
2560:       for (v = 0; v < numSubVertices; ++v) {
2561:         const PetscInt vertex    = subVertices[v];
2562:         const PetscInt subvertex = firstSubVertex+v;
2563:         PetscInt       dof, off, sdof, soff, d;

2565:         PetscSectionGetDof(coordSection, vertex, &dof);
2566:         PetscSectionGetOffset(coordSection, vertex, &off);
2567:         PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2568:         PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2569:         if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2570:         for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2571:       }
2572:       VecRestoreArray(coordinates,    &coords);
2573:       VecRestoreArray(subCoordinates, &subCoords);
2574:     }
2575:     DMSetCoordinatesLocal(subdm, subCoordinates);
2576:     VecDestroy(&subCoordinates);
2577:   }
2578:   /* Cleanup */
2579:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2580:   ISDestroy(&subvertexIS);
2581:   if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2582:   ISDestroy(&subcellIS);
2583:   return(0);
2584: }

2588: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2589: {
2590:   MPI_Comm         comm;
2591:   DMLabel          subpointMap;
2592:   IS              *subpointIS;
2593:   const PetscInt **subpoints;
2594:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *coneONew;
2595:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2596:   PetscErrorCode   ierr;

2599:   PetscObjectGetComm((PetscObject)dm,&comm);
2600:   /* Create subpointMap which marks the submesh */
2601:   DMLabelCreate("subpoint_map", &subpointMap);
2602:   DMPlexSetSubpointMap(subdm, subpointMap);
2603:   DMLabelDestroy(&subpointMap);
2604:   if (vertexLabel) {DMPlexMarkSubmesh_Interpolated(dm, vertexLabel, value, subpointMap, subdm);}
2605:   /* Setup chart */
2606:   DMGetDimension(dm, &dim);
2607:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2608:   for (d = 0; d <= dim; ++d) {
2609:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2610:     totSubPoints += numSubPoints[d];
2611:   }
2612:   DMPlexSetChart(subdm, 0, totSubPoints);
2613:   DMPlexSetVTKCellHeight(subdm, 1);
2614:   /* Set cone sizes */
2615:   firstSubPoint[dim] = 0;
2616:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
2617:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
2618:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2619:   for (d = 0; d <= dim; ++d) {
2620:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2621:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2622:   }
2623:   for (d = 0; d <= dim; ++d) {
2624:     for (p = 0; p < numSubPoints[d]; ++p) {
2625:       const PetscInt  point    = subpoints[d][p];
2626:       const PetscInt  subpoint = firstSubPoint[d] + p;
2627:       const PetscInt *cone;
2628:       PetscInt        coneSize, coneSizeNew, c, val;

2630:       DMPlexGetConeSize(dm, point, &coneSize);
2631:       DMPlexSetConeSize(subdm, subpoint, coneSize);
2632:       if (d == dim) {
2633:         DMPlexGetCone(dm, point, &cone);
2634:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2635:           DMLabelGetValue(subpointMap, cone[c], &val);
2636:           if (val >= 0) coneSizeNew++;
2637:         }
2638:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2639:       }
2640:     }
2641:   }
2642:   DMSetUp(subdm);
2643:   /* Set cones */
2644:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2645:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&coneONew);
2646:   for (d = 0; d <= dim; ++d) {
2647:     for (p = 0; p < numSubPoints[d]; ++p) {
2648:       const PetscInt  point    = subpoints[d][p];
2649:       const PetscInt  subpoint = firstSubPoint[d] + p;
2650:       const PetscInt *cone, *ornt;
2651:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

2653:       DMPlexGetConeSize(dm, point, &coneSize);
2654:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2655:       DMPlexGetCone(dm, point, &cone);
2656:       DMPlexGetConeOrientation(dm, point, &ornt);
2657:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2658:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2659:         if (subc >= 0) {
2660:           coneNew[coneSizeNew]  = firstSubPoint[d-1] + subc;
2661:           coneONew[coneSizeNew] = ornt[c];
2662:           ++coneSizeNew;
2663:         }
2664:       }
2665:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2666:       DMPlexSetCone(subdm, subpoint, coneNew);
2667:       DMPlexSetConeOrientation(subdm, subpoint, coneONew);
2668:     }
2669:   }
2670:   PetscFree2(coneNew,coneONew);
2671:   DMPlexSymmetrize(subdm);
2672:   DMPlexStratify(subdm);
2673:   /* Build coordinates */
2674:   {
2675:     PetscSection coordSection, subCoordSection;
2676:     Vec          coordinates, subCoordinates;
2677:     PetscScalar *coords, *subCoords;
2678:     PetscInt     numComp, coordSize;

2680:     DMGetCoordinateSection(dm, &coordSection);
2681:     DMGetCoordinatesLocal(dm, &coordinates);
2682:     DMGetCoordinateSection(subdm, &subCoordSection);
2683:     PetscSectionSetNumFields(subCoordSection, 1);
2684:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2685:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2686:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2687:     for (v = 0; v < numSubPoints[0]; ++v) {
2688:       const PetscInt vertex    = subpoints[0][v];
2689:       const PetscInt subvertex = firstSubPoint[0]+v;
2690:       PetscInt       dof;

2692:       PetscSectionGetDof(coordSection, vertex, &dof);
2693:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2694:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2695:     }
2696:     PetscSectionSetUp(subCoordSection);
2697:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2698:     VecCreate(comm, &subCoordinates);
2699:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2700:     VecSetType(subCoordinates,VECSTANDARD);
2701:     VecGetArray(coordinates,    &coords);
2702:     VecGetArray(subCoordinates, &subCoords);
2703:     for (v = 0; v < numSubPoints[0]; ++v) {
2704:       const PetscInt vertex    = subpoints[0][v];
2705:       const PetscInt subvertex = firstSubPoint[0]+v;
2706:       PetscInt dof, off, sdof, soff, d;

2708:       PetscSectionGetDof(coordSection, vertex, &dof);
2709:       PetscSectionGetOffset(coordSection, vertex, &off);
2710:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2711:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2712:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2713:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2714:     }
2715:     VecRestoreArray(coordinates,    &coords);
2716:     VecRestoreArray(subCoordinates, &subCoords);
2717:     DMSetCoordinatesLocal(subdm, subCoordinates);
2718:     VecDestroy(&subCoordinates);
2719:   }
2720:   /* Cleanup */
2721:   for (d = 0; d <= dim; ++d) {
2722:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
2723:     ISDestroy(&subpointIS[d]);
2724:   }
2725:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
2726:   return(0);
2727: }

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

2734:   Input Parameters:
2735: + dm           - The original mesh
2736: . vertexLabel  - The DMLabel marking vertices contained in the surface
2737: - value        - The label value to use

2739:   Output Parameter:
2740: . subdm - The surface mesh

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

2744:   Level: developer

2746: .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue()
2747: @*/
2748: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
2749: {
2750:   PetscInt       dim, depth;

2756:   DMGetDimension(dm, &dim);
2757:   DMPlexGetDepth(dm, &depth);
2758:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
2759:   DMSetType(*subdm, DMPLEX);
2760:   DMSetDimension(*subdm, dim-1);
2761:   if (depth == dim) {
2762:     DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
2763:   } else {
2764:     DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
2765:   }
2766:   return(0);
2767: }

2771: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2772: {
2773:   PetscInt       subPoint;

2776:   PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2777:   return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2778: }

2782: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
2783: {
2784:   MPI_Comm        comm;
2785:   DMLabel         subpointMap;
2786:   IS              subvertexIS;
2787:   const PetscInt *subVertices;
2788:   PetscInt        numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
2789:   PetscInt       *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
2790:   PetscInt        cMax, c, f;
2791:   PetscErrorCode  ierr;

2794:   PetscObjectGetComm((PetscObject)dm, &comm);
2795:   /* Create subpointMap which marks the submesh */
2796:   DMLabelCreate("subpoint_map", &subpointMap);
2797:   DMPlexSetSubpointMap(subdm, subpointMap);
2798:   DMLabelDestroy(&subpointMap);
2799:   DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
2800:   /* Setup chart */
2801:   DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2802:   DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2803:   DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2804:   DMPlexSetVTKCellHeight(subdm, 1);
2805:   /* Set cone sizes */
2806:   firstSubVertex = numSubCells;
2807:   firstSubFace   = numSubCells+numSubVertices;
2808:   newFacePoint   = firstSubFace;
2809:   DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2810:   if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2811:   for (c = 0; c < numSubCells; ++c) {
2812:     DMPlexSetConeSize(subdm, c, 1);
2813:   }
2814:   for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2815:     DMPlexSetConeSize(subdm, f, nFV);
2816:   }
2817:   DMSetUp(subdm);
2818:   /* Create face cones */
2819:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2820:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2821:   DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2822:   for (c = 0; c < numSubCells; ++c) {
2823:     const PetscInt  cell    = subCells[c];
2824:     const PetscInt  subcell = c;
2825:     const PetscInt *cone, *cells;
2826:     PetscInt        numCells, subVertex, p, v;

2828:     if (cell < cMax) continue;
2829:     DMPlexGetCone(dm, cell, &cone);
2830:     for (v = 0; v < nFV; ++v) {
2831:       PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
2832:       subface[v] = firstSubVertex+subVertex;
2833:     }
2834:     DMPlexSetCone(subdm, newFacePoint, subface);
2835:     DMPlexSetCone(subdm, subcell, &newFacePoint);
2836:     DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
2837:     /* Not true in parallel
2838:     if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2839:     for (p = 0; p < numCells; ++p) {
2840:       PetscInt negsubcell;

2842:       if (cells[p] >= cMax) continue;
2843:       /* I know this is a crap search */
2844:       for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
2845:         if (subCells[negsubcell] == cells[p]) break;
2846:       }
2847:       if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
2848:       DMPlexSetCone(subdm, negsubcell, &newFacePoint);
2849:     }
2850:     DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
2851:     ++newFacePoint;
2852:   }
2853:   DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2854:   DMPlexSymmetrize(subdm);
2855:   DMPlexStratify(subdm);
2856:   /* Build coordinates */
2857:   {
2858:     PetscSection coordSection, subCoordSection;
2859:     Vec          coordinates, subCoordinates;
2860:     PetscScalar *coords, *subCoords;
2861:     PetscInt     numComp, coordSize, v;

2863:     DMGetCoordinateSection(dm, &coordSection);
2864:     DMGetCoordinatesLocal(dm, &coordinates);
2865:     DMGetCoordinateSection(subdm, &subCoordSection);
2866:     PetscSectionSetNumFields(subCoordSection, 1);
2867:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2868:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2869:     PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2870:     for (v = 0; v < numSubVertices; ++v) {
2871:       const PetscInt vertex    = subVertices[v];
2872:       const PetscInt subvertex = firstSubVertex+v;
2873:       PetscInt       dof;

2875:       PetscSectionGetDof(coordSection, vertex, &dof);
2876:       PetscSectionSetDof(subCoordSection, subvertex, dof);
2877:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2878:     }
2879:     PetscSectionSetUp(subCoordSection);
2880:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
2881:     VecCreate(comm, &subCoordinates);
2882:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2883:     VecSetType(subCoordinates,VECSTANDARD);
2884:     VecGetArray(coordinates,    &coords);
2885:     VecGetArray(subCoordinates, &subCoords);
2886:     for (v = 0; v < numSubVertices; ++v) {
2887:       const PetscInt vertex    = subVertices[v];
2888:       const PetscInt subvertex = firstSubVertex+v;
2889:       PetscInt       dof, off, sdof, soff, d;

2891:       PetscSectionGetDof(coordSection, vertex, &dof);
2892:       PetscSectionGetOffset(coordSection, vertex, &off);
2893:       PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2894:       PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2895:       if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2896:       for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2897:     }
2898:     VecRestoreArray(coordinates,    &coords);
2899:     VecRestoreArray(subCoordinates, &subCoords);
2900:     DMSetCoordinatesLocal(subdm, subCoordinates);
2901:     VecDestroy(&subCoordinates);
2902:   }
2903:   /* Build SF */
2904:   CHKMEMQ;
2905:   {
2906:     PetscSF            sfPoint, sfPointSub;
2907:     const PetscSFNode *remotePoints;
2908:     PetscSFNode       *sremotePoints, *newLocalPoints, *newOwners;
2909:     const PetscInt    *localPoints;
2910:     PetscInt          *slocalPoints;
2911:     PetscInt           numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
2912:     PetscMPIInt        rank;

2914:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2915:     DMGetPointSF(dm, &sfPoint);
2916:     DMGetPointSF(subdm, &sfPointSub);
2917:     DMPlexGetChart(dm, &pStart, &pEnd);
2918:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2919:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2920:     if (numRoots >= 0) {
2921:       /* Only vertices should be shared */
2922:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2923:       for (p = 0; p < pEnd-pStart; ++p) {
2924:         newLocalPoints[p].rank  = -2;
2925:         newLocalPoints[p].index = -2;
2926:       }
2927:       /* Set subleaves */
2928:       for (l = 0; l < numLeaves; ++l) {
2929:         const PetscInt point    = localPoints[l];
2930:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2932:         if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
2933:         if (subPoint < 0) continue;
2934:         newLocalPoints[point-pStart].rank  = rank;
2935:         newLocalPoints[point-pStart].index = subPoint;
2936:         ++numSubLeaves;
2937:       }
2938:       /* Must put in owned subpoints */
2939:       for (p = pStart; p < pEnd; ++p) {
2940:         const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);

2942:         if (subPoint < 0) {
2943:           newOwners[p-pStart].rank  = -3;
2944:           newOwners[p-pStart].index = -3;
2945:         } else {
2946:           newOwners[p-pStart].rank  = rank;
2947:           newOwners[p-pStart].index = subPoint;
2948:         }
2949:       }
2950:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2951:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
2952:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2953:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
2954:       PetscMalloc1(numSubLeaves,    &slocalPoints);
2955:       PetscMalloc1(numSubLeaves, &sremotePoints);
2956:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
2957:         const PetscInt point    = localPoints[l];
2958:         const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);

2960:         if (subPoint < 0) continue;
2961:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
2962:         slocalPoints[sl]        = subPoint;
2963:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
2964:         sremotePoints[sl].index = newLocalPoints[point].index;
2965:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
2966:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
2967:         ++sl;
2968:       }
2969:       PetscFree2(newLocalPoints,newOwners);
2970:       if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
2971:       PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
2972:     }
2973:   }
2974:   CHKMEMQ;
2975:   /* Cleanup */
2976:   if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2977:   ISDestroy(&subvertexIS);
2978:   PetscFree(subCells);
2979:   return(0);
2980: }

2984: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char label[], PetscInt value, DM subdm)
2985: {
2986:   MPI_Comm         comm;
2987:   DMLabel          subpointMap;
2988:   IS              *subpointIS;
2989:   const PetscInt **subpoints;
2990:   PetscInt        *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2991:   PetscInt         totSubPoints = 0, maxConeSize, dim, p, d, v;
2992:   PetscErrorCode   ierr;

2995:   PetscObjectGetComm((PetscObject)dm,&comm);
2996:   /* Create subpointMap which marks the submesh */
2997:   DMLabelCreate("subpoint_map", &subpointMap);
2998:   DMPlexSetSubpointMap(subdm, subpointMap);
2999:   DMLabelDestroy(&subpointMap);
3000:   DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
3001:   /* Setup chart */
3002:   DMGetDimension(dm, &dim);
3003:   PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
3004:   for (d = 0; d <= dim; ++d) {
3005:     DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3006:     totSubPoints += numSubPoints[d];
3007:   }
3008:   DMPlexSetChart(subdm, 0, totSubPoints);
3009:   DMPlexSetVTKCellHeight(subdm, 1);
3010:   /* Set cone sizes */
3011:   firstSubPoint[dim] = 0;
3012:   firstSubPoint[0]   = firstSubPoint[dim] + numSubPoints[dim];
3013:   if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0]     + numSubPoints[0];}
3014:   if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
3015:   for (d = 0; d <= dim; ++d) {
3016:     DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3017:     if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
3018:   }
3019:   for (d = 0; d <= dim; ++d) {
3020:     for (p = 0; p < numSubPoints[d]; ++p) {
3021:       const PetscInt  point    = subpoints[d][p];
3022:       const PetscInt  subpoint = firstSubPoint[d] + p;
3023:       const PetscInt *cone;
3024:       PetscInt        coneSize, coneSizeNew, c, val;

3026:       DMPlexGetConeSize(dm, point, &coneSize);
3027:       DMPlexSetConeSize(subdm, subpoint, coneSize);
3028:       if (d == dim) {
3029:         DMPlexGetCone(dm, point, &cone);
3030:         for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3031:           DMLabelGetValue(subpointMap, cone[c], &val);
3032:           if (val >= 0) coneSizeNew++;
3033:         }
3034:         DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3035:       }
3036:     }
3037:   }
3038:   DMSetUp(subdm);
3039:   /* Set cones */
3040:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3041:   PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3042:   for (d = 0; d <= dim; ++d) {
3043:     for (p = 0; p < numSubPoints[d]; ++p) {
3044:       const PetscInt  point    = subpoints[d][p];
3045:       const PetscInt  subpoint = firstSubPoint[d] + p;
3046:       const PetscInt *cone, *ornt;
3047:       PetscInt        coneSize, subconeSize, coneSizeNew, c, subc;

3049:       DMPlexGetConeSize(dm, point, &coneSize);
3050:       DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3051:       DMPlexGetCone(dm, point, &cone);
3052:       DMPlexGetConeOrientation(dm, point, &ornt);
3053:       for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3054:         PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3055:         if (subc >= 0) {
3056:           coneNew[coneSizeNew]   = firstSubPoint[d-1] + subc;
3057:           orntNew[coneSizeNew++] = ornt[c];
3058:         }
3059:       }
3060:       if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3061:       DMPlexSetCone(subdm, subpoint, coneNew);
3062:       DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3063:     }
3064:   }
3065:   PetscFree2(coneNew,orntNew);
3066:   DMPlexSymmetrize(subdm);
3067:   DMPlexStratify(subdm);
3068:   /* Build coordinates */
3069:   {
3070:     PetscSection coordSection, subCoordSection;
3071:     Vec          coordinates, subCoordinates;
3072:     PetscScalar *coords, *subCoords;
3073:     PetscInt     numComp, coordSize;

3075:     DMGetCoordinateSection(dm, &coordSection);
3076:     DMGetCoordinatesLocal(dm, &coordinates);
3077:     DMGetCoordinateSection(subdm, &subCoordSection);
3078:     PetscSectionSetNumFields(subCoordSection, 1);
3079:     PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3080:     PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3081:     PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3082:     for (v = 0; v < numSubPoints[0]; ++v) {
3083:       const PetscInt vertex    = subpoints[0][v];
3084:       const PetscInt subvertex = firstSubPoint[0]+v;
3085:       PetscInt       dof;

3087:       PetscSectionGetDof(coordSection, vertex, &dof);
3088:       PetscSectionSetDof(subCoordSection, subvertex, dof);
3089:       PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3090:     }
3091:     PetscSectionSetUp(subCoordSection);
3092:     PetscSectionGetStorageSize(subCoordSection, &coordSize);
3093:     VecCreate(comm, &subCoordinates);
3094:     VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3095:     VecSetType(subCoordinates,VECSTANDARD);
3096:     VecGetArray(coordinates,    &coords);
3097:     VecGetArray(subCoordinates, &subCoords);
3098:     for (v = 0; v < numSubPoints[0]; ++v) {
3099:       const PetscInt vertex    = subpoints[0][v];
3100:       const PetscInt subvertex = firstSubPoint[0]+v;
3101:       PetscInt dof, off, sdof, soff, d;

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

3126:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3127:     DMGetPointSF(dm, &sfPoint);
3128:     DMGetPointSF(subdm, &sfPointSub);
3129:     DMPlexGetChart(dm, &pStart, &pEnd);
3130:     DMPlexGetChart(subdm, NULL, &numSubroots);
3131:     DMPlexCreateSubpointIS(subdm, &subpIS);
3132:     if (subpIS) {
3133:       ISGetIndices(subpIS, &subpoints);
3134:       ISGetLocalSize(subpIS, &numSubpoints);
3135:     }
3136:     PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3137:     if (numRoots >= 0) {
3138:       PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3139:       for (p = 0; p < pEnd-pStart; ++p) {
3140:         newLocalPoints[p].rank  = -2;
3141:         newLocalPoints[p].index = -2;
3142:       }
3143:       /* Set subleaves */
3144:       for (l = 0; l < numLeaves; ++l) {
3145:         const PetscInt point    = localPoints[l];
3146:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3148:         if (subpoint < 0) continue;
3149:         newLocalPoints[point-pStart].rank  = rank;
3150:         newLocalPoints[point-pStart].index = subpoint;
3151:         ++numSubleaves;
3152:       }
3153:       /* Must put in owned subpoints */
3154:       for (p = pStart; p < pEnd; ++p) {
3155:         const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);

3157:         if (subpoint < 0) {
3158:           newOwners[p-pStart].rank  = -3;
3159:           newOwners[p-pStart].index = -3;
3160:         } else {
3161:           newOwners[p-pStart].rank  = rank;
3162:           newOwners[p-pStart].index = subpoint;
3163:         }
3164:       }
3165:       PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3166:       PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3167:       PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3168:       PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3169:       PetscMalloc1(numSubleaves, &slocalPoints);
3170:       PetscMalloc1(numSubleaves, &sremotePoints);
3171:       for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3172:         const PetscInt point    = localPoints[l];
3173:         const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);

3175:         if (subpoint < 0) continue;
3176:         if (newLocalPoints[point].rank == rank) {++ll; continue;}
3177:         slocalPoints[sl]        = subpoint;
3178:         sremotePoints[sl].rank  = newLocalPoints[point].rank;
3179:         sremotePoints[sl].index = newLocalPoints[point].index;
3180:         if (sremotePoints[sl].rank  < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3181:         if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3182:         ++sl;
3183:       }
3184:       if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3185:       PetscFree2(newLocalPoints,newOwners);
3186:       PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3187:     }
3188:     if (subpIS) {
3189:       ISRestoreIndices(subpIS, &subpoints);
3190:       ISDestroy(&subpIS);
3191:     }
3192:   }
3193:   /* Cleanup */
3194:   for (d = 0; d <= dim; ++d) {
3195:     if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3196:     ISDestroy(&subpointIS[d]);
3197:   }
3198:   PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3199:   return(0);
3200: }

3204: /*
3205:   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.

3207:   Input Parameters:
3208: + dm          - The original mesh
3209: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3210: . label       - A label name, or NULL
3211: - value  - A label value

3213:   Output Parameter:
3214: . subdm - The surface mesh

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

3218:   Level: developer

3220: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3221: */
3222: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3223: {
3224:   PetscInt       dim, depth;

3230:   DMGetDimension(dm, &dim);
3231:   DMPlexGetDepth(dm, &depth);
3232:   DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3233:   DMSetType(*subdm, DMPLEX);
3234:   DMSetDimension(*subdm, dim-1);
3235:   if (depth == dim) {
3236:     DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3237:   } else {
3238:     DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3239:   }
3240:   return(0);
3241: }

3245: /*@
3246:   DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values

3248:   Input Parameter:
3249: . dm - The submesh DM

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

3254:   Level: developer

3256: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3257: @*/
3258: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3259: {
3260:   DM_Plex *mesh = (DM_Plex*) dm->data;

3265:   *subpointMap = mesh->subpointMap;
3266:   return(0);
3267: }

3271: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3272: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3273: {
3274:   DM_Plex       *mesh = (DM_Plex *) dm->data;
3275:   DMLabel        tmp;

3280:   tmp  = mesh->subpointMap;
3281:   mesh->subpointMap = subpointMap;
3282:   ++mesh->subpointMap->refct;
3283:   DMLabelDestroy(&tmp);
3284:   return(0);
3285: }

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

3292:   Input Parameter:
3293: . dm - The submesh DM

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

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

3300:   Level: developer

3302: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3303: @*/
3304: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3305: {
3306:   MPI_Comm        comm;
3307:   DMLabel         subpointMap;
3308:   IS              is;
3309:   const PetscInt *opoints;
3310:   PetscInt       *points, *depths;
3311:   PetscInt        depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3312:   PetscErrorCode  ierr;

3317:   PetscObjectGetComm((PetscObject)dm,&comm);
3318:   *subpointIS = NULL;
3319:   DMPlexGetSubpointMap(dm, &subpointMap);
3320:   DMPlexGetDepth(dm, &depth);
3321:   if (subpointMap && depth >= 0) {
3322:     DMPlexGetChart(dm, &pStart, &pEnd);
3323:     if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3324:     DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3325:     depths[0] = depth;
3326:     depths[1] = 0;
3327:     for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3328:     PetscMalloc1(pEnd, &points);
3329:     for(d = 0, off = 0; d <= depth; ++d) {
3330:       const PetscInt dep = depths[d];

3332:       DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3333:       DMLabelGetStratumSize(subpointMap, dep, &n);
3334:       if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3335:         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);
3336:       } else {
3337:         if (!n) {
3338:           if (d == 0) {
3339:             /* Missing cells */
3340:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3341:           } else {
3342:             /* Missing faces */
3343:             for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3344:           }
3345:         }
3346:       }
3347:       if (n) {
3348:         DMLabelGetStratumIS(subpointMap, dep, &is);
3349:         ISGetIndices(is, &opoints);
3350:         for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3351:         ISRestoreIndices(is, &opoints);
3352:         ISDestroy(&is);
3353:       }
3354:     }
3355:     DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3356:     if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3357:     ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3358:   }
3359:   return(0);
3360: }