Actual source code: plexsubmesh.c
petsc-3.11.4 2019-09-28
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: if (0) {
615: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
616: }
617: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
618: for (f = fStart; f < fEnd; ++f) {
619: PetscInt numCells;
621: DMPlexGetSupportSize(dmNew, f, &numCells);
622: if (numCells < 2) {
623: DMLabelSetValue(ghostLabel, f, 1);
624: } else {
625: const PetscInt *cells = NULL;
626: PetscInt vA, vB;
628: DMPlexGetSupport(dmNew, f, &cells);
629: DMLabelGetValue(vtkLabel, cells[0], &vA);
630: DMLabelGetValue(vtkLabel, cells[1], &vB);
631: if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
632: }
633: }
634: if (0) {
635: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
636: }
637: return(0);
638: }
640: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
641: {
642: DM refTree;
643: PetscSection pSec;
644: PetscInt *parents, *childIDs;
648: DMPlexGetReferenceTree(dm,&refTree);
649: DMPlexSetReferenceTree(dmNew,refTree);
650: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
651: if (pSec) {
652: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
653: PetscInt *childIDsShifted;
654: PetscSection pSecShifted;
656: PetscSectionGetChart(pSec,&pStart,&pEnd);
657: DMPlexGetDepth(dm,&depth);
658: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
659: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
660: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
661: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
662: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
663: for (p = pStartShifted; p < pEndShifted; p++) {
664: /* start off assuming no children */
665: PetscSectionSetDof(pSecShifted,p,0);
666: }
667: for (p = pStart; p < pEnd; p++) {
668: PetscInt dof;
669: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
671: PetscSectionGetDof(pSec,p,&dof);
672: PetscSectionSetDof(pSecShifted,pNew,dof);
673: }
674: PetscSectionSetUp(pSecShifted);
675: for (p = pStart; p < pEnd; p++) {
676: PetscInt dof;
677: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
679: PetscSectionGetDof(pSec,p,&dof);
680: if (dof) {
681: PetscInt off, offNew;
683: PetscSectionGetOffset(pSec,p,&off);
684: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
685: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
686: childIDsShifted[offNew] = childIDs[off];
687: }
688: }
689: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
690: PetscFree2(parentsShifted,childIDsShifted);
691: PetscSectionDestroy(&pSecShifted);
692: }
693: return(0);
694: }
696: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
697: {
698: PetscSF sf;
699: IS valueIS;
700: const PetscInt *values, *leaves;
701: PetscInt *depthShift;
702: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
703: PetscBool isper;
704: const PetscReal *maxCell, *L;
705: const DMBoundaryType *bd;
706: PetscErrorCode ierr;
709: DMGetPointSF(dm, &sf);
710: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
711: nleaves = PetscMax(0, nleaves);
712: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
713: /* Count ghost cells */
714: DMLabelGetValueIS(label, &valueIS);
715: ISGetLocalSize(valueIS, &numFS);
716: ISGetIndices(valueIS, &values);
717: Ng = 0;
718: for (fs = 0; fs < numFS; ++fs) {
719: IS faceIS;
720: const PetscInt *faces;
721: PetscInt numFaces, f, numBdFaces = 0;
723: DMLabelGetStratumIS(label, values[fs], &faceIS);
724: ISGetLocalSize(faceIS, &numFaces);
725: ISGetIndices(faceIS, &faces);
726: for (f = 0; f < numFaces; ++f) {
727: PetscInt numChildren;
729: PetscFindInt(faces[f], nleaves, leaves, &loc);
730: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
731: /* non-local and ancestors points don't get to register ghosts */
732: if (loc >= 0 || numChildren) continue;
733: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
734: }
735: Ng += numBdFaces;
736: ISDestroy(&faceIS);
737: }
738: DMPlexGetDepth(dm, &depth);
739: PetscMalloc1(2*(depth+1), &depthShift);
740: for (d = 0; d <= depth; d++) {
741: PetscInt dEnd;
743: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
744: depthShift[2*d] = dEnd;
745: depthShift[2*d+1] = 0;
746: }
747: if (depth >= 0) depthShift[2*depth+1] = Ng;
748: DMPlexShiftPointSetUp_Internal(depth,depthShift);
749: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
750: /* Step 3: Set cone/support sizes for new points */
751: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
752: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
753: for (c = cEnd; c < cEnd + Ng; ++c) {
754: DMPlexSetConeSize(gdm, c, 1);
755: }
756: for (fs = 0; fs < numFS; ++fs) {
757: IS faceIS;
758: const PetscInt *faces;
759: PetscInt numFaces, f;
761: DMLabelGetStratumIS(label, values[fs], &faceIS);
762: ISGetLocalSize(faceIS, &numFaces);
763: ISGetIndices(faceIS, &faces);
764: for (f = 0; f < numFaces; ++f) {
765: PetscInt size, numChildren;
767: PetscFindInt(faces[f], nleaves, leaves, &loc);
768: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
769: if (loc >= 0 || numChildren) continue;
770: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
771: DMPlexGetSupportSize(dm, faces[f], &size);
772: if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
773: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
774: }
775: ISRestoreIndices(faceIS, &faces);
776: ISDestroy(&faceIS);
777: }
778: /* Step 4: Setup ghosted DM */
779: DMSetUp(gdm);
780: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
781: /* Step 6: Set cones and supports for new points */
782: ghostCell = cEnd;
783: for (fs = 0; fs < numFS; ++fs) {
784: IS faceIS;
785: const PetscInt *faces;
786: PetscInt numFaces, f;
788: DMLabelGetStratumIS(label, values[fs], &faceIS);
789: ISGetLocalSize(faceIS, &numFaces);
790: ISGetIndices(faceIS, &faces);
791: for (f = 0; f < numFaces; ++f) {
792: PetscInt newFace = faces[f] + Ng, numChildren;
794: PetscFindInt(faces[f], nleaves, leaves, &loc);
795: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
796: if (loc >= 0 || numChildren) continue;
797: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
798: DMPlexSetCone(gdm, ghostCell, &newFace);
799: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
800: ++ghostCell;
801: }
802: ISRestoreIndices(faceIS, &faces);
803: ISDestroy(&faceIS);
804: }
805: ISRestoreIndices(valueIS, &values);
806: ISDestroy(&valueIS);
807: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
808: DMPlexShiftSF_Internal(dm, depthShift, gdm);
809: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
810: DMPlexShiftTree_Internal(dm, depthShift, gdm);
811: PetscFree(depthShift);
812: /* Step 7: Periodicity */
813: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
814: DMSetPeriodicity(gdm, isper, maxCell, L, bd);
815: if (numGhostCells) *numGhostCells = Ng;
816: return(0);
817: }
819: /*@C
820: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
822: Collective on dm
824: Input Parameters:
825: + dm - The original DM
826: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
828: Output Parameters:
829: + numGhostCells - The number of ghost cells added to the DM
830: - dmGhosted - The new DM
832: Note: If no label exists of that name, one will be created marking all boundary faces
834: Level: developer
836: .seealso: DMCreate()
837: @*/
838: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
839: {
840: DM gdm;
841: DMLabel label;
842: const char *name = labelName ? labelName : "Face Sets";
843: PetscInt dim;
844: PetscBool useCone, useClosure;
851: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
852: DMSetType(gdm, DMPLEX);
853: DMGetDimension(dm, &dim);
854: DMSetDimension(gdm, dim);
855: DMGetBasicAdjacency(dm, &useCone, &useClosure);
856: DMSetBasicAdjacency(gdm, useCone, useClosure);
857: DMGetLabel(dm, name, &label);
858: if (!label) {
859: /* Get label for boundary faces */
860: DMCreateLabel(dm, name);
861: DMGetLabel(dm, name, &label);
862: DMPlexMarkBoundaryFaces(dm, 1, label);
863: }
864: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
865: DMCopyBoundary(dm, gdm);
866: DMCopyDisc(dm, gdm);
867: gdm->setfromoptionscalled = dm->setfromoptionscalled;
868: *dmGhosted = gdm;
869: return(0);
870: }
872: /*
873: We are adding three kinds of points here:
874: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
875: Non-replicated: Points which exist on the fault, but are not replicated
876: Hybrid: Entirely new points, such as cohesive cells
878: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
879: each depth so that the new split/hybrid points can be inserted as a block.
880: */
881: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
882: {
883: MPI_Comm comm;
884: IS valueIS;
885: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
886: const PetscInt *values; /* List of depths for which we have replicated points */
887: IS *splitIS;
888: IS *unsplitIS;
889: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
890: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
891: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
892: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
893: const PetscInt **splitPoints; /* Replicated points for each depth */
894: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
895: PetscSection coordSection;
896: Vec coordinates;
897: PetscScalar *coords;
898: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
899: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
900: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
901: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
902: PetscInt *coneNew, *coneONew, *supportNew;
903: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
904: PetscErrorCode ierr;
907: PetscObjectGetComm((PetscObject)dm,&comm);
908: DMGetDimension(dm, &dim);
909: DMPlexGetDepth(dm, &depth);
910: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
911: /* Count split points and add cohesive cells */
912: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
913: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
914: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
915: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
916: for (d = 0; d <= depth; ++d) {
917: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
918: depthEnd[d] = pMaxNew[d];
919: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
920: numSplitPoints[d] = 0;
921: numUnsplitPoints[d] = 0;
922: numHybridPoints[d] = 0;
923: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
924: splitPoints[d] = NULL;
925: unsplitPoints[d] = NULL;
926: splitIS[d] = NULL;
927: unsplitIS[d] = NULL;
928: /* we are shifting the existing hybrid points with the stratum behind them, so
929: * the split comes at the end of the normal points, i.e., at depthMax[d] */
930: depthShift[2*d] = depthMax[d];
931: depthShift[2*d+1] = 0;
932: }
933: if (label) {
934: DMLabelGetValueIS(label, &valueIS);
935: ISGetLocalSize(valueIS, &numSP);
936: ISGetIndices(valueIS, &values);
937: }
938: for (sp = 0; sp < numSP; ++sp) {
939: const PetscInt dep = values[sp];
941: if ((dep < 0) || (dep > depth)) continue;
942: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
943: if (splitIS[dep]) {
944: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
945: ISGetIndices(splitIS[dep], &splitPoints[dep]);
946: }
947: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
948: if (unsplitIS[dep]) {
949: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
950: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
951: }
952: }
953: /* Calculate number of hybrid points */
954: 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 */
955: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
956: DMPlexShiftPointSetUp_Internal(depth,depthShift);
957: /* the end of the points in this stratum that come before the new points:
958: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
959: * added points */
960: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
961: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
962: /* Step 3: Set cone/support sizes for new points */
963: for (dep = 0; dep <= depth; ++dep) {
964: for (p = 0; p < numSplitPoints[dep]; ++p) {
965: const PetscInt oldp = splitPoints[dep][p];
966: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
967: const PetscInt splitp = p + pMaxNew[dep];
968: const PetscInt *support;
969: PetscInt coneSize, supportSize, qf, qn, qp, e;
971: DMPlexGetConeSize(dm, oldp, &coneSize);
972: DMPlexSetConeSize(sdm, splitp, coneSize);
973: DMPlexGetSupportSize(dm, oldp, &supportSize);
974: DMPlexSetSupportSize(sdm, splitp, supportSize);
975: if (dep == depth-1) {
976: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
978: /* Add cohesive cells, they are prisms */
979: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
980: } else if (dep == 0) {
981: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
983: DMPlexGetSupport(dm, oldp, &support);
984: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
985: PetscInt val;
987: DMLabelGetValue(label, support[e], &val);
988: if (val == 1) ++qf;
989: if ((val == 1) || (val == (shift + 1))) ++qn;
990: if ((val == 1) || (val == -(shift + 1))) ++qp;
991: }
992: /* Split old vertex: Edges into original vertex and new cohesive edge */
993: DMPlexSetSupportSize(sdm, newp, qn+1);
994: /* Split new vertex: Edges into split vertex and new cohesive edge */
995: DMPlexSetSupportSize(sdm, splitp, qp+1);
996: /* Add hybrid edge */
997: DMPlexSetConeSize(sdm, hybedge, 2);
998: DMPlexSetSupportSize(sdm, hybedge, qf);
999: } else if (dep == dim-2) {
1000: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1002: DMPlexGetSupport(dm, oldp, &support);
1003: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1004: PetscInt val;
1006: DMLabelGetValue(label, support[e], &val);
1007: if (val == dim-1) ++qf;
1008: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
1009: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
1010: }
1011: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1012: DMPlexSetSupportSize(sdm, newp, qn+1);
1013: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1014: DMPlexSetSupportSize(sdm, splitp, qp+1);
1015: /* Add hybrid face */
1016: DMPlexSetConeSize(sdm, hybface, 4);
1017: DMPlexSetSupportSize(sdm, hybface, qf);
1018: }
1019: }
1020: }
1021: for (dep = 0; dep <= depth; ++dep) {
1022: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1023: const PetscInt oldp = unsplitPoints[dep][p];
1024: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1025: const PetscInt *support;
1026: PetscInt coneSize, supportSize, qf, e, s;
1028: DMPlexGetConeSize(dm, oldp, &coneSize);
1029: DMPlexGetSupportSize(dm, oldp, &supportSize);
1030: DMPlexGetSupport(dm, oldp, &support);
1031: if (dep == 0) {
1032: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1034: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1035: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1036: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1037: if (e >= 0) ++qf;
1038: }
1039: DMPlexSetSupportSize(sdm, newp, qf+2);
1040: /* Add hybrid edge */
1041: DMPlexSetConeSize(sdm, hybedge, 2);
1042: for (e = 0, qf = 0; e < supportSize; ++e) {
1043: PetscInt val;
1045: DMLabelGetValue(label, support[e], &val);
1046: /* Split and unsplit edges produce hybrid faces */
1047: if (val == 1) ++qf;
1048: if (val == (shift2 + 1)) ++qf;
1049: }
1050: DMPlexSetSupportSize(sdm, hybedge, qf);
1051: } else if (dep == dim-2) {
1052: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1053: PetscInt val;
1055: for (e = 0, qf = 0; e < supportSize; ++e) {
1056: DMLabelGetValue(label, support[e], &val);
1057: if (val == dim-1) qf += 2;
1058: else ++qf;
1059: }
1060: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1061: DMPlexSetSupportSize(sdm, newp, qf+2);
1062: /* Add hybrid face */
1063: for (e = 0, qf = 0; e < supportSize; ++e) {
1064: DMLabelGetValue(label, support[e], &val);
1065: if (val == dim-1) ++qf;
1066: }
1067: DMPlexSetConeSize(sdm, hybface, 4);
1068: DMPlexSetSupportSize(sdm, hybface, qf);
1069: }
1070: }
1071: }
1072: /* Step 4: Setup split DM */
1073: DMSetUp(sdm);
1074: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1075: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1076: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1077: /* Step 6: Set cones and supports for new points */
1078: for (dep = 0; dep <= depth; ++dep) {
1079: for (p = 0; p < numSplitPoints[dep]; ++p) {
1080: const PetscInt oldp = splitPoints[dep][p];
1081: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1082: const PetscInt splitp = p + pMaxNew[dep];
1083: const PetscInt *cone, *support, *ornt;
1084: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1086: DMPlexGetConeSize(dm, oldp, &coneSize);
1087: DMPlexGetCone(dm, oldp, &cone);
1088: DMPlexGetConeOrientation(dm, oldp, &ornt);
1089: DMPlexGetSupportSize(dm, oldp, &supportSize);
1090: DMPlexGetSupport(dm, oldp, &support);
1091: if (dep == depth-1) {
1092: PetscBool hasUnsplit = PETSC_FALSE;
1093: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1094: const PetscInt *supportF;
1096: /* Split face: copy in old face to new face to start */
1097: DMPlexGetSupport(sdm, newp, &supportF);
1098: DMPlexSetSupport(sdm, splitp, supportF);
1099: /* Split old face: old vertices/edges in cone so no change */
1100: /* Split new face: new vertices/edges in cone */
1101: for (q = 0; q < coneSize; ++q) {
1102: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1103: if (v < 0) {
1104: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1105: 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);
1106: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1107: hasUnsplit = PETSC_TRUE;
1108: } else {
1109: coneNew[2+q] = v + pMaxNew[dep-1];
1110: if (dep > 1) {
1111: const PetscInt *econe;
1112: PetscInt econeSize, r, vs, vu;
1114: DMPlexGetConeSize(dm, cone[q], &econeSize);
1115: DMPlexGetCone(dm, cone[q], &econe);
1116: for (r = 0; r < econeSize; ++r) {
1117: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1118: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1119: if (vs >= 0) continue;
1120: 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);
1121: hasUnsplit = PETSC_TRUE;
1122: }
1123: }
1124: }
1125: }
1126: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1127: DMPlexSetConeOrientation(sdm, splitp, ornt);
1128: /* Face support */
1129: for (s = 0; s < supportSize; ++s) {
1130: PetscInt val;
1132: DMLabelGetValue(label, support[s], &val);
1133: if (val < 0) {
1134: /* Split old face: Replace negative side cell with cohesive cell */
1135: DMPlexInsertSupport(sdm, newp, s, hybcell);
1136: } else {
1137: /* Split new face: Replace positive side cell with cohesive cell */
1138: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1139: /* Get orientation for cohesive face */
1140: {
1141: const PetscInt *ncone, *nconeO;
1142: PetscInt nconeSize, nc;
1144: DMPlexGetConeSize(dm, support[s], &nconeSize);
1145: DMPlexGetCone(dm, support[s], &ncone);
1146: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1147: for (nc = 0; nc < nconeSize; ++nc) {
1148: if (ncone[nc] == oldp) {
1149: coneONew[0] = nconeO[nc];
1150: break;
1151: }
1152: }
1153: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1154: }
1155: }
1156: }
1157: /* Cohesive cell: Old and new split face, then new cohesive faces */
1158: coneNew[0] = newp; /* Extracted negative side orientation above */
1159: coneNew[1] = splitp;
1160: coneONew[1] = coneONew[0];
1161: for (q = 0; q < coneSize; ++q) {
1162: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1163: if (v < 0) {
1164: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1165: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1166: coneONew[2+q] = 0;
1167: } else {
1168: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1169: }
1170: coneONew[2+q] = 0;
1171: }
1172: DMPlexSetCone(sdm, hybcell, coneNew);
1173: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1174: /* Label the hybrid cells on the boundary of the split */
1175: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1176: } else if (dep == 0) {
1177: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1179: /* Split old vertex: Edges in old split faces and new cohesive edge */
1180: for (e = 0, qn = 0; e < supportSize; ++e) {
1181: PetscInt val;
1183: DMLabelGetValue(label, support[e], &val);
1184: if ((val == 1) || (val == (shift + 1))) {
1185: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1186: }
1187: }
1188: supportNew[qn] = hybedge;
1189: DMPlexSetSupport(sdm, newp, supportNew);
1190: /* Split new vertex: Edges in new split faces and new cohesive edge */
1191: for (e = 0, qp = 0; e < supportSize; ++e) {
1192: PetscInt val, edge;
1194: DMLabelGetValue(label, support[e], &val);
1195: if (val == 1) {
1196: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1197: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1198: supportNew[qp++] = edge + pMaxNew[dep+1];
1199: } else if (val == -(shift + 1)) {
1200: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1201: }
1202: }
1203: supportNew[qp] = hybedge;
1204: DMPlexSetSupport(sdm, splitp, supportNew);
1205: /* Hybrid edge: Old and new split vertex */
1206: coneNew[0] = newp;
1207: coneNew[1] = splitp;
1208: DMPlexSetCone(sdm, hybedge, coneNew);
1209: for (e = 0, qf = 0; e < supportSize; ++e) {
1210: PetscInt val, edge;
1212: DMLabelGetValue(label, support[e], &val);
1213: if (val == 1) {
1214: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1215: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1216: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1217: }
1218: }
1219: DMPlexSetSupport(sdm, hybedge, supportNew);
1220: } else if (dep == dim-2) {
1221: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1223: /* Split old edge: old vertices in cone so no change */
1224: /* Split new edge: new vertices in cone */
1225: for (q = 0; q < coneSize; ++q) {
1226: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1227: if (v < 0) {
1228: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1229: 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);
1230: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1231: } else {
1232: coneNew[q] = v + pMaxNew[dep-1];
1233: }
1234: }
1235: DMPlexSetCone(sdm, splitp, coneNew);
1236: /* Split old edge: Faces in positive side cells and old split faces */
1237: for (e = 0, q = 0; e < supportSize; ++e) {
1238: PetscInt val;
1240: DMLabelGetValue(label, support[e], &val);
1241: if (val == dim-1) {
1242: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1243: } else if (val == (shift + dim-1)) {
1244: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1245: }
1246: }
1247: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1248: DMPlexSetSupport(sdm, newp, supportNew);
1249: /* Split new edge: Faces in negative side cells and new split faces */
1250: for (e = 0, q = 0; e < supportSize; ++e) {
1251: PetscInt val, face;
1253: DMLabelGetValue(label, support[e], &val);
1254: if (val == dim-1) {
1255: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1256: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1257: supportNew[q++] = face + pMaxNew[dep+1];
1258: } else if (val == -(shift + dim-1)) {
1259: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1260: }
1261: }
1262: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1263: DMPlexSetSupport(sdm, splitp, supportNew);
1264: /* Hybrid face */
1265: coneNew[0] = newp;
1266: coneNew[1] = splitp;
1267: for (v = 0; v < coneSize; ++v) {
1268: PetscInt vertex;
1269: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1270: if (vertex < 0) {
1271: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1272: 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);
1273: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1274: } else {
1275: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1276: }
1277: }
1278: DMPlexSetCone(sdm, hybface, coneNew);
1279: for (e = 0, qf = 0; e < supportSize; ++e) {
1280: PetscInt val, face;
1282: DMLabelGetValue(label, support[e], &val);
1283: if (val == dim-1) {
1284: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1285: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1286: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1287: }
1288: }
1289: DMPlexSetSupport(sdm, hybface, supportNew);
1290: }
1291: }
1292: }
1293: for (dep = 0; dep <= depth; ++dep) {
1294: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1295: const PetscInt oldp = unsplitPoints[dep][p];
1296: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1297: const PetscInt *cone, *support, *ornt;
1298: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1300: DMPlexGetConeSize(dm, oldp, &coneSize);
1301: DMPlexGetCone(dm, oldp, &cone);
1302: DMPlexGetConeOrientation(dm, oldp, &ornt);
1303: DMPlexGetSupportSize(dm, oldp, &supportSize);
1304: DMPlexGetSupport(dm, oldp, &support);
1305: if (dep == 0) {
1306: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1308: /* Unsplit vertex */
1309: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1310: for (s = 0, q = 0; s < supportSize; ++s) {
1311: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1312: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1313: if (e >= 0) {
1314: supportNew[q++] = e + pMaxNew[dep+1];
1315: }
1316: }
1317: supportNew[q++] = hybedge;
1318: supportNew[q++] = hybedge;
1319: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1320: DMPlexSetSupport(sdm, newp, supportNew);
1321: /* Hybrid edge */
1322: coneNew[0] = newp;
1323: coneNew[1] = newp;
1324: DMPlexSetCone(sdm, hybedge, coneNew);
1325: for (e = 0, qf = 0; e < supportSize; ++e) {
1326: PetscInt val, edge;
1328: DMLabelGetValue(label, support[e], &val);
1329: if (val == 1) {
1330: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1331: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1332: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1333: } else if (val == (shift2 + 1)) {
1334: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1335: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1336: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1337: }
1338: }
1339: DMPlexSetSupport(sdm, hybedge, supportNew);
1340: } else if (dep == dim-2) {
1341: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1343: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1344: for (f = 0, qf = 0; f < supportSize; ++f) {
1345: PetscInt val, face;
1347: DMLabelGetValue(label, support[f], &val);
1348: if (val == dim-1) {
1349: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1350: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1351: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1352: supportNew[qf++] = face + pMaxNew[dep+1];
1353: } else {
1354: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1355: }
1356: }
1357: supportNew[qf++] = hybface;
1358: supportNew[qf++] = hybface;
1359: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1360: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1361: DMPlexSetSupport(sdm, newp, supportNew);
1362: /* Add hybrid face */
1363: coneNew[0] = newp;
1364: coneNew[1] = newp;
1365: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1366: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1367: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1368: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1369: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1370: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1371: DMPlexSetCone(sdm, hybface, coneNew);
1372: for (f = 0, qf = 0; f < supportSize; ++f) {
1373: PetscInt val, face;
1375: DMLabelGetValue(label, support[f], &val);
1376: if (val == dim-1) {
1377: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1378: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1379: }
1380: }
1381: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1382: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1383: DMPlexSetSupport(sdm, hybface, supportNew);
1384: }
1385: }
1386: }
1387: /* Step 6b: Replace split points in negative side cones */
1388: for (sp = 0; sp < numSP; ++sp) {
1389: PetscInt dep = values[sp];
1390: IS pIS;
1391: PetscInt numPoints;
1392: const PetscInt *points;
1394: if (dep >= 0) continue;
1395: DMLabelGetStratumIS(label, dep, &pIS);
1396: if (!pIS) continue;
1397: dep = -dep - shift;
1398: ISGetLocalSize(pIS, &numPoints);
1399: ISGetIndices(pIS, &points);
1400: for (p = 0; p < numPoints; ++p) {
1401: const PetscInt oldp = points[p];
1402: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1403: const PetscInt *cone;
1404: PetscInt coneSize, c;
1405: /* PetscBool replaced = PETSC_FALSE; */
1407: /* Negative edge: replace split vertex */
1408: /* Negative cell: replace split face */
1409: DMPlexGetConeSize(sdm, newp, &coneSize);
1410: DMPlexGetCone(sdm, newp, &cone);
1411: for (c = 0; c < coneSize; ++c) {
1412: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1413: PetscInt csplitp, cp, val;
1415: DMLabelGetValue(label, coldp, &val);
1416: if (val == dep-1) {
1417: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1418: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1419: csplitp = pMaxNew[dep-1] + cp;
1420: DMPlexInsertCone(sdm, newp, c, csplitp);
1421: /* replaced = PETSC_TRUE; */
1422: }
1423: }
1424: /* Cells with only a vertex or edge on the submesh have no replacement */
1425: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1426: }
1427: ISRestoreIndices(pIS, &points);
1428: ISDestroy(&pIS);
1429: }
1430: /* Step 7: Coordinates */
1431: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1432: DMGetCoordinateSection(sdm, &coordSection);
1433: DMGetCoordinatesLocal(sdm, &coordinates);
1434: VecGetArray(coordinates, &coords);
1435: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1436: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1437: const PetscInt splitp = pMaxNew[0] + v;
1438: PetscInt dof, off, soff, d;
1440: PetscSectionGetDof(coordSection, newp, &dof);
1441: PetscSectionGetOffset(coordSection, newp, &off);
1442: PetscSectionGetOffset(coordSection, splitp, &soff);
1443: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1444: }
1445: VecRestoreArray(coordinates, &coords);
1446: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1447: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1448: /* Step 9: Labels */
1449: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1450: DMGetNumLabels(sdm, &numLabels);
1451: for (dep = 0; dep <= depth; ++dep) {
1452: for (p = 0; p < numSplitPoints[dep]; ++p) {
1453: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1454: const PetscInt splitp = pMaxNew[dep] + p;
1455: PetscInt l;
1457: if (splitLabel) {
1458: const PetscInt val = 100 + dep;
1460: DMLabelSetValue(splitLabel, newp, val);
1461: DMLabelSetValue(splitLabel, splitp, -val);
1462: }
1463: for (l = 0; l < numLabels; ++l) {
1464: DMLabel mlabel;
1465: const char *lname;
1466: PetscInt val;
1467: PetscBool isDepth;
1469: DMGetLabelName(sdm, l, &lname);
1470: PetscStrcmp(lname, "depth", &isDepth);
1471: if (isDepth) continue;
1472: DMGetLabel(sdm, lname, &mlabel);
1473: DMLabelGetValue(mlabel, newp, &val);
1474: if (val >= 0) {
1475: DMLabelSetValue(mlabel, splitp, val);
1476: }
1477: }
1478: }
1479: }
1480: for (sp = 0; sp < numSP; ++sp) {
1481: const PetscInt dep = values[sp];
1483: if ((dep < 0) || (dep > depth)) continue;
1484: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1485: ISDestroy(&splitIS[dep]);
1486: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1487: ISDestroy(&unsplitIS[dep]);
1488: }
1489: if (label) {
1490: ISRestoreIndices(valueIS, &values);
1491: ISDestroy(&valueIS);
1492: }
1493: for (d = 0; d <= depth; ++d) {
1494: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1495: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1496: }
1497: 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);
1498: PetscFree3(coneNew, coneONew, supportNew);
1499: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1500: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1501: return(0);
1502: }
1504: /*@C
1505: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1507: Collective on dm
1509: Input Parameters:
1510: + dm - The original DM
1511: - label - The label specifying the boundary faces (this could be auto-generated)
1513: Output Parameters:
1514: + splitLabel - The label containing the split points, or NULL if no output is desired
1515: - dmSplit - The new DM
1517: Level: developer
1519: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1520: @*/
1521: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1522: {
1523: DM sdm;
1524: PetscInt dim;
1530: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1531: DMSetType(sdm, DMPLEX);
1532: DMGetDimension(dm, &dim);
1533: DMSetDimension(sdm, dim);
1534: switch (dim) {
1535: case 2:
1536: case 3:
1537: DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1538: break;
1539: default:
1540: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1541: }
1542: *dmSplit = sdm;
1543: return(0);
1544: }
1546: /* Returns the side of the surface for a given cell with a face on the surface */
1547: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1548: {
1549: const PetscInt *cone, *ornt;
1550: PetscInt dim, coneSize, c;
1551: PetscErrorCode ierr;
1554: *pos = PETSC_TRUE;
1555: DMGetDimension(dm, &dim);
1556: DMPlexGetConeSize(dm, cell, &coneSize);
1557: DMPlexGetCone(dm, cell, &cone);
1558: DMPlexGetConeOrientation(dm, cell, &ornt);
1559: for (c = 0; c < coneSize; ++c) {
1560: if (cone[c] == face) {
1561: PetscInt o = ornt[c];
1563: if (subdm) {
1564: const PetscInt *subcone, *subornt;
1565: PetscInt subpoint, subface, subconeSize, sc;
1567: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1568: PetscFindInt(face, numSubpoints, subpoints, &subface);
1569: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1570: DMPlexGetCone(subdm, subpoint, &subcone);
1571: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1572: for (sc = 0; sc < subconeSize; ++sc) {
1573: if (subcone[sc] == subface) {
1574: o = subornt[0];
1575: break;
1576: }
1577: }
1578: 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);
1579: }
1580: if (o >= 0) *pos = PETSC_TRUE;
1581: else *pos = PETSC_FALSE;
1582: break;
1583: }
1584: }
1585: 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);
1586: return(0);
1587: }
1589: /*@
1590: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1591: to complete the surface
1593: Input Parameters:
1594: + dm - The DM
1595: . label - A DMLabel marking the surface
1596: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1597: . flip - Flag to flip the submesh normal and replace points on the other side
1598: - subdm - The subDM associated with the label, or NULL
1600: Output Parameter:
1601: . label - A DMLabel marking all surface points
1603: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1605: Level: developer
1607: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1608: @*/
1609: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1610: {
1611: DMLabel depthLabel;
1612: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1613: const PetscInt *points, *subpoints;
1614: const PetscInt rev = flip ? -1 : 1;
1615: PetscInt *pMax;
1616: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1617: PetscErrorCode ierr;
1620: DMPlexGetDepth(dm, &depth);
1621: DMGetDimension(dm, &dim);
1622: pSize = PetscMax(depth, dim) + 1;
1623: PetscMalloc1(pSize,&pMax);
1624: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1625: DMPlexGetDepthLabel(dm, &depthLabel);
1626: DMGetDimension(dm, &dim);
1627: if (subdm) {
1628: DMPlexCreateSubpointIS(subdm, &subpointIS);
1629: if (subpointIS) {
1630: ISGetLocalSize(subpointIS, &numSubpoints);
1631: ISGetIndices(subpointIS, &subpoints);
1632: }
1633: }
1634: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1635: DMLabelGetStratumIS(label, dim-1, &dimIS);
1636: if (!dimIS) {
1637: PetscFree(pMax);
1638: ISDestroy(&subpointIS);
1639: return(0);
1640: }
1641: ISGetLocalSize(dimIS, &numPoints);
1642: ISGetIndices(dimIS, &points);
1643: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1644: const PetscInt *support;
1645: PetscInt supportSize, s;
1647: DMPlexGetSupportSize(dm, points[p], &supportSize);
1648: #if 0
1649: if (supportSize != 2) {
1650: const PetscInt *lp;
1651: PetscInt Nlp, pind;
1653: /* Check that for a cell with a single support face, that face is in the SF */
1654: /* THis check only works for the remote side. We would need root side information */
1655: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
1656: PetscFindInt(points[p], Nlp, lp, &pind);
1657: 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);
1658: }
1659: #endif
1660: DMPlexGetSupport(dm, points[p], &support);
1661: for (s = 0; s < supportSize; ++s) {
1662: const PetscInt *cone;
1663: PetscInt coneSize, c;
1664: PetscBool pos;
1666: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1667: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1668: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1669: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1670: /* Put faces touching the fault in the label */
1671: DMPlexGetConeSize(dm, support[s], &coneSize);
1672: DMPlexGetCone(dm, support[s], &cone);
1673: for (c = 0; c < coneSize; ++c) {
1674: const PetscInt point = cone[c];
1676: DMLabelGetValue(label, point, &val);
1677: if (val == -1) {
1678: PetscInt *closure = NULL;
1679: PetscInt closureSize, cl;
1681: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1682: for (cl = 0; cl < closureSize*2; cl += 2) {
1683: const PetscInt clp = closure[cl];
1684: PetscInt bval = -1;
1686: DMLabelGetValue(label, clp, &val);
1687: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1688: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1689: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1690: break;
1691: }
1692: }
1693: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1694: }
1695: }
1696: }
1697: }
1698: ISRestoreIndices(dimIS, &points);
1699: ISDestroy(&dimIS);
1700: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1701: ISDestroy(&subpointIS);
1702: /* Mark boundary points as unsplit */
1703: if (blabel) {
1704: DMLabelGetStratumIS(blabel, 1, &dimIS);
1705: ISGetLocalSize(dimIS, &numPoints);
1706: ISGetIndices(dimIS, &points);
1707: for (p = 0; p < numPoints; ++p) {
1708: const PetscInt point = points[p];
1709: PetscInt val, bval;
1711: DMLabelGetValue(blabel, point, &bval);
1712: if (bval >= 0) {
1713: DMLabelGetValue(label, point, &val);
1714: if ((val < 0) || (val > dim)) {
1715: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1716: DMLabelClearValue(blabel, point, bval);
1717: }
1718: }
1719: }
1720: for (p = 0; p < numPoints; ++p) {
1721: const PetscInt point = points[p];
1722: PetscInt val, bval;
1724: DMLabelGetValue(blabel, point, &bval);
1725: if (bval >= 0) {
1726: const PetscInt *cone, *support;
1727: PetscInt coneSize, supportSize, s, valA, valB, valE;
1729: /* Mark as unsplit */
1730: DMLabelGetValue(label, point, &val);
1731: 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);
1732: DMLabelClearValue(label, point, val);
1733: DMLabelSetValue(label, point, shift2+val);
1734: /* Check for cross-edge
1735: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1736: if (val != 0) continue;
1737: DMPlexGetSupport(dm, point, &support);
1738: DMPlexGetSupportSize(dm, point, &supportSize);
1739: for (s = 0; s < supportSize; ++s) {
1740: DMPlexGetCone(dm, support[s], &cone);
1741: DMPlexGetConeSize(dm, support[s], &coneSize);
1742: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1743: DMLabelGetValue(blabel, cone[0], &valA);
1744: DMLabelGetValue(blabel, cone[1], &valB);
1745: DMLabelGetValue(blabel, support[s], &valE);
1746: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1747: }
1748: }
1749: }
1750: ISRestoreIndices(dimIS, &points);
1751: ISDestroy(&dimIS);
1752: }
1753: /* Search for other cells/faces/edges connected to the fault by a vertex */
1754: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1755: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1756: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1757: cMax = cMax < 0 ? cEnd : cMax;
1758: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1759: DMLabelGetStratumIS(label, 0, &dimIS);
1760: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1761: if (dimIS && crossEdgeIS) {
1762: IS vertIS = dimIS;
1764: ISExpand(vertIS, crossEdgeIS, &dimIS);
1765: ISDestroy(&crossEdgeIS);
1766: ISDestroy(&vertIS);
1767: }
1768: if (!dimIS) {
1769: PetscFree(pMax);
1770: return(0);
1771: }
1772: ISGetLocalSize(dimIS, &numPoints);
1773: ISGetIndices(dimIS, &points);
1774: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1775: PetscInt *star = NULL;
1776: PetscInt starSize, s;
1777: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1779: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1780: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1781: while (again) {
1782: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1783: again = 0;
1784: for (s = 0; s < starSize*2; s += 2) {
1785: const PetscInt point = star[s];
1786: const PetscInt *cone;
1787: PetscInt coneSize, c;
1789: if ((point < cStart) || (point >= cMax)) continue;
1790: DMLabelGetValue(label, point, &val);
1791: if (val != -1) continue;
1792: again = again == 1 ? 1 : 2;
1793: DMPlexGetConeSize(dm, point, &coneSize);
1794: DMPlexGetCone(dm, point, &cone);
1795: for (c = 0; c < coneSize; ++c) {
1796: DMLabelGetValue(label, cone[c], &val);
1797: if (val != -1) {
1798: const PetscInt *ccone;
1799: PetscInt cconeSize, cc, side;
1801: 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);
1802: if (val > 0) side = 1;
1803: else side = -1;
1804: DMLabelSetValue(label, point, side*(shift+dim));
1805: /* Mark cell faces which touch the fault */
1806: DMPlexGetConeSize(dm, point, &cconeSize);
1807: DMPlexGetCone(dm, point, &ccone);
1808: for (cc = 0; cc < cconeSize; ++cc) {
1809: PetscInt *closure = NULL;
1810: PetscInt closureSize, cl;
1812: DMLabelGetValue(label, ccone[cc], &val);
1813: if (val != -1) continue;
1814: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1815: for (cl = 0; cl < closureSize*2; cl += 2) {
1816: const PetscInt clp = closure[cl];
1818: DMLabelGetValue(label, clp, &val);
1819: if (val == -1) continue;
1820: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1821: break;
1822: }
1823: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1824: }
1825: again = 1;
1826: break;
1827: }
1828: }
1829: }
1830: }
1831: /* Classify the rest by cell membership */
1832: for (s = 0; s < starSize*2; s += 2) {
1833: const PetscInt point = star[s];
1835: DMLabelGetValue(label, point, &val);
1836: if (val == -1) {
1837: PetscInt *sstar = NULL;
1838: PetscInt sstarSize, ss;
1839: PetscBool marked = PETSC_FALSE;
1841: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1842: for (ss = 0; ss < sstarSize*2; ss += 2) {
1843: const PetscInt spoint = sstar[ss];
1845: if ((spoint < cStart) || (spoint >= cMax)) continue;
1846: DMLabelGetValue(label, spoint, &val);
1847: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1848: DMLabelGetValue(depthLabel, point, &dep);
1849: if (val > 0) {
1850: DMLabelSetValue(label, point, shift+dep);
1851: } else {
1852: DMLabelSetValue(label, point, -(shift+dep));
1853: }
1854: marked = PETSC_TRUE;
1855: break;
1856: }
1857: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1858: DMLabelGetValue(depthLabel, point, &dep);
1859: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1860: }
1861: }
1862: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1863: }
1864: ISRestoreIndices(dimIS, &points);
1865: ISDestroy(&dimIS);
1866: /* If any faces touching the fault divide cells on either side, split them */
1867: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1868: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1869: ISExpand(facePosIS, faceNegIS, &dimIS);
1870: ISDestroy(&facePosIS);
1871: ISDestroy(&faceNegIS);
1872: ISGetLocalSize(dimIS, &numPoints);
1873: ISGetIndices(dimIS, &points);
1874: for (p = 0; p < numPoints; ++p) {
1875: const PetscInt point = points[p];
1876: const PetscInt *support;
1877: PetscInt supportSize, valA, valB;
1879: DMPlexGetSupportSize(dm, point, &supportSize);
1880: if (supportSize != 2) continue;
1881: DMPlexGetSupport(dm, point, &support);
1882: DMLabelGetValue(label, support[0], &valA);
1883: DMLabelGetValue(label, support[1], &valB);
1884: if ((valA == -1) || (valB == -1)) continue;
1885: if (valA*valB > 0) continue;
1886: /* Split the face */
1887: DMLabelGetValue(label, point, &valA);
1888: DMLabelClearValue(label, point, valA);
1889: DMLabelSetValue(label, point, dim-1);
1890: /* Label its closure:
1891: unmarked: label as unsplit
1892: incident: relabel as split
1893: split: do nothing
1894: */
1895: {
1896: PetscInt *closure = NULL;
1897: PetscInt closureSize, cl;
1899: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1900: for (cl = 0; cl < closureSize*2; cl += 2) {
1901: DMLabelGetValue(label, closure[cl], &valA);
1902: if (valA == -1) { /* Mark as unsplit */
1903: DMLabelGetValue(depthLabel, closure[cl], &dep);
1904: DMLabelSetValue(label, closure[cl], shift2+dep);
1905: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1906: DMLabelGetValue(depthLabel, closure[cl], &dep);
1907: DMLabelClearValue(label, closure[cl], valA);
1908: DMLabelSetValue(label, closure[cl], dep);
1909: }
1910: }
1911: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1912: }
1913: }
1914: ISRestoreIndices(dimIS, &points);
1915: ISDestroy(&dimIS);
1916: PetscFree(pMax);
1917: return(0);
1918: }
1920: /* Check that no cell have all vertices on the fault */
1921: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
1922: {
1923: IS subpointIS;
1924: const PetscInt *dmpoints;
1925: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
1926: PetscErrorCode ierr;
1929: if (!label) return(0);
1930: DMLabelGetDefaultValue(label, &defaultValue);
1931: DMPlexCreateSubpointIS(subdm, &subpointIS);
1932: if (!subpointIS) return(0);
1933: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
1934: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1935: ISGetIndices(subpointIS, &dmpoints);
1936: for (c = cStart; c < cEnd; ++c) {
1937: PetscBool invalidCell = PETSC_TRUE;
1938: PetscInt *closure = NULL;
1939: PetscInt closureSize, cl;
1941: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
1942: for (cl = 0; cl < closureSize*2; cl += 2) {
1943: PetscInt value = 0;
1945: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
1946: DMLabelGetValue(label, closure[cl], &value);
1947: if (value == defaultValue) {invalidCell = PETSC_FALSE; break;}
1948: }
1949: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
1950: if (invalidCell) {
1951: ISRestoreIndices(subpointIS, &dmpoints);
1952: ISDestroy(&subpointIS);
1953: DMDestroy(&subdm);
1954: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]);
1955: }
1956: }
1957: ISRestoreIndices(subpointIS, &dmpoints);
1958: ISDestroy(&subpointIS);
1959: return(0);
1960: }
1963: /*@
1964: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1966: Collective on dm
1968: Input Parameters:
1969: + dm - The original DM
1970: . label - The label specifying the interface vertices
1971: - bdlabel - The optional label specifying the interface boundary vertices
1973: Output Parameters:
1974: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
1975: . splitLabel - The label containing the split points, or NULL if no output is desired
1976: . dmInterface - The new interface DM, or NULL
1977: - dmHybrid - The new DM with cohesive cells
1979: Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points
1980: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
1981: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
1982: 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
1983: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
1985: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
1986: mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from
1987: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
1989: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
1990: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
1992: Level: developer
1994: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate()
1995: @*/
1996: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
1997: {
1998: DM idm;
1999: DMLabel subpointMap, hlabel, slabel = NULL;
2000: PetscInt dim;
2010: DMGetDimension(dm, &dim);
2011: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2012: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2013: DMPlexOrient(idm);
2014: DMPlexGetSubpointMap(idm, &subpointMap);
2015: DMLabelDuplicate(subpointMap, &hlabel);
2016: DMLabelClearStratum(hlabel, dim);
2017: if (splitLabel) {
2018: const char *name;
2019: char sname[PETSC_MAX_PATH_LEN];
2021: PetscObjectGetName((PetscObject) hlabel, &name);
2022: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2023: PetscStrcat(sname, " split");
2024: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2025: }
2026: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2027: if (dmInterface) {*dmInterface = idm;}
2028: else {DMDestroy(&idm);}
2029: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2030: if (hybridLabel) *hybridLabel = hlabel;
2031: else {DMLabelDestroy(&hlabel);}
2032: if (splitLabel) *splitLabel = slabel;
2033: return(0);
2034: }
2036: /* Here we need the explicit assumption that:
2038: For any marked cell, the marked vertices constitute a single face
2039: */
2040: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2041: {
2042: IS subvertexIS = NULL;
2043: const PetscInt *subvertices;
2044: PetscInt *pStart, *pEnd, *pMax, pSize;
2045: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2046: PetscErrorCode ierr;
2049: *numFaces = 0;
2050: *nFV = 0;
2051: DMPlexGetDepth(dm, &depth);
2052: DMGetDimension(dm, &dim);
2053: pSize = PetscMax(depth, dim) + 1;
2054: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
2055: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2056: for (d = 0; d <= depth; ++d) {
2057: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2058: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2059: }
2060: /* Loop over initial vertices and mark all faces in the collective star() */
2061: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2062: if (subvertexIS) {
2063: ISGetSize(subvertexIS, &numSubVerticesInitial);
2064: ISGetIndices(subvertexIS, &subvertices);
2065: }
2066: for (v = 0; v < numSubVerticesInitial; ++v) {
2067: const PetscInt vertex = subvertices[v];
2068: PetscInt *star = NULL;
2069: PetscInt starSize, s, numCells = 0, c;
2071: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2072: for (s = 0; s < starSize*2; s += 2) {
2073: const PetscInt point = star[s];
2074: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2075: }
2076: for (c = 0; c < numCells; ++c) {
2077: const PetscInt cell = star[c];
2078: PetscInt *closure = NULL;
2079: PetscInt closureSize, cl;
2080: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2082: DMLabelGetValue(subpointMap, cell, &cellLoc);
2083: if (cellLoc == 2) continue;
2084: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2085: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2086: for (cl = 0; cl < closureSize*2; cl += 2) {
2087: const PetscInt point = closure[cl];
2088: PetscInt vertexLoc;
2090: if ((point >= pStart[0]) && (point < pEnd[0])) {
2091: ++numCorners;
2092: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2093: if (vertexLoc == value) closure[faceSize++] = point;
2094: }
2095: }
2096: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2097: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2098: if (faceSize == *nFV) {
2099: const PetscInt *cells = NULL;
2100: PetscInt numCells, nc;
2102: ++(*numFaces);
2103: for (cl = 0; cl < faceSize; ++cl) {
2104: DMLabelSetValue(subpointMap, closure[cl], 0);
2105: }
2106: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2107: for (nc = 0; nc < numCells; ++nc) {
2108: DMLabelSetValue(subpointMap, cells[nc], 2);
2109: }
2110: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2111: }
2112: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2113: }
2114: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2115: }
2116: if (subvertexIS) {
2117: ISRestoreIndices(subvertexIS, &subvertices);
2118: }
2119: ISDestroy(&subvertexIS);
2120: PetscFree3(pStart,pEnd,pMax);
2121: return(0);
2122: }
2124: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2125: {
2126: IS subvertexIS = NULL;
2127: const PetscInt *subvertices;
2128: PetscInt *pStart, *pEnd, *pMax;
2129: PetscInt dim, d, numSubVerticesInitial = 0, v;
2130: PetscErrorCode ierr;
2133: DMGetDimension(dm, &dim);
2134: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2135: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2136: for (d = 0; d <= dim; ++d) {
2137: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2138: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2139: }
2140: /* Loop over initial vertices and mark all faces in the collective star() */
2141: if (vertexLabel) {
2142: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2143: if (subvertexIS) {
2144: ISGetSize(subvertexIS, &numSubVerticesInitial);
2145: ISGetIndices(subvertexIS, &subvertices);
2146: }
2147: }
2148: for (v = 0; v < numSubVerticesInitial; ++v) {
2149: const PetscInt vertex = subvertices[v];
2150: PetscInt *star = NULL;
2151: PetscInt starSize, s, numFaces = 0, f;
2153: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2154: for (s = 0; s < starSize*2; s += 2) {
2155: const PetscInt point = star[s];
2156: PetscInt faceLoc;
2158: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2159: if (markedFaces) {
2160: DMLabelGetValue(vertexLabel, point, &faceLoc);
2161: if (faceLoc < 0) continue;
2162: }
2163: star[numFaces++] = point;
2164: }
2165: }
2166: for (f = 0; f < numFaces; ++f) {
2167: const PetscInt face = star[f];
2168: PetscInt *closure = NULL;
2169: PetscInt closureSize, c;
2170: PetscInt faceLoc;
2172: DMLabelGetValue(subpointMap, face, &faceLoc);
2173: if (faceLoc == dim-1) continue;
2174: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2175: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2176: for (c = 0; c < closureSize*2; c += 2) {
2177: const PetscInt point = closure[c];
2178: PetscInt vertexLoc;
2180: if ((point >= pStart[0]) && (point < pEnd[0])) {
2181: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2182: if (vertexLoc != value) break;
2183: }
2184: }
2185: if (c == closureSize*2) {
2186: const PetscInt *support;
2187: PetscInt supportSize, s;
2189: for (c = 0; c < closureSize*2; c += 2) {
2190: const PetscInt point = closure[c];
2192: for (d = 0; d < dim; ++d) {
2193: if ((point >= pStart[d]) && (point < pEnd[d])) {
2194: DMLabelSetValue(subpointMap, point, d);
2195: break;
2196: }
2197: }
2198: }
2199: DMPlexGetSupportSize(dm, face, &supportSize);
2200: DMPlexGetSupport(dm, face, &support);
2201: for (s = 0; s < supportSize; ++s) {
2202: DMLabelSetValue(subpointMap, support[s], dim);
2203: }
2204: }
2205: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2206: }
2207: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2208: }
2209: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2210: ISDestroy(&subvertexIS);
2211: PetscFree3(pStart,pEnd,pMax);
2212: return(0);
2213: }
2215: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2216: {
2217: DMLabel label = NULL;
2218: const PetscInt *cone;
2219: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2220: PetscErrorCode ierr;
2223: *numFaces = 0;
2224: *nFV = 0;
2225: if (labelname) {DMGetLabel(dm, labelname, &label);}
2226: *subCells = NULL;
2227: DMGetDimension(dm, &dim);
2228: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2229: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2230: if (cMax < 0) return(0);
2231: if (label) {
2232: for (c = cMax; c < cEnd; ++c) {
2233: PetscInt val;
2235: DMLabelGetValue(label, c, &val);
2236: if (val == value) {
2237: ++(*numFaces);
2238: DMPlexGetConeSize(dm, c, &coneSize);
2239: }
2240: }
2241: } else {
2242: *numFaces = cEnd - cMax;
2243: DMPlexGetConeSize(dm, cMax, &coneSize);
2244: }
2245: PetscMalloc1(*numFaces *2, subCells);
2246: if (!(*numFaces)) return(0);
2247: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2248: for (c = cMax; c < cEnd; ++c) {
2249: const PetscInt *cells;
2250: PetscInt numCells;
2252: if (label) {
2253: PetscInt val;
2255: DMLabelGetValue(label, c, &val);
2256: if (val != value) continue;
2257: }
2258: DMPlexGetCone(dm, c, &cone);
2259: for (p = 0; p < *nFV; ++p) {
2260: DMLabelSetValue(subpointMap, cone[p], 0);
2261: }
2262: /* Negative face */
2263: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2264: /* Not true in parallel
2265: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2266: for (p = 0; p < numCells; ++p) {
2267: DMLabelSetValue(subpointMap, cells[p], 2);
2268: (*subCells)[subc++] = cells[p];
2269: }
2270: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2271: /* Positive face is not included */
2272: }
2273: return(0);
2274: }
2276: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2277: {
2278: PetscInt *pStart, *pEnd;
2279: PetscInt dim, cMax, cEnd, c, d;
2283: DMGetDimension(dm, &dim);
2284: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2285: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2286: if (cMax < 0) return(0);
2287: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2288: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2289: for (c = cMax; c < cEnd; ++c) {
2290: const PetscInt *cone;
2291: PetscInt *closure = NULL;
2292: PetscInt fconeSize, coneSize, closureSize, cl, val;
2294: if (label) {
2295: DMLabelGetValue(label, c, &val);
2296: if (val != value) continue;
2297: }
2298: DMPlexGetConeSize(dm, c, &coneSize);
2299: DMPlexGetCone(dm, c, &cone);
2300: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2301: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2302: /* Negative face */
2303: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2304: for (cl = 0; cl < closureSize*2; cl += 2) {
2305: const PetscInt point = closure[cl];
2307: for (d = 0; d <= dim; ++d) {
2308: if ((point >= pStart[d]) && (point < pEnd[d])) {
2309: DMLabelSetValue(subpointMap, point, d);
2310: break;
2311: }
2312: }
2313: }
2314: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2315: /* Cells -- positive face is not included */
2316: for (cl = 0; cl < 1; ++cl) {
2317: const PetscInt *support;
2318: PetscInt supportSize, s;
2320: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2321: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2322: DMPlexGetSupport(dm, cone[cl], &support);
2323: for (s = 0; s < supportSize; ++s) {
2324: DMLabelSetValue(subpointMap, support[s], dim);
2325: }
2326: }
2327: }
2328: PetscFree2(pStart, pEnd);
2329: return(0);
2330: }
2332: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2333: {
2334: MPI_Comm comm;
2335: PetscBool posOrient = PETSC_FALSE;
2336: const PetscInt debug = 0;
2337: PetscInt cellDim, faceSize, f;
2341: PetscObjectGetComm((PetscObject)dm,&comm);
2342: DMGetDimension(dm, &cellDim);
2343: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2345: if (cellDim == 1 && numCorners == 2) {
2346: /* Triangle */
2347: faceSize = numCorners-1;
2348: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2349: } else if (cellDim == 2 && numCorners == 3) {
2350: /* Triangle */
2351: faceSize = numCorners-1;
2352: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2353: } else if (cellDim == 3 && numCorners == 4) {
2354: /* Tetrahedron */
2355: faceSize = numCorners-1;
2356: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2357: } else if (cellDim == 1 && numCorners == 3) {
2358: /* Quadratic line */
2359: faceSize = 1;
2360: posOrient = PETSC_TRUE;
2361: } else if (cellDim == 2 && numCorners == 4) {
2362: /* Quads */
2363: faceSize = 2;
2364: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2365: posOrient = PETSC_TRUE;
2366: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2367: posOrient = PETSC_TRUE;
2368: } else {
2369: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2370: posOrient = PETSC_FALSE;
2371: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2372: }
2373: } else if (cellDim == 2 && numCorners == 6) {
2374: /* Quadratic triangle (I hate this) */
2375: /* Edges are determined by the first 2 vertices (corners of edges) */
2376: const PetscInt faceSizeTri = 3;
2377: PetscInt sortedIndices[3], i, iFace;
2378: PetscBool found = PETSC_FALSE;
2379: PetscInt faceVerticesTriSorted[9] = {
2380: 0, 3, 4, /* bottom */
2381: 1, 4, 5, /* right */
2382: 2, 3, 5, /* left */
2383: };
2384: PetscInt faceVerticesTri[9] = {
2385: 0, 3, 4, /* bottom */
2386: 1, 4, 5, /* right */
2387: 2, 5, 3, /* left */
2388: };
2390: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2391: PetscSortInt(faceSizeTri, sortedIndices);
2392: for (iFace = 0; iFace < 3; ++iFace) {
2393: const PetscInt ii = iFace*faceSizeTri;
2394: PetscInt fVertex, cVertex;
2396: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2397: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2398: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2399: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2400: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2401: faceVertices[fVertex] = origVertices[cVertex];
2402: break;
2403: }
2404: }
2405: }
2406: found = PETSC_TRUE;
2407: break;
2408: }
2409: }
2410: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2411: if (posOriented) *posOriented = PETSC_TRUE;
2412: return(0);
2413: } else if (cellDim == 2 && numCorners == 9) {
2414: /* Quadratic quad (I hate this) */
2415: /* Edges are determined by the first 2 vertices (corners of edges) */
2416: const PetscInt faceSizeQuad = 3;
2417: PetscInt sortedIndices[3], i, iFace;
2418: PetscBool found = PETSC_FALSE;
2419: PetscInt faceVerticesQuadSorted[12] = {
2420: 0, 1, 4, /* bottom */
2421: 1, 2, 5, /* right */
2422: 2, 3, 6, /* top */
2423: 0, 3, 7, /* left */
2424: };
2425: PetscInt faceVerticesQuad[12] = {
2426: 0, 1, 4, /* bottom */
2427: 1, 2, 5, /* right */
2428: 2, 3, 6, /* top */
2429: 3, 0, 7, /* left */
2430: };
2432: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2433: PetscSortInt(faceSizeQuad, sortedIndices);
2434: for (iFace = 0; iFace < 4; ++iFace) {
2435: const PetscInt ii = iFace*faceSizeQuad;
2436: PetscInt fVertex, cVertex;
2438: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2439: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2440: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2441: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2442: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2443: faceVertices[fVertex] = origVertices[cVertex];
2444: break;
2445: }
2446: }
2447: }
2448: found = PETSC_TRUE;
2449: break;
2450: }
2451: }
2452: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2453: if (posOriented) *posOriented = PETSC_TRUE;
2454: return(0);
2455: } else if (cellDim == 3 && numCorners == 8) {
2456: /* Hexes
2457: A hex is two oriented quads with the normal of the first
2458: pointing up at the second.
2460: 7---6
2461: /| /|
2462: 4---5 |
2463: | 1-|-2
2464: |/ |/
2465: 0---3
2467: Faces are determined by the first 4 vertices (corners of faces) */
2468: const PetscInt faceSizeHex = 4;
2469: PetscInt sortedIndices[4], i, iFace;
2470: PetscBool found = PETSC_FALSE;
2471: PetscInt faceVerticesHexSorted[24] = {
2472: 0, 1, 2, 3, /* bottom */
2473: 4, 5, 6, 7, /* top */
2474: 0, 3, 4, 5, /* front */
2475: 2, 3, 5, 6, /* right */
2476: 1, 2, 6, 7, /* back */
2477: 0, 1, 4, 7, /* left */
2478: };
2479: PetscInt faceVerticesHex[24] = {
2480: 1, 2, 3, 0, /* bottom */
2481: 4, 5, 6, 7, /* top */
2482: 0, 3, 5, 4, /* front */
2483: 3, 2, 6, 5, /* right */
2484: 2, 1, 7, 6, /* back */
2485: 1, 0, 4, 7, /* left */
2486: };
2488: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2489: PetscSortInt(faceSizeHex, sortedIndices);
2490: for (iFace = 0; iFace < 6; ++iFace) {
2491: const PetscInt ii = iFace*faceSizeHex;
2492: PetscInt fVertex, cVertex;
2494: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2495: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2496: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2497: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2498: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2499: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2500: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2501: faceVertices[fVertex] = origVertices[cVertex];
2502: break;
2503: }
2504: }
2505: }
2506: found = PETSC_TRUE;
2507: break;
2508: }
2509: }
2510: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2511: if (posOriented) *posOriented = PETSC_TRUE;
2512: return(0);
2513: } else if (cellDim == 3 && numCorners == 10) {
2514: /* Quadratic tet */
2515: /* Faces are determined by the first 3 vertices (corners of faces) */
2516: const PetscInt faceSizeTet = 6;
2517: PetscInt sortedIndices[6], i, iFace;
2518: PetscBool found = PETSC_FALSE;
2519: PetscInt faceVerticesTetSorted[24] = {
2520: 0, 1, 2, 6, 7, 8, /* bottom */
2521: 0, 3, 4, 6, 7, 9, /* front */
2522: 1, 4, 5, 7, 8, 9, /* right */
2523: 2, 3, 5, 6, 8, 9, /* left */
2524: };
2525: PetscInt faceVerticesTet[24] = {
2526: 0, 1, 2, 6, 7, 8, /* bottom */
2527: 0, 4, 3, 6, 7, 9, /* front */
2528: 1, 5, 4, 7, 8, 9, /* right */
2529: 2, 3, 5, 8, 6, 9, /* left */
2530: };
2532: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2533: PetscSortInt(faceSizeTet, sortedIndices);
2534: for (iFace=0; iFace < 4; ++iFace) {
2535: const PetscInt ii = iFace*faceSizeTet;
2536: PetscInt fVertex, cVertex;
2538: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2539: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2540: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2541: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2542: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2543: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2544: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2545: faceVertices[fVertex] = origVertices[cVertex];
2546: break;
2547: }
2548: }
2549: }
2550: found = PETSC_TRUE;
2551: break;
2552: }
2553: }
2554: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2555: if (posOriented) *posOriented = PETSC_TRUE;
2556: return(0);
2557: } else if (cellDim == 3 && numCorners == 27) {
2558: /* Quadratic hexes (I hate this)
2559: A hex is two oriented quads with the normal of the first
2560: pointing up at the second.
2562: 7---6
2563: /| /|
2564: 4---5 |
2565: | 3-|-2
2566: |/ |/
2567: 0---1
2569: Faces are determined by the first 4 vertices (corners of faces) */
2570: const PetscInt faceSizeQuadHex = 9;
2571: PetscInt sortedIndices[9], i, iFace;
2572: PetscBool found = PETSC_FALSE;
2573: PetscInt faceVerticesQuadHexSorted[54] = {
2574: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2575: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2576: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2577: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2578: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2579: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2580: };
2581: PetscInt faceVerticesQuadHex[54] = {
2582: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2583: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2584: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2585: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2586: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2587: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2588: };
2590: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2591: PetscSortInt(faceSizeQuadHex, sortedIndices);
2592: for (iFace = 0; iFace < 6; ++iFace) {
2593: const PetscInt ii = iFace*faceSizeQuadHex;
2594: PetscInt fVertex, cVertex;
2596: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2597: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2598: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2599: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2600: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2601: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2602: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2603: faceVertices[fVertex] = origVertices[cVertex];
2604: break;
2605: }
2606: }
2607: }
2608: found = PETSC_TRUE;
2609: break;
2610: }
2611: }
2612: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2613: if (posOriented) *posOriented = PETSC_TRUE;
2614: return(0);
2615: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2616: if (!posOrient) {
2617: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2618: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2619: } else {
2620: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2621: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2622: }
2623: if (posOriented) *posOriented = posOrient;
2624: return(0);
2625: }
2627: /*@
2628: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2629: in faceVertices. The orientation is such that the face normal points out of the cell
2631: Not collective
2633: Input Parameters:
2634: + dm - The original mesh
2635: . cell - The cell mesh point
2636: . faceSize - The number of vertices on the face
2637: . face - The face vertices
2638: . numCorners - The number of vertices on the cell
2639: . indices - Local numbering of face vertices in cell cone
2640: - origVertices - Original face vertices
2642: Output Parameter:
2643: + faceVertices - The face vertices properly oriented
2644: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2646: Level: developer
2648: .seealso: DMPlexGetCone()
2649: @*/
2650: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2651: {
2652: const PetscInt *cone = NULL;
2653: PetscInt coneSize, v, f, v2;
2654: PetscInt oppositeVertex = -1;
2655: PetscErrorCode ierr;
2658: DMPlexGetConeSize(dm, cell, &coneSize);
2659: DMPlexGetCone(dm, cell, &cone);
2660: for (v = 0, v2 = 0; v < coneSize; ++v) {
2661: PetscBool found = PETSC_FALSE;
2663: for (f = 0; f < faceSize; ++f) {
2664: if (face[f] == cone[v]) {
2665: found = PETSC_TRUE; break;
2666: }
2667: }
2668: if (found) {
2669: indices[v2] = v;
2670: origVertices[v2] = cone[v];
2671: ++v2;
2672: } else {
2673: oppositeVertex = v;
2674: }
2675: }
2676: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2677: return(0);
2678: }
2680: /*
2681: DMPlexInsertFace_Internal - Puts a face into the mesh
2683: Not collective
2685: Input Parameters:
2686: + dm - The DMPlex
2687: . numFaceVertex - The number of vertices in the face
2688: . faceVertices - The vertices in the face for dm
2689: . subfaceVertices - The vertices in the face for subdm
2690: . numCorners - The number of vertices in the cell
2691: . cell - A cell in dm containing the face
2692: . subcell - A cell in subdm containing the face
2693: . firstFace - First face in the mesh
2694: - newFacePoint - Next face in the mesh
2696: Output Parameters:
2697: . newFacePoint - Contains next face point number on input, updated on output
2699: Level: developer
2700: */
2701: 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)
2702: {
2703: MPI_Comm comm;
2704: DM_Plex *submesh = (DM_Plex*) subdm->data;
2705: const PetscInt *faces;
2706: PetscInt numFaces, coneSize;
2707: PetscErrorCode ierr;
2710: PetscObjectGetComm((PetscObject)dm,&comm);
2711: DMPlexGetConeSize(subdm, subcell, &coneSize);
2712: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2713: #if 0
2714: /* Cannot use this because support() has not been constructed yet */
2715: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2716: #else
2717: {
2718: PetscInt f;
2720: numFaces = 0;
2721: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2722: for (f = firstFace; f < *newFacePoint; ++f) {
2723: PetscInt dof, off, d;
2725: PetscSectionGetDof(submesh->coneSection, f, &dof);
2726: PetscSectionGetOffset(submesh->coneSection, f, &off);
2727: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2728: for (d = 0; d < dof; ++d) {
2729: const PetscInt p = submesh->cones[off+d];
2730: PetscInt v;
2732: for (v = 0; v < numFaceVertices; ++v) {
2733: if (subfaceVertices[v] == p) break;
2734: }
2735: if (v == numFaceVertices) break;
2736: }
2737: if (d == dof) {
2738: numFaces = 1;
2739: ((PetscInt*) faces)[0] = f;
2740: }
2741: }
2742: }
2743: #endif
2744: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2745: else if (numFaces == 1) {
2746: /* Add the other cell neighbor for this face */
2747: DMPlexSetCone(subdm, subcell, faces);
2748: } else {
2749: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2750: PetscBool posOriented;
2752: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2753: origVertices = &orientedVertices[numFaceVertices];
2754: indices = &orientedVertices[numFaceVertices*2];
2755: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2756: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2757: /* TODO: I know that routine should return a permutation, not the indices */
2758: for (v = 0; v < numFaceVertices; ++v) {
2759: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2760: for (ov = 0; ov < numFaceVertices; ++ov) {
2761: if (orientedVertices[ov] == vertex) {
2762: orientedSubVertices[ov] = subvertex;
2763: break;
2764: }
2765: }
2766: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2767: }
2768: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2769: DMPlexSetCone(subdm, subcell, newFacePoint);
2770: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2771: ++(*newFacePoint);
2772: }
2773: #if 0
2774: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2775: #else
2776: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2777: #endif
2778: return(0);
2779: }
2781: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2782: {
2783: MPI_Comm comm;
2784: DMLabel subpointMap;
2785: IS subvertexIS, subcellIS;
2786: const PetscInt *subVertices, *subCells;
2787: PetscInt numSubVertices, firstSubVertex, numSubCells;
2788: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2789: PetscInt vStart, vEnd, c, f;
2790: PetscErrorCode ierr;
2793: PetscObjectGetComm((PetscObject)dm,&comm);
2794: /* Create subpointMap which marks the submesh */
2795: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2796: DMPlexSetSubpointMap(subdm, subpointMap);
2797: DMLabelDestroy(&subpointMap);
2798: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2799: /* Setup chart */
2800: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2801: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2802: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2803: DMPlexSetVTKCellHeight(subdm, 1);
2804: /* Set cone sizes */
2805: firstSubVertex = numSubCells;
2806: firstSubFace = numSubCells+numSubVertices;
2807: newFacePoint = firstSubFace;
2808: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2809: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2810: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2811: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2812: for (c = 0; c < numSubCells; ++c) {
2813: DMPlexSetConeSize(subdm, c, 1);
2814: }
2815: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2816: DMPlexSetConeSize(subdm, f, nFV);
2817: }
2818: DMSetUp(subdm);
2819: /* Create face cones */
2820: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2821: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2822: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2823: for (c = 0; c < numSubCells; ++c) {
2824: const PetscInt cell = subCells[c];
2825: const PetscInt subcell = c;
2826: PetscInt *closure = NULL;
2827: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2829: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2830: for (cl = 0; cl < closureSize*2; cl += 2) {
2831: const PetscInt point = closure[cl];
2832: PetscInt subVertex;
2834: if ((point >= vStart) && (point < vEnd)) {
2835: ++numCorners;
2836: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2837: if (subVertex >= 0) {
2838: closure[faceSize] = point;
2839: subface[faceSize] = firstSubVertex+subVertex;
2840: ++faceSize;
2841: }
2842: }
2843: }
2844: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2845: if (faceSize == nFV) {
2846: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2847: }
2848: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2849: }
2850: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2851: DMPlexSymmetrize(subdm);
2852: DMPlexStratify(subdm);
2853: /* Build coordinates */
2854: {
2855: PetscSection coordSection, subCoordSection;
2856: Vec coordinates, subCoordinates;
2857: PetscScalar *coords, *subCoords;
2858: PetscInt numComp, coordSize, v;
2859: const char *name;
2861: DMGetCoordinateSection(dm, &coordSection);
2862: DMGetCoordinatesLocal(dm, &coordinates);
2863: DMGetCoordinateSection(subdm, &subCoordSection);
2864: PetscSectionSetNumFields(subCoordSection, 1);
2865: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2866: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2867: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2868: for (v = 0; v < numSubVertices; ++v) {
2869: const PetscInt vertex = subVertices[v];
2870: const PetscInt subvertex = firstSubVertex+v;
2871: PetscInt dof;
2873: PetscSectionGetDof(coordSection, vertex, &dof);
2874: PetscSectionSetDof(subCoordSection, subvertex, dof);
2875: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2876: }
2877: PetscSectionSetUp(subCoordSection);
2878: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2879: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2880: PetscObjectGetName((PetscObject)coordinates,&name);
2881: PetscObjectSetName((PetscObject)subCoordinates,name);
2882: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2883: VecSetType(subCoordinates,VECSTANDARD);
2884: if (coordSize) {
2885: VecGetArray(coordinates, &coords);
2886: VecGetArray(subCoordinates, &subCoords);
2887: for (v = 0; v < numSubVertices; ++v) {
2888: const PetscInt vertex = subVertices[v];
2889: const PetscInt subvertex = firstSubVertex+v;
2890: PetscInt dof, off, sdof, soff, d;
2892: PetscSectionGetDof(coordSection, vertex, &dof);
2893: PetscSectionGetOffset(coordSection, vertex, &off);
2894: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2895: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2896: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2897: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2898: }
2899: VecRestoreArray(coordinates, &coords);
2900: VecRestoreArray(subCoordinates, &subCoords);
2901: }
2902: DMSetCoordinatesLocal(subdm, subCoordinates);
2903: VecDestroy(&subCoordinates);
2904: }
2905: /* Cleanup */
2906: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2907: ISDestroy(&subvertexIS);
2908: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2909: ISDestroy(&subcellIS);
2910: return(0);
2911: }
2913: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2914: {
2915: PetscInt subPoint;
2918: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2919: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2920: }
2922: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2923: {
2924: MPI_Comm comm;
2925: DMLabel subpointMap;
2926: IS *subpointIS;
2927: const PetscInt **subpoints;
2928: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2929: PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2930: PetscMPIInt rank;
2931: PetscErrorCode ierr;
2934: PetscObjectGetComm((PetscObject)dm,&comm);
2935: MPI_Comm_rank(comm, &rank);
2936: /* Create subpointMap which marks the submesh */
2937: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2938: DMPlexSetSubpointMap(subdm, subpointMap);
2939: if (cellHeight) {
2940: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2941: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
2942: } else {
2943: DMLabel depth;
2944: IS pointIS;
2945: const PetscInt *points;
2946: PetscInt numPoints;
2948: DMPlexGetDepthLabel(dm, &depth);
2949: DMLabelGetStratumSize(label, value, &numPoints);
2950: DMLabelGetStratumIS(label, value, &pointIS);
2951: ISGetIndices(pointIS, &points);
2952: for (p = 0; p < numPoints; ++p) {
2953: PetscInt *closure = NULL;
2954: PetscInt closureSize, c, pdim;
2956: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2957: for (c = 0; c < closureSize*2; c += 2) {
2958: DMLabelGetValue(depth, closure[c], &pdim);
2959: DMLabelSetValue(subpointMap, closure[c], pdim);
2960: }
2961: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2962: }
2963: ISRestoreIndices(pointIS, &points);
2964: ISDestroy(&pointIS);
2965: }
2966: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2967: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2968: cMax = (cMax < 0) ? cEnd : cMax;
2969: /* Setup chart */
2970: DMGetDimension(dm, &dim);
2971: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2972: for (d = 0; d <= dim; ++d) {
2973: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2974: totSubPoints += numSubPoints[d];
2975: }
2976: DMPlexSetChart(subdm, 0, totSubPoints);
2977: DMPlexSetVTKCellHeight(subdm, cellHeight);
2978: /* Set cone sizes */
2979: firstSubPoint[dim] = 0;
2980: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2981: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2982: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2983: for (d = 0; d <= dim; ++d) {
2984: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2985: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2986: }
2987: for (d = 0; d <= dim; ++d) {
2988: for (p = 0; p < numSubPoints[d]; ++p) {
2989: const PetscInt point = subpoints[d][p];
2990: const PetscInt subpoint = firstSubPoint[d] + p;
2991: const PetscInt *cone;
2992: PetscInt coneSize, coneSizeNew, c, val;
2994: DMPlexGetConeSize(dm, point, &coneSize);
2995: DMPlexSetConeSize(subdm, subpoint, coneSize);
2996: if (cellHeight && (d == dim)) {
2997: DMPlexGetCone(dm, point, &cone);
2998: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2999: DMLabelGetValue(subpointMap, cone[c], &val);
3000: if (val >= 0) coneSizeNew++;
3001: }
3002: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3003: }
3004: }
3005: }
3006: DMLabelDestroy(&subpointMap);
3007: DMSetUp(subdm);
3008: /* Set cones */
3009: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3010: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
3011: for (d = 0; d <= dim; ++d) {
3012: for (p = 0; p < numSubPoints[d]; ++p) {
3013: const PetscInt point = subpoints[d][p];
3014: const PetscInt subpoint = firstSubPoint[d] + p;
3015: const PetscInt *cone, *ornt;
3016: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3018: if (d == dim-1) {
3019: const PetscInt *support, *cone, *ornt;
3020: PetscInt supportSize, coneSize, s, subc;
3022: DMPlexGetSupport(dm, point, &support);
3023: DMPlexGetSupportSize(dm, point, &supportSize);
3024: for (s = 0; s < supportSize; ++s) {
3025: if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
3026: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3027: if (subc >= 0) {
3028: const PetscInt ccell = subpoints[d+1][subc];
3030: DMPlexGetCone(dm, ccell, &cone);
3031: DMPlexGetConeSize(dm, ccell, &coneSize);
3032: DMPlexGetConeOrientation(dm, ccell, &ornt);
3033: for (c = 0; c < coneSize; ++c) {
3034: if (cone[c] == point) {
3035: fornt = ornt[c];
3036: break;
3037: }
3038: }
3039: break;
3040: }
3041: }
3042: }
3043: DMPlexGetConeSize(dm, point, &coneSize);
3044: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3045: DMPlexGetCone(dm, point, &cone);
3046: DMPlexGetConeOrientation(dm, point, &ornt);
3047: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3048: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3049: if (subc >= 0) {
3050: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3051: orntNew[coneSizeNew] = ornt[c];
3052: ++coneSizeNew;
3053: }
3054: }
3055: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3056: if (fornt < 0) {
3057: /* This should be replaced by a call to DMPlexReverseCell() */
3058: #if 0
3059: DMPlexReverseCell(subdm, subpoint);
3060: #else
3061: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3062: PetscInt faceSize, tmp;
3064: tmp = coneNew[c];
3065: coneNew[c] = coneNew[coneSizeNew-1-c];
3066: coneNew[coneSizeNew-1-c] = tmp;
3067: DMPlexGetConeSize(dm, cone[c], &faceSize);
3068: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3069: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3070: orntNew[coneSizeNew-1-c] = tmp;
3071: }
3072: }
3073: DMPlexSetCone(subdm, subpoint, coneNew);
3074: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3075: #endif
3076: }
3077: }
3078: PetscFree2(coneNew,orntNew);
3079: DMPlexSymmetrize(subdm);
3080: DMPlexStratify(subdm);
3081: /* Build coordinates */
3082: {
3083: PetscSection coordSection, subCoordSection;
3084: Vec coordinates, subCoordinates;
3085: PetscScalar *coords, *subCoords;
3086: PetscInt cdim, numComp, coordSize;
3087: const char *name;
3089: DMGetCoordinateDim(dm, &cdim);
3090: DMGetCoordinateSection(dm, &coordSection);
3091: DMGetCoordinatesLocal(dm, &coordinates);
3092: DMGetCoordinateSection(subdm, &subCoordSection);
3093: PetscSectionSetNumFields(subCoordSection, 1);
3094: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3095: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3096: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3097: for (v = 0; v < numSubPoints[0]; ++v) {
3098: const PetscInt vertex = subpoints[0][v];
3099: const PetscInt subvertex = firstSubPoint[0]+v;
3100: PetscInt dof;
3102: PetscSectionGetDof(coordSection, vertex, &dof);
3103: PetscSectionSetDof(subCoordSection, subvertex, dof);
3104: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3105: }
3106: PetscSectionSetUp(subCoordSection);
3107: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3108: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3109: PetscObjectGetName((PetscObject)coordinates,&name);
3110: PetscObjectSetName((PetscObject)subCoordinates,name);
3111: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3112: VecSetBlockSize(subCoordinates, cdim);
3113: VecSetType(subCoordinates,VECSTANDARD);
3114: VecGetArray(coordinates, &coords);
3115: VecGetArray(subCoordinates, &subCoords);
3116: for (v = 0; v < numSubPoints[0]; ++v) {
3117: const PetscInt vertex = subpoints[0][v];
3118: const PetscInt subvertex = firstSubPoint[0]+v;
3119: PetscInt dof, off, sdof, soff, d;
3121: PetscSectionGetDof(coordSection, vertex, &dof);
3122: PetscSectionGetOffset(coordSection, vertex, &off);
3123: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3124: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3125: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3126: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3127: }
3128: VecRestoreArray(coordinates, &coords);
3129: VecRestoreArray(subCoordinates, &subCoords);
3130: DMSetCoordinatesLocal(subdm, subCoordinates);
3131: VecDestroy(&subCoordinates);
3132: }
3133: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3134: {
3135: PetscSF sfPoint, sfPointSub;
3136: IS subpIS;
3137: const PetscSFNode *remotePoints;
3138: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3139: const PetscInt *localPoints, *subpoints;
3140: PetscInt *slocalPoints;
3141: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3142: PetscMPIInt rank;
3144: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3145: DMGetPointSF(dm, &sfPoint);
3146: DMGetPointSF(subdm, &sfPointSub);
3147: DMPlexGetChart(dm, &pStart, &pEnd);
3148: DMPlexGetChart(subdm, NULL, &numSubroots);
3149: DMPlexCreateSubpointIS(subdm, &subpIS);
3150: if (subpIS) {
3151: ISGetIndices(subpIS, &subpoints);
3152: ISGetLocalSize(subpIS, &numSubpoints);
3153: }
3154: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3155: if (numRoots >= 0) {
3156: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3157: for (p = 0; p < pEnd-pStart; ++p) {
3158: newLocalPoints[p].rank = -2;
3159: newLocalPoints[p].index = -2;
3160: }
3161: /* Set subleaves */
3162: for (l = 0; l < numLeaves; ++l) {
3163: const PetscInt point = localPoints[l];
3164: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3166: if (subpoint < 0) continue;
3167: newLocalPoints[point-pStart].rank = rank;
3168: newLocalPoints[point-pStart].index = subpoint;
3169: ++numSubleaves;
3170: }
3171: /* Must put in owned subpoints */
3172: for (p = pStart; p < pEnd; ++p) {
3173: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3175: if (subpoint < 0) {
3176: newOwners[p-pStart].rank = -3;
3177: newOwners[p-pStart].index = -3;
3178: } else {
3179: newOwners[p-pStart].rank = rank;
3180: newOwners[p-pStart].index = subpoint;
3181: }
3182: }
3183: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3184: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3185: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3186: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3187: PetscMalloc1(numSubleaves, &slocalPoints);
3188: PetscMalloc1(numSubleaves, &sremotePoints);
3189: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3190: const PetscInt point = localPoints[l];
3191: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3193: if (subpoint < 0) continue;
3194: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3195: slocalPoints[sl] = subpoint;
3196: sremotePoints[sl].rank = newLocalPoints[point].rank;
3197: sremotePoints[sl].index = newLocalPoints[point].index;
3198: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3199: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3200: ++sl;
3201: }
3202: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3203: PetscFree2(newLocalPoints,newOwners);
3204: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3205: }
3206: if (subpIS) {
3207: ISRestoreIndices(subpIS, &subpoints);
3208: ISDestroy(&subpIS);
3209: }
3210: }
3211: /* Cleanup */
3212: for (d = 0; d <= dim; ++d) {
3213: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3214: ISDestroy(&subpointIS[d]);
3215: }
3216: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3217: return(0);
3218: }
3220: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3221: {
3225: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3226: return(0);
3227: }
3229: /*@
3230: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3232: Input Parameters:
3233: + dm - The original mesh
3234: . vertexLabel - The DMLabel marking points contained in the surface
3235: . value - The label value to use
3236: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3238: Output Parameter:
3239: . subdm - The surface mesh
3241: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3243: Level: developer
3245: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3246: @*/
3247: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3248: {
3249: PetscInt dim, cdim, depth;
3255: DMGetDimension(dm, &dim);
3256: DMPlexGetDepth(dm, &depth);
3257: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3258: DMSetType(*subdm, DMPLEX);
3259: DMSetDimension(*subdm, dim-1);
3260: DMGetCoordinateDim(dm, &cdim);
3261: DMSetCoordinateDim(*subdm, cdim);
3262: if (depth == dim) {
3263: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3264: } else {
3265: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3266: }
3267: return(0);
3268: }
3270: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3271: {
3272: MPI_Comm comm;
3273: DMLabel subpointMap;
3274: IS subvertexIS;
3275: const PetscInt *subVertices;
3276: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3277: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3278: PetscInt cMax, c, f;
3279: PetscErrorCode ierr;
3282: PetscObjectGetComm((PetscObject)dm, &comm);
3283: /* Create subpointMap which marks the submesh */
3284: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3285: DMPlexSetSubpointMap(subdm, subpointMap);
3286: DMLabelDestroy(&subpointMap);
3287: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3288: /* Setup chart */
3289: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3290: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3291: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3292: DMPlexSetVTKCellHeight(subdm, 1);
3293: /* Set cone sizes */
3294: firstSubVertex = numSubCells;
3295: firstSubFace = numSubCells+numSubVertices;
3296: newFacePoint = firstSubFace;
3297: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3298: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3299: for (c = 0; c < numSubCells; ++c) {
3300: DMPlexSetConeSize(subdm, c, 1);
3301: }
3302: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3303: DMPlexSetConeSize(subdm, f, nFV);
3304: }
3305: DMSetUp(subdm);
3306: /* Create face cones */
3307: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3308: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3309: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3310: for (c = 0; c < numSubCells; ++c) {
3311: const PetscInt cell = subCells[c];
3312: const PetscInt subcell = c;
3313: const PetscInt *cone, *cells;
3314: PetscInt numCells, subVertex, p, v;
3316: if (cell < cMax) continue;
3317: DMPlexGetCone(dm, cell, &cone);
3318: for (v = 0; v < nFV; ++v) {
3319: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3320: subface[v] = firstSubVertex+subVertex;
3321: }
3322: DMPlexSetCone(subdm, newFacePoint, subface);
3323: DMPlexSetCone(subdm, subcell, &newFacePoint);
3324: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3325: /* Not true in parallel
3326: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3327: for (p = 0; p < numCells; ++p) {
3328: PetscInt negsubcell;
3330: if (cells[p] >= cMax) continue;
3331: /* I know this is a crap search */
3332: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3333: if (subCells[negsubcell] == cells[p]) break;
3334: }
3335: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3336: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3337: }
3338: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3339: ++newFacePoint;
3340: }
3341: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3342: DMPlexSymmetrize(subdm);
3343: DMPlexStratify(subdm);
3344: /* Build coordinates */
3345: {
3346: PetscSection coordSection, subCoordSection;
3347: Vec coordinates, subCoordinates;
3348: PetscScalar *coords, *subCoords;
3349: PetscInt cdim, numComp, coordSize, v;
3350: const char *name;
3352: DMGetCoordinateDim(dm, &cdim);
3353: DMGetCoordinateSection(dm, &coordSection);
3354: DMGetCoordinatesLocal(dm, &coordinates);
3355: DMGetCoordinateSection(subdm, &subCoordSection);
3356: PetscSectionSetNumFields(subCoordSection, 1);
3357: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3358: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3359: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3360: for (v = 0; v < numSubVertices; ++v) {
3361: const PetscInt vertex = subVertices[v];
3362: const PetscInt subvertex = firstSubVertex+v;
3363: PetscInt dof;
3365: PetscSectionGetDof(coordSection, vertex, &dof);
3366: PetscSectionSetDof(subCoordSection, subvertex, dof);
3367: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3368: }
3369: PetscSectionSetUp(subCoordSection);
3370: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3371: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3372: PetscObjectGetName((PetscObject)coordinates,&name);
3373: PetscObjectSetName((PetscObject)subCoordinates,name);
3374: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3375: VecSetBlockSize(subCoordinates, cdim);
3376: VecSetType(subCoordinates,VECSTANDARD);
3377: VecGetArray(coordinates, &coords);
3378: VecGetArray(subCoordinates, &subCoords);
3379: for (v = 0; v < numSubVertices; ++v) {
3380: const PetscInt vertex = subVertices[v];
3381: const PetscInt subvertex = firstSubVertex+v;
3382: PetscInt dof, off, sdof, soff, d;
3384: PetscSectionGetDof(coordSection, vertex, &dof);
3385: PetscSectionGetOffset(coordSection, vertex, &off);
3386: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3387: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3388: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3389: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3390: }
3391: VecRestoreArray(coordinates, &coords);
3392: VecRestoreArray(subCoordinates, &subCoords);
3393: DMSetCoordinatesLocal(subdm, subCoordinates);
3394: VecDestroy(&subCoordinates);
3395: }
3396: /* Build SF */
3397: CHKMEMQ;
3398: {
3399: PetscSF sfPoint, sfPointSub;
3400: const PetscSFNode *remotePoints;
3401: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3402: const PetscInt *localPoints;
3403: PetscInt *slocalPoints;
3404: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3405: PetscMPIInt rank;
3407: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3408: DMGetPointSF(dm, &sfPoint);
3409: DMGetPointSF(subdm, &sfPointSub);
3410: DMPlexGetChart(dm, &pStart, &pEnd);
3411: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3412: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3413: if (numRoots >= 0) {
3414: /* Only vertices should be shared */
3415: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3416: for (p = 0; p < pEnd-pStart; ++p) {
3417: newLocalPoints[p].rank = -2;
3418: newLocalPoints[p].index = -2;
3419: }
3420: /* Set subleaves */
3421: for (l = 0; l < numLeaves; ++l) {
3422: const PetscInt point = localPoints[l];
3423: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3425: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3426: if (subPoint < 0) continue;
3427: newLocalPoints[point-pStart].rank = rank;
3428: newLocalPoints[point-pStart].index = subPoint;
3429: ++numSubLeaves;
3430: }
3431: /* Must put in owned subpoints */
3432: for (p = pStart; p < pEnd; ++p) {
3433: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3435: if (subPoint < 0) {
3436: newOwners[p-pStart].rank = -3;
3437: newOwners[p-pStart].index = -3;
3438: } else {
3439: newOwners[p-pStart].rank = rank;
3440: newOwners[p-pStart].index = subPoint;
3441: }
3442: }
3443: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3444: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3445: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3446: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3447: PetscMalloc1(numSubLeaves, &slocalPoints);
3448: PetscMalloc1(numSubLeaves, &sremotePoints);
3449: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3450: const PetscInt point = localPoints[l];
3451: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3453: if (subPoint < 0) continue;
3454: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3455: slocalPoints[sl] = subPoint;
3456: sremotePoints[sl].rank = newLocalPoints[point].rank;
3457: sremotePoints[sl].index = newLocalPoints[point].index;
3458: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3459: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3460: ++sl;
3461: }
3462: PetscFree2(newLocalPoints,newOwners);
3463: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3464: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3465: }
3466: }
3467: CHKMEMQ;
3468: /* Cleanup */
3469: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3470: ISDestroy(&subvertexIS);
3471: PetscFree(subCells);
3472: return(0);
3473: }
3475: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3476: {
3477: DMLabel label = NULL;
3481: if (labelname) {DMGetLabel(dm, labelname, &label);}
3482: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3483: return(0);
3484: }
3486: /*@C
3487: 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.
3489: Input Parameters:
3490: + dm - The original mesh
3491: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3492: . label - A label name, or NULL
3493: - value - A label value
3495: Output Parameter:
3496: . subdm - The surface mesh
3498: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3500: Level: developer
3502: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3503: @*/
3504: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3505: {
3506: PetscInt dim, cdim, depth;
3512: DMGetDimension(dm, &dim);
3513: DMPlexGetDepth(dm, &depth);
3514: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3515: DMSetType(*subdm, DMPLEX);
3516: DMSetDimension(*subdm, dim-1);
3517: DMGetCoordinateDim(dm, &cdim);
3518: DMSetCoordinateDim(*subdm, cdim);
3519: if (depth == dim) {
3520: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3521: } else {
3522: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3523: }
3524: return(0);
3525: }
3527: /*@
3528: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3530: Input Parameters:
3531: + dm - The original mesh
3532: . cellLabel - The DMLabel marking cells contained in the new mesh
3533: - value - The label value to use
3535: Output Parameter:
3536: . subdm - The new mesh
3538: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3540: Level: developer
3542: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3543: @*/
3544: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3545: {
3546: PetscInt dim;
3552: DMGetDimension(dm, &dim);
3553: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3554: DMSetType(*subdm, DMPLEX);
3555: DMSetDimension(*subdm, dim);
3556: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3557: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3558: return(0);
3559: }
3561: /*@
3562: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3564: Input Parameter:
3565: . dm - The submesh DM
3567: Output Parameter:
3568: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3570: Level: developer
3572: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3573: @*/
3574: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3575: {
3579: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3580: return(0);
3581: }
3583: /*@
3584: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3586: Input Parameters:
3587: + dm - The submesh DM
3588: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3590: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3592: Level: developer
3594: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3595: @*/
3596: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3597: {
3598: DM_Plex *mesh = (DM_Plex *) dm->data;
3599: DMLabel tmp;
3604: tmp = mesh->subpointMap;
3605: mesh->subpointMap = subpointMap;
3606: PetscObjectReference((PetscObject) mesh->subpointMap);
3607: DMLabelDestroy(&tmp);
3608: return(0);
3609: }
3611: /*@
3612: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3614: Input Parameter:
3615: . dm - The submesh DM
3617: Output Parameter:
3618: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3620: Note: This IS is guaranteed to be sorted by the construction of the submesh
3622: Level: developer
3624: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3625: @*/
3626: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3627: {
3628: MPI_Comm comm;
3629: DMLabel subpointMap;
3630: IS is;
3631: const PetscInt *opoints;
3632: PetscInt *points, *depths;
3633: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3634: PetscErrorCode ierr;
3639: PetscObjectGetComm((PetscObject)dm,&comm);
3640: *subpointIS = NULL;
3641: DMPlexGetSubpointMap(dm, &subpointMap);
3642: DMPlexGetDepth(dm, &depth);
3643: if (subpointMap && depth >= 0) {
3644: DMPlexGetChart(dm, &pStart, &pEnd);
3645: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3646: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3647: depths[0] = depth;
3648: depths[1] = 0;
3649: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3650: PetscMalloc1(pEnd, &points);
3651: for(d = 0, off = 0; d <= depth; ++d) {
3652: const PetscInt dep = depths[d];
3654: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3655: DMLabelGetStratumSize(subpointMap, dep, &n);
3656: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3657: 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);
3658: } else {
3659: if (!n) {
3660: if (d == 0) {
3661: /* Missing cells */
3662: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3663: } else {
3664: /* Missing faces */
3665: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3666: }
3667: }
3668: }
3669: if (n) {
3670: DMLabelGetStratumIS(subpointMap, dep, &is);
3671: ISGetIndices(is, &opoints);
3672: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3673: ISRestoreIndices(is, &opoints);
3674: ISDestroy(&is);
3675: }
3676: }
3677: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3678: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3679: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3680: }
3681: return(0);
3682: }
3684: /*@
3685: DMPlexGetSubpoint - Return the subpoint corresponding to a point in the original mesh. If the DM
3686: is not a submesh, just return the input point.
3688: Note collective
3690: Input Parameters:
3691: + dm - The submesh DM
3692: - p - The point in the original, from which the submesh was created
3694: Output Parameter:
3695: . subp - The point in the submesh
3697: Level: developer
3699: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3700: @*/
3701: PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp)
3702: {
3703: DMLabel spmap;
3707: *subp = p;
3708: DMPlexGetSubpointMap(dm, &spmap);
3709: if (spmap) {
3710: IS subpointIS;
3711: const PetscInt *subpoints;
3712: PetscInt numSubpoints;
3714: /* TODO Cache the IS, making it look like an index */
3715: DMPlexCreateSubpointIS(dm, &subpointIS);
3716: ISGetLocalSize(subpointIS, &numSubpoints);
3717: ISGetIndices(subpointIS, &subpoints);
3718: PetscFindInt(p, numSubpoints, subpoints, subp);
3719: if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p);
3720: ISRestoreIndices(subpointIS, &subpoints);
3721: ISDestroy(&subpointIS);
3722: }
3723: return(0);
3724: }
3726: /*@
3727: DMPlexGetAuxiliaryPoint - For a given point in the DM, return the matching point in the auxiliary DM.
3729: Note collective
3731: Input Parameters:
3732: + dm - The DM
3733: . dmAux - The related auxiliary DM
3734: - p - The point in the original DM
3736: Output Parameter:
3737: . subp - The point in the auxiliary DM
3739: 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,
3740: then we map the point back to the original space.
3742: Level: developer
3744: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3745: @*/
3746: PetscErrorCode DMPlexGetAuxiliaryPoint(DM dm, DM dmAux, PetscInt p, PetscInt *subp)
3747: {
3748: DMLabel spmap;
3752: *subp = p;
3753: /* If dm is a submesh, do not get subpoint */
3754: DMPlexGetSubpointMap(dm, &spmap);
3755: if (dmAux && !spmap) {
3756: PetscInt h;
3758: DMPlexGetVTKCellHeight(dmAux, &h);
3759: DMPlexGetSubpointMap(dmAux, &spmap);
3760: if (spmap && !h) {DMLabelGetValue(spmap, p, subp);}
3761: else {DMPlexGetSubpoint(dmAux, p, subp);}
3762: }
3763: return(0);
3764: }