Actual source code: plexsubmesh.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
3: #include <petscsf.h>
7: PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt cellHeight, DMLabel label)
8: {
9: PetscInt fStart, fEnd, f;
14: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
15: for (f = fStart; f < fEnd; ++f) {
16: PetscInt supportSize;
18: DMPlexGetSupportSize(dm, f, &supportSize);
19: if (supportSize == 1) {DMLabelSetValue(label, f, 1);}
20: }
21: return(0);
22: }
26: /*@
27: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
29: Not Collective
31: Input Parameter:
32: . dm - The original DM
34: Output Parameter:
35: . label - The DMLabel marking boundary faces with value 1
37: Level: developer
39: .seealso: DMLabelCreate(), DMCreateLabel()
40: @*/
41: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, DMLabel label)
42: {
46: DMPlexMarkBoundaryFaces_Internal(dm, 0, label);
47: return(0);
48: }
52: PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
53: {
54: IS valueIS;
55: const PetscInt *values;
56: PetscInt numValues, v, cStart, cEnd;
57: PetscErrorCode ierr;
60: DMLabelGetNumValues(label, &numValues);
61: DMLabelGetValueIS(label, &valueIS);
62: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
63: ISGetIndices(valueIS, &values);
64: for (v = 0; v < numValues; ++v) {
65: IS pointIS;
66: const PetscInt *points;
67: PetscInt numPoints, p;
69: DMLabelGetStratumSize(label, values[v], &numPoints);
70: DMLabelGetStratumIS(label, values[v], &pointIS);
71: ISGetIndices(pointIS, &points);
72: for (p = 0; p < numPoints; ++p) {
73: PetscInt q = points[p];
74: PetscInt *closure = NULL;
75: PetscInt closureSize, c;
77: if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
78: continue;
79: }
80: DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
81: for (c = 0; c < closureSize*2; c += 2) {
82: DMLabelSetValue(label, closure[c], values[v]);
83: }
84: DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
85: }
86: ISRestoreIndices(pointIS, &points);
87: ISDestroy(&pointIS);
88: }
89: ISRestoreIndices(valueIS, &values);
90: ISDestroy(&valueIS);
91: return(0);
92: }
96: /*@
97: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
99: Input Parameters:
100: + dm - The DM
101: - label - A DMLabel marking the surface points
103: Output Parameter:
104: . label - A DMLabel marking all surface points in the transitive closure
106: Level: developer
108: .seealso: DMPlexLabelCohesiveComplete()
109: @*/
110: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
111: {
115: DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
116: return(0);
117: }
121: /*@
122: DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face
124: Input Parameters:
125: + dm - The DM
126: - label - A DMLabel marking the surface points
128: Output Parameter:
129: . label - A DMLabel incorporating cells
131: Level: developer
133: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
135: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
136: @*/
137: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
138: {
139: IS valueIS;
140: const PetscInt *values;
141: PetscInt numValues, v, cStart, cEnd, cEndInterior;
142: PetscErrorCode ierr;
145: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
146: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
147: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
148: DMLabelGetNumValues(label, &numValues);
149: DMLabelGetValueIS(label, &valueIS);
150: ISGetIndices(valueIS, &values);
151: for (v = 0; v < numValues; ++v) {
152: IS pointIS;
153: const PetscInt *points;
154: PetscInt numPoints, p;
156: DMLabelGetStratumSize(label, values[v], &numPoints);
157: DMLabelGetStratumIS(label, values[v], &pointIS);
158: ISGetIndices(pointIS, &points);
159: for (p = 0; p < numPoints; ++p) {
160: PetscInt *closure = NULL;
161: PetscInt closureSize, point, cl;
163: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
164: for (cl = closureSize-1; cl > 0; --cl) {
165: point = closure[cl*2];
166: if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
167: }
168: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
169: }
170: ISRestoreIndices(pointIS, &points);
171: ISDestroy(&pointIS);
172: }
173: ISRestoreIndices(valueIS, &values);
174: ISDestroy(&valueIS);
175: return(0);
176: }
180: /*@
181: DMPlexLabelClearCells - Remove cells from a label
183: Input Parameters:
184: + dm - The DM
185: - label - A DMLabel marking surface points and their adjacent cells
187: Output Parameter:
188: . label - A DMLabel without cells
190: Level: developer
192: Note: This undoes DMPlexLabelAddCells()
194: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
195: @*/
196: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
197: {
198: IS valueIS;
199: const PetscInt *values;
200: PetscInt numValues, v, cStart, cEnd, cEndInterior;
201: PetscErrorCode ierr;
204: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
205: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
206: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
207: DMLabelGetNumValues(label, &numValues);
208: DMLabelGetValueIS(label, &valueIS);
209: ISGetIndices(valueIS, &values);
210: for (v = 0; v < numValues; ++v) {
211: IS pointIS;
212: const PetscInt *points;
213: PetscInt numPoints, p;
215: DMLabelGetStratumSize(label, values[v], &numPoints);
216: DMLabelGetStratumIS(label, values[v], &pointIS);
217: ISGetIndices(pointIS, &points);
218: for (p = 0; p < numPoints; ++p) {
219: PetscInt point = points[p];
221: if (point >= cStart && point < cEnd) {
222: DMLabelClearValue(label,point,values[v]);
223: }
224: }
225: ISRestoreIndices(pointIS, &points);
226: ISDestroy(&pointIS);
227: }
228: ISRestoreIndices(valueIS, &values);
229: ISDestroy(&valueIS);
230: return(0);
231: }
235: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
236: * index (skipping first, which is (0,0)) */
237: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
238: {
239: PetscInt d, off = 0;
242: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
243: for (d = 0; d < depth; d++) {
244: PetscInt firstd = d;
245: PetscInt firstStart = depthShift[2*d];
246: PetscInt e;
248: for (e = d+1; e <= depth; e++) {
249: if (depthShift[2*e] < firstStart) {
250: firstd = e;
251: firstStart = depthShift[2*d];
252: }
253: }
254: if (firstd != d) {
255: PetscInt swap[2];
257: e = firstd;
258: swap[0] = depthShift[2*d];
259: swap[1] = depthShift[2*d+1];
260: depthShift[2*d] = depthShift[2*e];
261: depthShift[2*d+1] = depthShift[2*e+1];
262: depthShift[2*e] = swap[0];
263: depthShift[2*e+1] = swap[1];
264: }
265: }
266: /* convert (oldstart, added) to (oldstart, newstart) */
267: for (d = 0; d <= depth; d++) {
268: off += depthShift[2*d+1];
269: depthShift[2*d+1] = depthShift[2*d] + off;
270: }
271: return(0);
272: }
276: /* depthShift is a list of (old, new) pairs */
277: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
278: {
279: PetscInt d;
280: PetscInt newOff = 0;
282: for (d = 0; d <= depth; d++) {
283: if (p < depthShift[2*d]) return p + newOff;
284: else newOff = depthShift[2*d+1] - depthShift[2*d];
285: }
286: return p + newOff;
287: }
291: /* depthShift is a list of (old, new) pairs */
292: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
293: {
294: PetscInt d;
295: PetscInt newOff = 0;
297: for (d = 0; d <= depth; d++) {
298: if (p < depthShift[2*d+1]) return p + newOff;
299: else newOff = depthShift[2*d] - depthShift[2*d+1];
300: }
301: return p + newOff;
302: }
306: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
307: {
308: PetscInt depth = 0, d, pStart, pEnd, p;
312: DMPlexGetDepth(dm, &depth);
313: if (depth < 0) return(0);
314: /* Step 1: Expand chart */
315: DMPlexGetChart(dm, &pStart, &pEnd);
316: pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
317: DMPlexSetChart(dmNew, pStart, pEnd);
318: /* Step 2: Set cone and support sizes */
319: for (d = 0; d <= depth; ++d) {
320: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
321: for (p = pStart; p < pEnd; ++p) {
322: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
323: PetscInt size;
325: DMPlexGetConeSize(dm, p, &size);
326: DMPlexSetConeSize(dmNew, newp, size);
327: DMPlexGetSupportSize(dm, p, &size);
328: DMPlexSetSupportSize(dmNew, newp, size);
329: }
330: }
331: return(0);
332: }
336: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
337: {
338: PetscInt *newpoints;
339: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
343: DMPlexGetDepth(dm, &depth);
344: if (depth < 0) return(0);
345: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
346: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
347: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
348: /* Step 5: Set cones and supports */
349: DMPlexGetChart(dm, &pStart, &pEnd);
350: for (p = pStart; p < pEnd; ++p) {
351: const PetscInt *points = NULL, *orientations = NULL;
352: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
354: DMPlexGetConeSize(dm, p, &size);
355: DMPlexGetCone(dm, p, &points);
356: DMPlexGetConeOrientation(dm, p, &orientations);
357: for (i = 0; i < size; ++i) {
358: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
359: }
360: DMPlexSetCone(dmNew, newp, newpoints);
361: DMPlexSetConeOrientation(dmNew, newp, orientations);
362: DMPlexGetSupportSize(dm, p, &size);
363: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
364: DMPlexGetSupport(dm, p, &points);
365: for (i = 0; i < size; ++i) {
366: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
367: }
368: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
369: DMPlexSetSupport(dmNew, newp, newpoints);
370: }
371: PetscFree(newpoints);
372: return(0);
373: }
377: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
378: {
379: PetscSection coordSection, newCoordSection;
380: Vec coordinates, newCoordinates;
381: PetscScalar *coords, *newCoords;
382: PetscInt coordSize, sStart, sEnd;
383: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
384: PetscBool hasCells;
388: DMGetCoordinateDim(dm, &dim);
389: DMPlexGetDepth(dm, &depth);
390: /* Step 8: Convert coordinates */
391: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
392: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
393: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
394: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
395: DMGetCoordinateSection(dm, &coordSection);
396: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
397: PetscSectionSetNumFields(newCoordSection, 1);
398: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
399: PetscSectionGetChart(coordSection, &sStart, &sEnd);
400: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
401: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
402: if (hasCells) {
403: for (c = cStart; c < cEnd; ++c) {
404: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
406: PetscSectionGetDof(coordSection, c, &dof);
407: PetscSectionSetDof(newCoordSection, cNew, dof);
408: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
409: }
410: }
411: for (v = vStartNew; v < vEndNew; ++v) {
412: PetscSectionSetDof(newCoordSection, v, dim);
413: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
414: }
415: PetscSectionSetUp(newCoordSection);
416: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
417: PetscSectionGetStorageSize(newCoordSection, &coordSize);
418: VecCreate(PETSC_COMM_SELF, &newCoordinates);
419: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
420: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
421: VecSetBlockSize(newCoordinates, dim);
422: VecSetType(newCoordinates,VECSTANDARD);
423: DMSetCoordinatesLocal(dmNew, newCoordinates);
424: DMGetCoordinatesLocal(dm, &coordinates);
425: VecGetArray(coordinates, &coords);
426: VecGetArray(newCoordinates, &newCoords);
427: if (hasCells) {
428: for (c = cStart; c < cEnd; ++c) {
429: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
431: PetscSectionGetDof(coordSection, c, &dof);
432: PetscSectionGetOffset(coordSection, c, &off);
433: PetscSectionGetOffset(newCoordSection, cNew, &noff);
434: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
435: }
436: }
437: for (v = vStart; v < vEnd; ++v) {
438: PetscInt dof, off, noff, d;
440: PetscSectionGetDof(coordSection, v, &dof);
441: PetscSectionGetOffset(coordSection, v, &off);
442: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
443: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
444: }
445: VecRestoreArray(coordinates, &coords);
446: VecRestoreArray(newCoordinates, &newCoords);
447: VecDestroy(&newCoordinates);
448: PetscSectionDestroy(&newCoordSection);
449: return(0);
450: }
454: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
455: {
456: PetscInt depth = 0;
457: PetscSF sfPoint, sfPointNew;
458: const PetscSFNode *remotePoints;
459: PetscSFNode *gremotePoints;
460: const PetscInt *localPoints;
461: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
462: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
463: PetscErrorCode ierr;
466: DMPlexGetDepth(dm, &depth);
467: /* Step 9: Convert pointSF */
468: DMGetPointSF(dm, &sfPoint);
469: DMGetPointSF(dmNew, &sfPointNew);
470: DMPlexGetChart(dm, &pStart, &pEnd);
471: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
472: totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
473: if (numRoots >= 0) {
474: PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
475: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
476: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
477: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
478: PetscMalloc1(numLeaves, &glocalPoints);
479: PetscMalloc1(numLeaves, &gremotePoints);
480: for (l = 0; l < numLeaves; ++l) {
481: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
482: gremotePoints[l].rank = remotePoints[l].rank;
483: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
484: }
485: PetscFree2(newLocation,newRemoteLocation);
486: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
487: }
488: return(0);
489: }
493: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
494: {
495: PetscSF sfPoint;
496: DMLabel vtkLabel, ghostLabel;
497: const PetscSFNode *leafRemote;
498: const PetscInt *leafLocal;
499: PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
500: PetscMPIInt rank;
501: PetscErrorCode ierr;
504: DMPlexGetDepth(dm, &depth);
505: /* Step 10: Convert labels */
506: DMGetNumLabels(dm, &numLabels);
507: for (l = 0; l < numLabels; ++l) {
508: DMLabel label, newlabel;
509: const char *lname;
510: PetscBool isDepth;
511: IS valueIS;
512: const PetscInt *values;
513: PetscInt numValues, val;
515: DMGetLabelName(dm, l, &lname);
516: PetscStrcmp(lname, "depth", &isDepth);
517: if (isDepth) continue;
518: DMCreateLabel(dmNew, lname);
519: DMGetLabel(dm, lname, &label);
520: DMGetLabel(dmNew, lname, &newlabel);
521: DMLabelGetDefaultValue(label,&val);
522: DMLabelSetDefaultValue(newlabel,val);
523: DMLabelGetValueIS(label, &valueIS);
524: ISGetLocalSize(valueIS, &numValues);
525: ISGetIndices(valueIS, &values);
526: for (val = 0; val < numValues; ++val) {
527: IS pointIS;
528: const PetscInt *points;
529: PetscInt numPoints, p;
531: DMLabelGetStratumIS(label, values[val], &pointIS);
532: ISGetLocalSize(pointIS, &numPoints);
533: ISGetIndices(pointIS, &points);
534: for (p = 0; p < numPoints; ++p) {
535: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
537: DMLabelSetValue(newlabel, newpoint, values[val]);
538: }
539: ISRestoreIndices(pointIS, &points);
540: ISDestroy(&pointIS);
541: }
542: ISRestoreIndices(valueIS, &values);
543: ISDestroy(&valueIS);
544: }
545: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
546: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
547: DMGetPointSF(dm, &sfPoint);
548: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
549: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
550: DMCreateLabel(dmNew, "vtk");
551: DMCreateLabel(dmNew, "ghost");
552: DMGetLabel(dmNew, "vtk", &vtkLabel);
553: DMGetLabel(dmNew, "ghost", &ghostLabel);
554: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
555: for (; c < leafLocal[l] && c < cEnd; ++c) {
556: DMLabelSetValue(vtkLabel, c, 1);
557: }
558: if (leafLocal[l] >= cEnd) break;
559: if (leafRemote[l].rank == rank) {
560: DMLabelSetValue(vtkLabel, c, 1);
561: } else {
562: DMLabelSetValue(ghostLabel, c, 2);
563: }
564: }
565: for (; c < cEnd; ++c) {
566: DMLabelSetValue(vtkLabel, c, 1);
567: }
568: if (0) {
569: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
570: }
571: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
572: for (f = fStart; f < fEnd; ++f) {
573: PetscInt numCells;
575: DMPlexGetSupportSize(dmNew, f, &numCells);
576: if (numCells < 2) {
577: DMLabelSetValue(ghostLabel, f, 1);
578: } else {
579: const PetscInt *cells = NULL;
580: PetscInt vA, vB;
582: DMPlexGetSupport(dmNew, f, &cells);
583: DMLabelGetValue(vtkLabel, cells[0], &vA);
584: DMLabelGetValue(vtkLabel, cells[1], &vB);
585: if (!vA && !vB) {DMLabelSetValue(ghostLabel, f, 1);}
586: }
587: }
588: if (0) {
589: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
590: }
591: return(0);
592: }
596: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
597: {
598: DM refTree;
599: PetscSection pSec;
600: PetscInt *parents, *childIDs;
604: DMPlexGetReferenceTree(dm,&refTree);
605: DMPlexSetReferenceTree(dmNew,refTree);
606: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
607: if (pSec) {
608: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
609: PetscInt *childIDsShifted;
610: PetscSection pSecShifted;
612: PetscSectionGetChart(pSec,&pStart,&pEnd);
613: DMPlexGetDepth(dm,&depth);
614: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
615: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
616: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
617: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
618: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
619: for (p = pStartShifted; p < pEndShifted; p++) {
620: /* start off assuming no children */
621: PetscSectionSetDof(pSecShifted,p,0);
622: }
623: for (p = pStart; p < pEnd; p++) {
624: PetscInt dof;
625: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
627: PetscSectionGetDof(pSec,p,&dof);
628: PetscSectionSetDof(pSecShifted,pNew,dof);
629: }
630: PetscSectionSetUp(pSecShifted);
631: for (p = pStart; p < pEnd; p++) {
632: PetscInt dof;
633: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
635: PetscSectionGetDof(pSec,p,&dof);
636: if (dof) {
637: PetscInt off, offNew;
639: PetscSectionGetOffset(pSec,p,&off);
640: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
641: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
642: childIDsShifted[offNew] = childIDs[off];
643: }
644: }
645: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
646: PetscFree2(parentsShifted,childIDsShifted);
647: PetscSectionDestroy(&pSecShifted);
648: }
649: return(0);
650: }
654: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
655: {
656: PetscSF sf;
657: IS valueIS;
658: const PetscInt *values, *leaves;
659: PetscInt *depthShift;
660: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
661: PetscErrorCode ierr;
664: DMGetPointSF(dm, &sf);
665: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
666: nleaves = PetscMax(0, nleaves);
667: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
668: /* Count ghost cells */
669: DMLabelGetValueIS(label, &valueIS);
670: ISGetLocalSize(valueIS, &numFS);
671: ISGetIndices(valueIS, &values);
672: Ng = 0;
673: for (fs = 0; fs < numFS; ++fs) {
674: IS faceIS;
675: const PetscInt *faces;
676: PetscInt numFaces, f, numBdFaces = 0;
678: DMLabelGetStratumIS(label, values[fs], &faceIS);
679: ISGetLocalSize(faceIS, &numFaces);
680: ISGetIndices(faceIS, &faces);
681: for (f = 0; f < numFaces; ++f) {
682: PetscInt numChildren;
684: PetscFindInt(faces[f], nleaves, leaves, &loc);
685: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
686: /* non-local and ancestors points don't get to register ghosts */
687: if (loc >= 0 || numChildren) continue;
688: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
689: }
690: Ng += numBdFaces;
691: ISDestroy(&faceIS);
692: }
693: DMPlexGetDepth(dm, &depth);
694: PetscMalloc1(2*(depth+1), &depthShift);
695: for (d = 0; d <= depth; d++) {
696: PetscInt dEnd;
698: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
699: depthShift[2*d] = dEnd;
700: depthShift[2*d+1] = 0;
701: }
702: if (depth >= 0) depthShift[2*depth+1] = Ng;
703: DMPlexShiftPointSetUp_Internal(depth,depthShift);
704: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
705: /* Step 3: Set cone/support sizes for new points */
706: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
707: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
708: for (c = cEnd; c < cEnd + Ng; ++c) {
709: DMPlexSetConeSize(gdm, c, 1);
710: }
711: for (fs = 0; fs < numFS; ++fs) {
712: IS faceIS;
713: const PetscInt *faces;
714: PetscInt numFaces, f;
716: DMLabelGetStratumIS(label, values[fs], &faceIS);
717: ISGetLocalSize(faceIS, &numFaces);
718: ISGetIndices(faceIS, &faces);
719: for (f = 0; f < numFaces; ++f) {
720: PetscInt size, numChildren;
722: PetscFindInt(faces[f], nleaves, leaves, &loc);
723: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
724: if (loc >= 0 || numChildren) continue;
725: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
726: DMPlexGetSupportSize(dm, faces[f], &size);
727: if (size != 1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
728: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
729: }
730: ISRestoreIndices(faceIS, &faces);
731: ISDestroy(&faceIS);
732: }
733: /* Step 4: Setup ghosted DM */
734: DMSetUp(gdm);
735: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
736: /* Step 6: Set cones and supports for new points */
737: ghostCell = cEnd;
738: for (fs = 0; fs < numFS; ++fs) {
739: IS faceIS;
740: const PetscInt *faces;
741: PetscInt numFaces, f;
743: DMLabelGetStratumIS(label, values[fs], &faceIS);
744: ISGetLocalSize(faceIS, &numFaces);
745: ISGetIndices(faceIS, &faces);
746: for (f = 0; f < numFaces; ++f) {
747: PetscInt newFace = faces[f] + Ng, numChildren;
749: PetscFindInt(faces[f], nleaves, leaves, &loc);
750: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
751: if (loc >= 0 || numChildren) continue;
752: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
753: DMPlexSetCone(gdm, ghostCell, &newFace);
754: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
755: ++ghostCell;
756: }
757: ISRestoreIndices(faceIS, &faces);
758: ISDestroy(&faceIS);
759: }
760: ISRestoreIndices(valueIS, &values);
761: ISDestroy(&valueIS);
762: /* Step 7: Stratify */
763: DMPlexStratify(gdm);
764: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
765: DMPlexShiftSF_Internal(dm, depthShift, gdm);
766: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
767: DMPlexShiftTree_Internal(dm, depthShift, gdm);
768: PetscFree(depthShift);
769: /* Step 7: Periodicity */
770: if (dm->maxCell) {
771: const PetscReal *maxCell, *L;
772: const DMBoundaryType *bd;
773: DMGetPeriodicity(dm, &maxCell, &L, &bd);
774: DMSetPeriodicity(gdm, maxCell, L, bd);
775: }
776: if (numGhostCells) *numGhostCells = Ng;
777: return(0);
778: }
782: /*@C
783: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
785: Collective on dm
787: Input Parameters:
788: + dm - The original DM
789: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
791: Output Parameters:
792: + numGhostCells - The number of ghost cells added to the DM
793: - dmGhosted - The new DM
795: Note: If no label exists of that name, one will be created marking all boundary faces
797: Level: developer
799: .seealso: DMCreate()
800: @*/
801: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
802: {
803: DM gdm;
804: DMLabel label;
805: const char *name = labelName ? labelName : "Face Sets";
806: PetscInt dim;
807: PetscBool flag;
814: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
815: DMSetType(gdm, DMPLEX);
816: DMGetDimension(dm, &dim);
817: DMSetDimension(gdm, dim);
818: DMPlexGetAdjacencyUseCone(dm, &flag);
819: DMPlexSetAdjacencyUseCone(gdm, flag);
820: DMPlexGetAdjacencyUseClosure(dm, &flag);
821: DMPlexSetAdjacencyUseClosure(gdm, flag);
822: DMGetLabel(dm, name, &label);
823: if (!label) {
824: /* Get label for boundary faces */
825: DMCreateLabel(dm, name);
826: DMGetLabel(dm, name, &label);
827: DMPlexMarkBoundaryFaces(dm, label);
828: }
829: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
830: DMCopyBoundary(dm, gdm);
831: *dmGhosted = gdm;
832: return(0);
833: }
837: /*
838: We are adding three kinds of points here:
839: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
840: Non-replicated: Points which exist on the fault, but are not replicated
841: Hybrid: Entirely new points, such as cohesive cells
843: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
844: each depth so that the new split/hybrid points can be inserted as a block.
845: */
846: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
847: {
848: MPI_Comm comm;
849: IS valueIS;
850: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
851: const PetscInt *values; /* List of depths for which we have replicated points */
852: IS *splitIS;
853: IS *unsplitIS;
854: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
855: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
856: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
857: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
858: const PetscInt **splitPoints; /* Replicated points for each depth */
859: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
860: PetscSection coordSection;
861: Vec coordinates;
862: PetscScalar *coords;
863: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
864: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
865: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
866: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
867: PetscInt *coneNew, *coneONew, *supportNew;
868: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
869: PetscErrorCode ierr;
872: PetscObjectGetComm((PetscObject)dm,&comm);
873: DMGetDimension(dm, &dim);
874: DMPlexGetDepth(dm, &depth);
875: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
876: /* Count split points and add cohesive cells */
877: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
878: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
879: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
880: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
881: for (d = 0; d <= depth; ++d) {
882: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
883: depthEnd[d] = pMaxNew[d];
884: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
885: numSplitPoints[d] = 0;
886: numUnsplitPoints[d] = 0;
887: numHybridPoints[d] = 0;
888: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
889: splitPoints[d] = NULL;
890: unsplitPoints[d] = NULL;
891: splitIS[d] = NULL;
892: unsplitIS[d] = NULL;
893: /* we are shifting the existing hybrid points with the stratum behind them, so
894: * the split comes at the end of the normal points, i.e., at depthMax[d] */
895: depthShift[2*d] = depthMax[d];
896: depthShift[2*d+1] = 0;
897: }
898: if (label) {
899: DMLabelGetValueIS(label, &valueIS);
900: ISGetLocalSize(valueIS, &numSP);
901: ISGetIndices(valueIS, &values);
902: }
903: for (sp = 0; sp < numSP; ++sp) {
904: const PetscInt dep = values[sp];
906: if ((dep < 0) || (dep > depth)) continue;
907: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
908: if (splitIS[dep]) {
909: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
910: ISGetIndices(splitIS[dep], &splitPoints[dep]);
911: }
912: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
913: if (unsplitIS[dep]) {
914: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
915: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
916: }
917: }
918: /* Calculate number of hybrid points */
919: 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 */
920: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
921: DMPlexShiftPointSetUp_Internal(depth,depthShift);
922: /* the end of the points in this stratum that come before the new points:
923: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
924: * added points */
925: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
926: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
927: /* Step 3: Set cone/support sizes for new points */
928: for (dep = 0; dep <= depth; ++dep) {
929: for (p = 0; p < numSplitPoints[dep]; ++p) {
930: const PetscInt oldp = splitPoints[dep][p];
931: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
932: const PetscInt splitp = p + pMaxNew[dep];
933: const PetscInt *support;
934: PetscInt coneSize, supportSize, qf, qn, qp, e;
936: DMPlexGetConeSize(dm, oldp, &coneSize);
937: DMPlexSetConeSize(sdm, splitp, coneSize);
938: DMPlexGetSupportSize(dm, oldp, &supportSize);
939: DMPlexSetSupportSize(sdm, splitp, supportSize);
940: if (dep == depth-1) {
941: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
943: /* Add cohesive cells, they are prisms */
944: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
945: } else if (dep == 0) {
946: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
948: DMPlexGetSupport(dm, oldp, &support);
949: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
950: PetscInt val;
952: DMLabelGetValue(label, support[e], &val);
953: if (val == 1) ++qf;
954: if ((val == 1) || (val == (shift + 1))) ++qn;
955: if ((val == 1) || (val == -(shift + 1))) ++qp;
956: }
957: /* Split old vertex: Edges into original vertex and new cohesive edge */
958: DMPlexSetSupportSize(sdm, newp, qn+1);
959: /* Split new vertex: Edges into split vertex and new cohesive edge */
960: DMPlexSetSupportSize(sdm, splitp, qp+1);
961: /* Add hybrid edge */
962: DMPlexSetConeSize(sdm, hybedge, 2);
963: DMPlexSetSupportSize(sdm, hybedge, qf);
964: } else if (dep == dim-2) {
965: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
967: DMPlexGetSupport(dm, oldp, &support);
968: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
969: PetscInt val;
971: DMLabelGetValue(label, support[e], &val);
972: if (val == dim-1) ++qf;
973: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
974: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
975: }
976: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
977: DMPlexSetSupportSize(sdm, newp, qn+1);
978: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
979: DMPlexSetSupportSize(sdm, splitp, qp+1);
980: /* Add hybrid face */
981: DMPlexSetConeSize(sdm, hybface, 4);
982: DMPlexSetSupportSize(sdm, hybface, qf);
983: }
984: }
985: }
986: for (dep = 0; dep <= depth; ++dep) {
987: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
988: const PetscInt oldp = unsplitPoints[dep][p];
989: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
990: const PetscInt *support;
991: PetscInt coneSize, supportSize, qf, e, s;
993: DMPlexGetConeSize(dm, oldp, &coneSize);
994: DMPlexGetSupportSize(dm, oldp, &supportSize);
995: DMPlexGetSupport(dm, oldp, &support);
996: if (dep == 0) {
997: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
999: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1000: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1001: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1002: if (e >= 0) ++qf;
1003: }
1004: DMPlexSetSupportSize(sdm, newp, qf+2);
1005: /* Add hybrid edge */
1006: DMPlexSetConeSize(sdm, hybedge, 2);
1007: for (e = 0, qf = 0; e < supportSize; ++e) {
1008: PetscInt val;
1010: DMLabelGetValue(label, support[e], &val);
1011: /* Split and unsplit edges produce hybrid faces */
1012: if (val == 1) ++qf;
1013: if (val == (shift2 + 1)) ++qf;
1014: }
1015: DMPlexSetSupportSize(sdm, hybedge, qf);
1016: } else if (dep == dim-2) {
1017: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1018: PetscInt val;
1020: for (e = 0, qf = 0; e < supportSize; ++e) {
1021: DMLabelGetValue(label, support[e], &val);
1022: if (val == dim-1) qf += 2;
1023: else ++qf;
1024: }
1025: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1026: DMPlexSetSupportSize(sdm, newp, qf+2);
1027: /* Add hybrid face */
1028: for (e = 0, qf = 0; e < supportSize; ++e) {
1029: DMLabelGetValue(label, support[e], &val);
1030: if (val == dim-1) ++qf;
1031: }
1032: DMPlexSetConeSize(sdm, hybface, 4);
1033: DMPlexSetSupportSize(sdm, hybface, qf);
1034: }
1035: }
1036: }
1037: /* Step 4: Setup split DM */
1038: DMSetUp(sdm);
1039: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1040: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1041: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1042: /* Step 6: Set cones and supports for new points */
1043: for (dep = 0; dep <= depth; ++dep) {
1044: for (p = 0; p < numSplitPoints[dep]; ++p) {
1045: const PetscInt oldp = splitPoints[dep][p];
1046: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1047: const PetscInt splitp = p + pMaxNew[dep];
1048: const PetscInt *cone, *support, *ornt;
1049: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1051: DMPlexGetConeSize(dm, oldp, &coneSize);
1052: DMPlexGetCone(dm, oldp, &cone);
1053: DMPlexGetConeOrientation(dm, oldp, &ornt);
1054: DMPlexGetSupportSize(dm, oldp, &supportSize);
1055: DMPlexGetSupport(dm, oldp, &support);
1056: if (dep == depth-1) {
1057: PetscBool hasUnsplit = PETSC_FALSE;
1058: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1059: const PetscInt *supportF;
1061: /* Split face: copy in old face to new face to start */
1062: DMPlexGetSupport(sdm, newp, &supportF);
1063: DMPlexSetSupport(sdm, splitp, supportF);
1064: /* Split old face: old vertices/edges in cone so no change */
1065: /* Split new face: new vertices/edges in cone */
1066: for (q = 0; q < coneSize; ++q) {
1067: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1068: if (v < 0) {
1069: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1070: 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);
1071: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1072: hasUnsplit = PETSC_TRUE;
1073: } else {
1074: coneNew[2+q] = v + pMaxNew[dep-1];
1075: if (dep > 1) {
1076: const PetscInt *econe;
1077: PetscInt econeSize, r, vs, vu;
1079: DMPlexGetConeSize(dm, cone[q], &econeSize);
1080: DMPlexGetCone(dm, cone[q], &econe);
1081: for (r = 0; r < econeSize; ++r) {
1082: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1083: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1084: if (vs >= 0) continue;
1085: 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);
1086: hasUnsplit = PETSC_TRUE;
1087: }
1088: }
1089: }
1090: }
1091: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1092: DMPlexSetConeOrientation(sdm, splitp, ornt);
1093: /* Face support */
1094: for (s = 0; s < supportSize; ++s) {
1095: PetscInt val;
1097: DMLabelGetValue(label, support[s], &val);
1098: if (val < 0) {
1099: /* Split old face: Replace negative side cell with cohesive cell */
1100: DMPlexInsertSupport(sdm, newp, s, hybcell);
1101: } else {
1102: /* Split new face: Replace positive side cell with cohesive cell */
1103: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1104: /* Get orientation for cohesive face */
1105: {
1106: const PetscInt *ncone, *nconeO;
1107: PetscInt nconeSize, nc;
1109: DMPlexGetConeSize(dm, support[s], &nconeSize);
1110: DMPlexGetCone(dm, support[s], &ncone);
1111: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1112: for (nc = 0; nc < nconeSize; ++nc) {
1113: if (ncone[nc] == oldp) {
1114: coneONew[0] = nconeO[nc];
1115: break;
1116: }
1117: }
1118: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1119: }
1120: }
1121: }
1122: /* Cohesive cell: Old and new split face, then new cohesive faces */
1123: coneNew[0] = newp; /* Extracted negative side orientation above */
1124: coneNew[1] = splitp;
1125: coneONew[1] = coneONew[0];
1126: for (q = 0; q < coneSize; ++q) {
1127: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1128: if (v < 0) {
1129: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1130: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1131: coneONew[2+q] = 0;
1132: } else {
1133: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1134: }
1135: coneONew[2+q] = 0;
1136: }
1137: DMPlexSetCone(sdm, hybcell, coneNew);
1138: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1139: /* Label the hybrid cells on the boundary of the split */
1140: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1141: } else if (dep == 0) {
1142: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1144: /* Split old vertex: Edges in old split faces and new cohesive edge */
1145: for (e = 0, qn = 0; e < supportSize; ++e) {
1146: PetscInt val;
1148: DMLabelGetValue(label, support[e], &val);
1149: if ((val == 1) || (val == (shift + 1))) {
1150: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1151: }
1152: }
1153: supportNew[qn] = hybedge;
1154: DMPlexSetSupport(sdm, newp, supportNew);
1155: /* Split new vertex: Edges in new split faces and new cohesive edge */
1156: for (e = 0, qp = 0; e < supportSize; ++e) {
1157: PetscInt val, edge;
1159: DMLabelGetValue(label, support[e], &val);
1160: if (val == 1) {
1161: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1162: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1163: supportNew[qp++] = edge + pMaxNew[dep+1];
1164: } else if (val == -(shift + 1)) {
1165: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1166: }
1167: }
1168: supportNew[qp] = hybedge;
1169: DMPlexSetSupport(sdm, splitp, supportNew);
1170: /* Hybrid edge: Old and new split vertex */
1171: coneNew[0] = newp;
1172: coneNew[1] = splitp;
1173: DMPlexSetCone(sdm, hybedge, coneNew);
1174: for (e = 0, qf = 0; e < supportSize; ++e) {
1175: PetscInt val, edge;
1177: DMLabelGetValue(label, support[e], &val);
1178: if (val == 1) {
1179: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1180: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1181: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1182: }
1183: }
1184: DMPlexSetSupport(sdm, hybedge, supportNew);
1185: } else if (dep == dim-2) {
1186: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1188: /* Split old edge: old vertices in cone so no change */
1189: /* Split new edge: new vertices in cone */
1190: for (q = 0; q < coneSize; ++q) {
1191: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1192: if (v < 0) {
1193: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1194: 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);
1195: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1196: } else {
1197: coneNew[q] = v + pMaxNew[dep-1];
1198: }
1199: }
1200: DMPlexSetCone(sdm, splitp, coneNew);
1201: /* Split old edge: Faces in positive side cells and old split faces */
1202: for (e = 0, q = 0; e < supportSize; ++e) {
1203: PetscInt val;
1205: DMLabelGetValue(label, support[e], &val);
1206: if (val == dim-1) {
1207: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1208: } else if (val == (shift + dim-1)) {
1209: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1210: }
1211: }
1212: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1213: DMPlexSetSupport(sdm, newp, supportNew);
1214: /* Split new edge: Faces in negative side cells and new split faces */
1215: for (e = 0, q = 0; e < supportSize; ++e) {
1216: PetscInt val, face;
1218: DMLabelGetValue(label, support[e], &val);
1219: if (val == dim-1) {
1220: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1221: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1222: supportNew[q++] = face + pMaxNew[dep+1];
1223: } else if (val == -(shift + dim-1)) {
1224: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1225: }
1226: }
1227: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1228: DMPlexSetSupport(sdm, splitp, supportNew);
1229: /* Hybrid face */
1230: coneNew[0] = newp;
1231: coneNew[1] = splitp;
1232: for (v = 0; v < coneSize; ++v) {
1233: PetscInt vertex;
1234: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1235: if (vertex < 0) {
1236: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1237: 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);
1238: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1239: } else {
1240: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1241: }
1242: }
1243: DMPlexSetCone(sdm, hybface, coneNew);
1244: for (e = 0, qf = 0; e < supportSize; ++e) {
1245: PetscInt val, face;
1247: DMLabelGetValue(label, support[e], &val);
1248: if (val == dim-1) {
1249: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1250: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1251: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1252: }
1253: }
1254: DMPlexSetSupport(sdm, hybface, supportNew);
1255: }
1256: }
1257: }
1258: for (dep = 0; dep <= depth; ++dep) {
1259: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1260: const PetscInt oldp = unsplitPoints[dep][p];
1261: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1262: const PetscInt *cone, *support, *ornt;
1263: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1265: DMPlexGetConeSize(dm, oldp, &coneSize);
1266: DMPlexGetCone(dm, oldp, &cone);
1267: DMPlexGetConeOrientation(dm, oldp, &ornt);
1268: DMPlexGetSupportSize(dm, oldp, &supportSize);
1269: DMPlexGetSupport(dm, oldp, &support);
1270: if (dep == 0) {
1271: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1273: /* Unsplit vertex */
1274: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1275: for (s = 0, q = 0; s < supportSize; ++s) {
1276: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1277: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1278: if (e >= 0) {
1279: supportNew[q++] = e + pMaxNew[dep+1];
1280: }
1281: }
1282: supportNew[q++] = hybedge;
1283: supportNew[q++] = hybedge;
1284: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1285: DMPlexSetSupport(sdm, newp, supportNew);
1286: /* Hybrid edge */
1287: coneNew[0] = newp;
1288: coneNew[1] = newp;
1289: DMPlexSetCone(sdm, hybedge, coneNew);
1290: for (e = 0, qf = 0; e < supportSize; ++e) {
1291: PetscInt val, edge;
1293: DMLabelGetValue(label, support[e], &val);
1294: if (val == 1) {
1295: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1296: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1297: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1298: } else if (val == (shift2 + 1)) {
1299: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1300: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1301: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1302: }
1303: }
1304: DMPlexSetSupport(sdm, hybedge, supportNew);
1305: } else if (dep == dim-2) {
1306: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1308: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1309: for (f = 0, qf = 0; f < supportSize; ++f) {
1310: PetscInt val, face;
1312: DMLabelGetValue(label, support[f], &val);
1313: if (val == dim-1) {
1314: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1315: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1316: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1317: supportNew[qf++] = face + pMaxNew[dep+1];
1318: } else {
1319: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1320: }
1321: }
1322: supportNew[qf++] = hybface;
1323: supportNew[qf++] = hybface;
1324: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1325: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1326: DMPlexSetSupport(sdm, newp, supportNew);
1327: /* Add hybrid face */
1328: coneNew[0] = newp;
1329: coneNew[1] = newp;
1330: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1331: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1332: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1333: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1334: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1335: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1336: DMPlexSetCone(sdm, hybface, coneNew);
1337: for (f = 0, qf = 0; f < supportSize; ++f) {
1338: PetscInt val, face;
1340: DMLabelGetValue(label, support[f], &val);
1341: if (val == dim-1) {
1342: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1343: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1344: }
1345: }
1346: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1347: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1348: DMPlexSetSupport(sdm, hybface, supportNew);
1349: }
1350: }
1351: }
1352: /* Step 6b: Replace split points in negative side cones */
1353: for (sp = 0; sp < numSP; ++sp) {
1354: PetscInt dep = values[sp];
1355: IS pIS;
1356: PetscInt numPoints;
1357: const PetscInt *points;
1359: if (dep >= 0) continue;
1360: DMLabelGetStratumIS(label, dep, &pIS);
1361: if (!pIS) continue;
1362: dep = -dep - shift;
1363: ISGetLocalSize(pIS, &numPoints);
1364: ISGetIndices(pIS, &points);
1365: for (p = 0; p < numPoints; ++p) {
1366: const PetscInt oldp = points[p];
1367: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1368: const PetscInt *cone;
1369: PetscInt coneSize, c;
1370: /* PetscBool replaced = PETSC_FALSE; */
1372: /* Negative edge: replace split vertex */
1373: /* Negative cell: replace split face */
1374: DMPlexGetConeSize(sdm, newp, &coneSize);
1375: DMPlexGetCone(sdm, newp, &cone);
1376: for (c = 0; c < coneSize; ++c) {
1377: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1378: PetscInt csplitp, cp, val;
1380: DMLabelGetValue(label, coldp, &val);
1381: if (val == dep-1) {
1382: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1383: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1384: csplitp = pMaxNew[dep-1] + cp;
1385: DMPlexInsertCone(sdm, newp, c, csplitp);
1386: /* replaced = PETSC_TRUE; */
1387: }
1388: }
1389: /* Cells with only a vertex or edge on the submesh have no replacement */
1390: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1391: }
1392: ISRestoreIndices(pIS, &points);
1393: ISDestroy(&pIS);
1394: }
1395: /* Step 7: Stratify */
1396: DMPlexStratify(sdm);
1397: /* Step 8: Coordinates */
1398: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1399: DMGetCoordinateSection(sdm, &coordSection);
1400: DMGetCoordinatesLocal(sdm, &coordinates);
1401: VecGetArray(coordinates, &coords);
1402: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1403: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1404: const PetscInt splitp = pMaxNew[0] + v;
1405: PetscInt dof, off, soff, d;
1407: PetscSectionGetDof(coordSection, newp, &dof);
1408: PetscSectionGetOffset(coordSection, newp, &off);
1409: PetscSectionGetOffset(coordSection, splitp, &soff);
1410: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1411: }
1412: VecRestoreArray(coordinates, &coords);
1413: /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1414: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1415: /* Step 10: Labels */
1416: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1417: DMGetNumLabels(sdm, &numLabels);
1418: for (dep = 0; dep <= depth; ++dep) {
1419: for (p = 0; p < numSplitPoints[dep]; ++p) {
1420: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1421: const PetscInt splitp = pMaxNew[dep] + p;
1422: PetscInt l;
1424: for (l = 0; l < numLabels; ++l) {
1425: DMLabel mlabel;
1426: const char *lname;
1427: PetscInt val;
1428: PetscBool isDepth;
1430: DMGetLabelName(sdm, l, &lname);
1431: PetscStrcmp(lname, "depth", &isDepth);
1432: if (isDepth) continue;
1433: DMGetLabel(sdm, lname, &mlabel);
1434: DMLabelGetValue(mlabel, newp, &val);
1435: if (val >= 0) {
1436: DMLabelSetValue(mlabel, splitp, val);
1437: #if 0
1438: /* Do not put cohesive edges into the label */
1439: if (dep == 0) {
1440: const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1441: DMLabelSetValue(mlabel, cedge, val);
1442: } else if (dep == dim-2) {
1443: const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1444: DMLabelSetValue(mlabel, cface, val);
1445: }
1446: /* Do not put cohesive faces into the label */
1447: #endif
1448: }
1449: }
1450: }
1451: }
1452: for (sp = 0; sp < numSP; ++sp) {
1453: const PetscInt dep = values[sp];
1455: if ((dep < 0) || (dep > depth)) continue;
1456: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1457: ISDestroy(&splitIS[dep]);
1458: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1459: ISDestroy(&unsplitIS[dep]);
1460: }
1461: if (label) {
1462: ISRestoreIndices(valueIS, &values);
1463: ISDestroy(&valueIS);
1464: }
1465: for (d = 0; d <= depth; ++d) {
1466: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1467: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1468: }
1469: 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);
1470: PetscFree3(coneNew, coneONew, supportNew);
1471: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1472: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1473: return(0);
1474: }
1478: /*@C
1479: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1481: Collective on dm
1483: Input Parameters:
1484: + dm - The original DM
1485: - label - The label specifying the boundary faces (this could be auto-generated)
1487: Output Parameters:
1488: - dmSplit - The new DM
1490: Level: developer
1492: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1493: @*/
1494: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1495: {
1496: DM sdm;
1497: PetscInt dim;
1503: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1504: DMSetType(sdm, DMPLEX);
1505: DMGetDimension(dm, &dim);
1506: DMSetDimension(sdm, dim);
1507: switch (dim) {
1508: case 2:
1509: case 3:
1510: DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1511: break;
1512: default:
1513: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1514: }
1515: *dmSplit = sdm;
1516: return(0);
1517: }
1521: /* Returns the side of the surface for a given cell with a face on the surface */
1522: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1523: {
1524: const PetscInt *cone, *ornt;
1525: PetscInt dim, coneSize, c;
1526: PetscErrorCode ierr;
1529: *pos = PETSC_TRUE;
1530: DMGetDimension(dm, &dim);
1531: DMPlexGetConeSize(dm, cell, &coneSize);
1532: DMPlexGetCone(dm, cell, &cone);
1533: DMPlexGetConeOrientation(dm, cell, &ornt);
1534: for (c = 0; c < coneSize; ++c) {
1535: if (cone[c] == face) {
1536: PetscInt o = ornt[c];
1538: if (subdm) {
1539: const PetscInt *subcone, *subornt;
1540: PetscInt subpoint, subface, subconeSize, sc;
1542: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1543: PetscFindInt(face, numSubpoints, subpoints, &subface);
1544: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1545: DMPlexGetCone(subdm, subpoint, &subcone);
1546: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1547: for (sc = 0; sc < subconeSize; ++sc) {
1548: if (subcone[sc] == subface) {
1549: o = subornt[0];
1550: break;
1551: }
1552: }
1553: 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);
1554: }
1555: if (o >= 0) *pos = PETSC_TRUE;
1556: else *pos = PETSC_FALSE;
1557: break;
1558: }
1559: }
1560: 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);
1561: return(0);
1562: }
1566: /*@
1567: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1568: to complete the surface
1570: Input Parameters:
1571: + dm - The DM
1572: . label - A DMLabel marking the surface
1573: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1574: . flip - Flag to flip the submesh normal and replace points on the other side
1575: - subdm - The subDM associated with the label, or NULL
1577: Output Parameter:
1578: . label - A DMLabel marking all surface points
1580: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1582: Level: developer
1584: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1585: @*/
1586: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1587: {
1588: DMLabel depthLabel;
1589: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1590: const PetscInt *points, *subpoints;
1591: const PetscInt rev = flip ? -1 : 1;
1592: PetscInt *pMax;
1593: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1594: PetscErrorCode ierr;
1597: DMPlexGetDepth(dm, &depth);
1598: DMGetDimension(dm, &dim);
1599: pSize = PetscMax(depth, dim) + 1;
1600: PetscMalloc1(pSize,&pMax);
1601: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1602: DMPlexGetDepthLabel(dm, &depthLabel);
1603: DMGetDimension(dm, &dim);
1604: if (subdm) {
1605: DMPlexCreateSubpointIS(subdm, &subpointIS);
1606: if (subpointIS) {
1607: ISGetLocalSize(subpointIS, &numSubpoints);
1608: ISGetIndices(subpointIS, &subpoints);
1609: }
1610: }
1611: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1612: DMLabelGetStratumIS(label, dim-1, &dimIS);
1613: if (!dimIS) {
1614: PetscFree(pMax);
1615: ISDestroy(&subpointIS);
1616: return(0);
1617: }
1618: ISGetLocalSize(dimIS, &numPoints);
1619: ISGetIndices(dimIS, &points);
1620: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1621: const PetscInt *support;
1622: PetscInt supportSize, s;
1624: DMPlexGetSupportSize(dm, points[p], &supportSize);
1625: if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1626: DMPlexGetSupport(dm, points[p], &support);
1627: for (s = 0; s < supportSize; ++s) {
1628: const PetscInt *cone;
1629: PetscInt coneSize, c;
1630: PetscBool pos;
1632: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1633: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1634: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1635: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1636: /* Put faces touching the fault in the label */
1637: DMPlexGetConeSize(dm, support[s], &coneSize);
1638: DMPlexGetCone(dm, support[s], &cone);
1639: for (c = 0; c < coneSize; ++c) {
1640: const PetscInt point = cone[c];
1642: DMLabelGetValue(label, point, &val);
1643: if (val == -1) {
1644: PetscInt *closure = NULL;
1645: PetscInt closureSize, cl;
1647: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1648: for (cl = 0; cl < closureSize*2; cl += 2) {
1649: const PetscInt clp = closure[cl];
1650: PetscInt bval = -1;
1652: DMLabelGetValue(label, clp, &val);
1653: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1654: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1655: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1656: break;
1657: }
1658: }
1659: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1660: }
1661: }
1662: }
1663: }
1664: ISRestoreIndices(dimIS, &points);
1665: ISDestroy(&dimIS);
1666: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1667: ISDestroy(&subpointIS);
1668: /* Mark boundary points as unsplit */
1669: if (blabel) {
1670: DMLabelGetStratumIS(blabel, 1, &dimIS);
1671: ISGetLocalSize(dimIS, &numPoints);
1672: ISGetIndices(dimIS, &points);
1673: for (p = 0; p < numPoints; ++p) {
1674: const PetscInt point = points[p];
1675: PetscInt val, bval;
1677: DMLabelGetValue(blabel, point, &bval);
1678: if (bval >= 0) {
1679: DMLabelGetValue(label, point, &val);
1680: if ((val < 0) || (val > dim)) {
1681: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1682: DMLabelClearValue(blabel, point, bval);
1683: }
1684: }
1685: }
1686: for (p = 0; p < numPoints; ++p) {
1687: const PetscInt point = points[p];
1688: PetscInt val, bval;
1690: DMLabelGetValue(blabel, point, &bval);
1691: if (bval >= 0) {
1692: const PetscInt *cone, *support;
1693: PetscInt coneSize, supportSize, s, valA, valB, valE;
1695: /* Mark as unsplit */
1696: DMLabelGetValue(label, point, &val);
1697: 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);
1698: DMLabelClearValue(label, point, val);
1699: DMLabelSetValue(label, point, shift2+val);
1700: /* Check for cross-edge
1701: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1702: if (val != 0) continue;
1703: DMPlexGetSupport(dm, point, &support);
1704: DMPlexGetSupportSize(dm, point, &supportSize);
1705: for (s = 0; s < supportSize; ++s) {
1706: DMPlexGetCone(dm, support[s], &cone);
1707: DMPlexGetConeSize(dm, support[s], &coneSize);
1708: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1709: DMLabelGetValue(blabel, cone[0], &valA);
1710: DMLabelGetValue(blabel, cone[1], &valB);
1711: DMLabelGetValue(blabel, support[s], &valE);
1712: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1713: }
1714: }
1715: }
1716: ISRestoreIndices(dimIS, &points);
1717: ISDestroy(&dimIS);
1718: }
1719: /* Search for other cells/faces/edges connected to the fault by a vertex */
1720: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1721: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1722: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1723: cMax = cMax < 0 ? cEnd : cMax;
1724: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1725: DMLabelGetStratumIS(label, 0, &dimIS);
1726: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1727: if (dimIS && crossEdgeIS) {
1728: IS vertIS = dimIS;
1730: ISExpand(vertIS, crossEdgeIS, &dimIS);
1731: ISDestroy(&crossEdgeIS);
1732: ISDestroy(&vertIS);
1733: }
1734: if (!dimIS) {
1735: PetscFree(pMax);
1736: return(0);
1737: }
1738: ISGetLocalSize(dimIS, &numPoints);
1739: ISGetIndices(dimIS, &points);
1740: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1741: PetscInt *star = NULL;
1742: PetscInt starSize, s;
1743: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1745: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1746: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1747: while (again) {
1748: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1749: again = 0;
1750: for (s = 0; s < starSize*2; s += 2) {
1751: const PetscInt point = star[s];
1752: const PetscInt *cone;
1753: PetscInt coneSize, c;
1755: if ((point < cStart) || (point >= cMax)) continue;
1756: DMLabelGetValue(label, point, &val);
1757: if (val != -1) continue;
1758: again = again == 1 ? 1 : 2;
1759: DMPlexGetConeSize(dm, point, &coneSize);
1760: DMPlexGetCone(dm, point, &cone);
1761: for (c = 0; c < coneSize; ++c) {
1762: DMLabelGetValue(label, cone[c], &val);
1763: if (val != -1) {
1764: const PetscInt *ccone;
1765: PetscInt cconeSize, cc, side;
1767: if (PetscAbs(val) < shift) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d on cell %d has an invalid label %d", cone[c], point, val);
1768: if (val > 0) side = 1;
1769: else side = -1;
1770: DMLabelSetValue(label, point, side*(shift+dim));
1771: /* Mark cell faces which touch the fault */
1772: DMPlexGetConeSize(dm, point, &cconeSize);
1773: DMPlexGetCone(dm, point, &ccone);
1774: for (cc = 0; cc < cconeSize; ++cc) {
1775: PetscInt *closure = NULL;
1776: PetscInt closureSize, cl;
1778: DMLabelGetValue(label, ccone[cc], &val);
1779: if (val != -1) continue;
1780: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1781: for (cl = 0; cl < closureSize*2; cl += 2) {
1782: const PetscInt clp = closure[cl];
1784: DMLabelGetValue(label, clp, &val);
1785: if (val == -1) continue;
1786: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1787: break;
1788: }
1789: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1790: }
1791: again = 1;
1792: break;
1793: }
1794: }
1795: }
1796: }
1797: /* Classify the rest by cell membership */
1798: for (s = 0; s < starSize*2; s += 2) {
1799: const PetscInt point = star[s];
1801: DMLabelGetValue(label, point, &val);
1802: if (val == -1) {
1803: PetscInt *sstar = NULL;
1804: PetscInt sstarSize, ss;
1805: PetscBool marked = PETSC_FALSE;
1807: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1808: for (ss = 0; ss < sstarSize*2; ss += 2) {
1809: const PetscInt spoint = sstar[ss];
1811: if ((spoint < cStart) || (spoint >= cMax)) continue;
1812: DMLabelGetValue(label, spoint, &val);
1813: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1814: DMLabelGetValue(depthLabel, point, &dep);
1815: if (val > 0) {
1816: DMLabelSetValue(label, point, shift+dep);
1817: } else {
1818: DMLabelSetValue(label, point, -(shift+dep));
1819: }
1820: marked = PETSC_TRUE;
1821: break;
1822: }
1823: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1824: DMLabelGetValue(depthLabel, point, &dep);
1825: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1826: }
1827: }
1828: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1829: }
1830: ISRestoreIndices(dimIS, &points);
1831: ISDestroy(&dimIS);
1832: /* If any faces touching the fault divide cells on either side, split them */
1833: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1834: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1835: ISExpand(facePosIS, faceNegIS, &dimIS);
1836: ISDestroy(&facePosIS);
1837: ISDestroy(&faceNegIS);
1838: ISGetLocalSize(dimIS, &numPoints);
1839: ISGetIndices(dimIS, &points);
1840: for (p = 0; p < numPoints; ++p) {
1841: const PetscInt point = points[p];
1842: const PetscInt *support;
1843: PetscInt supportSize, valA, valB;
1845: DMPlexGetSupportSize(dm, point, &supportSize);
1846: if (supportSize != 2) continue;
1847: DMPlexGetSupport(dm, point, &support);
1848: DMLabelGetValue(label, support[0], &valA);
1849: DMLabelGetValue(label, support[1], &valB);
1850: if ((valA == -1) || (valB == -1)) continue;
1851: if (valA*valB > 0) continue;
1852: /* Split the face */
1853: DMLabelGetValue(label, point, &valA);
1854: DMLabelClearValue(label, point, valA);
1855: DMLabelSetValue(label, point, dim-1);
1856: /* Label its closure:
1857: unmarked: label as unsplit
1858: incident: relabel as split
1859: split: do nothing
1860: */
1861: {
1862: PetscInt *closure = NULL;
1863: PetscInt closureSize, cl;
1865: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1866: for (cl = 0; cl < closureSize*2; cl += 2) {
1867: DMLabelGetValue(label, closure[cl], &valA);
1868: if (valA == -1) { /* Mark as unsplit */
1869: DMLabelGetValue(depthLabel, closure[cl], &dep);
1870: DMLabelSetValue(label, closure[cl], shift2+dep);
1871: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1872: DMLabelGetValue(depthLabel, closure[cl], &dep);
1873: DMLabelClearValue(label, closure[cl], valA);
1874: DMLabelSetValue(label, closure[cl], dep);
1875: }
1876: }
1877: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1878: }
1879: }
1880: ISRestoreIndices(dimIS, &points);
1881: ISDestroy(&dimIS);
1882: PetscFree(pMax);
1883: return(0);
1884: }
1888: /*@C
1889: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1891: Collective on dm
1893: Input Parameters:
1894: + dm - The original DM
1895: - labelName - The label specifying the interface vertices
1897: Output Parameters:
1898: + hybridLabel - The label fully marking the interface
1899: - dmHybrid - The new DM
1901: Level: developer
1903: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1904: @*/
1905: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1906: {
1907: DM idm;
1908: DMLabel subpointMap, hlabel;
1909: PetscInt dim;
1916: DMGetDimension(dm, &dim);
1917: DMPlexCreateSubmesh(dm, label, 1, &idm);
1918: DMPlexOrient(idm);
1919: DMPlexGetSubpointMap(idm, &subpointMap);
1920: DMLabelDuplicate(subpointMap, &hlabel);
1921: DMLabelClearStratum(hlabel, dim);
1922: DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1923: DMDestroy(&idm);
1924: DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1925: if (hybridLabel) *hybridLabel = hlabel;
1926: else {DMLabelDestroy(&hlabel);}
1927: return(0);
1928: }
1932: /* Here we need the explicit assumption that:
1934: For any marked cell, the marked vertices constitute a single face
1935: */
1936: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1937: {
1938: IS subvertexIS = NULL;
1939: const PetscInt *subvertices;
1940: PetscInt *pStart, *pEnd, *pMax, pSize;
1941: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
1942: PetscErrorCode ierr;
1945: *numFaces = 0;
1946: *nFV = 0;
1947: DMPlexGetDepth(dm, &depth);
1948: DMGetDimension(dm, &dim);
1949: pSize = PetscMax(depth, dim) + 1;
1950: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1951: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1952: for (d = 0; d <= depth; ++d) {
1953: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1954: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1955: }
1956: /* Loop over initial vertices and mark all faces in the collective star() */
1957: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1958: if (subvertexIS) {
1959: ISGetSize(subvertexIS, &numSubVerticesInitial);
1960: ISGetIndices(subvertexIS, &subvertices);
1961: }
1962: for (v = 0; v < numSubVerticesInitial; ++v) {
1963: const PetscInt vertex = subvertices[v];
1964: PetscInt *star = NULL;
1965: PetscInt starSize, s, numCells = 0, c;
1967: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1968: for (s = 0; s < starSize*2; s += 2) {
1969: const PetscInt point = star[s];
1970: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1971: }
1972: for (c = 0; c < numCells; ++c) {
1973: const PetscInt cell = star[c];
1974: PetscInt *closure = NULL;
1975: PetscInt closureSize, cl;
1976: PetscInt cellLoc, numCorners = 0, faceSize = 0;
1978: DMLabelGetValue(subpointMap, cell, &cellLoc);
1979: if (cellLoc == 2) continue;
1980: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1981: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1982: for (cl = 0; cl < closureSize*2; cl += 2) {
1983: const PetscInt point = closure[cl];
1984: PetscInt vertexLoc;
1986: if ((point >= pStart[0]) && (point < pEnd[0])) {
1987: ++numCorners;
1988: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1989: if (vertexLoc == value) closure[faceSize++] = point;
1990: }
1991: }
1992: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1993: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1994: if (faceSize == *nFV) {
1995: const PetscInt *cells = NULL;
1996: PetscInt numCells, nc;
1998: ++(*numFaces);
1999: for (cl = 0; cl < faceSize; ++cl) {
2000: DMLabelSetValue(subpointMap, closure[cl], 0);
2001: }
2002: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2003: for (nc = 0; nc < numCells; ++nc) {
2004: DMLabelSetValue(subpointMap, cells[nc], 2);
2005: }
2006: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2007: }
2008: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2009: }
2010: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2011: }
2012: if (subvertexIS) {
2013: ISRestoreIndices(subvertexIS, &subvertices);
2014: }
2015: ISDestroy(&subvertexIS);
2016: PetscFree3(pStart,pEnd,pMax);
2017: return(0);
2018: }
2022: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2023: {
2024: IS subvertexIS = NULL;
2025: const PetscInt *subvertices;
2026: PetscInt *pStart, *pEnd, *pMax;
2027: PetscInt dim, d, numSubVerticesInitial = 0, v;
2028: PetscErrorCode ierr;
2031: DMGetDimension(dm, &dim);
2032: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2033: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2034: for (d = 0; d <= dim; ++d) {
2035: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2036: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2037: }
2038: /* Loop over initial vertices and mark all faces in the collective star() */
2039: if (vertexLabel) {
2040: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2041: if (subvertexIS) {
2042: ISGetSize(subvertexIS, &numSubVerticesInitial);
2043: ISGetIndices(subvertexIS, &subvertices);
2044: }
2045: }
2046: for (v = 0; v < numSubVerticesInitial; ++v) {
2047: const PetscInt vertex = subvertices[v];
2048: PetscInt *star = NULL;
2049: PetscInt starSize, s, numFaces = 0, f;
2051: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2052: for (s = 0; s < starSize*2; s += 2) {
2053: const PetscInt point = star[s];
2054: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2055: }
2056: for (f = 0; f < numFaces; ++f) {
2057: const PetscInt face = star[f];
2058: PetscInt *closure = NULL;
2059: PetscInt closureSize, c;
2060: PetscInt faceLoc;
2062: DMLabelGetValue(subpointMap, face, &faceLoc);
2063: if (faceLoc == dim-1) continue;
2064: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2065: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2066: for (c = 0; c < closureSize*2; c += 2) {
2067: const PetscInt point = closure[c];
2068: PetscInt vertexLoc;
2070: if ((point >= pStart[0]) && (point < pEnd[0])) {
2071: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2072: if (vertexLoc != value) break;
2073: }
2074: }
2075: if (c == closureSize*2) {
2076: const PetscInt *support;
2077: PetscInt supportSize, s;
2079: for (c = 0; c < closureSize*2; c += 2) {
2080: const PetscInt point = closure[c];
2082: for (d = 0; d < dim; ++d) {
2083: if ((point >= pStart[d]) && (point < pEnd[d])) {
2084: DMLabelSetValue(subpointMap, point, d);
2085: break;
2086: }
2087: }
2088: }
2089: DMPlexGetSupportSize(dm, face, &supportSize);
2090: DMPlexGetSupport(dm, face, &support);
2091: for (s = 0; s < supportSize; ++s) {
2092: DMLabelSetValue(subpointMap, support[s], dim);
2093: }
2094: }
2095: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2096: }
2097: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2098: }
2099: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2100: ISDestroy(&subvertexIS);
2101: PetscFree3(pStart,pEnd,pMax);
2102: return(0);
2103: }
2107: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2108: {
2109: DMLabel label = NULL;
2110: const PetscInt *cone;
2111: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2112: PetscErrorCode ierr;
2115: *numFaces = 0;
2116: *nFV = 0;
2117: if (labelname) {DMGetLabel(dm, labelname, &label);}
2118: *subCells = NULL;
2119: DMGetDimension(dm, &dim);
2120: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2121: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2122: if (cMax < 0) return(0);
2123: if (label) {
2124: for (c = cMax; c < cEnd; ++c) {
2125: PetscInt val;
2127: DMLabelGetValue(label, c, &val);
2128: if (val == value) {
2129: ++(*numFaces);
2130: DMPlexGetConeSize(dm, c, &coneSize);
2131: }
2132: }
2133: } else {
2134: *numFaces = cEnd - cMax;
2135: DMPlexGetConeSize(dm, cMax, &coneSize);
2136: }
2137: PetscMalloc1(*numFaces *2, subCells);
2138: if (!(*numFaces)) return(0);
2139: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2140: for (c = cMax; c < cEnd; ++c) {
2141: const PetscInt *cells;
2142: PetscInt numCells;
2144: if (label) {
2145: PetscInt val;
2147: DMLabelGetValue(label, c, &val);
2148: if (val != value) continue;
2149: }
2150: DMPlexGetCone(dm, c, &cone);
2151: for (p = 0; p < *nFV; ++p) {
2152: DMLabelSetValue(subpointMap, cone[p], 0);
2153: }
2154: /* Negative face */
2155: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2156: /* Not true in parallel
2157: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2158: for (p = 0; p < numCells; ++p) {
2159: DMLabelSetValue(subpointMap, cells[p], 2);
2160: (*subCells)[subc++] = cells[p];
2161: }
2162: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2163: /* Positive face is not included */
2164: }
2165: return(0);
2166: }
2170: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2171: {
2172: PetscInt *pStart, *pEnd;
2173: PetscInt dim, cMax, cEnd, c, d;
2177: DMGetDimension(dm, &dim);
2178: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2179: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2180: if (cMax < 0) return(0);
2181: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2182: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2183: for (c = cMax; c < cEnd; ++c) {
2184: const PetscInt *cone;
2185: PetscInt *closure = NULL;
2186: PetscInt fconeSize, coneSize, closureSize, cl, val;
2188: if (label) {
2189: DMLabelGetValue(label, c, &val);
2190: if (val != value) continue;
2191: }
2192: DMPlexGetConeSize(dm, c, &coneSize);
2193: DMPlexGetCone(dm, c, &cone);
2194: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2195: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2196: /* Negative face */
2197: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2198: for (cl = 0; cl < closureSize*2; cl += 2) {
2199: const PetscInt point = closure[cl];
2201: for (d = 0; d <= dim; ++d) {
2202: if ((point >= pStart[d]) && (point < pEnd[d])) {
2203: DMLabelSetValue(subpointMap, point, d);
2204: break;
2205: }
2206: }
2207: }
2208: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2209: /* Cells -- positive face is not included */
2210: for (cl = 0; cl < 1; ++cl) {
2211: const PetscInt *support;
2212: PetscInt supportSize, s;
2214: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2215: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2216: DMPlexGetSupport(dm, cone[cl], &support);
2217: for (s = 0; s < supportSize; ++s) {
2218: DMLabelSetValue(subpointMap, support[s], dim);
2219: }
2220: }
2221: }
2222: PetscFree2(pStart, pEnd);
2223: return(0);
2224: }
2228: PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2229: {
2230: MPI_Comm comm;
2231: PetscBool posOrient = PETSC_FALSE;
2232: const PetscInt debug = 0;
2233: PetscInt cellDim, faceSize, f;
2237: PetscObjectGetComm((PetscObject)dm,&comm);
2238: DMGetDimension(dm, &cellDim);
2239: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2241: if (cellDim == 1 && numCorners == 2) {
2242: /* Triangle */
2243: faceSize = numCorners-1;
2244: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2245: } else if (cellDim == 2 && numCorners == 3) {
2246: /* Triangle */
2247: faceSize = numCorners-1;
2248: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2249: } else if (cellDim == 3 && numCorners == 4) {
2250: /* Tetrahedron */
2251: faceSize = numCorners-1;
2252: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2253: } else if (cellDim == 1 && numCorners == 3) {
2254: /* Quadratic line */
2255: faceSize = 1;
2256: posOrient = PETSC_TRUE;
2257: } else if (cellDim == 2 && numCorners == 4) {
2258: /* Quads */
2259: faceSize = 2;
2260: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2261: posOrient = PETSC_TRUE;
2262: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2263: posOrient = PETSC_TRUE;
2264: } else {
2265: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2266: posOrient = PETSC_FALSE;
2267: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2268: }
2269: } else if (cellDim == 2 && numCorners == 6) {
2270: /* Quadratic triangle (I hate this) */
2271: /* Edges are determined by the first 2 vertices (corners of edges) */
2272: const PetscInt faceSizeTri = 3;
2273: PetscInt sortedIndices[3], i, iFace;
2274: PetscBool found = PETSC_FALSE;
2275: PetscInt faceVerticesTriSorted[9] = {
2276: 0, 3, 4, /* bottom */
2277: 1, 4, 5, /* right */
2278: 2, 3, 5, /* left */
2279: };
2280: PetscInt faceVerticesTri[9] = {
2281: 0, 3, 4, /* bottom */
2282: 1, 4, 5, /* right */
2283: 2, 5, 3, /* left */
2284: };
2286: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2287: PetscSortInt(faceSizeTri, sortedIndices);
2288: for (iFace = 0; iFace < 3; ++iFace) {
2289: const PetscInt ii = iFace*faceSizeTri;
2290: PetscInt fVertex, cVertex;
2292: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2293: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2294: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2295: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2296: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2297: faceVertices[fVertex] = origVertices[cVertex];
2298: break;
2299: }
2300: }
2301: }
2302: found = PETSC_TRUE;
2303: break;
2304: }
2305: }
2306: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2307: if (posOriented) *posOriented = PETSC_TRUE;
2308: return(0);
2309: } else if (cellDim == 2 && numCorners == 9) {
2310: /* Quadratic quad (I hate this) */
2311: /* Edges are determined by the first 2 vertices (corners of edges) */
2312: const PetscInt faceSizeQuad = 3;
2313: PetscInt sortedIndices[3], i, iFace;
2314: PetscBool found = PETSC_FALSE;
2315: PetscInt faceVerticesQuadSorted[12] = {
2316: 0, 1, 4, /* bottom */
2317: 1, 2, 5, /* right */
2318: 2, 3, 6, /* top */
2319: 0, 3, 7, /* left */
2320: };
2321: PetscInt faceVerticesQuad[12] = {
2322: 0, 1, 4, /* bottom */
2323: 1, 2, 5, /* right */
2324: 2, 3, 6, /* top */
2325: 3, 0, 7, /* left */
2326: };
2328: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2329: PetscSortInt(faceSizeQuad, sortedIndices);
2330: for (iFace = 0; iFace < 4; ++iFace) {
2331: const PetscInt ii = iFace*faceSizeQuad;
2332: PetscInt fVertex, cVertex;
2334: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2335: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2336: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2337: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2338: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2339: faceVertices[fVertex] = origVertices[cVertex];
2340: break;
2341: }
2342: }
2343: }
2344: found = PETSC_TRUE;
2345: break;
2346: }
2347: }
2348: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2349: if (posOriented) *posOriented = PETSC_TRUE;
2350: return(0);
2351: } else if (cellDim == 3 && numCorners == 8) {
2352: /* Hexes
2353: A hex is two oriented quads with the normal of the first
2354: pointing up at the second.
2356: 7---6
2357: /| /|
2358: 4---5 |
2359: | 1-|-2
2360: |/ |/
2361: 0---3
2363: Faces are determined by the first 4 vertices (corners of faces) */
2364: const PetscInt faceSizeHex = 4;
2365: PetscInt sortedIndices[4], i, iFace;
2366: PetscBool found = PETSC_FALSE;
2367: PetscInt faceVerticesHexSorted[24] = {
2368: 0, 1, 2, 3, /* bottom */
2369: 4, 5, 6, 7, /* top */
2370: 0, 3, 4, 5, /* front */
2371: 2, 3, 5, 6, /* right */
2372: 1, 2, 6, 7, /* back */
2373: 0, 1, 4, 7, /* left */
2374: };
2375: PetscInt faceVerticesHex[24] = {
2376: 1, 2, 3, 0, /* bottom */
2377: 4, 5, 6, 7, /* top */
2378: 0, 3, 5, 4, /* front */
2379: 3, 2, 6, 5, /* right */
2380: 2, 1, 7, 6, /* back */
2381: 1, 0, 4, 7, /* left */
2382: };
2384: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2385: PetscSortInt(faceSizeHex, sortedIndices);
2386: for (iFace = 0; iFace < 6; ++iFace) {
2387: const PetscInt ii = iFace*faceSizeHex;
2388: PetscInt fVertex, cVertex;
2390: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2391: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2392: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2393: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2394: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2395: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2396: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2397: faceVertices[fVertex] = origVertices[cVertex];
2398: break;
2399: }
2400: }
2401: }
2402: found = PETSC_TRUE;
2403: break;
2404: }
2405: }
2406: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2407: if (posOriented) *posOriented = PETSC_TRUE;
2408: return(0);
2409: } else if (cellDim == 3 && numCorners == 10) {
2410: /* Quadratic tet */
2411: /* Faces are determined by the first 3 vertices (corners of faces) */
2412: const PetscInt faceSizeTet = 6;
2413: PetscInt sortedIndices[6], i, iFace;
2414: PetscBool found = PETSC_FALSE;
2415: PetscInt faceVerticesTetSorted[24] = {
2416: 0, 1, 2, 6, 7, 8, /* bottom */
2417: 0, 3, 4, 6, 7, 9, /* front */
2418: 1, 4, 5, 7, 8, 9, /* right */
2419: 2, 3, 5, 6, 8, 9, /* left */
2420: };
2421: PetscInt faceVerticesTet[24] = {
2422: 0, 1, 2, 6, 7, 8, /* bottom */
2423: 0, 4, 3, 6, 7, 9, /* front */
2424: 1, 5, 4, 7, 8, 9, /* right */
2425: 2, 3, 5, 8, 6, 9, /* left */
2426: };
2428: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2429: PetscSortInt(faceSizeTet, sortedIndices);
2430: for (iFace=0; iFace < 4; ++iFace) {
2431: const PetscInt ii = iFace*faceSizeTet;
2432: PetscInt fVertex, cVertex;
2434: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2435: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2436: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2437: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2438: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2439: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2440: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2441: faceVertices[fVertex] = origVertices[cVertex];
2442: break;
2443: }
2444: }
2445: }
2446: found = PETSC_TRUE;
2447: break;
2448: }
2449: }
2450: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2451: if (posOriented) *posOriented = PETSC_TRUE;
2452: return(0);
2453: } else if (cellDim == 3 && numCorners == 27) {
2454: /* Quadratic hexes (I hate this)
2455: A hex is two oriented quads with the normal of the first
2456: pointing up at the second.
2458: 7---6
2459: /| /|
2460: 4---5 |
2461: | 3-|-2
2462: |/ |/
2463: 0---1
2465: Faces are determined by the first 4 vertices (corners of faces) */
2466: const PetscInt faceSizeQuadHex = 9;
2467: PetscInt sortedIndices[9], i, iFace;
2468: PetscBool found = PETSC_FALSE;
2469: PetscInt faceVerticesQuadHexSorted[54] = {
2470: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2471: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2472: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2473: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2474: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2475: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2476: };
2477: PetscInt faceVerticesQuadHex[54] = {
2478: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2479: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2480: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2481: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2482: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2483: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2484: };
2486: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2487: PetscSortInt(faceSizeQuadHex, sortedIndices);
2488: for (iFace = 0; iFace < 6; ++iFace) {
2489: const PetscInt ii = iFace*faceSizeQuadHex;
2490: PetscInt fVertex, cVertex;
2492: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2493: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2494: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2495: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2496: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2497: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2498: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2499: faceVertices[fVertex] = origVertices[cVertex];
2500: break;
2501: }
2502: }
2503: }
2504: found = PETSC_TRUE;
2505: break;
2506: }
2507: }
2508: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2509: if (posOriented) *posOriented = PETSC_TRUE;
2510: return(0);
2511: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2512: if (!posOrient) {
2513: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2514: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2515: } else {
2516: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2517: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2518: }
2519: if (posOriented) *posOriented = posOrient;
2520: return(0);
2521: }
2525: /*
2526: Given a cell and a face, as a set of vertices,
2527: return the oriented face, as a set of vertices, in faceVertices
2528: The orientation is such that the face normal points out of the cell
2529: */
2530: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2531: {
2532: const PetscInt *cone = NULL;
2533: PetscInt coneSize, v, f, v2;
2534: PetscInt oppositeVertex = -1;
2535: PetscErrorCode ierr;
2538: DMPlexGetConeSize(dm, cell, &coneSize);
2539: DMPlexGetCone(dm, cell, &cone);
2540: for (v = 0, v2 = 0; v < coneSize; ++v) {
2541: PetscBool found = PETSC_FALSE;
2543: for (f = 0; f < faceSize; ++f) {
2544: if (face[f] == cone[v]) {
2545: found = PETSC_TRUE; break;
2546: }
2547: }
2548: if (found) {
2549: indices[v2] = v;
2550: origVertices[v2] = cone[v];
2551: ++v2;
2552: } else {
2553: oppositeVertex = v;
2554: }
2555: }
2556: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2557: return(0);
2558: }
2562: /*
2563: DMPlexInsertFace_Internal - Puts a face into the mesh
2565: Not collective
2567: Input Parameters:
2568: + dm - The DMPlex
2569: . numFaceVertex - The number of vertices in the face
2570: . faceVertices - The vertices in the face for dm
2571: . subfaceVertices - The vertices in the face for subdm
2572: . numCorners - The number of vertices in the cell
2573: . cell - A cell in dm containing the face
2574: . subcell - A cell in subdm containing the face
2575: . firstFace - First face in the mesh
2576: - newFacePoint - Next face in the mesh
2578: Output Parameters:
2579: . newFacePoint - Contains next face point number on input, updated on output
2581: Level: developer
2582: */
2583: 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)
2584: {
2585: MPI_Comm comm;
2586: DM_Plex *submesh = (DM_Plex*) subdm->data;
2587: const PetscInt *faces;
2588: PetscInt numFaces, coneSize;
2589: PetscErrorCode ierr;
2592: PetscObjectGetComm((PetscObject)dm,&comm);
2593: DMPlexGetConeSize(subdm, subcell, &coneSize);
2594: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2595: #if 0
2596: /* Cannot use this because support() has not been constructed yet */
2597: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2598: #else
2599: {
2600: PetscInt f;
2602: numFaces = 0;
2603: DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2604: for (f = firstFace; f < *newFacePoint; ++f) {
2605: PetscInt dof, off, d;
2607: PetscSectionGetDof(submesh->coneSection, f, &dof);
2608: PetscSectionGetOffset(submesh->coneSection, f, &off);
2609: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2610: for (d = 0; d < dof; ++d) {
2611: const PetscInt p = submesh->cones[off+d];
2612: PetscInt v;
2614: for (v = 0; v < numFaceVertices; ++v) {
2615: if (subfaceVertices[v] == p) break;
2616: }
2617: if (v == numFaceVertices) break;
2618: }
2619: if (d == dof) {
2620: numFaces = 1;
2621: ((PetscInt*) faces)[0] = f;
2622: }
2623: }
2624: }
2625: #endif
2626: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2627: else if (numFaces == 1) {
2628: /* Add the other cell neighbor for this face */
2629: DMPlexSetCone(subdm, subcell, faces);
2630: } else {
2631: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2632: PetscBool posOriented;
2634: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2635: origVertices = &orientedVertices[numFaceVertices];
2636: indices = &orientedVertices[numFaceVertices*2];
2637: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2638: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2639: /* TODO: I know that routine should return a permutation, not the indices */
2640: for (v = 0; v < numFaceVertices; ++v) {
2641: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2642: for (ov = 0; ov < numFaceVertices; ++ov) {
2643: if (orientedVertices[ov] == vertex) {
2644: orientedSubVertices[ov] = subvertex;
2645: break;
2646: }
2647: }
2648: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2649: }
2650: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2651: DMPlexSetCone(subdm, subcell, newFacePoint);
2652: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);
2653: ++(*newFacePoint);
2654: }
2655: #if 0
2656: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2657: #else
2658: DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);
2659: #endif
2660: return(0);
2661: }
2665: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2666: {
2667: MPI_Comm comm;
2668: DMLabel subpointMap;
2669: IS subvertexIS, subcellIS;
2670: const PetscInt *subVertices, *subCells;
2671: PetscInt numSubVertices, firstSubVertex, numSubCells;
2672: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2673: PetscInt vStart, vEnd, c, f;
2674: PetscErrorCode ierr;
2677: PetscObjectGetComm((PetscObject)dm,&comm);
2678: /* Create subpointMap which marks the submesh */
2679: DMLabelCreate("subpoint_map", &subpointMap);
2680: DMPlexSetSubpointMap(subdm, subpointMap);
2681: DMLabelDestroy(&subpointMap);
2682: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2683: /* Setup chart */
2684: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2685: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2686: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2687: DMPlexSetVTKCellHeight(subdm, 1);
2688: /* Set cone sizes */
2689: firstSubVertex = numSubCells;
2690: firstSubFace = numSubCells+numSubVertices;
2691: newFacePoint = firstSubFace;
2692: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2693: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2694: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2695: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2696: for (c = 0; c < numSubCells; ++c) {
2697: DMPlexSetConeSize(subdm, c, 1);
2698: }
2699: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2700: DMPlexSetConeSize(subdm, f, nFV);
2701: }
2702: DMSetUp(subdm);
2703: /* Create face cones */
2704: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2705: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2706: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2707: for (c = 0; c < numSubCells; ++c) {
2708: const PetscInt cell = subCells[c];
2709: const PetscInt subcell = c;
2710: PetscInt *closure = NULL;
2711: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2713: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2714: for (cl = 0; cl < closureSize*2; cl += 2) {
2715: const PetscInt point = closure[cl];
2716: PetscInt subVertex;
2718: if ((point >= vStart) && (point < vEnd)) {
2719: ++numCorners;
2720: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2721: if (subVertex >= 0) {
2722: closure[faceSize] = point;
2723: subface[faceSize] = firstSubVertex+subVertex;
2724: ++faceSize;
2725: }
2726: }
2727: }
2728: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2729: if (faceSize == nFV) {
2730: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2731: }
2732: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2733: }
2734: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
2735: DMPlexSymmetrize(subdm);
2736: DMPlexStratify(subdm);
2737: /* Build coordinates */
2738: {
2739: PetscSection coordSection, subCoordSection;
2740: Vec coordinates, subCoordinates;
2741: PetscScalar *coords, *subCoords;
2742: PetscInt numComp, coordSize, v;
2743: const char *name;
2745: DMGetCoordinateSection(dm, &coordSection);
2746: DMGetCoordinatesLocal(dm, &coordinates);
2747: DMGetCoordinateSection(subdm, &subCoordSection);
2748: PetscSectionSetNumFields(subCoordSection, 1);
2749: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2750: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2751: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2752: for (v = 0; v < numSubVertices; ++v) {
2753: const PetscInt vertex = subVertices[v];
2754: const PetscInt subvertex = firstSubVertex+v;
2755: PetscInt dof;
2757: PetscSectionGetDof(coordSection, vertex, &dof);
2758: PetscSectionSetDof(subCoordSection, subvertex, dof);
2759: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2760: }
2761: PetscSectionSetUp(subCoordSection);
2762: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2763: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2764: PetscObjectGetName((PetscObject)coordinates,&name);
2765: PetscObjectSetName((PetscObject)subCoordinates,name);
2766: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2767: VecSetType(subCoordinates,VECSTANDARD);
2768: if (coordSize) {
2769: VecGetArray(coordinates, &coords);
2770: VecGetArray(subCoordinates, &subCoords);
2771: for (v = 0; v < numSubVertices; ++v) {
2772: const PetscInt vertex = subVertices[v];
2773: const PetscInt subvertex = firstSubVertex+v;
2774: PetscInt dof, off, sdof, soff, d;
2776: PetscSectionGetDof(coordSection, vertex, &dof);
2777: PetscSectionGetOffset(coordSection, vertex, &off);
2778: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2779: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2780: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2781: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2782: }
2783: VecRestoreArray(coordinates, &coords);
2784: VecRestoreArray(subCoordinates, &subCoords);
2785: }
2786: DMSetCoordinatesLocal(subdm, subCoordinates);
2787: VecDestroy(&subCoordinates);
2788: }
2789: /* Cleanup */
2790: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2791: ISDestroy(&subvertexIS);
2792: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2793: ISDestroy(&subcellIS);
2794: return(0);
2795: }
2799: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2800: {
2801: PetscInt subPoint;
2804: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2805: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2806: }
2810: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2811: {
2812: MPI_Comm comm;
2813: DMLabel subpointMap;
2814: IS *subpointIS;
2815: const PetscInt **subpoints;
2816: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2817: PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v;
2818: PetscErrorCode ierr;
2821: PetscObjectGetComm((PetscObject)dm,&comm);
2822: /* Create subpointMap which marks the submesh */
2823: DMLabelCreate("subpoint_map", &subpointMap);
2824: DMPlexSetSubpointMap(subdm, subpointMap);
2825: if (cellHeight) {
2826: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2827: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2828: } else {
2829: DMLabel depth;
2830: IS pointIS;
2831: const PetscInt *points;
2832: PetscInt numPoints;
2834: DMPlexGetDepthLabel(dm, &depth);
2835: DMLabelGetStratumSize(label, value, &numPoints);
2836: DMLabelGetStratumIS(label, value, &pointIS);
2837: ISGetIndices(pointIS, &points);
2838: for (p = 0; p < numPoints; ++p) {
2839: PetscInt *closure = NULL;
2840: PetscInt closureSize, c, pdim;
2842: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2843: for (c = 0; c < closureSize*2; c += 2) {
2844: DMLabelGetValue(depth, closure[c], &pdim);
2845: DMLabelSetValue(subpointMap, closure[c], pdim);
2846: }
2847: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2848: }
2849: ISRestoreIndices(pointIS, &points);
2850: ISDestroy(&pointIS);
2851: }
2852: DMLabelDestroy(&subpointMap);
2853: /* Setup chart */
2854: DMGetDimension(dm, &dim);
2855: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2856: for (d = 0; d <= dim; ++d) {
2857: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2858: totSubPoints += numSubPoints[d];
2859: }
2860: DMPlexSetChart(subdm, 0, totSubPoints);
2861: DMPlexSetVTKCellHeight(subdm, cellHeight);
2862: /* Set cone sizes */
2863: firstSubPoint[dim] = 0;
2864: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2865: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2866: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2867: for (d = 0; d <= dim; ++d) {
2868: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2869: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2870: }
2871: for (d = 0; d <= dim; ++d) {
2872: for (p = 0; p < numSubPoints[d]; ++p) {
2873: const PetscInt point = subpoints[d][p];
2874: const PetscInt subpoint = firstSubPoint[d] + p;
2875: const PetscInt *cone;
2876: PetscInt coneSize, coneSizeNew, c, val;
2878: DMPlexGetConeSize(dm, point, &coneSize);
2879: DMPlexSetConeSize(subdm, subpoint, coneSize);
2880: if (cellHeight && (d == dim)) {
2881: DMPlexGetCone(dm, point, &cone);
2882: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2883: DMLabelGetValue(subpointMap, cone[c], &val);
2884: if (val >= 0) coneSizeNew++;
2885: }
2886: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2887: }
2888: }
2889: }
2890: DMSetUp(subdm);
2891: /* Set cones */
2892: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2893: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2894: for (d = 0; d <= dim; ++d) {
2895: for (p = 0; p < numSubPoints[d]; ++p) {
2896: const PetscInt point = subpoints[d][p];
2897: const PetscInt subpoint = firstSubPoint[d] + p;
2898: const PetscInt *cone, *ornt;
2899: PetscInt coneSize, subconeSize, coneSizeNew, c, subc;
2901: DMPlexGetConeSize(dm, point, &coneSize);
2902: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2903: DMPlexGetCone(dm, point, &cone);
2904: DMPlexGetConeOrientation(dm, point, &ornt);
2905: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2906: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2907: if (subc >= 0) {
2908: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2909: orntNew[coneSizeNew] = ornt[c];
2910: ++coneSizeNew;
2911: }
2912: }
2913: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2914: DMPlexSetCone(subdm, subpoint, coneNew);
2915: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2916: }
2917: }
2918: PetscFree2(coneNew,orntNew);
2919: DMPlexSymmetrize(subdm);
2920: DMPlexStratify(subdm);
2921: /* Build coordinates */
2922: {
2923: PetscSection coordSection, subCoordSection;
2924: Vec coordinates, subCoordinates;
2925: PetscScalar *coords, *subCoords;
2926: PetscInt numComp, coordSize;
2927: const char *name;
2929: DMGetCoordinateSection(dm, &coordSection);
2930: DMGetCoordinatesLocal(dm, &coordinates);
2931: DMGetCoordinateSection(subdm, &subCoordSection);
2932: PetscSectionSetNumFields(subCoordSection, 1);
2933: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2934: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2935: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2936: for (v = 0; v < numSubPoints[0]; ++v) {
2937: const PetscInt vertex = subpoints[0][v];
2938: const PetscInt subvertex = firstSubPoint[0]+v;
2939: PetscInt dof;
2941: PetscSectionGetDof(coordSection, vertex, &dof);
2942: PetscSectionSetDof(subCoordSection, subvertex, dof);
2943: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2944: }
2945: PetscSectionSetUp(subCoordSection);
2946: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2947: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2948: PetscObjectGetName((PetscObject)coordinates,&name);
2949: PetscObjectSetName((PetscObject)subCoordinates,name);
2950: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2951: VecSetType(subCoordinates,VECSTANDARD);
2952: VecGetArray(coordinates, &coords);
2953: VecGetArray(subCoordinates, &subCoords);
2954: for (v = 0; v < numSubPoints[0]; ++v) {
2955: const PetscInt vertex = subpoints[0][v];
2956: const PetscInt subvertex = firstSubPoint[0]+v;
2957: PetscInt dof, off, sdof, soff, d;
2959: PetscSectionGetDof(coordSection, vertex, &dof);
2960: PetscSectionGetOffset(coordSection, vertex, &off);
2961: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2962: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2963: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2964: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2965: }
2966: VecRestoreArray(coordinates, &coords);
2967: VecRestoreArray(subCoordinates, &subCoords);
2968: DMSetCoordinatesLocal(subdm, subCoordinates);
2969: VecDestroy(&subCoordinates);
2970: }
2971: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
2972: {
2973: PetscSF sfPoint, sfPointSub;
2974: IS subpIS;
2975: const PetscSFNode *remotePoints;
2976: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
2977: const PetscInt *localPoints, *subpoints;
2978: PetscInt *slocalPoints;
2979: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
2980: PetscMPIInt rank;
2982: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2983: DMGetPointSF(dm, &sfPoint);
2984: DMGetPointSF(subdm, &sfPointSub);
2985: DMPlexGetChart(dm, &pStart, &pEnd);
2986: DMPlexGetChart(subdm, NULL, &numSubroots);
2987: DMPlexCreateSubpointIS(subdm, &subpIS);
2988: if (subpIS) {
2989: ISGetIndices(subpIS, &subpoints);
2990: ISGetLocalSize(subpIS, &numSubpoints);
2991: }
2992: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
2993: if (numRoots >= 0) {
2994: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
2995: for (p = 0; p < pEnd-pStart; ++p) {
2996: newLocalPoints[p].rank = -2;
2997: newLocalPoints[p].index = -2;
2998: }
2999: /* Set subleaves */
3000: for (l = 0; l < numLeaves; ++l) {
3001: const PetscInt point = localPoints[l];
3002: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3004: if (subpoint < 0) continue;
3005: newLocalPoints[point-pStart].rank = rank;
3006: newLocalPoints[point-pStart].index = subpoint;
3007: ++numSubleaves;
3008: }
3009: /* Must put in owned subpoints */
3010: for (p = pStart; p < pEnd; ++p) {
3011: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3013: if (subpoint < 0) {
3014: newOwners[p-pStart].rank = -3;
3015: newOwners[p-pStart].index = -3;
3016: } else {
3017: newOwners[p-pStart].rank = rank;
3018: newOwners[p-pStart].index = subpoint;
3019: }
3020: }
3021: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3022: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3023: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3024: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3025: PetscMalloc1(numSubleaves, &slocalPoints);
3026: PetscMalloc1(numSubleaves, &sremotePoints);
3027: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3028: const PetscInt point = localPoints[l];
3029: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3031: if (subpoint < 0) continue;
3032: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3033: slocalPoints[sl] = subpoint;
3034: sremotePoints[sl].rank = newLocalPoints[point].rank;
3035: sremotePoints[sl].index = newLocalPoints[point].index;
3036: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3037: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3038: ++sl;
3039: }
3040: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3041: PetscFree2(newLocalPoints,newOwners);
3042: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3043: }
3044: if (subpIS) {
3045: ISRestoreIndices(subpIS, &subpoints);
3046: ISDestroy(&subpIS);
3047: }
3048: }
3049: /* Cleanup */
3050: for (d = 0; d <= dim; ++d) {
3051: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3052: ISDestroy(&subpointIS[d]);
3053: }
3054: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3055: return(0);
3056: }
3060: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3061: {
3065: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3066: return(0);
3067: }
3071: /*@
3072: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3074: Input Parameters:
3075: + dm - The original mesh
3076: . vertexLabel - The DMLabel marking vertices contained in the surface
3077: - value - The label value to use
3079: Output Parameter:
3080: . subdm - The surface mesh
3082: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3084: Level: developer
3086: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3087: @*/
3088: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3089: {
3090: PetscInt dim, depth;
3096: DMGetDimension(dm, &dim);
3097: DMPlexGetDepth(dm, &depth);
3098: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3099: DMSetType(*subdm, DMPLEX);
3100: DMSetDimension(*subdm, dim-1);
3101: if (depth == dim) {
3102: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3103: } else {
3104: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3105: }
3106: return(0);
3107: }
3111: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3112: {
3113: MPI_Comm comm;
3114: DMLabel subpointMap;
3115: IS subvertexIS;
3116: const PetscInt *subVertices;
3117: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3118: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3119: PetscInt cMax, c, f;
3120: PetscErrorCode ierr;
3123: PetscObjectGetComm((PetscObject)dm, &comm);
3124: /* Create subpointMap which marks the submesh */
3125: DMLabelCreate("subpoint_map", &subpointMap);
3126: DMPlexSetSubpointMap(subdm, subpointMap);
3127: DMLabelDestroy(&subpointMap);
3128: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3129: /* Setup chart */
3130: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3131: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3132: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3133: DMPlexSetVTKCellHeight(subdm, 1);
3134: /* Set cone sizes */
3135: firstSubVertex = numSubCells;
3136: firstSubFace = numSubCells+numSubVertices;
3137: newFacePoint = firstSubFace;
3138: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3139: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3140: for (c = 0; c < numSubCells; ++c) {
3141: DMPlexSetConeSize(subdm, c, 1);
3142: }
3143: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3144: DMPlexSetConeSize(subdm, f, nFV);
3145: }
3146: DMSetUp(subdm);
3147: /* Create face cones */
3148: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3149: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3150: DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3151: for (c = 0; c < numSubCells; ++c) {
3152: const PetscInt cell = subCells[c];
3153: const PetscInt subcell = c;
3154: const PetscInt *cone, *cells;
3155: PetscInt numCells, subVertex, p, v;
3157: if (cell < cMax) continue;
3158: DMPlexGetCone(dm, cell, &cone);
3159: for (v = 0; v < nFV; ++v) {
3160: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3161: subface[v] = firstSubVertex+subVertex;
3162: }
3163: DMPlexSetCone(subdm, newFacePoint, subface);
3164: DMPlexSetCone(subdm, subcell, &newFacePoint);
3165: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3166: /* Not true in parallel
3167: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3168: for (p = 0; p < numCells; ++p) {
3169: PetscInt negsubcell;
3171: if (cells[p] >= cMax) continue;
3172: /* I know this is a crap search */
3173: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3174: if (subCells[negsubcell] == cells[p]) break;
3175: }
3176: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3177: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3178: }
3179: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3180: ++newFacePoint;
3181: }
3182: DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);
3183: DMPlexSymmetrize(subdm);
3184: DMPlexStratify(subdm);
3185: /* Build coordinates */
3186: {
3187: PetscSection coordSection, subCoordSection;
3188: Vec coordinates, subCoordinates;
3189: PetscScalar *coords, *subCoords;
3190: PetscInt numComp, coordSize, v;
3191: const char *name;
3193: DMGetCoordinateSection(dm, &coordSection);
3194: DMGetCoordinatesLocal(dm, &coordinates);
3195: DMGetCoordinateSection(subdm, &subCoordSection);
3196: PetscSectionSetNumFields(subCoordSection, 1);
3197: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3198: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3199: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3200: for (v = 0; v < numSubVertices; ++v) {
3201: const PetscInt vertex = subVertices[v];
3202: const PetscInt subvertex = firstSubVertex+v;
3203: PetscInt dof;
3205: PetscSectionGetDof(coordSection, vertex, &dof);
3206: PetscSectionSetDof(subCoordSection, subvertex, dof);
3207: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3208: }
3209: PetscSectionSetUp(subCoordSection);
3210: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3211: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3212: PetscObjectGetName((PetscObject)coordinates,&name);
3213: PetscObjectSetName((PetscObject)subCoordinates,name);
3214: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3215: VecSetType(subCoordinates,VECSTANDARD);
3216: VecGetArray(coordinates, &coords);
3217: VecGetArray(subCoordinates, &subCoords);
3218: for (v = 0; v < numSubVertices; ++v) {
3219: const PetscInt vertex = subVertices[v];
3220: const PetscInt subvertex = firstSubVertex+v;
3221: PetscInt dof, off, sdof, soff, d;
3223: PetscSectionGetDof(coordSection, vertex, &dof);
3224: PetscSectionGetOffset(coordSection, vertex, &off);
3225: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3226: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3227: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3228: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3229: }
3230: VecRestoreArray(coordinates, &coords);
3231: VecRestoreArray(subCoordinates, &subCoords);
3232: DMSetCoordinatesLocal(subdm, subCoordinates);
3233: VecDestroy(&subCoordinates);
3234: }
3235: /* Build SF */
3236: CHKMEMQ;
3237: {
3238: PetscSF sfPoint, sfPointSub;
3239: const PetscSFNode *remotePoints;
3240: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3241: const PetscInt *localPoints;
3242: PetscInt *slocalPoints;
3243: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3244: PetscMPIInt rank;
3246: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3247: DMGetPointSF(dm, &sfPoint);
3248: DMGetPointSF(subdm, &sfPointSub);
3249: DMPlexGetChart(dm, &pStart, &pEnd);
3250: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3251: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3252: if (numRoots >= 0) {
3253: /* Only vertices should be shared */
3254: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3255: for (p = 0; p < pEnd-pStart; ++p) {
3256: newLocalPoints[p].rank = -2;
3257: newLocalPoints[p].index = -2;
3258: }
3259: /* Set subleaves */
3260: for (l = 0; l < numLeaves; ++l) {
3261: const PetscInt point = localPoints[l];
3262: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3264: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3265: if (subPoint < 0) continue;
3266: newLocalPoints[point-pStart].rank = rank;
3267: newLocalPoints[point-pStart].index = subPoint;
3268: ++numSubLeaves;
3269: }
3270: /* Must put in owned subpoints */
3271: for (p = pStart; p < pEnd; ++p) {
3272: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3274: if (subPoint < 0) {
3275: newOwners[p-pStart].rank = -3;
3276: newOwners[p-pStart].index = -3;
3277: } else {
3278: newOwners[p-pStart].rank = rank;
3279: newOwners[p-pStart].index = subPoint;
3280: }
3281: }
3282: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3283: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3284: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3285: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3286: PetscMalloc1(numSubLeaves, &slocalPoints);
3287: PetscMalloc1(numSubLeaves, &sremotePoints);
3288: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3289: const PetscInt point = localPoints[l];
3290: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3292: if (subPoint < 0) continue;
3293: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3294: slocalPoints[sl] = subPoint;
3295: sremotePoints[sl].rank = newLocalPoints[point].rank;
3296: sremotePoints[sl].index = newLocalPoints[point].index;
3297: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3298: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3299: ++sl;
3300: }
3301: PetscFree2(newLocalPoints,newOwners);
3302: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3303: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3304: }
3305: }
3306: CHKMEMQ;
3307: /* Cleanup */
3308: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3309: ISDestroy(&subvertexIS);
3310: PetscFree(subCells);
3311: return(0);
3312: }
3316: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3317: {
3318: DMLabel label = NULL;
3322: if (labelname) {DMGetLabel(dm, labelname, &label);}
3323: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3324: return(0);
3325: }
3329: /*
3330: 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.
3332: Input Parameters:
3333: + dm - The original mesh
3334: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3335: . label - A label name, or NULL
3336: - value - A label value
3338: Output Parameter:
3339: . subdm - The surface mesh
3341: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3343: Level: developer
3345: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3346: */
3347: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3348: {
3349: PetscInt dim, depth;
3355: DMGetDimension(dm, &dim);
3356: DMPlexGetDepth(dm, &depth);
3357: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3358: DMSetType(*subdm, DMPLEX);
3359: DMSetDimension(*subdm, dim-1);
3360: if (depth == dim) {
3361: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3362: } else {
3363: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3364: }
3365: return(0);
3366: }
3370: /*@
3371: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3373: Input Parameters:
3374: + dm - The original mesh
3375: . cellLabel - The DMLabel marking cells contained in the new mesh
3376: - value - The label value to use
3378: Output Parameter:
3379: . subdm - The new mesh
3381: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3383: Level: developer
3385: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3386: @*/
3387: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3388: {
3389: PetscInt dim;
3395: DMGetDimension(dm, &dim);
3396: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3397: DMSetType(*subdm, DMPLEX);
3398: DMSetDimension(*subdm, dim);
3399: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3400: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3401: return(0);
3402: }
3406: /*@
3407: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3409: Input Parameter:
3410: . dm - The submesh DM
3412: Output Parameter:
3413: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3415: Level: developer
3417: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3418: @*/
3419: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3420: {
3424: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3425: return(0);
3426: }
3430: /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */
3431: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3432: {
3433: DM_Plex *mesh = (DM_Plex *) dm->data;
3434: DMLabel tmp;
3439: tmp = mesh->subpointMap;
3440: mesh->subpointMap = subpointMap;
3441: ++mesh->subpointMap->refct;
3442: DMLabelDestroy(&tmp);
3443: return(0);
3444: }
3448: /*@
3449: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3451: Input Parameter:
3452: . dm - The submesh DM
3454: Output Parameter:
3455: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3457: Note: This IS is guaranteed to be sorted by the construction of the submesh
3459: Level: developer
3461: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3462: @*/
3463: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3464: {
3465: MPI_Comm comm;
3466: DMLabel subpointMap;
3467: IS is;
3468: const PetscInt *opoints;
3469: PetscInt *points, *depths;
3470: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3471: PetscErrorCode ierr;
3476: PetscObjectGetComm((PetscObject)dm,&comm);
3477: *subpointIS = NULL;
3478: DMPlexGetSubpointMap(dm, &subpointMap);
3479: DMPlexGetDepth(dm, &depth);
3480: if (subpointMap && depth >= 0) {
3481: DMPlexGetChart(dm, &pStart, &pEnd);
3482: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3483: DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);
3484: depths[0] = depth;
3485: depths[1] = 0;
3486: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3487: PetscMalloc1(pEnd, &points);
3488: for(d = 0, off = 0; d <= depth; ++d) {
3489: const PetscInt dep = depths[d];
3491: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3492: DMLabelGetStratumSize(subpointMap, dep, &n);
3493: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3494: 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);
3495: } else {
3496: if (!n) {
3497: if (d == 0) {
3498: /* Missing cells */
3499: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3500: } else {
3501: /* Missing faces */
3502: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3503: }
3504: }
3505: }
3506: if (n) {
3507: DMLabelGetStratumIS(subpointMap, dep, &is);
3508: ISGetIndices(is, &opoints);
3509: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3510: ISRestoreIndices(is, &opoints);
3511: ISDestroy(&is);
3512: }
3513: }
3514: DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);
3515: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3516: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3517: }
3518: return(0);
3519: }