Actual source code: plexsubmesh.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petscsf.h>
5: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
6: {
7: PetscInt fStart, fEnd, f;
12: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
13: for (f = fStart; f < fEnd; ++f) {
14: PetscInt supportSize;
16: DMPlexGetSupportSize(dm, f, &supportSize);
17: if (supportSize == 1) {
18: if (val < 0) {
19: PetscInt *closure = NULL;
20: PetscInt clSize, cl, cval;
22: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
23: for (cl = 0; cl < clSize*2; cl += 2) {
24: DMLabelGetValue(label, closure[cl], &cval);
25: if (cval < 0) continue;
26: DMLabelSetValue(label, f, cval);
27: break;
28: }
29: if (cl == clSize*2) {DMLabelSetValue(label, f, 1);}
30: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
31: } else {
32: DMLabelSetValue(label, f, val);
33: }
34: }
35: }
36: return(0);
37: }
39: /*@
40: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
42: Not Collective
44: Input Parameter:
45: + dm - The original DM
46: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)
48: Output Parameter:
49: . label - The DMLabel marking boundary faces with the given value
51: Level: developer
53: .seealso: DMLabelCreate(), DMCreateLabel()
54: @*/
55: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
56: {
60: DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
61: return(0);
62: }
64: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
65: {
66: IS valueIS;
67: PetscSF sfPoint;
68: const PetscInt *values;
69: PetscInt numValues, v, cStart, cEnd, nroots;
70: PetscErrorCode ierr;
73: DMLabelGetNumValues(label, &numValues);
74: DMLabelGetValueIS(label, &valueIS);
75: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
76: ISGetIndices(valueIS, &values);
77: for (v = 0; v < numValues; ++v) {
78: IS pointIS;
79: const PetscInt *points;
80: PetscInt numPoints, p;
82: DMLabelGetStratumSize(label, values[v], &numPoints);
83: DMLabelGetStratumIS(label, values[v], &pointIS);
84: ISGetIndices(pointIS, &points);
85: for (p = 0; p < numPoints; ++p) {
86: PetscInt q = points[p];
87: PetscInt *closure = NULL;
88: PetscInt closureSize, c;
90: if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
91: continue;
92: }
93: DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
94: for (c = 0; c < closureSize*2; c += 2) {
95: DMLabelSetValue(label, closure[c], values[v]);
96: }
97: DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
98: }
99: ISRestoreIndices(pointIS, &points);
100: ISDestroy(&pointIS);
101: }
102: ISRestoreIndices(valueIS, &values);
103: ISDestroy(&valueIS);
104: DMGetPointSF(dm, &sfPoint);
105: PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
106: if (nroots >= 0) {
107: DMLabel lblRoots, lblLeaves;
108: IS valueIS, pointIS;
109: const PetscInt *values;
110: PetscInt numValues, v;
111: PetscErrorCode ierr;
113: DMGetPointSF(dm, &sfPoint);
114: /* Pull point contributions from remote leaves into local roots */
115: DMLabelGather(label, sfPoint, &lblLeaves);
116: DMLabelGetValueIS(lblLeaves, &valueIS);
117: ISGetLocalSize(valueIS, &numValues);
118: ISGetIndices(valueIS, &values);
119: for (v = 0; v < numValues; ++v) {
120: const PetscInt value = values[v];
121:
122: DMLabelGetStratumIS(lblLeaves, value, &pointIS);
123: DMLabelInsertIS(label, pointIS, value);
124: ISDestroy(&pointIS);
125: }
126: ISRestoreIndices(valueIS, &values);
127: ISDestroy(&valueIS);
128: DMLabelDestroy(&lblLeaves);
129: /* Push point contributions from roots into remote leaves */
130: DMLabelDistribute(label, sfPoint, &lblRoots);
131: DMLabelGetValueIS(lblRoots, &valueIS);
132: ISGetLocalSize(valueIS, &numValues);
133: ISGetIndices(valueIS, &values);
134: for (v = 0; v < numValues; ++v) {
135: const PetscInt value = values[v];
136:
137: DMLabelGetStratumIS(lblRoots, value, &pointIS);
138: DMLabelInsertIS(label, pointIS, value);
139: ISDestroy(&pointIS);
140: }
141: ISRestoreIndices(valueIS, &values);
142: ISDestroy(&valueIS);
143: DMLabelDestroy(&lblRoots);
144: }
145: return(0);
146: }
148: /*@
149: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
151: Input Parameters:
152: + dm - The DM
153: - label - A DMLabel marking the surface points
155: Output Parameter:
156: . label - A DMLabel marking all surface points in the transitive closure
158: Level: developer
160: .seealso: DMPlexLabelCohesiveComplete()
161: @*/
162: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
163: {
167: DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
168: return(0);
169: }
171: /*@
172: DMPlexLabelAddCells - Starting with a label marking faces on a surface, we add a cell for each face
174: Input Parameters:
175: + dm - The DM
176: - label - A DMLabel marking the surface points
178: Output Parameter:
179: . label - A DMLabel incorporating cells
181: Level: developer
183: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
185: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete()
186: @*/
187: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
188: {
189: IS valueIS;
190: const PetscInt *values;
191: PetscInt numValues, v, cStart, cEnd, cEndInterior;
192: PetscErrorCode ierr;
195: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
196: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
197: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
198: DMLabelGetNumValues(label, &numValues);
199: DMLabelGetValueIS(label, &valueIS);
200: ISGetIndices(valueIS, &values);
201: for (v = 0; v < numValues; ++v) {
202: IS pointIS;
203: const PetscInt *points;
204: PetscInt numPoints, p;
206: DMLabelGetStratumSize(label, values[v], &numPoints);
207: DMLabelGetStratumIS(label, values[v], &pointIS);
208: ISGetIndices(pointIS, &points);
209: for (p = 0; p < numPoints; ++p) {
210: PetscInt *closure = NULL;
211: PetscInt closureSize, point, cl;
213: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
214: for (cl = closureSize-1; cl > 0; --cl) {
215: point = closure[cl*2];
216: if ((point >= cStart) && (point < cEnd)) {DMLabelSetValue(label, point, values[v]); break;}
217: }
218: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
219: }
220: ISRestoreIndices(pointIS, &points);
221: ISDestroy(&pointIS);
222: }
223: ISRestoreIndices(valueIS, &values);
224: ISDestroy(&valueIS);
225: return(0);
226: }
228: /*@
229: DMPlexLabelClearCells - Remove cells from a label
231: Input Parameters:
232: + dm - The DM
233: - label - A DMLabel marking surface points and their adjacent cells
235: Output Parameter:
236: . label - A DMLabel without cells
238: Level: developer
240: Note: This undoes DMPlexLabelAddCells()
242: .seealso: DMPlexLabelComplete(), DMPlexLabelCohesiveComplete(), DMPlexLabelAddCells()
243: @*/
244: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
245: {
246: IS valueIS;
247: const PetscInt *values;
248: PetscInt numValues, v, cStart, cEnd, cEndInterior;
249: PetscErrorCode ierr;
252: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
253: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
254: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
255: DMLabelGetNumValues(label, &numValues);
256: DMLabelGetValueIS(label, &valueIS);
257: ISGetIndices(valueIS, &values);
258: for (v = 0; v < numValues; ++v) {
259: IS pointIS;
260: const PetscInt *points;
261: PetscInt numPoints, p;
263: DMLabelGetStratumSize(label, values[v], &numPoints);
264: DMLabelGetStratumIS(label, values[v], &pointIS);
265: ISGetIndices(pointIS, &points);
266: for (p = 0; p < numPoints; ++p) {
267: PetscInt point = points[p];
269: if (point >= cStart && point < cEnd) {
270: DMLabelClearValue(label,point,values[v]);
271: }
272: }
273: ISRestoreIndices(pointIS, &points);
274: ISDestroy(&pointIS);
275: }
276: ISRestoreIndices(valueIS, &values);
277: ISDestroy(&valueIS);
278: return(0);
279: }
281: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
282: * index (skipping first, which is (0,0)) */
283: PETSC_STATIC_INLINE PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
284: {
285: PetscInt d, off = 0;
288: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
289: for (d = 0; d < depth; d++) {
290: PetscInt firstd = d;
291: PetscInt firstStart = depthShift[2*d];
292: PetscInt e;
294: for (e = d+1; e <= depth; e++) {
295: if (depthShift[2*e] < firstStart) {
296: firstd = e;
297: firstStart = depthShift[2*d];
298: }
299: }
300: if (firstd != d) {
301: PetscInt swap[2];
303: e = firstd;
304: swap[0] = depthShift[2*d];
305: swap[1] = depthShift[2*d+1];
306: depthShift[2*d] = depthShift[2*e];
307: depthShift[2*d+1] = depthShift[2*e+1];
308: depthShift[2*e] = swap[0];
309: depthShift[2*e+1] = swap[1];
310: }
311: }
312: /* convert (oldstart, added) to (oldstart, newstart) */
313: for (d = 0; d <= depth; d++) {
314: off += depthShift[2*d+1];
315: depthShift[2*d+1] = depthShift[2*d] + off;
316: }
317: return(0);
318: }
320: /* depthShift is a list of (old, new) pairs */
321: PETSC_STATIC_INLINE PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
322: {
323: PetscInt d;
324: PetscInt newOff = 0;
326: for (d = 0; d <= depth; d++) {
327: if (p < depthShift[2*d]) return p + newOff;
328: else newOff = depthShift[2*d+1] - depthShift[2*d];
329: }
330: return p + newOff;
331: }
333: /* depthShift is a list of (old, new) pairs */
334: PETSC_STATIC_INLINE PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
335: {
336: PetscInt d;
337: PetscInt newOff = 0;
339: for (d = 0; d <= depth; d++) {
340: if (p < depthShift[2*d+1]) return p + newOff;
341: else newOff = depthShift[2*d] - depthShift[2*d+1];
342: }
343: return p + newOff;
344: }
346: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
347: {
348: PetscInt depth = 0, d, pStart, pEnd, p;
349: DMLabel depthLabel;
353: DMPlexGetDepth(dm, &depth);
354: if (depth < 0) return(0);
355: /* Step 1: Expand chart */
356: DMPlexGetChart(dm, &pStart, &pEnd);
357: pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
358: DMPlexSetChart(dmNew, pStart, pEnd);
359: DMCreateLabel(dmNew,"depth");
360: DMPlexGetDepthLabel(dmNew,&depthLabel);
361: /* Step 2: Set cone and support sizes */
362: for (d = 0; d <= depth; ++d) {
363: PetscInt pStartNew, pEndNew;
364: IS pIS;
366: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
367: pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
368: pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
369: ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
370: DMLabelSetStratumIS(depthLabel, d, pIS);
371: ISDestroy(&pIS);
372: for (p = pStart; p < pEnd; ++p) {
373: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
374: PetscInt size;
376: DMPlexGetConeSize(dm, p, &size);
377: DMPlexSetConeSize(dmNew, newp, size);
378: DMPlexGetSupportSize(dm, p, &size);
379: DMPlexSetSupportSize(dmNew, newp, size);
380: }
381: }
382: return(0);
383: }
385: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
386: {
387: PetscInt *newpoints;
388: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
392: DMPlexGetDepth(dm, &depth);
393: if (depth < 0) return(0);
394: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
395: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
396: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
397: /* Step 5: Set cones and supports */
398: DMPlexGetChart(dm, &pStart, &pEnd);
399: for (p = pStart; p < pEnd; ++p) {
400: const PetscInt *points = NULL, *orientations = NULL;
401: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
403: DMPlexGetConeSize(dm, p, &size);
404: DMPlexGetCone(dm, p, &points);
405: DMPlexGetConeOrientation(dm, p, &orientations);
406: for (i = 0; i < size; ++i) {
407: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
408: }
409: DMPlexSetCone(dmNew, newp, newpoints);
410: DMPlexSetConeOrientation(dmNew, newp, orientations);
411: DMPlexGetSupportSize(dm, p, &size);
412: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
413: DMPlexGetSupport(dm, p, &points);
414: for (i = 0; i < size; ++i) {
415: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
416: }
417: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
418: DMPlexSetSupport(dmNew, newp, newpoints);
419: }
420: PetscFree(newpoints);
421: return(0);
422: }
424: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
425: {
426: PetscSection coordSection, newCoordSection;
427: Vec coordinates, newCoordinates;
428: PetscScalar *coords, *newCoords;
429: PetscInt coordSize, sStart, sEnd;
430: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
431: PetscBool hasCells;
435: DMGetCoordinateDim(dm, &dim);
436: DMSetCoordinateDim(dmNew, dim);
437: DMPlexGetDepth(dm, &depth);
438: /* Step 8: Convert coordinates */
439: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
440: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
441: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
442: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
443: DMGetCoordinateSection(dm, &coordSection);
444: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
445: PetscSectionSetNumFields(newCoordSection, 1);
446: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
447: PetscSectionGetChart(coordSection, &sStart, &sEnd);
448: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
449: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
450: if (hasCells) {
451: for (c = cStart; c < cEnd; ++c) {
452: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
454: PetscSectionGetDof(coordSection, c, &dof);
455: PetscSectionSetDof(newCoordSection, cNew, dof);
456: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
457: }
458: }
459: for (v = vStartNew; v < vEndNew; ++v) {
460: PetscSectionSetDof(newCoordSection, v, dim);
461: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
462: }
463: PetscSectionSetUp(newCoordSection);
464: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
465: PetscSectionGetStorageSize(newCoordSection, &coordSize);
466: VecCreate(PETSC_COMM_SELF, &newCoordinates);
467: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
468: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
469: VecSetBlockSize(newCoordinates, dim);
470: VecSetType(newCoordinates,VECSTANDARD);
471: DMSetCoordinatesLocal(dmNew, newCoordinates);
472: DMGetCoordinatesLocal(dm, &coordinates);
473: VecGetArray(coordinates, &coords);
474: VecGetArray(newCoordinates, &newCoords);
475: if (hasCells) {
476: for (c = cStart; c < cEnd; ++c) {
477: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
479: PetscSectionGetDof(coordSection, c, &dof);
480: PetscSectionGetOffset(coordSection, c, &off);
481: PetscSectionGetOffset(newCoordSection, cNew, &noff);
482: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
483: }
484: }
485: for (v = vStart; v < vEnd; ++v) {
486: PetscInt dof, off, noff, d;
488: PetscSectionGetDof(coordSection, v, &dof);
489: PetscSectionGetOffset(coordSection, v, &off);
490: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
491: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
492: }
493: VecRestoreArray(coordinates, &coords);
494: VecRestoreArray(newCoordinates, &newCoords);
495: VecDestroy(&newCoordinates);
496: PetscSectionDestroy(&newCoordSection);
497: return(0);
498: }
500: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
501: {
502: PetscInt depth = 0;
503: PetscSF sfPoint, sfPointNew;
504: const PetscSFNode *remotePoints;
505: PetscSFNode *gremotePoints;
506: const PetscInt *localPoints;
507: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
508: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
509: PetscErrorCode ierr;
512: DMPlexGetDepth(dm, &depth);
513: /* Step 9: Convert pointSF */
514: DMGetPointSF(dm, &sfPoint);
515: DMGetPointSF(dmNew, &sfPointNew);
516: DMPlexGetChart(dm, &pStart, &pEnd);
517: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
518: totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
519: if (numRoots >= 0) {
520: PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
521: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
522: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
523: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
524: PetscMalloc1(numLeaves, &glocalPoints);
525: PetscMalloc1(numLeaves, &gremotePoints);
526: for (l = 0; l < numLeaves; ++l) {
527: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
528: gremotePoints[l].rank = remotePoints[l].rank;
529: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
530: }
531: PetscFree2(newLocation,newRemoteLocation);
532: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
533: }
534: return(0);
535: }
537: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
538: {
539: PetscSF sfPoint;
540: DMLabel vtkLabel, ghostLabel;
541: const PetscSFNode *leafRemote;
542: const PetscInt *leafLocal;
543: PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
544: PetscMPIInt rank;
545: PetscErrorCode ierr;
548: DMPlexGetDepth(dm, &depth);
549: /* Step 10: Convert labels */
550: DMGetNumLabels(dm, &numLabels);
551: for (l = 0; l < numLabels; ++l) {
552: DMLabel label, newlabel;
553: const char *lname;
554: PetscBool isDepth, isDim;
555: IS valueIS;
556: const PetscInt *values;
557: PetscInt numValues, val;
559: DMGetLabelName(dm, l, &lname);
560: PetscStrcmp(lname, "depth", &isDepth);
561: if (isDepth) continue;
562: PetscStrcmp(lname, "dim", &isDim);
563: if (isDim) continue;
564: DMCreateLabel(dmNew, lname);
565: DMGetLabel(dm, lname, &label);
566: DMGetLabel(dmNew, lname, &newlabel);
567: DMLabelGetDefaultValue(label,&val);
568: DMLabelSetDefaultValue(newlabel,val);
569: DMLabelGetValueIS(label, &valueIS);
570: ISGetLocalSize(valueIS, &numValues);
571: ISGetIndices(valueIS, &values);
572: for (val = 0; val < numValues; ++val) {
573: IS pointIS;
574: const PetscInt *points;
575: PetscInt numPoints, p;
577: DMLabelGetStratumIS(label, values[val], &pointIS);
578: ISGetLocalSize(pointIS, &numPoints);
579: ISGetIndices(pointIS, &points);
580: for (p = 0; p < numPoints; ++p) {
581: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
583: DMLabelSetValue(newlabel, newpoint, values[val]);
584: }
585: ISRestoreIndices(pointIS, &points);
586: ISDestroy(&pointIS);
587: }
588: ISRestoreIndices(valueIS, &values);
589: ISDestroy(&valueIS);
590: }
591: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
592: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
593: DMGetPointSF(dm, &sfPoint);
594: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
595: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
596: DMCreateLabel(dmNew, "vtk");
597: DMCreateLabel(dmNew, "ghost");
598: DMGetLabel(dmNew, "vtk", &vtkLabel);
599: DMGetLabel(dmNew, "ghost", &ghostLabel);
600: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
601: for (; c < leafLocal[l] && c < cEnd; ++c) {
602: DMLabelSetValue(vtkLabel, c, 1);
603: }
604: if (leafLocal[l] >= cEnd) break;
605: if (leafRemote[l].rank == rank) {
606: DMLabelSetValue(vtkLabel, c, 1);
607: } else {
608: DMLabelSetValue(ghostLabel, c, 2);
609: }
610: }
611: for (; c < cEnd; ++c) {
612: DMLabelSetValue(vtkLabel, c, 1);
613: }
614: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
615: for (f = fStart; f < fEnd; ++f) {
616: PetscInt numCells;
618: DMPlexGetSupportSize(dmNew, f, &numCells);
619: if (numCells < 2) {
620: DMLabelSetValue(ghostLabel, f, 1);
621: } else {
622: const PetscInt *cells = NULL;
623: PetscInt vA, vB;
625: DMPlexGetSupport(dmNew, f, &cells);
626: DMLabelGetValue(vtkLabel, cells[0], &vA);
627: DMLabelGetValue(vtkLabel, cells[1], &vB);
628: if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
629: }
630: }
631: return(0);
632: }
634: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
635: {
636: DM refTree;
637: PetscSection pSec;
638: PetscInt *parents, *childIDs;
642: DMPlexGetReferenceTree(dm,&refTree);
643: DMPlexSetReferenceTree(dmNew,refTree);
644: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
645: if (pSec) {
646: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
647: PetscInt *childIDsShifted;
648: PetscSection pSecShifted;
650: PetscSectionGetChart(pSec,&pStart,&pEnd);
651: DMPlexGetDepth(dm,&depth);
652: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
653: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
654: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
655: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
656: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
657: for (p = pStartShifted; p < pEndShifted; p++) {
658: /* start off assuming no children */
659: PetscSectionSetDof(pSecShifted,p,0);
660: }
661: for (p = pStart; p < pEnd; p++) {
662: PetscInt dof;
663: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
665: PetscSectionGetDof(pSec,p,&dof);
666: PetscSectionSetDof(pSecShifted,pNew,dof);
667: }
668: PetscSectionSetUp(pSecShifted);
669: for (p = pStart; p < pEnd; p++) {
670: PetscInt dof;
671: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
673: PetscSectionGetDof(pSec,p,&dof);
674: if (dof) {
675: PetscInt off, offNew;
677: PetscSectionGetOffset(pSec,p,&off);
678: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
679: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
680: childIDsShifted[offNew] = childIDs[off];
681: }
682: }
683: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
684: PetscFree2(parentsShifted,childIDsShifted);
685: PetscSectionDestroy(&pSecShifted);
686: }
687: return(0);
688: }
690: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
691: {
692: PetscSF sf;
693: IS valueIS;
694: const PetscInt *values, *leaves;
695: PetscInt *depthShift;
696: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
697: PetscBool isper;
698: const PetscReal *maxCell, *L;
699: const DMBoundaryType *bd;
700: PetscErrorCode ierr;
703: DMGetPointSF(dm, &sf);
704: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
705: nleaves = PetscMax(0, nleaves);
706: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
707: /* Count ghost cells */
708: DMLabelGetValueIS(label, &valueIS);
709: ISGetLocalSize(valueIS, &numFS);
710: ISGetIndices(valueIS, &values);
711: Ng = 0;
712: for (fs = 0; fs < numFS; ++fs) {
713: IS faceIS;
714: const PetscInt *faces;
715: PetscInt numFaces, f, numBdFaces = 0;
717: DMLabelGetStratumIS(label, values[fs], &faceIS);
718: ISGetLocalSize(faceIS, &numFaces);
719: ISGetIndices(faceIS, &faces);
720: for (f = 0; f < numFaces; ++f) {
721: PetscInt numChildren;
723: PetscFindInt(faces[f], nleaves, leaves, &loc);
724: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
725: /* non-local and ancestors points don't get to register ghosts */
726: if (loc >= 0 || numChildren) continue;
727: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
728: }
729: Ng += numBdFaces;
730: ISDestroy(&faceIS);
731: }
732: DMPlexGetDepth(dm, &depth);
733: PetscMalloc1(2*(depth+1), &depthShift);
734: for (d = 0; d <= depth; d++) {
735: PetscInt dEnd;
737: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
738: depthShift[2*d] = dEnd;
739: depthShift[2*d+1] = 0;
740: }
741: if (depth >= 0) depthShift[2*depth+1] = Ng;
742: DMPlexShiftPointSetUp_Internal(depth,depthShift);
743: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
744: /* Step 3: Set cone/support sizes for new points */
745: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
746: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
747: for (c = cEnd; c < cEnd + Ng; ++c) {
748: DMPlexSetConeSize(gdm, c, 1);
749: }
750: for (fs = 0; fs < numFS; ++fs) {
751: IS faceIS;
752: const PetscInt *faces;
753: PetscInt numFaces, f;
755: DMLabelGetStratumIS(label, values[fs], &faceIS);
756: ISGetLocalSize(faceIS, &numFaces);
757: ISGetIndices(faceIS, &faces);
758: for (f = 0; f < numFaces; ++f) {
759: PetscInt size, numChildren;
761: PetscFindInt(faces[f], nleaves, leaves, &loc);
762: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
763: if (loc >= 0 || numChildren) continue;
764: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
765: DMPlexGetSupportSize(dm, faces[f], &size);
766: if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
767: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
768: }
769: ISRestoreIndices(faceIS, &faces);
770: ISDestroy(&faceIS);
771: }
772: /* Step 4: Setup ghosted DM */
773: DMSetUp(gdm);
774: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
775: /* Step 6: Set cones and supports for new points */
776: ghostCell = cEnd;
777: for (fs = 0; fs < numFS; ++fs) {
778: IS faceIS;
779: const PetscInt *faces;
780: PetscInt numFaces, f;
782: DMLabelGetStratumIS(label, values[fs], &faceIS);
783: ISGetLocalSize(faceIS, &numFaces);
784: ISGetIndices(faceIS, &faces);
785: for (f = 0; f < numFaces; ++f) {
786: PetscInt newFace = faces[f] + Ng, numChildren;
788: PetscFindInt(faces[f], nleaves, leaves, &loc);
789: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
790: if (loc >= 0 || numChildren) continue;
791: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
792: DMPlexSetCone(gdm, ghostCell, &newFace);
793: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
794: ++ghostCell;
795: }
796: ISRestoreIndices(faceIS, &faces);
797: ISDestroy(&faceIS);
798: }
799: ISRestoreIndices(valueIS, &values);
800: ISDestroy(&valueIS);
801: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
802: DMPlexShiftSF_Internal(dm, depthShift, gdm);
803: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
804: DMPlexShiftTree_Internal(dm, depthShift, gdm);
805: PetscFree(depthShift);
806: /* Step 7: Periodicity */
807: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
808: DMSetPeriodicity(gdm, isper, maxCell, L, bd);
809: if (numGhostCells) *numGhostCells = Ng;
810: return(0);
811: }
813: /*@C
814: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
816: Collective on dm
818: Input Parameters:
819: + dm - The original DM
820: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
822: Output Parameters:
823: + numGhostCells - The number of ghost cells added to the DM
824: - dmGhosted - The new DM
826: Note: If no label exists of that name, one will be created marking all boundary faces
828: Level: developer
830: .seealso: DMCreate()
831: @*/
832: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
833: {
834: DM gdm;
835: DMLabel label;
836: const char *name = labelName ? labelName : "Face Sets";
837: PetscInt dim;
838: PetscBool useCone, useClosure;
845: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
846: DMSetType(gdm, DMPLEX);
847: DMGetDimension(dm, &dim);
848: DMSetDimension(gdm, dim);
849: DMGetBasicAdjacency(dm, &useCone, &useClosure);
850: DMSetBasicAdjacency(gdm, useCone, useClosure);
851: DMGetLabel(dm, name, &label);
852: if (!label) {
853: /* Get label for boundary faces */
854: DMCreateLabel(dm, name);
855: DMGetLabel(dm, name, &label);
856: DMPlexMarkBoundaryFaces(dm, 1, label);
857: }
858: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
859: DMCopyBoundary(dm, gdm);
860: DMCopyDisc(dm, gdm);
861: gdm->setfromoptionscalled = dm->setfromoptionscalled;
862: *dmGhosted = gdm;
863: return(0);
864: }
866: /*
867: We are adding three kinds of points here:
868: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
869: Non-replicated: Points which exist on the fault, but are not replicated
870: Hybrid: Entirely new points, such as cohesive cells
872: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
873: each depth so that the new split/hybrid points can be inserted as a block.
874: */
875: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
876: {
877: MPI_Comm comm;
878: IS valueIS;
879: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
880: const PetscInt *values; /* List of depths for which we have replicated points */
881: IS *splitIS;
882: IS *unsplitIS;
883: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
884: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
885: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
886: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
887: const PetscInt **splitPoints; /* Replicated points for each depth */
888: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
889: PetscSection coordSection;
890: Vec coordinates;
891: PetscScalar *coords;
892: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
893: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
894: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
895: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
896: PetscInt *coneNew, *coneONew, *supportNew;
897: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
898: PetscErrorCode ierr;
901: PetscObjectGetComm((PetscObject)dm,&comm);
902: DMGetDimension(dm, &dim);
903: DMPlexGetDepth(dm, &depth);
904: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
905: /* Count split points and add cohesive cells */
906: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
907: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
908: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
909: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
910: for (d = 0; d <= depth; ++d) {
911: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
912: depthEnd[d] = pMaxNew[d];
913: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
914: numSplitPoints[d] = 0;
915: numUnsplitPoints[d] = 0;
916: numHybridPoints[d] = 0;
917: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
918: splitPoints[d] = NULL;
919: unsplitPoints[d] = NULL;
920: splitIS[d] = NULL;
921: unsplitIS[d] = NULL;
922: /* we are shifting the existing hybrid points with the stratum behind them, so
923: * the split comes at the end of the normal points, i.e., at depthMax[d] */
924: depthShift[2*d] = depthMax[d];
925: depthShift[2*d+1] = 0;
926: }
927: if (label) {
928: DMLabelGetValueIS(label, &valueIS);
929: ISGetLocalSize(valueIS, &numSP);
930: ISGetIndices(valueIS, &values);
931: }
932: for (sp = 0; sp < numSP; ++sp) {
933: const PetscInt dep = values[sp];
935: if ((dep < 0) || (dep > depth)) continue;
936: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
937: if (splitIS[dep]) {
938: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
939: ISGetIndices(splitIS[dep], &splitPoints[dep]);
940: }
941: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
942: if (unsplitIS[dep]) {
943: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
944: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
945: }
946: }
947: /* Calculate number of hybrid points */
948: 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 */
949: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
950: DMPlexShiftPointSetUp_Internal(depth,depthShift);
951: /* the end of the points in this stratum that come before the new points:
952: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
953: * added points */
954: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
955: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
956: /* Step 3: Set cone/support sizes for new points */
957: for (dep = 0; dep <= depth; ++dep) {
958: for (p = 0; p < numSplitPoints[dep]; ++p) {
959: const PetscInt oldp = splitPoints[dep][p];
960: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
961: const PetscInt splitp = p + pMaxNew[dep];
962: const PetscInt *support;
963: PetscInt coneSize, supportSize, qf, qn, qp, e;
965: DMPlexGetConeSize(dm, oldp, &coneSize);
966: DMPlexSetConeSize(sdm, splitp, coneSize);
967: DMPlexGetSupportSize(dm, oldp, &supportSize);
968: DMPlexSetSupportSize(sdm, splitp, supportSize);
969: if (dep == depth-1) {
970: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
972: /* Add cohesive cells, they are prisms */
973: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
974: } else if (dep == 0) {
975: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
977: DMPlexGetSupport(dm, oldp, &support);
978: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
979: PetscInt val;
981: DMLabelGetValue(label, support[e], &val);
982: if (val == 1) ++qf;
983: if ((val == 1) || (val == (shift + 1))) ++qn;
984: if ((val == 1) || (val == -(shift + 1))) ++qp;
985: }
986: /* Split old vertex: Edges into original vertex and new cohesive edge */
987: DMPlexSetSupportSize(sdm, newp, qn+1);
988: /* Split new vertex: Edges into split vertex and new cohesive edge */
989: DMPlexSetSupportSize(sdm, splitp, qp+1);
990: /* Add hybrid edge */
991: DMPlexSetConeSize(sdm, hybedge, 2);
992: DMPlexSetSupportSize(sdm, hybedge, qf);
993: } else if (dep == dim-2) {
994: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
996: DMPlexGetSupport(dm, oldp, &support);
997: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
998: PetscInt val;
1000: DMLabelGetValue(label, support[e], &val);
1001: if (val == dim-1) ++qf;
1002: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
1003: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1004: }
1005: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1006: DMPlexSetSupportSize(sdm, newp, qn+1);
1007: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1008: DMPlexSetSupportSize(sdm, splitp, qp+1);
1009: /* Add hybrid face */
1010: DMPlexSetConeSize(sdm, hybface, 4);
1011: DMPlexSetSupportSize(sdm, hybface, qf);
1012: }
1013: }
1014: }
1015: for (dep = 0; dep <= depth; ++dep) {
1016: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1017: const PetscInt oldp = unsplitPoints[dep][p];
1018: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1019: const PetscInt *support;
1020: PetscInt coneSize, supportSize, qf, e, s;
1022: DMPlexGetConeSize(dm, oldp, &coneSize);
1023: DMPlexGetSupportSize(dm, oldp, &supportSize);
1024: DMPlexGetSupport(dm, oldp, &support);
1025: if (dep == 0) {
1026: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1028: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1029: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1030: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1031: if (e >= 0) ++qf;
1032: }
1033: DMPlexSetSupportSize(sdm, newp, qf+2);
1034: /* Add hybrid edge */
1035: DMPlexSetConeSize(sdm, hybedge, 2);
1036: for (e = 0, qf = 0; e < supportSize; ++e) {
1037: PetscInt val;
1039: DMLabelGetValue(label, support[e], &val);
1040: /* Split and unsplit edges produce hybrid faces */
1041: if (val == 1) ++qf;
1042: if (val == (shift2 + 1)) ++qf;
1043: }
1044: DMPlexSetSupportSize(sdm, hybedge, qf);
1045: } else if (dep == dim-2) {
1046: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1047: PetscInt val;
1049: for (e = 0, qf = 0; e < supportSize; ++e) {
1050: DMLabelGetValue(label, support[e], &val);
1051: if (val == dim-1) qf += 2;
1052: else ++qf;
1053: }
1054: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1055: DMPlexSetSupportSize(sdm, newp, qf+2);
1056: /* Add hybrid face */
1057: for (e = 0, qf = 0; e < supportSize; ++e) {
1058: DMLabelGetValue(label, support[e], &val);
1059: if (val == dim-1) ++qf;
1060: }
1061: DMPlexSetConeSize(sdm, hybface, 4);
1062: DMPlexSetSupportSize(sdm, hybface, qf);
1063: }
1064: }
1065: }
1066: /* Step 4: Setup split DM */
1067: DMSetUp(sdm);
1068: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1069: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1070: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1071: /* Step 6: Set cones and supports for new points */
1072: for (dep = 0; dep <= depth; ++dep) {
1073: for (p = 0; p < numSplitPoints[dep]; ++p) {
1074: const PetscInt oldp = splitPoints[dep][p];
1075: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1076: const PetscInt splitp = p + pMaxNew[dep];
1077: const PetscInt *cone, *support, *ornt;
1078: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1080: DMPlexGetConeSize(dm, oldp, &coneSize);
1081: DMPlexGetCone(dm, oldp, &cone);
1082: DMPlexGetConeOrientation(dm, oldp, &ornt);
1083: DMPlexGetSupportSize(dm, oldp, &supportSize);
1084: DMPlexGetSupport(dm, oldp, &support);
1085: if (dep == depth-1) {
1086: PetscBool hasUnsplit = PETSC_FALSE;
1087: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1088: const PetscInt *supportF;
1090: /* Split face: copy in old face to new face to start */
1091: DMPlexGetSupport(sdm, newp, &supportF);
1092: DMPlexSetSupport(sdm, splitp, supportF);
1093: /* Split old face: old vertices/edges in cone so no change */
1094: /* Split new face: new vertices/edges in cone */
1095: for (q = 0; q < coneSize; ++q) {
1096: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1097: if (v < 0) {
1098: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1099: 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);
1100: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1101: hasUnsplit = PETSC_TRUE;
1102: } else {
1103: coneNew[2+q] = v + pMaxNew[dep-1];
1104: if (dep > 1) {
1105: const PetscInt *econe;
1106: PetscInt econeSize, r, vs, vu;
1108: DMPlexGetConeSize(dm, cone[q], &econeSize);
1109: DMPlexGetCone(dm, cone[q], &econe);
1110: for (r = 0; r < econeSize; ++r) {
1111: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1112: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1113: if (vs >= 0) continue;
1114: 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);
1115: hasUnsplit = PETSC_TRUE;
1116: }
1117: }
1118: }
1119: }
1120: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1121: DMPlexSetConeOrientation(sdm, splitp, ornt);
1122: /* Face support */
1123: for (s = 0; s < supportSize; ++s) {
1124: PetscInt val;
1126: DMLabelGetValue(label, support[s], &val);
1127: if (val < 0) {
1128: /* Split old face: Replace negative side cell with cohesive cell */
1129: DMPlexInsertSupport(sdm, newp, s, hybcell);
1130: } else {
1131: /* Split new face: Replace positive side cell with cohesive cell */
1132: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1133: /* Get orientation for cohesive face */
1134: {
1135: const PetscInt *ncone, *nconeO;
1136: PetscInt nconeSize, nc;
1138: DMPlexGetConeSize(dm, support[s], &nconeSize);
1139: DMPlexGetCone(dm, support[s], &ncone);
1140: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1141: for (nc = 0; nc < nconeSize; ++nc) {
1142: if (ncone[nc] == oldp) {
1143: coneONew[0] = nconeO[nc];
1144: break;
1145: }
1146: }
1147: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1148: }
1149: }
1150: }
1151: /* Cohesive cell: Old and new split face, then new cohesive faces */
1152: coneNew[0] = newp; /* Extracted negative side orientation above */
1153: coneNew[1] = splitp;
1154: coneONew[1] = coneONew[0];
1155: for (q = 0; q < coneSize; ++q) {
1156: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1157: if (v < 0) {
1158: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1159: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1160: coneONew[2+q] = 0;
1161: } else {
1162: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1163: }
1164: coneONew[2+q] = 0;
1165: }
1166: DMPlexSetCone(sdm, hybcell, coneNew);
1167: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1168: /* Label the hybrid cells on the boundary of the split */
1169: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1170: } else if (dep == 0) {
1171: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1173: /* Split old vertex: Edges in old split faces and new cohesive edge */
1174: for (e = 0, qn = 0; e < supportSize; ++e) {
1175: PetscInt val;
1177: DMLabelGetValue(label, support[e], &val);
1178: if ((val == 1) || (val == (shift + 1))) {
1179: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1180: }
1181: }
1182: supportNew[qn] = hybedge;
1183: DMPlexSetSupport(sdm, newp, supportNew);
1184: /* Split new vertex: Edges in new split faces and new cohesive edge */
1185: for (e = 0, qp = 0; e < supportSize; ++e) {
1186: PetscInt val, edge;
1188: DMLabelGetValue(label, support[e], &val);
1189: if (val == 1) {
1190: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1191: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1192: supportNew[qp++] = edge + pMaxNew[dep+1];
1193: } else if (val == -(shift + 1)) {
1194: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1195: }
1196: }
1197: supportNew[qp] = hybedge;
1198: DMPlexSetSupport(sdm, splitp, supportNew);
1199: /* Hybrid edge: Old and new split vertex */
1200: coneNew[0] = newp;
1201: coneNew[1] = splitp;
1202: DMPlexSetCone(sdm, hybedge, coneNew);
1203: for (e = 0, qf = 0; e < supportSize; ++e) {
1204: PetscInt val, edge;
1206: DMLabelGetValue(label, support[e], &val);
1207: if (val == 1) {
1208: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1209: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1210: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1211: }
1212: }
1213: DMPlexSetSupport(sdm, hybedge, supportNew);
1214: } else if (dep == dim-2) {
1215: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1217: /* Split old edge: old vertices in cone so no change */
1218: /* Split new edge: new vertices in cone */
1219: for (q = 0; q < coneSize; ++q) {
1220: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1221: if (v < 0) {
1222: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1223: 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);
1224: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1225: } else {
1226: coneNew[q] = v + pMaxNew[dep-1];
1227: }
1228: }
1229: DMPlexSetCone(sdm, splitp, coneNew);
1230: /* Split old edge: Faces in positive side cells and old split faces */
1231: for (e = 0, q = 0; e < supportSize; ++e) {
1232: PetscInt val;
1234: DMLabelGetValue(label, support[e], &val);
1235: if (val == dim-1) {
1236: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1237: } else if (val == (shift + dim-1)) {
1238: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1239: }
1240: }
1241: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1242: DMPlexSetSupport(sdm, newp, supportNew);
1243: /* Split new edge: Faces in negative side cells and new split faces */
1244: for (e = 0, q = 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[q++] = face + pMaxNew[dep+1];
1252: } else if (val == -(shift + dim-1)) {
1253: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1254: }
1255: }
1256: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1257: DMPlexSetSupport(sdm, splitp, supportNew);
1258: /* Hybrid face */
1259: coneNew[0] = newp;
1260: coneNew[1] = splitp;
1261: for (v = 0; v < coneSize; ++v) {
1262: PetscInt vertex;
1263: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1264: if (vertex < 0) {
1265: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1266: 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);
1267: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1268: } else {
1269: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1270: }
1271: }
1272: DMPlexSetCone(sdm, hybface, coneNew);
1273: for (e = 0, qf = 0; e < supportSize; ++e) {
1274: PetscInt val, face;
1276: DMLabelGetValue(label, support[e], &val);
1277: if (val == dim-1) {
1278: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1279: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1280: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1281: }
1282: }
1283: DMPlexSetSupport(sdm, hybface, supportNew);
1284: }
1285: }
1286: }
1287: for (dep = 0; dep <= depth; ++dep) {
1288: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1289: const PetscInt oldp = unsplitPoints[dep][p];
1290: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1291: const PetscInt *cone, *support, *ornt;
1292: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1294: DMPlexGetConeSize(dm, oldp, &coneSize);
1295: DMPlexGetCone(dm, oldp, &cone);
1296: DMPlexGetConeOrientation(dm, oldp, &ornt);
1297: DMPlexGetSupportSize(dm, oldp, &supportSize);
1298: DMPlexGetSupport(dm, oldp, &support);
1299: if (dep == 0) {
1300: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1302: /* Unsplit vertex */
1303: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1304: for (s = 0, q = 0; s < supportSize; ++s) {
1305: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1306: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1307: if (e >= 0) {
1308: supportNew[q++] = e + pMaxNew[dep+1];
1309: }
1310: }
1311: supportNew[q++] = hybedge;
1312: supportNew[q++] = hybedge;
1313: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1314: DMPlexSetSupport(sdm, newp, supportNew);
1315: /* Hybrid edge */
1316: coneNew[0] = newp;
1317: coneNew[1] = newp;
1318: DMPlexSetCone(sdm, hybedge, coneNew);
1319: for (e = 0, qf = 0; e < supportSize; ++e) {
1320: PetscInt val, edge;
1322: DMLabelGetValue(label, support[e], &val);
1323: if (val == 1) {
1324: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1325: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1326: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1327: } else if (val == (shift2 + 1)) {
1328: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1329: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1330: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1331: }
1332: }
1333: DMPlexSetSupport(sdm, hybedge, supportNew);
1334: } else if (dep == dim-2) {
1335: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1337: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1338: for (f = 0, qf = 0; f < supportSize; ++f) {
1339: PetscInt val, face;
1341: DMLabelGetValue(label, support[f], &val);
1342: if (val == dim-1) {
1343: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1344: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1345: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1346: supportNew[qf++] = face + pMaxNew[dep+1];
1347: } else {
1348: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1349: }
1350: }
1351: supportNew[qf++] = hybface;
1352: supportNew[qf++] = hybface;
1353: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1354: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1355: DMPlexSetSupport(sdm, newp, supportNew);
1356: /* Add hybrid face */
1357: coneNew[0] = newp;
1358: coneNew[1] = newp;
1359: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1360: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1361: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1362: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1363: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1364: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1365: DMPlexSetCone(sdm, hybface, coneNew);
1366: for (f = 0, qf = 0; f < supportSize; ++f) {
1367: PetscInt val, face;
1369: DMLabelGetValue(label, support[f], &val);
1370: if (val == dim-1) {
1371: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1372: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1373: }
1374: }
1375: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1376: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1377: DMPlexSetSupport(sdm, hybface, supportNew);
1378: }
1379: }
1380: }
1381: /* Step 6b: Replace split points in negative side cones */
1382: for (sp = 0; sp < numSP; ++sp) {
1383: PetscInt dep = values[sp];
1384: IS pIS;
1385: PetscInt numPoints;
1386: const PetscInt *points;
1388: if (dep >= 0) continue;
1389: DMLabelGetStratumIS(label, dep, &pIS);
1390: if (!pIS) continue;
1391: dep = -dep - shift;
1392: ISGetLocalSize(pIS, &numPoints);
1393: ISGetIndices(pIS, &points);
1394: for (p = 0; p < numPoints; ++p) {
1395: const PetscInt oldp = points[p];
1396: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1397: const PetscInt *cone;
1398: PetscInt coneSize, c;
1399: /* PetscBool replaced = PETSC_FALSE; */
1401: /* Negative edge: replace split vertex */
1402: /* Negative cell: replace split face */
1403: DMPlexGetConeSize(sdm, newp, &coneSize);
1404: DMPlexGetCone(sdm, newp, &cone);
1405: for (c = 0; c < coneSize; ++c) {
1406: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1407: PetscInt csplitp, cp, val;
1409: DMLabelGetValue(label, coldp, &val);
1410: if (val == dep-1) {
1411: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1412: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1413: csplitp = pMaxNew[dep-1] + cp;
1414: DMPlexInsertCone(sdm, newp, c, csplitp);
1415: /* replaced = PETSC_TRUE; */
1416: }
1417: }
1418: /* Cells with only a vertex or edge on the submesh have no replacement */
1419: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1420: }
1421: ISRestoreIndices(pIS, &points);
1422: ISDestroy(&pIS);
1423: }
1424: /* Step 7: Coordinates */
1425: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1426: DMGetCoordinateSection(sdm, &coordSection);
1427: DMGetCoordinatesLocal(sdm, &coordinates);
1428: VecGetArray(coordinates, &coords);
1429: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1430: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1431: const PetscInt splitp = pMaxNew[0] + v;
1432: PetscInt dof, off, soff, d;
1434: PetscSectionGetDof(coordSection, newp, &dof);
1435: PetscSectionGetOffset(coordSection, newp, &off);
1436: PetscSectionGetOffset(coordSection, splitp, &soff);
1437: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1438: }
1439: VecRestoreArray(coordinates, &coords);
1440: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1441: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1442: /* Step 9: Labels */
1443: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1444: DMGetNumLabels(sdm, &numLabels);
1445: for (dep = 0; dep <= depth; ++dep) {
1446: for (p = 0; p < numSplitPoints[dep]; ++p) {
1447: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1448: const PetscInt splitp = pMaxNew[dep] + p;
1449: PetscInt l;
1451: if (splitLabel) {
1452: const PetscInt val = 100 + dep;
1454: DMLabelSetValue(splitLabel, newp, val);
1455: DMLabelSetValue(splitLabel, splitp, -val);
1456: }
1457: for (l = 0; l < numLabels; ++l) {
1458: DMLabel mlabel;
1459: const char *lname;
1460: PetscInt val;
1461: PetscBool isDepth;
1463: DMGetLabelName(sdm, l, &lname);
1464: PetscStrcmp(lname, "depth", &isDepth);
1465: if (isDepth) continue;
1466: DMGetLabel(sdm, lname, &mlabel);
1467: DMLabelGetValue(mlabel, newp, &val);
1468: if (val >= 0) {
1469: DMLabelSetValue(mlabel, splitp, val);
1470: }
1471: }
1472: }
1473: }
1474: for (sp = 0; sp < numSP; ++sp) {
1475: const PetscInt dep = values[sp];
1477: if ((dep < 0) || (dep > depth)) continue;
1478: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1479: ISDestroy(&splitIS[dep]);
1480: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1481: ISDestroy(&unsplitIS[dep]);
1482: }
1483: if (label) {
1484: ISRestoreIndices(valueIS, &values);
1485: ISDestroy(&valueIS);
1486: }
1487: for (d = 0; d <= depth; ++d) {
1488: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1489: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1490: }
1491: 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);
1492: PetscFree3(coneNew, coneONew, supportNew);
1493: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1494: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1495: return(0);
1496: }
1498: /*@C
1499: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1501: Collective on dm
1503: Input Parameters:
1504: + dm - The original DM
1505: - label - The label specifying the boundary faces (this could be auto-generated)
1507: Output Parameters:
1508: + splitLabel - The label containing the split points, or NULL if no output is desired
1509: - dmSplit - The new DM
1511: Level: developer
1513: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1514: @*/
1515: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1516: {
1517: DM sdm;
1518: PetscInt dim;
1524: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1525: DMSetType(sdm, DMPLEX);
1526: DMGetDimension(dm, &dim);
1527: DMSetDimension(sdm, dim);
1528: switch (dim) {
1529: case 2:
1530: case 3:
1531: DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1532: break;
1533: default:
1534: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1535: }
1536: *dmSplit = sdm;
1537: return(0);
1538: }
1540: /* Returns the side of the surface for a given cell with a face on the surface */
1541: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1542: {
1543: const PetscInt *cone, *ornt;
1544: PetscInt dim, coneSize, c;
1545: PetscErrorCode ierr;
1548: *pos = PETSC_TRUE;
1549: DMGetDimension(dm, &dim);
1550: DMPlexGetConeSize(dm, cell, &coneSize);
1551: DMPlexGetCone(dm, cell, &cone);
1552: DMPlexGetConeOrientation(dm, cell, &ornt);
1553: for (c = 0; c < coneSize; ++c) {
1554: if (cone[c] == face) {
1555: PetscInt o = ornt[c];
1557: if (subdm) {
1558: const PetscInt *subcone, *subornt;
1559: PetscInt subpoint, subface, subconeSize, sc;
1561: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1562: PetscFindInt(face, numSubpoints, subpoints, &subface);
1563: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1564: DMPlexGetCone(subdm, subpoint, &subcone);
1565: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1566: for (sc = 0; sc < subconeSize; ++sc) {
1567: if (subcone[sc] == subface) {
1568: o = subornt[0];
1569: break;
1570: }
1571: }
1572: 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);
1573: }
1574: if (o >= 0) *pos = PETSC_TRUE;
1575: else *pos = PETSC_FALSE;
1576: break;
1577: }
1578: }
1579: 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);
1580: return(0);
1581: }
1583: /*@
1584: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1585: to complete the surface
1587: Input Parameters:
1588: + dm - The DM
1589: . label - A DMLabel marking the surface
1590: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1591: . flip - Flag to flip the submesh normal and replace points on the other side
1592: - subdm - The subDM associated with the label, or NULL
1594: Output Parameter:
1595: . label - A DMLabel marking all surface points
1597: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1599: Level: developer
1601: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1602: @*/
1603: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1604: {
1605: DMLabel depthLabel;
1606: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1607: const PetscInt *points, *subpoints;
1608: const PetscInt rev = flip ? -1 : 1;
1609: PetscInt *pMax;
1610: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1611: PetscErrorCode ierr;
1614: DMPlexGetDepth(dm, &depth);
1615: DMGetDimension(dm, &dim);
1616: pSize = PetscMax(depth, dim) + 1;
1617: PetscMalloc1(pSize,&pMax);
1618: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1619: DMPlexGetDepthLabel(dm, &depthLabel);
1620: DMGetDimension(dm, &dim);
1621: if (subdm) {
1622: DMPlexCreateSubpointIS(subdm, &subpointIS);
1623: if (subpointIS) {
1624: ISGetLocalSize(subpointIS, &numSubpoints);
1625: ISGetIndices(subpointIS, &subpoints);
1626: }
1627: }
1628: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1629: DMLabelGetStratumIS(label, dim-1, &dimIS);
1630: if (!dimIS) {
1631: PetscFree(pMax);
1632: ISDestroy(&subpointIS);
1633: return(0);
1634: }
1635: ISGetLocalSize(dimIS, &numPoints);
1636: ISGetIndices(dimIS, &points);
1637: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1638: const PetscInt *support;
1639: PetscInt supportSize, s;
1641: DMPlexGetSupportSize(dm, points[p], &supportSize);
1642: #if 0
1643: if (supportSize != 2) {
1644: const PetscInt *lp;
1645: PetscInt Nlp, pind;
1647: /* Check that for a cell with a single support face, that face is in the SF */
1648: /* THis check only works for the remote side. We would need root side information */
1649: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1650: PetscFindInt(points[p], Nlp, lp, &pind);
1651: if (pind < 0) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports, and the face is not shared with another process", points[p], supportSize);
1652: }
1653: #endif
1654: DMPlexGetSupport(dm, points[p], &support);
1655: for (s = 0; s < supportSize; ++s) {
1656: const PetscInt *cone;
1657: PetscInt coneSize, c;
1658: PetscBool pos;
1660: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1661: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1662: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1663: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1664: /* Put faces touching the fault in the label */
1665: DMPlexGetConeSize(dm, support[s], &coneSize);
1666: DMPlexGetCone(dm, support[s], &cone);
1667: for (c = 0; c < coneSize; ++c) {
1668: const PetscInt point = cone[c];
1670: DMLabelGetValue(label, point, &val);
1671: if (val == -1) {
1672: PetscInt *closure = NULL;
1673: PetscInt closureSize, cl;
1675: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1676: for (cl = 0; cl < closureSize*2; cl += 2) {
1677: const PetscInt clp = closure[cl];
1678: PetscInt bval = -1;
1680: DMLabelGetValue(label, clp, &val);
1681: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1682: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1683: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1684: break;
1685: }
1686: }
1687: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1688: }
1689: }
1690: }
1691: }
1692: ISRestoreIndices(dimIS, &points);
1693: ISDestroy(&dimIS);
1694: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1695: ISDestroy(&subpointIS);
1696: /* Mark boundary points as unsplit */
1697: if (blabel) {
1698: DMLabelGetStratumIS(blabel, 1, &dimIS);
1699: ISGetLocalSize(dimIS, &numPoints);
1700: ISGetIndices(dimIS, &points);
1701: for (p = 0; p < numPoints; ++p) {
1702: const PetscInt point = points[p];
1703: PetscInt val, bval;
1705: DMLabelGetValue(blabel, point, &bval);
1706: if (bval >= 0) {
1707: DMLabelGetValue(label, point, &val);
1708: if ((val < 0) || (val > dim)) {
1709: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1710: DMLabelClearValue(blabel, point, bval);
1711: }
1712: }
1713: }
1714: for (p = 0; p < numPoints; ++p) {
1715: const PetscInt point = points[p];
1716: PetscInt val, bval;
1718: DMLabelGetValue(blabel, point, &bval);
1719: if (bval >= 0) {
1720: const PetscInt *cone, *support;
1721: PetscInt coneSize, supportSize, s, valA, valB, valE;
1723: /* Mark as unsplit */
1724: DMLabelGetValue(label, point, &val);
1725: 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);
1726: DMLabelClearValue(label, point, val);
1727: DMLabelSetValue(label, point, shift2+val);
1728: /* Check for cross-edge
1729: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1730: if (val != 0) continue;
1731: DMPlexGetSupport(dm, point, &support);
1732: DMPlexGetSupportSize(dm, point, &supportSize);
1733: for (s = 0; s < supportSize; ++s) {
1734: DMPlexGetCone(dm, support[s], &cone);
1735: DMPlexGetConeSize(dm, support[s], &coneSize);
1736: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1737: DMLabelGetValue(blabel, cone[0], &valA);
1738: DMLabelGetValue(blabel, cone[1], &valB);
1739: DMLabelGetValue(blabel, support[s], &valE);
1740: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1741: }
1742: }
1743: }
1744: ISRestoreIndices(dimIS, &points);
1745: ISDestroy(&dimIS);
1746: }
1747: /* Search for other cells/faces/edges connected to the fault by a vertex */
1748: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1749: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1750: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1751: cMax = cMax < 0 ? cEnd : cMax;
1752: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1753: DMLabelGetStratumIS(label, 0, &dimIS);
1754: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1755: if (dimIS && crossEdgeIS) {
1756: IS vertIS = dimIS;
1758: ISExpand(vertIS, crossEdgeIS, &dimIS);
1759: ISDestroy(&crossEdgeIS);
1760: ISDestroy(&vertIS);
1761: }
1762: if (!dimIS) {
1763: PetscFree(pMax);
1764: return(0);
1765: }
1766: ISGetLocalSize(dimIS, &numPoints);
1767: ISGetIndices(dimIS, &points);
1768: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1769: PetscInt *star = NULL;
1770: PetscInt starSize, s;
1771: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1773: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1774: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1775: while (again) {
1776: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1777: again = 0;
1778: for (s = 0; s < starSize*2; s += 2) {
1779: const PetscInt point = star[s];
1780: const PetscInt *cone;
1781: PetscInt coneSize, c;
1783: if ((point < cStart) || (point >= cMax)) continue;
1784: DMLabelGetValue(label, point, &val);
1785: if (val != -1) continue;
1786: again = again == 1 ? 1 : 2;
1787: DMPlexGetConeSize(dm, point, &coneSize);
1788: DMPlexGetCone(dm, point, &cone);
1789: for (c = 0; c < coneSize; ++c) {
1790: DMLabelGetValue(label, cone[c], &val);
1791: if (val != -1) {
1792: const PetscInt *ccone;
1793: PetscInt cconeSize, cc, side;
1795: 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);
1796: if (val > 0) side = 1;
1797: else side = -1;
1798: DMLabelSetValue(label, point, side*(shift+dim));
1799: /* Mark cell faces which touch the fault */
1800: DMPlexGetConeSize(dm, point, &cconeSize);
1801: DMPlexGetCone(dm, point, &ccone);
1802: for (cc = 0; cc < cconeSize; ++cc) {
1803: PetscInt *closure = NULL;
1804: PetscInt closureSize, cl;
1806: DMLabelGetValue(label, ccone[cc], &val);
1807: if (val != -1) continue;
1808: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1809: for (cl = 0; cl < closureSize*2; cl += 2) {
1810: const PetscInt clp = closure[cl];
1812: DMLabelGetValue(label, clp, &val);
1813: if (val == -1) continue;
1814: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1815: break;
1816: }
1817: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1818: }
1819: again = 1;
1820: break;
1821: }
1822: }
1823: }
1824: }
1825: /* Classify the rest by cell membership */
1826: for (s = 0; s < starSize*2; s += 2) {
1827: const PetscInt point = star[s];
1829: DMLabelGetValue(label, point, &val);
1830: if (val == -1) {
1831: PetscInt *sstar = NULL;
1832: PetscInt sstarSize, ss;
1833: PetscBool marked = PETSC_FALSE;
1835: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1836: for (ss = 0; ss < sstarSize*2; ss += 2) {
1837: const PetscInt spoint = sstar[ss];
1839: if ((spoint < cStart) || (spoint >= cMax)) continue;
1840: DMLabelGetValue(label, spoint, &val);
1841: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1842: DMLabelGetValue(depthLabel, point, &dep);
1843: if (val > 0) {
1844: DMLabelSetValue(label, point, shift+dep);
1845: } else {
1846: DMLabelSetValue(label, point, -(shift+dep));
1847: }
1848: marked = PETSC_TRUE;
1849: break;
1850: }
1851: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1852: DMLabelGetValue(depthLabel, point, &dep);
1853: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1854: }
1855: }
1856: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1857: }
1858: ISRestoreIndices(dimIS, &points);
1859: ISDestroy(&dimIS);
1860: /* If any faces touching the fault divide cells on either side, split them */
1861: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1862: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1863: ISExpand(facePosIS, faceNegIS, &dimIS);
1864: ISDestroy(&facePosIS);
1865: ISDestroy(&faceNegIS);
1866: ISGetLocalSize(dimIS, &numPoints);
1867: ISGetIndices(dimIS, &points);
1868: for (p = 0; p < numPoints; ++p) {
1869: const PetscInt point = points[p];
1870: const PetscInt *support;
1871: PetscInt supportSize, valA, valB;
1873: DMPlexGetSupportSize(dm, point, &supportSize);
1874: if (supportSize != 2) continue;
1875: DMPlexGetSupport(dm, point, &support);
1876: DMLabelGetValue(label, support[0], &valA);
1877: DMLabelGetValue(label, support[1], &valB);
1878: if ((valA == -1) || (valB == -1)) continue;
1879: if (valA*valB > 0) continue;
1880: /* Split the face */
1881: DMLabelGetValue(label, point, &valA);
1882: DMLabelClearValue(label, point, valA);
1883: DMLabelSetValue(label, point, dim-1);
1884: /* Label its closure:
1885: unmarked: label as unsplit
1886: incident: relabel as split
1887: split: do nothing
1888: */
1889: {
1890: PetscInt *closure = NULL;
1891: PetscInt closureSize, cl;
1893: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1894: for (cl = 0; cl < closureSize*2; cl += 2) {
1895: DMLabelGetValue(label, closure[cl], &valA);
1896: if (valA == -1) { /* Mark as unsplit */
1897: DMLabelGetValue(depthLabel, closure[cl], &dep);
1898: DMLabelSetValue(label, closure[cl], shift2+dep);
1899: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1900: DMLabelGetValue(depthLabel, closure[cl], &dep);
1901: DMLabelClearValue(label, closure[cl], valA);
1902: DMLabelSetValue(label, closure[cl], dep);
1903: }
1904: }
1905: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1906: }
1907: }
1908: ISRestoreIndices(dimIS, &points);
1909: ISDestroy(&dimIS);
1910: PetscFree(pMax);
1911: return(0);
1912: }
1914: /* Check that no cell have all vertices on the fault */
1915: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
1916: {
1917: IS subpointIS;
1918: const PetscInt *dmpoints;
1919: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
1920: PetscErrorCode ierr;
1923: if (!label) return(0);
1924: DMLabelGetDefaultValue(label, &defaultValue);
1925: DMPlexCreateSubpointIS(subdm, &subpointIS);
1926: if (!subpointIS) return(0);
1927: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
1928: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1929: ISGetIndices(subpointIS, &dmpoints);
1930: for (c = cStart; c < cEnd; ++c) {
1931: PetscBool invalidCell = PETSC_TRUE;
1932: PetscInt *closure = NULL;
1933: PetscInt closureSize, cl;
1935: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
1936: for (cl = 0; cl < closureSize*2; cl += 2) {
1937: PetscInt value = 0;
1939: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
1940: DMLabelGetValue(label, closure[cl], &value);
1941: if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
1942: }
1943: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
1944: if (invalidCell) {
1945: ISRestoreIndices(subpointIS, &dmpoints);
1946: ISDestroy(&subpointIS);
1947: DMDestroy(&subdm);
1948: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
1949: }
1950: }
1951: ISRestoreIndices(subpointIS, &dmpoints);
1952: ISDestroy(&subpointIS);
1953: return(0);
1954: }
1957: /*@
1958: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1960: Collective on dm
1962: Input Parameters:
1963: + dm - The original DM
1964: . label - The label specifying the interface vertices
1965: - bdlabel - The optional label specifying the interface boundary vertices
1967: Output Parameters:
1968: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
1969: . splitLabel - The label containing the split points, or NULL if no output is desired
1970: . dmInterface - The new interface DM, or NULL
1971: - dmHybrid - The new DM with cohesive cells
1973: Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
1974: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
1975: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
1976: one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
1977: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
1979: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
1980: mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
1981: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
1983: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
1984: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
1986: Level: developer
1988: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
1989: @*/
1990: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
1991: {
1992: DM idm;
1993: DMLabel subpointMap, hlabel, slabel = NULL;
1994: PetscInt dim;
2004: DMGetDimension(dm, &dim);
2005: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2006: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2007: DMPlexOrient(idm);
2008: DMPlexGetSubpointMap(idm, &subpointMap);
2009: DMLabelDuplicate(subpointMap, &hlabel);
2010: DMLabelClearStratum(hlabel, dim);
2011: if (splitLabel) {
2012: const char *name;
2013: char sname[PETSC_MAX_PATH_LEN];
2015: PetscObjectGetName((PetscObject) hlabel, &name);
2016: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2017: PetscStrcat(sname, " split");
2018: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2019: }
2020: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2021: if (dmInterface) {*dmInterface = idm;}
2022: else {DMDestroy(&idm);}
2023: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2024: if (hybridLabel) *hybridLabel = hlabel;
2025: else {DMLabelDestroy(&hlabel);}
2026: if (splitLabel) *splitLabel = slabel;
2027: return(0);
2028: }
2030: /* Here we need the explicit assumption that:
2032: For any marked cell, the marked vertices constitute a single face
2033: */
2034: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2035: {
2036: IS subvertexIS = NULL;
2037: const PetscInt *subvertices;
2038: PetscInt *pStart, *pEnd, *pMax, pSize;
2039: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2040: PetscErrorCode ierr;
2043: *numFaces = 0;
2044: *nFV = 0;
2045: DMPlexGetDepth(dm, &depth);
2046: DMGetDimension(dm, &dim);
2047: pSize = PetscMax(depth, dim) + 1;
2048: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
2049: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2050: for (d = 0; d <= depth; ++d) {
2051: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2052: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2053: }
2054: /* Loop over initial vertices and mark all faces in the collective star() */
2055: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2056: if (subvertexIS) {
2057: ISGetSize(subvertexIS, &numSubVerticesInitial);
2058: ISGetIndices(subvertexIS, &subvertices);
2059: }
2060: for (v = 0; v < numSubVerticesInitial; ++v) {
2061: const PetscInt vertex = subvertices[v];
2062: PetscInt *star = NULL;
2063: PetscInt starSize, s, numCells = 0, c;
2065: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2066: for (s = 0; s < starSize*2; s += 2) {
2067: const PetscInt point = star[s];
2068: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2069: }
2070: for (c = 0; c < numCells; ++c) {
2071: const PetscInt cell = star[c];
2072: PetscInt *closure = NULL;
2073: PetscInt closureSize, cl;
2074: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2076: DMLabelGetValue(subpointMap, cell, &cellLoc);
2077: if (cellLoc == 2) continue;
2078: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2079: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2080: for (cl = 0; cl < closureSize*2; cl += 2) {
2081: const PetscInt point = closure[cl];
2082: PetscInt vertexLoc;
2084: if ((point >= pStart[0]) && (point < pEnd[0])) {
2085: ++numCorners;
2086: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2087: if (vertexLoc == value) closure[faceSize++] = point;
2088: }
2089: }
2090: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2091: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2092: if (faceSize == *nFV) {
2093: const PetscInt *cells = NULL;
2094: PetscInt numCells, nc;
2096: ++(*numFaces);
2097: for (cl = 0; cl < faceSize; ++cl) {
2098: DMLabelSetValue(subpointMap, closure[cl], 0);
2099: }
2100: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2101: for (nc = 0; nc < numCells; ++nc) {
2102: DMLabelSetValue(subpointMap, cells[nc], 2);
2103: }
2104: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2105: }
2106: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2107: }
2108: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2109: }
2110: if (subvertexIS) {
2111: ISRestoreIndices(subvertexIS, &subvertices);
2112: }
2113: ISDestroy(&subvertexIS);
2114: PetscFree3(pStart,pEnd,pMax);
2115: return(0);
2116: }
2118: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2119: {
2120: IS subvertexIS = NULL;
2121: const PetscInt *subvertices;
2122: PetscInt *pStart, *pEnd, *pMax;
2123: PetscInt dim, d, numSubVerticesInitial = 0, v;
2124: PetscErrorCode ierr;
2127: DMGetDimension(dm, &dim);
2128: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2129: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2130: for (d = 0; d <= dim; ++d) {
2131: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2132: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2133: }
2134: /* Loop over initial vertices and mark all faces in the collective star() */
2135: if (vertexLabel) {
2136: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2137: if (subvertexIS) {
2138: ISGetSize(subvertexIS, &numSubVerticesInitial);
2139: ISGetIndices(subvertexIS, &subvertices);
2140: }
2141: }
2142: for (v = 0; v < numSubVerticesInitial; ++v) {
2143: const PetscInt vertex = subvertices[v];
2144: PetscInt *star = NULL;
2145: PetscInt starSize, s, numFaces = 0, f;
2147: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2148: for (s = 0; s < starSize*2; s += 2) {
2149: const PetscInt point = star[s];
2150: PetscInt faceLoc;
2152: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2153: if (markedFaces) {
2154: DMLabelGetValue(vertexLabel, point, &faceLoc);
2155: if (faceLoc < 0) continue;
2156: }
2157: star[numFaces++] = point;
2158: }
2159: }
2160: for (f = 0; f < numFaces; ++f) {
2161: const PetscInt face = star[f];
2162: PetscInt *closure = NULL;
2163: PetscInt closureSize, c;
2164: PetscInt faceLoc;
2166: DMLabelGetValue(subpointMap, face, &faceLoc);
2167: if (faceLoc == dim-1) continue;
2168: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2169: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2170: for (c = 0; c < closureSize*2; c += 2) {
2171: const PetscInt point = closure[c];
2172: PetscInt vertexLoc;
2174: if ((point >= pStart[0]) && (point < pEnd[0])) {
2175: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2176: if (vertexLoc != value) break;
2177: }
2178: }
2179: if (c == closureSize*2) {
2180: const PetscInt *support;
2181: PetscInt supportSize, s;
2183: for (c = 0; c < closureSize*2; c += 2) {
2184: const PetscInt point = closure[c];
2186: for (d = 0; d < dim; ++d) {
2187: if ((point >= pStart[d]) && (point < pEnd[d])) {
2188: DMLabelSetValue(subpointMap, point, d);
2189: break;
2190: }
2191: }
2192: }
2193: DMPlexGetSupportSize(dm, face, &supportSize);
2194: DMPlexGetSupport(dm, face, &support);
2195: for (s = 0; s < supportSize; ++s) {
2196: DMLabelSetValue(subpointMap, support[s], dim);
2197: }
2198: }
2199: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2200: }
2201: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2202: }
2203: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2204: ISDestroy(&subvertexIS);
2205: PetscFree3(pStart,pEnd,pMax);
2206: return(0);
2207: }
2209: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2210: {
2211: DMLabel label = NULL;
2212: const PetscInt *cone;
2213: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2214: PetscErrorCode ierr;
2217: *numFaces = 0;
2218: *nFV = 0;
2219: if (labelname) {DMGetLabel(dm, labelname, &label);}
2220: *subCells = NULL;
2221: DMGetDimension(dm, &dim);
2222: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2223: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2224: if (cMax < 0) return(0);
2225: if (label) {
2226: for (c = cMax; c < cEnd; ++c) {
2227: PetscInt val;
2229: DMLabelGetValue(label, c, &val);
2230: if (val == value) {
2231: ++(*numFaces);
2232: DMPlexGetConeSize(dm, c, &coneSize);
2233: }
2234: }
2235: } else {
2236: *numFaces = cEnd - cMax;
2237: DMPlexGetConeSize(dm, cMax, &coneSize);
2238: }
2239: PetscMalloc1(*numFaces *2, subCells);
2240: if (!(*numFaces)) return(0);
2241: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2242: for (c = cMax; c < cEnd; ++c) {
2243: const PetscInt *cells;
2244: PetscInt numCells;
2246: if (label) {
2247: PetscInt val;
2249: DMLabelGetValue(label, c, &val);
2250: if (val != value) continue;
2251: }
2252: DMPlexGetCone(dm, c, &cone);
2253: for (p = 0; p < *nFV; ++p) {
2254: DMLabelSetValue(subpointMap, cone[p], 0);
2255: }
2256: /* Negative face */
2257: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2258: /* Not true in parallel
2259: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2260: for (p = 0; p < numCells; ++p) {
2261: DMLabelSetValue(subpointMap, cells[p], 2);
2262: (*subCells)[subc++] = cells[p];
2263: }
2264: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2265: /* Positive face is not included */
2266: }
2267: return(0);
2268: }
2270: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2271: {
2272: PetscInt *pStart, *pEnd;
2273: PetscInt dim, cMax, cEnd, c, d;
2277: DMGetDimension(dm, &dim);
2278: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2279: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2280: if (cMax < 0) return(0);
2281: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2282: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2283: for (c = cMax; c < cEnd; ++c) {
2284: const PetscInt *cone;
2285: PetscInt *closure = NULL;
2286: PetscInt fconeSize, coneSize, closureSize, cl, val;
2288: if (label) {
2289: DMLabelGetValue(label, c, &val);
2290: if (val != value) continue;
2291: }
2292: DMPlexGetConeSize(dm, c, &coneSize);
2293: DMPlexGetCone(dm, c, &cone);
2294: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2295: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2296: /* Negative face */
2297: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2298: for (cl = 0; cl < closureSize*2; cl += 2) {
2299: const PetscInt point = closure[cl];
2301: for (d = 0; d <= dim; ++d) {
2302: if ((point >= pStart[d]) && (point < pEnd[d])) {
2303: DMLabelSetValue(subpointMap, point, d);
2304: break;
2305: }
2306: }
2307: }
2308: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2309: /* Cells -- positive face is not included */
2310: for (cl = 0; cl < 1; ++cl) {
2311: const PetscInt *support;
2312: PetscInt supportSize, s;
2314: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2315: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2316: DMPlexGetSupport(dm, cone[cl], &support);
2317: for (s = 0; s < supportSize; ++s) {
2318: DMLabelSetValue(subpointMap, support[s], dim);
2319: }
2320: }
2321: }
2322: PetscFree2(pStart, pEnd);
2323: return(0);
2324: }
2326: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2327: {
2328: MPI_Comm comm;
2329: PetscBool posOrient = PETSC_FALSE;
2330: const PetscInt debug = 0;
2331: PetscInt cellDim, faceSize, f;
2335: PetscObjectGetComm((PetscObject)dm,&comm);
2336: DMGetDimension(dm, &cellDim);
2337: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2339: if (cellDim == 1 && numCorners == 2) {
2340: /* Triangle */
2341: faceSize = numCorners-1;
2342: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2343: } else if (cellDim == 2 && numCorners == 3) {
2344: /* Triangle */
2345: faceSize = numCorners-1;
2346: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2347: } else if (cellDim == 3 && numCorners == 4) {
2348: /* Tetrahedron */
2349: faceSize = numCorners-1;
2350: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2351: } else if (cellDim == 1 && numCorners == 3) {
2352: /* Quadratic line */
2353: faceSize = 1;
2354: posOrient = PETSC_TRUE;
2355: } else if (cellDim == 2 && numCorners == 4) {
2356: /* Quads */
2357: faceSize = 2;
2358: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2359: posOrient = PETSC_TRUE;
2360: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2361: posOrient = PETSC_TRUE;
2362: } else {
2363: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2364: posOrient = PETSC_FALSE;
2365: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2366: }
2367: } else if (cellDim == 2 && numCorners == 6) {
2368: /* Quadratic triangle (I hate this) */
2369: /* Edges are determined by the first 2 vertices (corners of edges) */
2370: const PetscInt faceSizeTri = 3;
2371: PetscInt sortedIndices[3], i, iFace;
2372: PetscBool found = PETSC_FALSE;
2373: PetscInt faceVerticesTriSorted[9] = {
2374: 0, 3, 4, /* bottom */
2375: 1, 4, 5, /* right */
2376: 2, 3, 5, /* left */
2377: };
2378: PetscInt faceVerticesTri[9] = {
2379: 0, 3, 4, /* bottom */
2380: 1, 4, 5, /* right */
2381: 2, 5, 3, /* left */
2382: };
2384: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2385: PetscSortInt(faceSizeTri, sortedIndices);
2386: for (iFace = 0; iFace < 3; ++iFace) {
2387: const PetscInt ii = iFace*faceSizeTri;
2388: PetscInt fVertex, cVertex;
2390: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2391: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2392: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2393: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2394: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2395: faceVertices[fVertex] = origVertices[cVertex];
2396: break;
2397: }
2398: }
2399: }
2400: found = PETSC_TRUE;
2401: break;
2402: }
2403: }
2404: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2405: if (posOriented) *posOriented = PETSC_TRUE;
2406: return(0);
2407: } else if (cellDim == 2 && numCorners == 9) {
2408: /* Quadratic quad (I hate this) */
2409: /* Edges are determined by the first 2 vertices (corners of edges) */
2410: const PetscInt faceSizeQuad = 3;
2411: PetscInt sortedIndices[3], i, iFace;
2412: PetscBool found = PETSC_FALSE;
2413: PetscInt faceVerticesQuadSorted[12] = {
2414: 0, 1, 4, /* bottom */
2415: 1, 2, 5, /* right */
2416: 2, 3, 6, /* top */
2417: 0, 3, 7, /* left */
2418: };
2419: PetscInt faceVerticesQuad[12] = {
2420: 0, 1, 4, /* bottom */
2421: 1, 2, 5, /* right */
2422: 2, 3, 6, /* top */
2423: 3, 0, 7, /* left */
2424: };
2426: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2427: PetscSortInt(faceSizeQuad, sortedIndices);
2428: for (iFace = 0; iFace < 4; ++iFace) {
2429: const PetscInt ii = iFace*faceSizeQuad;
2430: PetscInt fVertex, cVertex;
2432: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2433: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2434: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2435: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2436: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2437: faceVertices[fVertex] = origVertices[cVertex];
2438: break;
2439: }
2440: }
2441: }
2442: found = PETSC_TRUE;
2443: break;
2444: }
2445: }
2446: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2447: if (posOriented) *posOriented = PETSC_TRUE;
2448: return(0);
2449: } else if (cellDim == 3 && numCorners == 8) {
2450: /* Hexes
2451: A hex is two oriented quads with the normal of the first
2452: pointing up at the second.
2454: 7---6
2455: /| /|
2456: 4---5 |
2457: | 1-|-2
2458: |/ |/
2459: 0---3
2461: Faces are determined by the first 4 vertices (corners of faces) */
2462: const PetscInt faceSizeHex = 4;
2463: PetscInt sortedIndices[4], i, iFace;
2464: PetscBool found = PETSC_FALSE;
2465: PetscInt faceVerticesHexSorted[24] = {
2466: 0, 1, 2, 3, /* bottom */
2467: 4, 5, 6, 7, /* top */
2468: 0, 3, 4, 5, /* front */
2469: 2, 3, 5, 6, /* right */
2470: 1, 2, 6, 7, /* back */
2471: 0, 1, 4, 7, /* left */
2472: };
2473: PetscInt faceVerticesHex[24] = {
2474: 1, 2, 3, 0, /* bottom */
2475: 4, 5, 6, 7, /* top */
2476: 0, 3, 5, 4, /* front */
2477: 3, 2, 6, 5, /* right */
2478: 2, 1, 7, 6, /* back */
2479: 1, 0, 4, 7, /* left */
2480: };
2482: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2483: PetscSortInt(faceSizeHex, sortedIndices);
2484: for (iFace = 0; iFace < 6; ++iFace) {
2485: const PetscInt ii = iFace*faceSizeHex;
2486: PetscInt fVertex, cVertex;
2488: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2489: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2490: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2491: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2492: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2493: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2494: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2495: faceVertices[fVertex] = origVertices[cVertex];
2496: break;
2497: }
2498: }
2499: }
2500: found = PETSC_TRUE;
2501: break;
2502: }
2503: }
2504: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2505: if (posOriented) *posOriented = PETSC_TRUE;
2506: return(0);
2507: } else if (cellDim == 3 && numCorners == 10) {
2508: /* Quadratic tet */
2509: /* Faces are determined by the first 3 vertices (corners of faces) */
2510: const PetscInt faceSizeTet = 6;
2511: PetscInt sortedIndices[6], i, iFace;
2512: PetscBool found = PETSC_FALSE;
2513: PetscInt faceVerticesTetSorted[24] = {
2514: 0, 1, 2, 6, 7, 8, /* bottom */
2515: 0, 3, 4, 6, 7, 9, /* front */
2516: 1, 4, 5, 7, 8, 9, /* right */
2517: 2, 3, 5, 6, 8, 9, /* left */
2518: };
2519: PetscInt faceVerticesTet[24] = {
2520: 0, 1, 2, 6, 7, 8, /* bottom */
2521: 0, 4, 3, 6, 7, 9, /* front */
2522: 1, 5, 4, 7, 8, 9, /* right */
2523: 2, 3, 5, 8, 6, 9, /* left */
2524: };
2526: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2527: PetscSortInt(faceSizeTet, sortedIndices);
2528: for (iFace=0; iFace < 4; ++iFace) {
2529: const PetscInt ii = iFace*faceSizeTet;
2530: PetscInt fVertex, cVertex;
2532: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2533: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2534: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2535: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2536: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2537: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2538: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2539: faceVertices[fVertex] = origVertices[cVertex];
2540: break;
2541: }
2542: }
2543: }
2544: found = PETSC_TRUE;
2545: break;
2546: }
2547: }
2548: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2549: if (posOriented) *posOriented = PETSC_TRUE;
2550: return(0);
2551: } else if (cellDim == 3 && numCorners == 27) {
2552: /* Quadratic hexes (I hate this)
2553: A hex is two oriented quads with the normal of the first
2554: pointing up at the second.
2556: 7---6
2557: /| /|
2558: 4---5 |
2559: | 3-|-2
2560: |/ |/
2561: 0---1
2563: Faces are determined by the first 4 vertices (corners of faces) */
2564: const PetscInt faceSizeQuadHex = 9;
2565: PetscInt sortedIndices[9], i, iFace;
2566: PetscBool found = PETSC_FALSE;
2567: PetscInt faceVerticesQuadHexSorted[54] = {
2568: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2569: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2570: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2571: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2572: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2573: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2574: };
2575: PetscInt faceVerticesQuadHex[54] = {
2576: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2577: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2578: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2579: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2580: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2581: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2582: };
2584: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2585: PetscSortInt(faceSizeQuadHex, sortedIndices);
2586: for (iFace = 0; iFace < 6; ++iFace) {
2587: const PetscInt ii = iFace*faceSizeQuadHex;
2588: PetscInt fVertex, cVertex;
2590: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2591: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2592: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2593: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2594: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2595: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2596: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2597: faceVertices[fVertex] = origVertices[cVertex];
2598: break;
2599: }
2600: }
2601: }
2602: found = PETSC_TRUE;
2603: break;
2604: }
2605: }
2606: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2607: if (posOriented) *posOriented = PETSC_TRUE;
2608: return(0);
2609: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2610: if (!posOrient) {
2611: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2612: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2613: } else {
2614: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2615: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2616: }
2617: if (posOriented) *posOriented = posOrient;
2618: return(0);
2619: }
2621: /*@
2622: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2623: in faceVertices. The orientation is such that the face normal points out of the cell
2625: Not collective
2627: Input Parameters:
2628: + dm - The original mesh
2629: . cell - The cell mesh point
2630: . faceSize - The number of vertices on the face
2631: . face - The face vertices
2632: . numCorners - The number of vertices on the cell
2633: . indices - Local numbering of face vertices in cell cone
2634: - origVertices - Original face vertices
2636: Output Parameter:
2637: + faceVertices - The face vertices properly oriented
2638: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2640: Level: developer
2642: .seealso: DMPlexGetCone()
2643: @*/
2644: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2645: {
2646: const PetscInt *cone = NULL;
2647: PetscInt coneSize, v, f, v2;
2648: PetscInt oppositeVertex = -1;
2649: PetscErrorCode ierr;
2652: DMPlexGetConeSize(dm, cell, &coneSize);
2653: DMPlexGetCone(dm, cell, &cone);
2654: for (v = 0, v2 = 0; v < coneSize; ++v) {
2655: PetscBool found = PETSC_FALSE;
2657: for (f = 0; f < faceSize; ++f) {
2658: if (face[f] == cone[v]) {
2659: found = PETSC_TRUE; break;
2660: }
2661: }
2662: if (found) {
2663: indices[v2] = v;
2664: origVertices[v2] = cone[v];
2665: ++v2;
2666: } else {
2667: oppositeVertex = v;
2668: }
2669: }
2670: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2671: return(0);
2672: }
2674: /*
2675: DMPlexInsertFace_Internal - Puts a face into the mesh
2677: Not collective
2679: Input Parameters:
2680: + dm - The DMPlex
2681: . numFaceVertex - The number of vertices in the face
2682: . faceVertices - The vertices in the face for dm
2683: . subfaceVertices - The vertices in the face for subdm
2684: . numCorners - The number of vertices in the cell
2685: . cell - A cell in dm containing the face
2686: . subcell - A cell in subdm containing the face
2687: . firstFace - First face in the mesh
2688: - newFacePoint - Next face in the mesh
2690: Output Parameters:
2691: . newFacePoint - Contains next face point number on input, updated on output
2693: Level: developer
2694: */
2695: 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)
2696: {
2697: MPI_Comm comm;
2698: DM_Plex *submesh = (DM_Plex*) subdm->data;
2699: const PetscInt *faces;
2700: PetscInt numFaces, coneSize;
2701: PetscErrorCode ierr;
2704: PetscObjectGetComm((PetscObject)dm,&comm);
2705: DMPlexGetConeSize(subdm, subcell, &coneSize);
2706: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2707: #if 0
2708: /* Cannot use this because support() has not been constructed yet */
2709: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2710: #else
2711: {
2712: PetscInt f;
2714: numFaces = 0;
2715: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2716: for (f = firstFace; f < *newFacePoint; ++f) {
2717: PetscInt dof, off, d;
2719: PetscSectionGetDof(submesh->coneSection, f, &dof);
2720: PetscSectionGetOffset(submesh->coneSection, f, &off);
2721: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2722: for (d = 0; d < dof; ++d) {
2723: const PetscInt p = submesh->cones[off+d];
2724: PetscInt v;
2726: for (v = 0; v < numFaceVertices; ++v) {
2727: if (subfaceVertices[v] == p) break;
2728: }
2729: if (v == numFaceVertices) break;
2730: }
2731: if (d == dof) {
2732: numFaces = 1;
2733: ((PetscInt*) faces)[0] = f;
2734: }
2735: }
2736: }
2737: #endif
2738: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2739: else if (numFaces == 1) {
2740: /* Add the other cell neighbor for this face */
2741: DMPlexSetCone(subdm, subcell, faces);
2742: } else {
2743: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2744: PetscBool posOriented;
2746: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2747: origVertices = &orientedVertices[numFaceVertices];
2748: indices = &orientedVertices[numFaceVertices*2];
2749: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2750: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2751: /* TODO: I know that routine should return a permutation, not the indices */
2752: for (v = 0; v < numFaceVertices; ++v) {
2753: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2754: for (ov = 0; ov < numFaceVertices; ++ov) {
2755: if (orientedVertices[ov] == vertex) {
2756: orientedSubVertices[ov] = subvertex;
2757: break;
2758: }
2759: }
2760: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2761: }
2762: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2763: DMPlexSetCone(subdm, subcell, newFacePoint);
2764: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2765: ++(*newFacePoint);
2766: }
2767: #if 0
2768: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2769: #else
2770: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2771: #endif
2772: return(0);
2773: }
2775: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2776: {
2777: MPI_Comm comm;
2778: DMLabel subpointMap;
2779: IS subvertexIS, subcellIS;
2780: const PetscInt *subVertices, *subCells;
2781: PetscInt numSubVertices, firstSubVertex, numSubCells;
2782: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2783: PetscInt vStart, vEnd, c, f;
2784: PetscErrorCode ierr;
2787: PetscObjectGetComm((PetscObject)dm,&comm);
2788: /* Create subpointMap which marks the submesh */
2789: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2790: DMPlexSetSubpointMap(subdm, subpointMap);
2791: DMLabelDestroy(&subpointMap);
2792: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2793: /* Setup chart */
2794: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2795: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2796: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2797: DMPlexSetVTKCellHeight(subdm, 1);
2798: /* Set cone sizes */
2799: firstSubVertex = numSubCells;
2800: firstSubFace = numSubCells+numSubVertices;
2801: newFacePoint = firstSubFace;
2802: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2803: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2804: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2805: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2806: for (c = 0; c < numSubCells; ++c) {
2807: DMPlexSetConeSize(subdm, c, 1);
2808: }
2809: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2810: DMPlexSetConeSize(subdm, f, nFV);
2811: }
2812: DMSetUp(subdm);
2813: /* Create face cones */
2814: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2815: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2816: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2817: for (c = 0; c < numSubCells; ++c) {
2818: const PetscInt cell = subCells[c];
2819: const PetscInt subcell = c;
2820: PetscInt *closure = NULL;
2821: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2823: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2824: for (cl = 0; cl < closureSize*2; cl += 2) {
2825: const PetscInt point = closure[cl];
2826: PetscInt subVertex;
2828: if ((point >= vStart) && (point < vEnd)) {
2829: ++numCorners;
2830: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2831: if (subVertex >= 0) {
2832: closure[faceSize] = point;
2833: subface[faceSize] = firstSubVertex+subVertex;
2834: ++faceSize;
2835: }
2836: }
2837: }
2838: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2839: if (faceSize == nFV) {
2840: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2841: }
2842: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2843: }
2844: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2845: DMPlexSymmetrize(subdm);
2846: DMPlexStratify(subdm);
2847: /* Build coordinates */
2848: {
2849: PetscSection coordSection, subCoordSection;
2850: Vec coordinates, subCoordinates;
2851: PetscScalar *coords, *subCoords;
2852: PetscInt numComp, coordSize, v;
2853: const char *name;
2855: DMGetCoordinateSection(dm, &coordSection);
2856: DMGetCoordinatesLocal(dm, &coordinates);
2857: DMGetCoordinateSection(subdm, &subCoordSection);
2858: PetscSectionSetNumFields(subCoordSection, 1);
2859: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2860: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2861: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2862: for (v = 0; v < numSubVertices; ++v) {
2863: const PetscInt vertex = subVertices[v];
2864: const PetscInt subvertex = firstSubVertex+v;
2865: PetscInt dof;
2867: PetscSectionGetDof(coordSection, vertex, &dof);
2868: PetscSectionSetDof(subCoordSection, subvertex, dof);
2869: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2870: }
2871: PetscSectionSetUp(subCoordSection);
2872: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2873: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2874: PetscObjectGetName((PetscObject)coordinates,&name);
2875: PetscObjectSetName((PetscObject)subCoordinates,name);
2876: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2877: VecSetType(subCoordinates,VECSTANDARD);
2878: if (coordSize) {
2879: VecGetArray(coordinates, &coords);
2880: VecGetArray(subCoordinates, &subCoords);
2881: for (v = 0; v < numSubVertices; ++v) {
2882: const PetscInt vertex = subVertices[v];
2883: const PetscInt subvertex = firstSubVertex+v;
2884: PetscInt dof, off, sdof, soff, d;
2886: PetscSectionGetDof(coordSection, vertex, &dof);
2887: PetscSectionGetOffset(coordSection, vertex, &off);
2888: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2889: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2890: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2891: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2892: }
2893: VecRestoreArray(coordinates, &coords);
2894: VecRestoreArray(subCoordinates, &subCoords);
2895: }
2896: DMSetCoordinatesLocal(subdm, subCoordinates);
2897: VecDestroy(&subCoordinates);
2898: }
2899: /* Cleanup */
2900: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2901: ISDestroy(&subvertexIS);
2902: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2903: ISDestroy(&subcellIS);
2904: return(0);
2905: }
2907: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2908: {
2909: PetscInt subPoint;
2912: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2913: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2914: }
2916: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2917: {
2918: MPI_Comm comm;
2919: DMLabel subpointMap;
2920: IS *subpointIS;
2921: const PetscInt **subpoints;
2922: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2923: PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2924: PetscMPIInt rank;
2925: PetscErrorCode ierr;
2928: PetscObjectGetComm((PetscObject)dm,&comm);
2929: MPI_Comm_rank(comm, &rank);
2930: /* Create subpointMap which marks the submesh */
2931: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2932: DMPlexSetSubpointMap(subdm, subpointMap);
2933: if (cellHeight) {
2934: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2935: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
2936: } else {
2937: DMLabel depth;
2938: IS pointIS;
2939: const PetscInt *points;
2940: PetscInt numPoints;
2942: DMPlexGetDepthLabel(dm, &depth);
2943: DMLabelGetStratumSize(label, value, &numPoints);
2944: DMLabelGetStratumIS(label, value, &pointIS);
2945: ISGetIndices(pointIS, &points);
2946: for (p = 0; p < numPoints; ++p) {
2947: PetscInt *closure = NULL;
2948: PetscInt closureSize, c, pdim;
2950: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2951: for (c = 0; c < closureSize*2; c += 2) {
2952: DMLabelGetValue(depth, closure[c], &pdim);
2953: DMLabelSetValue(subpointMap, closure[c], pdim);
2954: }
2955: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2956: }
2957: ISRestoreIndices(pointIS, &points);
2958: ISDestroy(&pointIS);
2959: }
2960: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2961: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2962: cMax = (cMax < 0) ? cEnd : cMax;
2963: /* Setup chart */
2964: DMGetDimension(dm, &dim);
2965: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2966: for (d = 0; d <= dim; ++d) {
2967: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2968: totSubPoints += numSubPoints[d];
2969: }
2970: DMPlexSetChart(subdm, 0, totSubPoints);
2971: DMPlexSetVTKCellHeight(subdm, cellHeight);
2972: /* Set cone sizes */
2973: firstSubPoint[dim] = 0;
2974: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2975: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2976: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2977: for (d = 0; d <= dim; ++d) {
2978: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2979: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2980: }
2981: for (d = 0; d <= dim; ++d) {
2982: for (p = 0; p < numSubPoints[d]; ++p) {
2983: const PetscInt point = subpoints[d][p];
2984: const PetscInt subpoint = firstSubPoint[d] + p;
2985: const PetscInt *cone;
2986: PetscInt coneSize, coneSizeNew, c, val;
2988: DMPlexGetConeSize(dm, point, &coneSize);
2989: DMPlexSetConeSize(subdm, subpoint, coneSize);
2990: if (cellHeight && (d == dim)) {
2991: DMPlexGetCone(dm, point, &cone);
2992: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2993: DMLabelGetValue(subpointMap, cone[c], &val);
2994: if (val >= 0) coneSizeNew++;
2995: }
2996: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2997: }
2998: }
2999: }
3000: DMLabelDestroy(&subpointMap);
3001: DMSetUp(subdm);
3002: /* Set cones */
3003: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3004: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3005: for (d = 0; d <= dim; ++d) {
3006: for (p = 0; p < numSubPoints[d]; ++p) {
3007: const PetscInt point = subpoints[d][p];
3008: const PetscInt subpoint = firstSubPoint[d] + p;
3009: const PetscInt *cone, *ornt;
3010: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3012: if (d == dim-1) {
3013: const PetscInt *support, *cone, *ornt;
3014: PetscInt supportSize, coneSize, s, subc;
3016: DMPlexGetSupport(dm, point, &support);
3017: DMPlexGetSupportSize(dm, point, &supportSize);
3018: for (s = 0; s < supportSize; ++s) {
3019: if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
3020: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3021: if (subc >= 0) {
3022: const PetscInt ccell = subpoints[d+1][subc];
3024: DMPlexGetCone(dm, ccell, &cone);
3025: DMPlexGetConeSize(dm, ccell, &coneSize);
3026: DMPlexGetConeOrientation(dm, ccell, &ornt);
3027: for (c = 0; c < coneSize; ++c) {
3028: if (cone[c] == point) {
3029: fornt = ornt[c];
3030: break;
3031: }
3032: }
3033: break;
3034: }
3035: }
3036: }
3037: DMPlexGetConeSize(dm, point, &coneSize);
3038: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3039: DMPlexGetCone(dm, point, &cone);
3040: DMPlexGetConeOrientation(dm, point, &ornt);
3041: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3042: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3043: if (subc >= 0) {
3044: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3045: orntNew[coneSizeNew] = ornt[c];
3046: ++coneSizeNew;
3047: }
3048: }
3049: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3050: if (fornt < 0) {
3051: /* This should be replaced by a call to DMPlexReverseCell() */
3052: #if 0
3053: DMPlexReverseCell(subdm, subpoint);
3054: #else
3055: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3056: PetscInt faceSize, tmp;
3058: tmp = coneNew[c];
3059: coneNew[c] = coneNew[coneSizeNew-1-c];
3060: coneNew[coneSizeNew-1-c] = tmp;
3061: DMPlexGetConeSize(dm, cone[c], &faceSize);
3062: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3063: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3064: orntNew[coneSizeNew-1-c] = tmp;
3065: }
3066: }
3067: DMPlexSetCone(subdm, subpoint, coneNew);
3068: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3069: #endif
3070: }
3071: }
3072: PetscFree2(coneNew,orntNew);
3073: DMPlexSymmetrize(subdm);
3074: DMPlexStratify(subdm);
3075: /* Build coordinates */
3076: {
3077: PetscSection coordSection, subCoordSection;
3078: Vec coordinates, subCoordinates;
3079: PetscScalar *coords, *subCoords;
3080: PetscInt cdim, numComp, coordSize;
3081: const char *name;
3083: DMGetCoordinateDim(dm, &cdim);
3084: DMGetCoordinateSection(dm, &coordSection);
3085: DMGetCoordinatesLocal(dm, &coordinates);
3086: DMGetCoordinateSection(subdm, &subCoordSection);
3087: PetscSectionSetNumFields(subCoordSection, 1);
3088: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3089: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3090: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3091: for (v = 0; v < numSubPoints[0]; ++v) {
3092: const PetscInt vertex = subpoints[0][v];
3093: const PetscInt subvertex = firstSubPoint[0]+v;
3094: PetscInt dof;
3096: PetscSectionGetDof(coordSection, vertex, &dof);
3097: PetscSectionSetDof(subCoordSection, subvertex, dof);
3098: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3099: }
3100: PetscSectionSetUp(subCoordSection);
3101: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3102: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3103: PetscObjectGetName((PetscObject)coordinates,&name);
3104: PetscObjectSetName((PetscObject)subCoordinates,name);
3105: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3106: VecSetBlockSize(subCoordinates, cdim);
3107: VecSetType(subCoordinates,VECSTANDARD);
3108: VecGetArray(coordinates, &coords);
3109: VecGetArray(subCoordinates, &subCoords);
3110: for (v = 0; v < numSubPoints[0]; ++v) {
3111: const PetscInt vertex = subpoints[0][v];
3112: const PetscInt subvertex = firstSubPoint[0]+v;
3113: PetscInt dof, off, sdof, soff, d;
3115: PetscSectionGetDof(coordSection, vertex, &dof);
3116: PetscSectionGetOffset(coordSection, vertex, &off);
3117: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3118: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3119: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3120: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3121: }
3122: VecRestoreArray(coordinates, &coords);
3123: VecRestoreArray(subCoordinates, &subCoords);
3124: DMSetCoordinatesLocal(subdm, subCoordinates);
3125: VecDestroy(&subCoordinates);
3126: }
3127: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3128: {
3129: PetscSF sfPoint, sfPointSub;
3130: IS subpIS;
3131: const PetscSFNode *remotePoints;
3132: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3133: const PetscInt *localPoints, *subpoints;
3134: PetscInt *slocalPoints;
3135: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3136: PetscMPIInt rank;
3138: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3139: DMGetPointSF(dm, &sfPoint);
3140: DMGetPointSF(subdm, &sfPointSub);
3141: DMPlexGetChart(dm, &pStart, &pEnd);
3142: DMPlexGetChart(subdm, NULL, &numSubroots);
3143: DMPlexCreateSubpointIS(subdm, &subpIS);
3144: if (subpIS) {
3145: ISGetIndices(subpIS, &subpoints);
3146: ISGetLocalSize(subpIS, &numSubpoints);
3147: }
3148: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3149: if (numRoots >= 0) {
3150: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3151: for (p = 0; p < pEnd-pStart; ++p) {
3152: newLocalPoints[p].rank = -2;
3153: newLocalPoints[p].index = -2;
3154: }
3155: /* Set subleaves */
3156: for (l = 0; l < numLeaves; ++l) {
3157: const PetscInt point = localPoints[l];
3158: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3160: if (subpoint < 0) continue;
3161: newLocalPoints[point-pStart].rank = rank;
3162: newLocalPoints[point-pStart].index = subpoint;
3163: ++numSubleaves;
3164: }
3165: /* Must put in owned subpoints */
3166: for (p = pStart; p < pEnd; ++p) {
3167: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3169: if (subpoint < 0) {
3170: newOwners[p-pStart].rank = -3;
3171: newOwners[p-pStart].index = -3;
3172: } else {
3173: newOwners[p-pStart].rank = rank;
3174: newOwners[p-pStart].index = subpoint;
3175: }
3176: }
3177: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3178: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3179: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3180: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3181: PetscMalloc1(numSubleaves, &slocalPoints);
3182: PetscMalloc1(numSubleaves, &sremotePoints);
3183: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3184: const PetscInt point = localPoints[l];
3185: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3187: if (subpoint < 0) continue;
3188: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3189: slocalPoints[sl] = subpoint;
3190: sremotePoints[sl].rank = newLocalPoints[point].rank;
3191: sremotePoints[sl].index = newLocalPoints[point].index;
3192: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3193: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3194: ++sl;
3195: }
3196: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3197: PetscFree2(newLocalPoints,newOwners);
3198: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3199: }
3200: if (subpIS) {
3201: ISRestoreIndices(subpIS, &subpoints);
3202: ISDestroy(&subpIS);
3203: }
3204: }
3205: /* Cleanup */
3206: for (d = 0; d <= dim; ++d) {
3207: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3208: ISDestroy(&subpointIS[d]);
3209: }
3210: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3211: return(0);
3212: }
3214: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3215: {
3219: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3220: return(0);
3221: }
3223: /*@
3224: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3226: Input Parameters:
3227: + dm - The original mesh
3228: . vertexLabel - The DMLabel marking points contained in the surface
3229: . value - The label value to use
3230: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3232: Output Parameter:
3233: . subdm - The surface mesh
3235: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3237: Level: developer
3239: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3240: @*/
3241: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3242: {
3243: PetscInt dim, cdim, depth;
3249: DMGetDimension(dm, &dim);
3250: DMPlexGetDepth(dm, &depth);
3251: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3252: DMSetType(*subdm, DMPLEX);
3253: DMSetDimension(*subdm, dim-1);
3254: DMGetCoordinateDim(dm, &cdim);
3255: DMSetCoordinateDim(*subdm, cdim);
3256: if (depth == dim) {
3257: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3258: } else {
3259: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3260: }
3261: return(0);
3262: }
3264: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3265: {
3266: MPI_Comm comm;
3267: DMLabel subpointMap;
3268: IS subvertexIS;
3269: const PetscInt *subVertices;
3270: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3271: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3272: PetscInt cMax, c, f;
3273: PetscErrorCode ierr;
3276: PetscObjectGetComm((PetscObject)dm, &comm);
3277: /* Create subpointMap which marks the submesh */
3278: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3279: DMPlexSetSubpointMap(subdm, subpointMap);
3280: DMLabelDestroy(&subpointMap);
3281: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3282: /* Setup chart */
3283: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3284: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3285: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3286: DMPlexSetVTKCellHeight(subdm, 1);
3287: /* Set cone sizes */
3288: firstSubVertex = numSubCells;
3289: firstSubFace = numSubCells+numSubVertices;
3290: newFacePoint = firstSubFace;
3291: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3292: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3293: for (c = 0; c < numSubCells; ++c) {
3294: DMPlexSetConeSize(subdm, c, 1);
3295: }
3296: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3297: DMPlexSetConeSize(subdm, f, nFV);
3298: }
3299: DMSetUp(subdm);
3300: /* Create face cones */
3301: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3302: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3303: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3304: for (c = 0; c < numSubCells; ++c) {
3305: const PetscInt cell = subCells[c];
3306: const PetscInt subcell = c;
3307: const PetscInt *cone, *cells;
3308: PetscInt numCells, subVertex, p, v;
3310: if (cell < cMax) continue;
3311: DMPlexGetCone(dm, cell, &cone);
3312: for (v = 0; v < nFV; ++v) {
3313: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3314: subface[v] = firstSubVertex+subVertex;
3315: }
3316: DMPlexSetCone(subdm, newFacePoint, subface);
3317: DMPlexSetCone(subdm, subcell, &newFacePoint);
3318: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3319: /* Not true in parallel
3320: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3321: for (p = 0; p < numCells; ++p) {
3322: PetscInt negsubcell;
3324: if (cells[p] >= cMax) continue;
3325: /* I know this is a crap search */
3326: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3327: if (subCells[negsubcell] == cells[p]) break;
3328: }
3329: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3330: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3331: }
3332: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3333: ++newFacePoint;
3334: }
3335: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3336: DMPlexSymmetrize(subdm);
3337: DMPlexStratify(subdm);
3338: /* Build coordinates */
3339: {
3340: PetscSection coordSection, subCoordSection;
3341: Vec coordinates, subCoordinates;
3342: PetscScalar *coords, *subCoords;
3343: PetscInt cdim, numComp, coordSize, v;
3344: const char *name;
3346: DMGetCoordinateDim(dm, &cdim);
3347: DMGetCoordinateSection(dm, &coordSection);
3348: DMGetCoordinatesLocal(dm, &coordinates);
3349: DMGetCoordinateSection(subdm, &subCoordSection);
3350: PetscSectionSetNumFields(subCoordSection, 1);
3351: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3352: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3353: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3354: for (v = 0; v < numSubVertices; ++v) {
3355: const PetscInt vertex = subVertices[v];
3356: const PetscInt subvertex = firstSubVertex+v;
3357: PetscInt dof;
3359: PetscSectionGetDof(coordSection, vertex, &dof);
3360: PetscSectionSetDof(subCoordSection, subvertex, dof);
3361: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3362: }
3363: PetscSectionSetUp(subCoordSection);
3364: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3365: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3366: PetscObjectGetName((PetscObject)coordinates,&name);
3367: PetscObjectSetName((PetscObject)subCoordinates,name);
3368: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3369: VecSetBlockSize(subCoordinates, cdim);
3370: VecSetType(subCoordinates,VECSTANDARD);
3371: VecGetArray(coordinates, &coords);
3372: VecGetArray(subCoordinates, &subCoords);
3373: for (v = 0; v < numSubVertices; ++v) {
3374: const PetscInt vertex = subVertices[v];
3375: const PetscInt subvertex = firstSubVertex+v;
3376: PetscInt dof, off, sdof, soff, d;
3378: PetscSectionGetDof(coordSection, vertex, &dof);
3379: PetscSectionGetOffset(coordSection, vertex, &off);
3380: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3381: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3382: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3383: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3384: }
3385: VecRestoreArray(coordinates, &coords);
3386: VecRestoreArray(subCoordinates, &subCoords);
3387: DMSetCoordinatesLocal(subdm, subCoordinates);
3388: VecDestroy(&subCoordinates);
3389: }
3390: /* Build SF */
3391: CHKMEMQ;
3392: {
3393: PetscSF sfPoint, sfPointSub;
3394: const PetscSFNode *remotePoints;
3395: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3396: const PetscInt *localPoints;
3397: PetscInt *slocalPoints;
3398: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3399: PetscMPIInt rank;
3401: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3402: DMGetPointSF(dm, &sfPoint);
3403: DMGetPointSF(subdm, &sfPointSub);
3404: DMPlexGetChart(dm, &pStart, &pEnd);
3405: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3406: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3407: if (numRoots >= 0) {
3408: /* Only vertices should be shared */
3409: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3410: for (p = 0; p < pEnd-pStart; ++p) {
3411: newLocalPoints[p].rank = -2;
3412: newLocalPoints[p].index = -2;
3413: }
3414: /* Set subleaves */
3415: for (l = 0; l < numLeaves; ++l) {
3416: const PetscInt point = localPoints[l];
3417: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3419: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3420: if (subPoint < 0) continue;
3421: newLocalPoints[point-pStart].rank = rank;
3422: newLocalPoints[point-pStart].index = subPoint;
3423: ++numSubLeaves;
3424: }
3425: /* Must put in owned subpoints */
3426: for (p = pStart; p < pEnd; ++p) {
3427: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3429: if (subPoint < 0) {
3430: newOwners[p-pStart].rank = -3;
3431: newOwners[p-pStart].index = -3;
3432: } else {
3433: newOwners[p-pStart].rank = rank;
3434: newOwners[p-pStart].index = subPoint;
3435: }
3436: }
3437: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3438: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3439: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3440: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3441: PetscMalloc1(numSubLeaves, &slocalPoints);
3442: PetscMalloc1(numSubLeaves, &sremotePoints);
3443: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3444: const PetscInt point = localPoints[l];
3445: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3447: if (subPoint < 0) continue;
3448: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3449: slocalPoints[sl] = subPoint;
3450: sremotePoints[sl].rank = newLocalPoints[point].rank;
3451: sremotePoints[sl].index = newLocalPoints[point].index;
3452: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3453: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3454: ++sl;
3455: }
3456: PetscFree2(newLocalPoints,newOwners);
3457: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3458: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3459: }
3460: }
3461: CHKMEMQ;
3462: /* Cleanup */
3463: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3464: ISDestroy(&subvertexIS);
3465: PetscFree(subCells);
3466: return(0);
3467: }
3469: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3470: {
3471: DMLabel label = NULL;
3475: if (labelname) {DMGetLabel(dm, labelname, &label);}
3476: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3477: return(0);
3478: }
3480: /*@C
3481: 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.
3483: Input Parameters:
3484: + dm - The original mesh
3485: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3486: . label - A label name, or NULL
3487: - value - A label value
3489: Output Parameter:
3490: . subdm - The surface mesh
3492: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3494: Level: developer
3496: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3497: @*/
3498: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3499: {
3500: PetscInt dim, cdim, depth;
3506: DMGetDimension(dm, &dim);
3507: DMPlexGetDepth(dm, &depth);
3508: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3509: DMSetType(*subdm, DMPLEX);
3510: DMSetDimension(*subdm, dim-1);
3511: DMGetCoordinateDim(dm, &cdim);
3512: DMSetCoordinateDim(*subdm, cdim);
3513: if (depth == dim) {
3514: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3515: } else {
3516: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3517: }
3518: return(0);
3519: }
3521: /*@
3522: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3524: Input Parameters:
3525: + dm - The original mesh
3526: . cellLabel - The DMLabel marking cells contained in the new mesh
3527: - value - The label value to use
3529: Output Parameter:
3530: . subdm - The new mesh
3532: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3534: Level: developer
3536: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3537: @*/
3538: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3539: {
3540: PetscInt dim;
3546: DMGetDimension(dm, &dim);
3547: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3548: DMSetType(*subdm, DMPLEX);
3549: DMSetDimension(*subdm, dim);
3550: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3551: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3552: return(0);
3553: }
3555: /*@
3556: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3558: Input Parameter:
3559: . dm - The submesh DM
3561: Output Parameter:
3562: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3564: Level: developer
3566: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3567: @*/
3568: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3569: {
3573: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3574: return(0);
3575: }
3577: /*@
3578: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3580: Input Parameters:
3581: + dm - The submesh DM
3582: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3584: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3586: Level: developer
3588: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3589: @*/
3590: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3591: {
3592: DM_Plex *mesh = (DM_Plex *) dm->data;
3593: DMLabel tmp;
3598: tmp = mesh->subpointMap;
3599: mesh->subpointMap = subpointMap;
3600: PetscObjectReference((PetscObject) mesh->subpointMap);
3601: DMLabelDestroy(&tmp);
3602: return(0);
3603: }
3605: /*@
3606: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3608: Input Parameter:
3609: . dm - The submesh DM
3611: Output Parameter:
3612: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3614: Note: This IS is guaranteed to be sorted by the construction of the submesh
3616: Level: developer
3618: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3619: @*/
3620: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3621: {
3622: MPI_Comm comm;
3623: DMLabel subpointMap;
3624: IS is;
3625: const PetscInt *opoints;
3626: PetscInt *points, *depths;
3627: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3628: PetscErrorCode ierr;
3633: PetscObjectGetComm((PetscObject)dm,&comm);
3634: *subpointIS = NULL;
3635: DMPlexGetSubpointMap(dm, &subpointMap);
3636: DMPlexGetDepth(dm, &depth);
3637: if (subpointMap && depth >= 0) {
3638: DMPlexGetChart(dm, &pStart, &pEnd);
3639: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3640: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3641: depths[0] = depth;
3642: depths[1] = 0;
3643: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3644: PetscMalloc1(pEnd, &points);
3645: for(d = 0, off = 0; d <= depth; ++d) {
3646: const PetscInt dep = depths[d];
3648: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3649: DMLabelGetStratumSize(subpointMap, dep, &n);
3650: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3651: 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);
3652: } else {
3653: if (!n) {
3654: if (d == 0) {
3655: /* Missing cells */
3656: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3657: } else {
3658: /* Missing faces */
3659: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3660: }
3661: }
3662: }
3663: if (n) {
3664: DMLabelGetStratumIS(subpointMap, dep, &is);
3665: ISGetIndices(is, &opoints);
3666: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3667: ISRestoreIndices(is, &opoints);
3668: ISDestroy(&is);
3669: }
3670: }
3671: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3672: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3673: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3674: }
3675: return(0);
3676: }
3678: /*@
3679: DMPlexGetSubpoint - Return the subpoint corresponding to a point in the original mesh. If the DM
3680: is not a submesh, just return the input point.
3682: Note collective
3684: Input Parameters:
3685: + dm - The submesh DM
3686: - p - The point in the original, from which the submesh was created
3688: Output Parameter:
3689: . subp - The point in the submesh
3691: Level: developer
3693: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3694: @*/
3695: PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp)
3696: {
3697: DMLabel spmap;
3701: *subp = p;
3702: DMPlexGetSubpointMap(dm, &spmap);
3703: if (spmap) {
3704: IS subpointIS;
3705: const PetscInt *subpoints;
3706: PetscInt numSubpoints;
3708: /* TODO Cache the IS, making it look like an index */
3709: DMPlexCreateSubpointIS(dm, &subpointIS);
3710: ISGetLocalSize(subpointIS, &numSubpoints);
3711: ISGetIndices(subpointIS, &subpoints);
3712: PetscFindInt(p, numSubpoints, subpoints, subp);
3713: if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p);
3714: ISRestoreIndices(subpointIS, &subpoints);
3715: ISDestroy(&subpointIS);
3716: }
3717: return(0);
3718: }
3720: /*@
3721: DMPlexGetAuxiliaryPoint - For a given point in the DM, return the matching point in the auxiliary DM.
3723: Note collective
3725: Input Parameters:
3726: + dm - The DM
3727: . dmAux - The related auxiliary DM
3728: - p - The point in the original DM
3730: Output Parameter:
3731: . subp - The point in the auxiliary DM
3733: Notes: If the DM is a submesh, we assume the dmAux is as well and just return the point. If only dmAux is a submesh,
3734: then we map the point back to the original space.
3736: Level: developer
3738: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3739: @*/
3740: PetscErrorCode DMPlexGetAuxiliaryPoint(DM dm, DM dmAux, PetscInt p, PetscInt *subp)
3741: {
3742: DMLabel spmap;
3746: *subp = p;
3747: /* If dm is a submesh, do not get subpoint */
3748: DMPlexGetSubpointMap(dm, &spmap);
3749: if (dmAux && !spmap) {
3750: PetscInt h;
3752: DMPlexGetVTKCellHeight(dmAux, &h);
3753: DMPlexGetSubpointMap(dmAux, &spmap);
3754: if (spmap && !h) {DMLabelGetValue(spmap, p, subp);}
3755: else {DMPlexGetSubpoint(dmAux, p, subp);}
3756: }
3757: return(0);
3758: }