Actual source code: plexsubmesh.c
petsc-3.10.5 2019-03-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 flag;
851: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
852: DMSetType(gdm, DMPLEX);
853: DMGetDimension(dm, &dim);
854: DMSetDimension(gdm, dim);
855: DMPlexGetAdjacencyUseCone(dm, &flag);
856: DMPlexSetAdjacencyUseCone(gdm, flag);
857: DMPlexGetAdjacencyUseClosure(dm, &flag);
858: DMPlexSetAdjacencyUseClosure(gdm, flag);
859: DMGetLabel(dm, name, &label);
860: if (!label) {
861: /* Get label for boundary faces */
862: DMCreateLabel(dm, name);
863: DMGetLabel(dm, name, &label);
864: DMPlexMarkBoundaryFaces(dm, 1, label);
865: }
866: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
867: DMCopyBoundary(dm, gdm);
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: Level: developer
1981: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1982: @*/
1983: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
1984: {
1985: DM idm;
1986: DMLabel subpointMap, hlabel, slabel = NULL;
1987: PetscInt dim;
1997: DMGetDimension(dm, &dim);
1998: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
1999: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2000: DMPlexOrient(idm);
2001: DMPlexGetSubpointMap(idm, &subpointMap);
2002: DMLabelDuplicate(subpointMap, &hlabel);
2003: DMLabelClearStratum(hlabel, dim);
2004: if (splitLabel) {
2005: const char *name;
2006: char sname[PETSC_MAX_PATH_LEN];
2008: DMLabelGetName(hlabel, &name);
2009: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2010: PetscStrcat(sname, " split");
2011: DMLabelCreate(sname, &slabel);
2012: }
2013: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);
2014: if (dmInterface) {*dmInterface = idm;}
2015: else {DMDestroy(&idm);}
2016: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2017: if (hybridLabel) *hybridLabel = hlabel;
2018: else {DMLabelDestroy(&hlabel);}
2019: if (splitLabel) *splitLabel = slabel;
2020: return(0);
2021: }
2023: /* Here we need the explicit assumption that:
2025: For any marked cell, the marked vertices constitute a single face
2026: */
2027: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2028: {
2029: IS subvertexIS = NULL;
2030: const PetscInt *subvertices;
2031: PetscInt *pStart, *pEnd, *pMax, pSize;
2032: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2033: PetscErrorCode ierr;
2036: *numFaces = 0;
2037: *nFV = 0;
2038: DMPlexGetDepth(dm, &depth);
2039: DMGetDimension(dm, &dim);
2040: pSize = PetscMax(depth, dim) + 1;
2041: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
2042: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2043: for (d = 0; d <= depth; ++d) {
2044: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2045: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2046: }
2047: /* Loop over initial vertices and mark all faces in the collective star() */
2048: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
2049: if (subvertexIS) {
2050: ISGetSize(subvertexIS, &numSubVerticesInitial);
2051: ISGetIndices(subvertexIS, &subvertices);
2052: }
2053: for (v = 0; v < numSubVerticesInitial; ++v) {
2054: const PetscInt vertex = subvertices[v];
2055: PetscInt *star = NULL;
2056: PetscInt starSize, s, numCells = 0, c;
2058: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2059: for (s = 0; s < starSize*2; s += 2) {
2060: const PetscInt point = star[s];
2061: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2062: }
2063: for (c = 0; c < numCells; ++c) {
2064: const PetscInt cell = star[c];
2065: PetscInt *closure = NULL;
2066: PetscInt closureSize, cl;
2067: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2069: DMLabelGetValue(subpointMap, cell, &cellLoc);
2070: if (cellLoc == 2) continue;
2071: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
2072: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2073: for (cl = 0; cl < closureSize*2; cl += 2) {
2074: const PetscInt point = closure[cl];
2075: PetscInt vertexLoc;
2077: if ((point >= pStart[0]) && (point < pEnd[0])) {
2078: ++numCorners;
2079: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2080: if (vertexLoc == value) closure[faceSize++] = point;
2081: }
2082: }
2083: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
2084: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2085: if (faceSize == *nFV) {
2086: const PetscInt *cells = NULL;
2087: PetscInt numCells, nc;
2089: ++(*numFaces);
2090: for (cl = 0; cl < faceSize; ++cl) {
2091: DMLabelSetValue(subpointMap, closure[cl], 0);
2092: }
2093: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2094: for (nc = 0; nc < numCells; ++nc) {
2095: DMLabelSetValue(subpointMap, cells[nc], 2);
2096: }
2097: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2098: }
2099: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2100: }
2101: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2102: }
2103: if (subvertexIS) {
2104: ISRestoreIndices(subvertexIS, &subvertices);
2105: }
2106: ISDestroy(&subvertexIS);
2107: PetscFree3(pStart,pEnd,pMax);
2108: return(0);
2109: }
2111: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2112: {
2113: IS subvertexIS = NULL;
2114: const PetscInt *subvertices;
2115: PetscInt *pStart, *pEnd, *pMax;
2116: PetscInt dim, d, numSubVerticesInitial = 0, v;
2117: PetscErrorCode ierr;
2120: DMGetDimension(dm, &dim);
2121: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2122: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2123: for (d = 0; d <= dim; ++d) {
2124: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2125: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2126: }
2127: /* Loop over initial vertices and mark all faces in the collective star() */
2128: if (vertexLabel) {
2129: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2130: if (subvertexIS) {
2131: ISGetSize(subvertexIS, &numSubVerticesInitial);
2132: ISGetIndices(subvertexIS, &subvertices);
2133: }
2134: }
2135: for (v = 0; v < numSubVerticesInitial; ++v) {
2136: const PetscInt vertex = subvertices[v];
2137: PetscInt *star = NULL;
2138: PetscInt starSize, s, numFaces = 0, f;
2140: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2141: for (s = 0; s < starSize*2; s += 2) {
2142: const PetscInt point = star[s];
2143: PetscInt faceLoc;
2145: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) {
2146: if (markedFaces) {
2147: DMLabelGetValue(vertexLabel, point, &faceLoc);
2148: if (faceLoc < 0) continue;
2149: }
2150: star[numFaces++] = point;
2151: }
2152: }
2153: for (f = 0; f < numFaces; ++f) {
2154: const PetscInt face = star[f];
2155: PetscInt *closure = NULL;
2156: PetscInt closureSize, c;
2157: PetscInt faceLoc;
2159: DMLabelGetValue(subpointMap, face, &faceLoc);
2160: if (faceLoc == dim-1) continue;
2161: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2162: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2163: for (c = 0; c < closureSize*2; c += 2) {
2164: const PetscInt point = closure[c];
2165: PetscInt vertexLoc;
2167: if ((point >= pStart[0]) && (point < pEnd[0])) {
2168: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2169: if (vertexLoc != value) break;
2170: }
2171: }
2172: if (c == closureSize*2) {
2173: const PetscInt *support;
2174: PetscInt supportSize, s;
2176: for (c = 0; c < closureSize*2; c += 2) {
2177: const PetscInt point = closure[c];
2179: for (d = 0; d < dim; ++d) {
2180: if ((point >= pStart[d]) && (point < pEnd[d])) {
2181: DMLabelSetValue(subpointMap, point, d);
2182: break;
2183: }
2184: }
2185: }
2186: DMPlexGetSupportSize(dm, face, &supportSize);
2187: DMPlexGetSupport(dm, face, &support);
2188: for (s = 0; s < supportSize; ++s) {
2189: DMLabelSetValue(subpointMap, support[s], dim);
2190: }
2191: }
2192: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2193: }
2194: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2195: }
2196: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2197: ISDestroy(&subvertexIS);
2198: PetscFree3(pStart,pEnd,pMax);
2199: return(0);
2200: }
2202: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2203: {
2204: DMLabel label = NULL;
2205: const PetscInt *cone;
2206: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2207: PetscErrorCode ierr;
2210: *numFaces = 0;
2211: *nFV = 0;
2212: if (labelname) {DMGetLabel(dm, labelname, &label);}
2213: *subCells = NULL;
2214: DMGetDimension(dm, &dim);
2215: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2216: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2217: if (cMax < 0) return(0);
2218: if (label) {
2219: for (c = cMax; c < cEnd; ++c) {
2220: PetscInt val;
2222: DMLabelGetValue(label, c, &val);
2223: if (val == value) {
2224: ++(*numFaces);
2225: DMPlexGetConeSize(dm, c, &coneSize);
2226: }
2227: }
2228: } else {
2229: *numFaces = cEnd - cMax;
2230: DMPlexGetConeSize(dm, cMax, &coneSize);
2231: }
2232: PetscMalloc1(*numFaces *2, subCells);
2233: if (!(*numFaces)) return(0);
2234: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2235: for (c = cMax; c < cEnd; ++c) {
2236: const PetscInt *cells;
2237: PetscInt numCells;
2239: if (label) {
2240: PetscInt val;
2242: DMLabelGetValue(label, c, &val);
2243: if (val != value) continue;
2244: }
2245: DMPlexGetCone(dm, c, &cone);
2246: for (p = 0; p < *nFV; ++p) {
2247: DMLabelSetValue(subpointMap, cone[p], 0);
2248: }
2249: /* Negative face */
2250: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2251: /* Not true in parallel
2252: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2253: for (p = 0; p < numCells; ++p) {
2254: DMLabelSetValue(subpointMap, cells[p], 2);
2255: (*subCells)[subc++] = cells[p];
2256: }
2257: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2258: /* Positive face is not included */
2259: }
2260: return(0);
2261: }
2263: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2264: {
2265: PetscInt *pStart, *pEnd;
2266: PetscInt dim, cMax, cEnd, c, d;
2270: DMGetDimension(dm, &dim);
2271: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2272: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2273: if (cMax < 0) return(0);
2274: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2275: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2276: for (c = cMax; c < cEnd; ++c) {
2277: const PetscInt *cone;
2278: PetscInt *closure = NULL;
2279: PetscInt fconeSize, coneSize, closureSize, cl, val;
2281: if (label) {
2282: DMLabelGetValue(label, c, &val);
2283: if (val != value) continue;
2284: }
2285: DMPlexGetConeSize(dm, c, &coneSize);
2286: DMPlexGetCone(dm, c, &cone);
2287: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2288: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2289: /* Negative face */
2290: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2291: for (cl = 0; cl < closureSize*2; cl += 2) {
2292: const PetscInt point = closure[cl];
2294: for (d = 0; d <= dim; ++d) {
2295: if ((point >= pStart[d]) && (point < pEnd[d])) {
2296: DMLabelSetValue(subpointMap, point, d);
2297: break;
2298: }
2299: }
2300: }
2301: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2302: /* Cells -- positive face is not included */
2303: for (cl = 0; cl < 1; ++cl) {
2304: const PetscInt *support;
2305: PetscInt supportSize, s;
2307: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2308: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2309: DMPlexGetSupport(dm, cone[cl], &support);
2310: for (s = 0; s < supportSize; ++s) {
2311: DMLabelSetValue(subpointMap, support[s], dim);
2312: }
2313: }
2314: }
2315: PetscFree2(pStart, pEnd);
2316: return(0);
2317: }
2319: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2320: {
2321: MPI_Comm comm;
2322: PetscBool posOrient = PETSC_FALSE;
2323: const PetscInt debug = 0;
2324: PetscInt cellDim, faceSize, f;
2328: PetscObjectGetComm((PetscObject)dm,&comm);
2329: DMGetDimension(dm, &cellDim);
2330: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2332: if (cellDim == 1 && numCorners == 2) {
2333: /* Triangle */
2334: faceSize = numCorners-1;
2335: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2336: } else if (cellDim == 2 && numCorners == 3) {
2337: /* Triangle */
2338: faceSize = numCorners-1;
2339: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2340: } else if (cellDim == 3 && numCorners == 4) {
2341: /* Tetrahedron */
2342: faceSize = numCorners-1;
2343: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2344: } else if (cellDim == 1 && numCorners == 3) {
2345: /* Quadratic line */
2346: faceSize = 1;
2347: posOrient = PETSC_TRUE;
2348: } else if (cellDim == 2 && numCorners == 4) {
2349: /* Quads */
2350: faceSize = 2;
2351: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2352: posOrient = PETSC_TRUE;
2353: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2354: posOrient = PETSC_TRUE;
2355: } else {
2356: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2357: posOrient = PETSC_FALSE;
2358: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2359: }
2360: } else if (cellDim == 2 && numCorners == 6) {
2361: /* Quadratic triangle (I hate this) */
2362: /* Edges are determined by the first 2 vertices (corners of edges) */
2363: const PetscInt faceSizeTri = 3;
2364: PetscInt sortedIndices[3], i, iFace;
2365: PetscBool found = PETSC_FALSE;
2366: PetscInt faceVerticesTriSorted[9] = {
2367: 0, 3, 4, /* bottom */
2368: 1, 4, 5, /* right */
2369: 2, 3, 5, /* left */
2370: };
2371: PetscInt faceVerticesTri[9] = {
2372: 0, 3, 4, /* bottom */
2373: 1, 4, 5, /* right */
2374: 2, 5, 3, /* left */
2375: };
2377: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2378: PetscSortInt(faceSizeTri, sortedIndices);
2379: for (iFace = 0; iFace < 3; ++iFace) {
2380: const PetscInt ii = iFace*faceSizeTri;
2381: PetscInt fVertex, cVertex;
2383: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2384: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2385: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2386: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2387: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2388: faceVertices[fVertex] = origVertices[cVertex];
2389: break;
2390: }
2391: }
2392: }
2393: found = PETSC_TRUE;
2394: break;
2395: }
2396: }
2397: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2398: if (posOriented) *posOriented = PETSC_TRUE;
2399: return(0);
2400: } else if (cellDim == 2 && numCorners == 9) {
2401: /* Quadratic quad (I hate this) */
2402: /* Edges are determined by the first 2 vertices (corners of edges) */
2403: const PetscInt faceSizeQuad = 3;
2404: PetscInt sortedIndices[3], i, iFace;
2405: PetscBool found = PETSC_FALSE;
2406: PetscInt faceVerticesQuadSorted[12] = {
2407: 0, 1, 4, /* bottom */
2408: 1, 2, 5, /* right */
2409: 2, 3, 6, /* top */
2410: 0, 3, 7, /* left */
2411: };
2412: PetscInt faceVerticesQuad[12] = {
2413: 0, 1, 4, /* bottom */
2414: 1, 2, 5, /* right */
2415: 2, 3, 6, /* top */
2416: 3, 0, 7, /* left */
2417: };
2419: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2420: PetscSortInt(faceSizeQuad, sortedIndices);
2421: for (iFace = 0; iFace < 4; ++iFace) {
2422: const PetscInt ii = iFace*faceSizeQuad;
2423: PetscInt fVertex, cVertex;
2425: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2426: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2427: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2428: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2429: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2430: faceVertices[fVertex] = origVertices[cVertex];
2431: break;
2432: }
2433: }
2434: }
2435: found = PETSC_TRUE;
2436: break;
2437: }
2438: }
2439: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2440: if (posOriented) *posOriented = PETSC_TRUE;
2441: return(0);
2442: } else if (cellDim == 3 && numCorners == 8) {
2443: /* Hexes
2444: A hex is two oriented quads with the normal of the first
2445: pointing up at the second.
2447: 7---6
2448: /| /|
2449: 4---5 |
2450: | 1-|-2
2451: |/ |/
2452: 0---3
2454: Faces are determined by the first 4 vertices (corners of faces) */
2455: const PetscInt faceSizeHex = 4;
2456: PetscInt sortedIndices[4], i, iFace;
2457: PetscBool found = PETSC_FALSE;
2458: PetscInt faceVerticesHexSorted[24] = {
2459: 0, 1, 2, 3, /* bottom */
2460: 4, 5, 6, 7, /* top */
2461: 0, 3, 4, 5, /* front */
2462: 2, 3, 5, 6, /* right */
2463: 1, 2, 6, 7, /* back */
2464: 0, 1, 4, 7, /* left */
2465: };
2466: PetscInt faceVerticesHex[24] = {
2467: 1, 2, 3, 0, /* bottom */
2468: 4, 5, 6, 7, /* top */
2469: 0, 3, 5, 4, /* front */
2470: 3, 2, 6, 5, /* right */
2471: 2, 1, 7, 6, /* back */
2472: 1, 0, 4, 7, /* left */
2473: };
2475: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2476: PetscSortInt(faceSizeHex, sortedIndices);
2477: for (iFace = 0; iFace < 6; ++iFace) {
2478: const PetscInt ii = iFace*faceSizeHex;
2479: PetscInt fVertex, cVertex;
2481: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2482: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2483: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2484: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2485: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2486: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2487: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2488: faceVertices[fVertex] = origVertices[cVertex];
2489: break;
2490: }
2491: }
2492: }
2493: found = PETSC_TRUE;
2494: break;
2495: }
2496: }
2497: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2498: if (posOriented) *posOriented = PETSC_TRUE;
2499: return(0);
2500: } else if (cellDim == 3 && numCorners == 10) {
2501: /* Quadratic tet */
2502: /* Faces are determined by the first 3 vertices (corners of faces) */
2503: const PetscInt faceSizeTet = 6;
2504: PetscInt sortedIndices[6], i, iFace;
2505: PetscBool found = PETSC_FALSE;
2506: PetscInt faceVerticesTetSorted[24] = {
2507: 0, 1, 2, 6, 7, 8, /* bottom */
2508: 0, 3, 4, 6, 7, 9, /* front */
2509: 1, 4, 5, 7, 8, 9, /* right */
2510: 2, 3, 5, 6, 8, 9, /* left */
2511: };
2512: PetscInt faceVerticesTet[24] = {
2513: 0, 1, 2, 6, 7, 8, /* bottom */
2514: 0, 4, 3, 6, 7, 9, /* front */
2515: 1, 5, 4, 7, 8, 9, /* right */
2516: 2, 3, 5, 8, 6, 9, /* left */
2517: };
2519: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2520: PetscSortInt(faceSizeTet, sortedIndices);
2521: for (iFace=0; iFace < 4; ++iFace) {
2522: const PetscInt ii = iFace*faceSizeTet;
2523: PetscInt fVertex, cVertex;
2525: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2526: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2527: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2528: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2529: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2530: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2531: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2532: faceVertices[fVertex] = origVertices[cVertex];
2533: break;
2534: }
2535: }
2536: }
2537: found = PETSC_TRUE;
2538: break;
2539: }
2540: }
2541: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2542: if (posOriented) *posOriented = PETSC_TRUE;
2543: return(0);
2544: } else if (cellDim == 3 && numCorners == 27) {
2545: /* Quadratic hexes (I hate this)
2546: A hex is two oriented quads with the normal of the first
2547: pointing up at the second.
2549: 7---6
2550: /| /|
2551: 4---5 |
2552: | 3-|-2
2553: |/ |/
2554: 0---1
2556: Faces are determined by the first 4 vertices (corners of faces) */
2557: const PetscInt faceSizeQuadHex = 9;
2558: PetscInt sortedIndices[9], i, iFace;
2559: PetscBool found = PETSC_FALSE;
2560: PetscInt faceVerticesQuadHexSorted[54] = {
2561: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2562: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2563: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2564: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2565: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2566: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2567: };
2568: PetscInt faceVerticesQuadHex[54] = {
2569: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2570: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2571: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2572: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2573: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2574: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2575: };
2577: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2578: PetscSortInt(faceSizeQuadHex, sortedIndices);
2579: for (iFace = 0; iFace < 6; ++iFace) {
2580: const PetscInt ii = iFace*faceSizeQuadHex;
2581: PetscInt fVertex, cVertex;
2583: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2584: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2585: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2586: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2587: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2588: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2589: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2590: faceVertices[fVertex] = origVertices[cVertex];
2591: break;
2592: }
2593: }
2594: }
2595: found = PETSC_TRUE;
2596: break;
2597: }
2598: }
2599: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2600: if (posOriented) *posOriented = PETSC_TRUE;
2601: return(0);
2602: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2603: if (!posOrient) {
2604: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2605: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2606: } else {
2607: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2608: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2609: }
2610: if (posOriented) *posOriented = posOrient;
2611: return(0);
2612: }
2614: /*@
2615: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2616: in faceVertices. The orientation is such that the face normal points out of the cell
2618: Not collective
2620: Input Parameters:
2621: + dm - The original mesh
2622: . cell - The cell mesh point
2623: . faceSize - The number of vertices on the face
2624: . face - The face vertices
2625: . numCorners - The number of vertices on the cell
2626: . indices - Local numbering of face vertices in cell cone
2627: - origVertices - Original face vertices
2629: Output Parameter:
2630: + faceVertices - The face vertices properly oriented
2631: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2633: Level: developer
2635: .seealso: DMPlexGetCone()
2636: @*/
2637: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2638: {
2639: const PetscInt *cone = NULL;
2640: PetscInt coneSize, v, f, v2;
2641: PetscInt oppositeVertex = -1;
2642: PetscErrorCode ierr;
2645: DMPlexGetConeSize(dm, cell, &coneSize);
2646: DMPlexGetCone(dm, cell, &cone);
2647: for (v = 0, v2 = 0; v < coneSize; ++v) {
2648: PetscBool found = PETSC_FALSE;
2650: for (f = 0; f < faceSize; ++f) {
2651: if (face[f] == cone[v]) {
2652: found = PETSC_TRUE; break;
2653: }
2654: }
2655: if (found) {
2656: indices[v2] = v;
2657: origVertices[v2] = cone[v];
2658: ++v2;
2659: } else {
2660: oppositeVertex = v;
2661: }
2662: }
2663: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2664: return(0);
2665: }
2667: /*
2668: DMPlexInsertFace_Internal - Puts a face into the mesh
2670: Not collective
2672: Input Parameters:
2673: + dm - The DMPlex
2674: . numFaceVertex - The number of vertices in the face
2675: . faceVertices - The vertices in the face for dm
2676: . subfaceVertices - The vertices in the face for subdm
2677: . numCorners - The number of vertices in the cell
2678: . cell - A cell in dm containing the face
2679: . subcell - A cell in subdm containing the face
2680: . firstFace - First face in the mesh
2681: - newFacePoint - Next face in the mesh
2683: Output Parameters:
2684: . newFacePoint - Contains next face point number on input, updated on output
2686: Level: developer
2687: */
2688: 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)
2689: {
2690: MPI_Comm comm;
2691: DM_Plex *submesh = (DM_Plex*) subdm->data;
2692: const PetscInt *faces;
2693: PetscInt numFaces, coneSize;
2694: PetscErrorCode ierr;
2697: PetscObjectGetComm((PetscObject)dm,&comm);
2698: DMPlexGetConeSize(subdm, subcell, &coneSize);
2699: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2700: #if 0
2701: /* Cannot use this because support() has not been constructed yet */
2702: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2703: #else
2704: {
2705: PetscInt f;
2707: numFaces = 0;
2708: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2709: for (f = firstFace; f < *newFacePoint; ++f) {
2710: PetscInt dof, off, d;
2712: PetscSectionGetDof(submesh->coneSection, f, &dof);
2713: PetscSectionGetOffset(submesh->coneSection, f, &off);
2714: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2715: for (d = 0; d < dof; ++d) {
2716: const PetscInt p = submesh->cones[off+d];
2717: PetscInt v;
2719: for (v = 0; v < numFaceVertices; ++v) {
2720: if (subfaceVertices[v] == p) break;
2721: }
2722: if (v == numFaceVertices) break;
2723: }
2724: if (d == dof) {
2725: numFaces = 1;
2726: ((PetscInt*) faces)[0] = f;
2727: }
2728: }
2729: }
2730: #endif
2731: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2732: else if (numFaces == 1) {
2733: /* Add the other cell neighbor for this face */
2734: DMPlexSetCone(subdm, subcell, faces);
2735: } else {
2736: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2737: PetscBool posOriented;
2739: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2740: origVertices = &orientedVertices[numFaceVertices];
2741: indices = &orientedVertices[numFaceVertices*2];
2742: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2743: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2744: /* TODO: I know that routine should return a permutation, not the indices */
2745: for (v = 0; v < numFaceVertices; ++v) {
2746: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2747: for (ov = 0; ov < numFaceVertices; ++ov) {
2748: if (orientedVertices[ov] == vertex) {
2749: orientedSubVertices[ov] = subvertex;
2750: break;
2751: }
2752: }
2753: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2754: }
2755: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2756: DMPlexSetCone(subdm, subcell, newFacePoint);
2757: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2758: ++(*newFacePoint);
2759: }
2760: #if 0
2761: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2762: #else
2763: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2764: #endif
2765: return(0);
2766: }
2768: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2769: {
2770: MPI_Comm comm;
2771: DMLabel subpointMap;
2772: IS subvertexIS, subcellIS;
2773: const PetscInt *subVertices, *subCells;
2774: PetscInt numSubVertices, firstSubVertex, numSubCells;
2775: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2776: PetscInt vStart, vEnd, c, f;
2777: PetscErrorCode ierr;
2780: PetscObjectGetComm((PetscObject)dm,&comm);
2781: /* Create subpointMap which marks the submesh */
2782: DMLabelCreate("subpoint_map", &subpointMap);
2783: DMPlexSetSubpointMap(subdm, subpointMap);
2784: DMLabelDestroy(&subpointMap);
2785: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2786: /* Setup chart */
2787: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2788: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2789: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2790: DMPlexSetVTKCellHeight(subdm, 1);
2791: /* Set cone sizes */
2792: firstSubVertex = numSubCells;
2793: firstSubFace = numSubCells+numSubVertices;
2794: newFacePoint = firstSubFace;
2795: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2796: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2797: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2798: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2799: for (c = 0; c < numSubCells; ++c) {
2800: DMPlexSetConeSize(subdm, c, 1);
2801: }
2802: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2803: DMPlexSetConeSize(subdm, f, nFV);
2804: }
2805: DMSetUp(subdm);
2806: /* Create face cones */
2807: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2808: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2809: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2810: for (c = 0; c < numSubCells; ++c) {
2811: const PetscInt cell = subCells[c];
2812: const PetscInt subcell = c;
2813: PetscInt *closure = NULL;
2814: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2816: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2817: for (cl = 0; cl < closureSize*2; cl += 2) {
2818: const PetscInt point = closure[cl];
2819: PetscInt subVertex;
2821: if ((point >= vStart) && (point < vEnd)) {
2822: ++numCorners;
2823: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2824: if (subVertex >= 0) {
2825: closure[faceSize] = point;
2826: subface[faceSize] = firstSubVertex+subVertex;
2827: ++faceSize;
2828: }
2829: }
2830: }
2831: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2832: if (faceSize == nFV) {
2833: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2834: }
2835: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2836: }
2837: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2838: DMPlexSymmetrize(subdm);
2839: DMPlexStratify(subdm);
2840: /* Build coordinates */
2841: {
2842: PetscSection coordSection, subCoordSection;
2843: Vec coordinates, subCoordinates;
2844: PetscScalar *coords, *subCoords;
2845: PetscInt numComp, coordSize, v;
2846: const char *name;
2848: DMGetCoordinateSection(dm, &coordSection);
2849: DMGetCoordinatesLocal(dm, &coordinates);
2850: DMGetCoordinateSection(subdm, &subCoordSection);
2851: PetscSectionSetNumFields(subCoordSection, 1);
2852: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2853: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2854: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2855: for (v = 0; v < numSubVertices; ++v) {
2856: const PetscInt vertex = subVertices[v];
2857: const PetscInt subvertex = firstSubVertex+v;
2858: PetscInt dof;
2860: PetscSectionGetDof(coordSection, vertex, &dof);
2861: PetscSectionSetDof(subCoordSection, subvertex, dof);
2862: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2863: }
2864: PetscSectionSetUp(subCoordSection);
2865: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2866: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2867: PetscObjectGetName((PetscObject)coordinates,&name);
2868: PetscObjectSetName((PetscObject)subCoordinates,name);
2869: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2870: VecSetType(subCoordinates,VECSTANDARD);
2871: if (coordSize) {
2872: VecGetArray(coordinates, &coords);
2873: VecGetArray(subCoordinates, &subCoords);
2874: for (v = 0; v < numSubVertices; ++v) {
2875: const PetscInt vertex = subVertices[v];
2876: const PetscInt subvertex = firstSubVertex+v;
2877: PetscInt dof, off, sdof, soff, d;
2879: PetscSectionGetDof(coordSection, vertex, &dof);
2880: PetscSectionGetOffset(coordSection, vertex, &off);
2881: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2882: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2883: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2884: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2885: }
2886: VecRestoreArray(coordinates, &coords);
2887: VecRestoreArray(subCoordinates, &subCoords);
2888: }
2889: DMSetCoordinatesLocal(subdm, subCoordinates);
2890: VecDestroy(&subCoordinates);
2891: }
2892: /* Cleanup */
2893: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2894: ISDestroy(&subvertexIS);
2895: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2896: ISDestroy(&subcellIS);
2897: return(0);
2898: }
2900: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2901: {
2902: PetscInt subPoint;
2905: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2906: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2907: }
2909: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2910: {
2911: MPI_Comm comm;
2912: DMLabel subpointMap;
2913: IS *subpointIS;
2914: const PetscInt **subpoints;
2915: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2916: PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2917: PetscMPIInt rank;
2918: PetscErrorCode ierr;
2921: PetscObjectGetComm((PetscObject)dm,&comm);
2922: MPI_Comm_rank(comm, &rank);
2923: /* Create subpointMap which marks the submesh */
2924: DMLabelCreate("subpoint_map", &subpointMap);
2925: DMPlexSetSubpointMap(subdm, subpointMap);
2926: if (cellHeight) {
2927: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2928: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);}
2929: } else {
2930: DMLabel depth;
2931: IS pointIS;
2932: const PetscInt *points;
2933: PetscInt numPoints;
2935: DMPlexGetDepthLabel(dm, &depth);
2936: DMLabelGetStratumSize(label, value, &numPoints);
2937: DMLabelGetStratumIS(label, value, &pointIS);
2938: ISGetIndices(pointIS, &points);
2939: for (p = 0; p < numPoints; ++p) {
2940: PetscInt *closure = NULL;
2941: PetscInt closureSize, c, pdim;
2943: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2944: for (c = 0; c < closureSize*2; c += 2) {
2945: DMLabelGetValue(depth, closure[c], &pdim);
2946: DMLabelSetValue(subpointMap, closure[c], pdim);
2947: }
2948: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2949: }
2950: ISRestoreIndices(pointIS, &points);
2951: ISDestroy(&pointIS);
2952: }
2953: DMLabelDestroy(&subpointMap);
2954: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2955: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2956: cMax = (cMax < 0) ? cEnd : cMax;
2957: /* Setup chart */
2958: DMGetDimension(dm, &dim);
2959: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2960: for (d = 0; d <= dim; ++d) {
2961: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2962: totSubPoints += numSubPoints[d];
2963: }
2964: DMPlexSetChart(subdm, 0, totSubPoints);
2965: DMPlexSetVTKCellHeight(subdm, cellHeight);
2966: /* Set cone sizes */
2967: firstSubPoint[dim] = 0;
2968: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2969: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2970: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2971: for (d = 0; d <= dim; ++d) {
2972: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2973: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2974: }
2975: for (d = 0; d <= dim; ++d) {
2976: for (p = 0; p < numSubPoints[d]; ++p) {
2977: const PetscInt point = subpoints[d][p];
2978: const PetscInt subpoint = firstSubPoint[d] + p;
2979: const PetscInt *cone;
2980: PetscInt coneSize, coneSizeNew, c, val;
2982: DMPlexGetConeSize(dm, point, &coneSize);
2983: DMPlexSetConeSize(subdm, subpoint, coneSize);
2984: if (cellHeight && (d == dim)) {
2985: DMPlexGetCone(dm, point, &cone);
2986: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2987: DMLabelGetValue(subpointMap, cone[c], &val);
2988: if (val >= 0) coneSizeNew++;
2989: }
2990: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2991: }
2992: }
2993: }
2994: DMSetUp(subdm);
2995: /* Set cones */
2996: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2997: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2998: for (d = 0; d <= dim; ++d) {
2999: for (p = 0; p < numSubPoints[d]; ++p) {
3000: const PetscInt point = subpoints[d][p];
3001: const PetscInt subpoint = firstSubPoint[d] + p;
3002: const PetscInt *cone, *ornt;
3003: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3005: if (d == dim-1) {
3006: const PetscInt *support, *cone, *ornt;
3007: PetscInt supportSize, coneSize, s, subc;
3009: DMPlexGetSupport(dm, point, &support);
3010: DMPlexGetSupportSize(dm, point, &supportSize);
3011: for (s = 0; s < supportSize; ++s) {
3012: if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
3013: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
3014: if (subc >= 0) {
3015: const PetscInt ccell = subpoints[d+1][subc];
3017: DMPlexGetCone(dm, ccell, &cone);
3018: DMPlexGetConeSize(dm, ccell, &coneSize);
3019: DMPlexGetConeOrientation(dm, ccell, &ornt);
3020: for (c = 0; c < coneSize; ++c) {
3021: if (cone[c] == point) {
3022: fornt = ornt[c];
3023: break;
3024: }
3025: }
3026: break;
3027: }
3028: }
3029: }
3030: DMPlexGetConeSize(dm, point, &coneSize);
3031: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3032: DMPlexGetCone(dm, point, &cone);
3033: DMPlexGetConeOrientation(dm, point, &ornt);
3034: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3035: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
3036: if (subc >= 0) {
3037: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
3038: orntNew[coneSizeNew] = ornt[c];
3039: ++coneSizeNew;
3040: }
3041: }
3042: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
3043: if (fornt < 0) {
3044: /* This should be replaced by a call to DMPlexReverseCell() */
3045: #if 0
3046: DMPlexReverseCell(subdm, subpoint);
3047: #else
3048: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
3049: PetscInt faceSize, tmp;
3051: tmp = coneNew[c];
3052: coneNew[c] = coneNew[coneSizeNew-1-c];
3053: coneNew[coneSizeNew-1-c] = tmp;
3054: DMPlexGetConeSize(dm, cone[c], &faceSize);
3055: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
3056: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
3057: orntNew[coneSizeNew-1-c] = tmp;
3058: }
3059: }
3060: DMPlexSetCone(subdm, subpoint, coneNew);
3061: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3062: #endif
3063: }
3064: }
3065: PetscFree2(coneNew,orntNew);
3066: DMPlexSymmetrize(subdm);
3067: DMPlexStratify(subdm);
3068: /* Build coordinates */
3069: {
3070: PetscSection coordSection, subCoordSection;
3071: Vec coordinates, subCoordinates;
3072: PetscScalar *coords, *subCoords;
3073: PetscInt cdim, numComp, coordSize;
3074: const char *name;
3076: DMGetCoordinateDim(dm, &cdim);
3077: DMGetCoordinateSection(dm, &coordSection);
3078: DMGetCoordinatesLocal(dm, &coordinates);
3079: DMGetCoordinateSection(subdm, &subCoordSection);
3080: PetscSectionSetNumFields(subCoordSection, 1);
3081: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3082: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3083: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
3084: for (v = 0; v < numSubPoints[0]; ++v) {
3085: const PetscInt vertex = subpoints[0][v];
3086: const PetscInt subvertex = firstSubPoint[0]+v;
3087: PetscInt dof;
3089: PetscSectionGetDof(coordSection, vertex, &dof);
3090: PetscSectionSetDof(subCoordSection, subvertex, dof);
3091: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3092: }
3093: PetscSectionSetUp(subCoordSection);
3094: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3095: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3096: PetscObjectGetName((PetscObject)coordinates,&name);
3097: PetscObjectSetName((PetscObject)subCoordinates,name);
3098: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3099: VecSetBlockSize(subCoordinates, cdim);
3100: VecSetType(subCoordinates,VECSTANDARD);
3101: VecGetArray(coordinates, &coords);
3102: VecGetArray(subCoordinates, &subCoords);
3103: for (v = 0; v < numSubPoints[0]; ++v) {
3104: const PetscInt vertex = subpoints[0][v];
3105: const PetscInt subvertex = firstSubPoint[0]+v;
3106: PetscInt dof, off, sdof, soff, d;
3108: PetscSectionGetDof(coordSection, vertex, &dof);
3109: PetscSectionGetOffset(coordSection, vertex, &off);
3110: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3111: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3112: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3113: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3114: }
3115: VecRestoreArray(coordinates, &coords);
3116: VecRestoreArray(subCoordinates, &subCoords);
3117: DMSetCoordinatesLocal(subdm, subCoordinates);
3118: VecDestroy(&subCoordinates);
3119: }
3120: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3121: {
3122: PetscSF sfPoint, sfPointSub;
3123: IS subpIS;
3124: const PetscSFNode *remotePoints;
3125: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3126: const PetscInt *localPoints, *subpoints;
3127: PetscInt *slocalPoints;
3128: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3129: PetscMPIInt rank;
3131: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3132: DMGetPointSF(dm, &sfPoint);
3133: DMGetPointSF(subdm, &sfPointSub);
3134: DMPlexGetChart(dm, &pStart, &pEnd);
3135: DMPlexGetChart(subdm, NULL, &numSubroots);
3136: DMPlexCreateSubpointIS(subdm, &subpIS);
3137: if (subpIS) {
3138: ISGetIndices(subpIS, &subpoints);
3139: ISGetLocalSize(subpIS, &numSubpoints);
3140: }
3141: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3142: if (numRoots >= 0) {
3143: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3144: for (p = 0; p < pEnd-pStart; ++p) {
3145: newLocalPoints[p].rank = -2;
3146: newLocalPoints[p].index = -2;
3147: }
3148: /* Set subleaves */
3149: for (l = 0; l < numLeaves; ++l) {
3150: const PetscInt point = localPoints[l];
3151: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3153: if (subpoint < 0) continue;
3154: newLocalPoints[point-pStart].rank = rank;
3155: newLocalPoints[point-pStart].index = subpoint;
3156: ++numSubleaves;
3157: }
3158: /* Must put in owned subpoints */
3159: for (p = pStart; p < pEnd; ++p) {
3160: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3162: if (subpoint < 0) {
3163: newOwners[p-pStart].rank = -3;
3164: newOwners[p-pStart].index = -3;
3165: } else {
3166: newOwners[p-pStart].rank = rank;
3167: newOwners[p-pStart].index = subpoint;
3168: }
3169: }
3170: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3171: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3172: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3173: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3174: PetscMalloc1(numSubleaves, &slocalPoints);
3175: PetscMalloc1(numSubleaves, &sremotePoints);
3176: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3177: const PetscInt point = localPoints[l];
3178: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3180: if (subpoint < 0) continue;
3181: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3182: slocalPoints[sl] = subpoint;
3183: sremotePoints[sl].rank = newLocalPoints[point].rank;
3184: sremotePoints[sl].index = newLocalPoints[point].index;
3185: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3186: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3187: ++sl;
3188: }
3189: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3190: PetscFree2(newLocalPoints,newOwners);
3191: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3192: }
3193: if (subpIS) {
3194: ISRestoreIndices(subpIS, &subpoints);
3195: ISDestroy(&subpIS);
3196: }
3197: }
3198: /* Cleanup */
3199: for (d = 0; d <= dim; ++d) {
3200: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3201: ISDestroy(&subpointIS[d]);
3202: }
3203: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3204: return(0);
3205: }
3207: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3208: {
3212: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3213: return(0);
3214: }
3216: /*@
3217: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3219: Input Parameters:
3220: + dm - The original mesh
3221: . vertexLabel - The DMLabel marking points contained in the surface
3222: . value - The label value to use
3223: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3225: Output Parameter:
3226: . subdm - The surface mesh
3228: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3230: Level: developer
3232: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3233: @*/
3234: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3235: {
3236: PetscInt dim, cdim, depth;
3242: DMGetDimension(dm, &dim);
3243: DMPlexGetDepth(dm, &depth);
3244: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3245: DMSetType(*subdm, DMPLEX);
3246: DMSetDimension(*subdm, dim-1);
3247: DMGetCoordinateDim(dm, &cdim);
3248: DMSetCoordinateDim(*subdm, cdim);
3249: if (depth == dim) {
3250: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3251: } else {
3252: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3253: }
3254: return(0);
3255: }
3257: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3258: {
3259: MPI_Comm comm;
3260: DMLabel subpointMap;
3261: IS subvertexIS;
3262: const PetscInt *subVertices;
3263: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3264: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3265: PetscInt cMax, c, f;
3266: PetscErrorCode ierr;
3269: PetscObjectGetComm((PetscObject)dm, &comm);
3270: /* Create subpointMap which marks the submesh */
3271: DMLabelCreate("subpoint_map", &subpointMap);
3272: DMPlexSetSubpointMap(subdm, subpointMap);
3273: DMLabelDestroy(&subpointMap);
3274: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3275: /* Setup chart */
3276: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3277: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3278: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3279: DMPlexSetVTKCellHeight(subdm, 1);
3280: /* Set cone sizes */
3281: firstSubVertex = numSubCells;
3282: firstSubFace = numSubCells+numSubVertices;
3283: newFacePoint = firstSubFace;
3284: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3285: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3286: for (c = 0; c < numSubCells; ++c) {
3287: DMPlexSetConeSize(subdm, c, 1);
3288: }
3289: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3290: DMPlexSetConeSize(subdm, f, nFV);
3291: }
3292: DMSetUp(subdm);
3293: /* Create face cones */
3294: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3295: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3296: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3297: for (c = 0; c < numSubCells; ++c) {
3298: const PetscInt cell = subCells[c];
3299: const PetscInt subcell = c;
3300: const PetscInt *cone, *cells;
3301: PetscInt numCells, subVertex, p, v;
3303: if (cell < cMax) continue;
3304: DMPlexGetCone(dm, cell, &cone);
3305: for (v = 0; v < nFV; ++v) {
3306: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3307: subface[v] = firstSubVertex+subVertex;
3308: }
3309: DMPlexSetCone(subdm, newFacePoint, subface);
3310: DMPlexSetCone(subdm, subcell, &newFacePoint);
3311: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3312: /* Not true in parallel
3313: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3314: for (p = 0; p < numCells; ++p) {
3315: PetscInt negsubcell;
3317: if (cells[p] >= cMax) continue;
3318: /* I know this is a crap search */
3319: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3320: if (subCells[negsubcell] == cells[p]) break;
3321: }
3322: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3323: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3324: }
3325: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3326: ++newFacePoint;
3327: }
3328: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3329: DMPlexSymmetrize(subdm);
3330: DMPlexStratify(subdm);
3331: /* Build coordinates */
3332: {
3333: PetscSection coordSection, subCoordSection;
3334: Vec coordinates, subCoordinates;
3335: PetscScalar *coords, *subCoords;
3336: PetscInt cdim, numComp, coordSize, v;
3337: const char *name;
3339: DMGetCoordinateDim(dm, &cdim);
3340: DMGetCoordinateSection(dm, &coordSection);
3341: DMGetCoordinatesLocal(dm, &coordinates);
3342: DMGetCoordinateSection(subdm, &subCoordSection);
3343: PetscSectionSetNumFields(subCoordSection, 1);
3344: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3345: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3346: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3347: for (v = 0; v < numSubVertices; ++v) {
3348: const PetscInt vertex = subVertices[v];
3349: const PetscInt subvertex = firstSubVertex+v;
3350: PetscInt dof;
3352: PetscSectionGetDof(coordSection, vertex, &dof);
3353: PetscSectionSetDof(subCoordSection, subvertex, dof);
3354: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3355: }
3356: PetscSectionSetUp(subCoordSection);
3357: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3358: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3359: PetscObjectGetName((PetscObject)coordinates,&name);
3360: PetscObjectSetName((PetscObject)subCoordinates,name);
3361: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3362: VecSetBlockSize(subCoordinates, cdim);
3363: VecSetType(subCoordinates,VECSTANDARD);
3364: VecGetArray(coordinates, &coords);
3365: VecGetArray(subCoordinates, &subCoords);
3366: for (v = 0; v < numSubVertices; ++v) {
3367: const PetscInt vertex = subVertices[v];
3368: const PetscInt subvertex = firstSubVertex+v;
3369: PetscInt dof, off, sdof, soff, d;
3371: PetscSectionGetDof(coordSection, vertex, &dof);
3372: PetscSectionGetOffset(coordSection, vertex, &off);
3373: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3374: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3375: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3376: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3377: }
3378: VecRestoreArray(coordinates, &coords);
3379: VecRestoreArray(subCoordinates, &subCoords);
3380: DMSetCoordinatesLocal(subdm, subCoordinates);
3381: VecDestroy(&subCoordinates);
3382: }
3383: /* Build SF */
3384: CHKMEMQ;
3385: {
3386: PetscSF sfPoint, sfPointSub;
3387: const PetscSFNode *remotePoints;
3388: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3389: const PetscInt *localPoints;
3390: PetscInt *slocalPoints;
3391: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3392: PetscMPIInt rank;
3394: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3395: DMGetPointSF(dm, &sfPoint);
3396: DMGetPointSF(subdm, &sfPointSub);
3397: DMPlexGetChart(dm, &pStart, &pEnd);
3398: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3399: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3400: if (numRoots >= 0) {
3401: /* Only vertices should be shared */
3402: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3403: for (p = 0; p < pEnd-pStart; ++p) {
3404: newLocalPoints[p].rank = -2;
3405: newLocalPoints[p].index = -2;
3406: }
3407: /* Set subleaves */
3408: for (l = 0; l < numLeaves; ++l) {
3409: const PetscInt point = localPoints[l];
3410: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3412: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3413: if (subPoint < 0) continue;
3414: newLocalPoints[point-pStart].rank = rank;
3415: newLocalPoints[point-pStart].index = subPoint;
3416: ++numSubLeaves;
3417: }
3418: /* Must put in owned subpoints */
3419: for (p = pStart; p < pEnd; ++p) {
3420: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3422: if (subPoint < 0) {
3423: newOwners[p-pStart].rank = -3;
3424: newOwners[p-pStart].index = -3;
3425: } else {
3426: newOwners[p-pStart].rank = rank;
3427: newOwners[p-pStart].index = subPoint;
3428: }
3429: }
3430: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3431: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3432: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3433: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3434: PetscMalloc1(numSubLeaves, &slocalPoints);
3435: PetscMalloc1(numSubLeaves, &sremotePoints);
3436: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3437: const PetscInt point = localPoints[l];
3438: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3440: if (subPoint < 0) continue;
3441: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3442: slocalPoints[sl] = subPoint;
3443: sremotePoints[sl].rank = newLocalPoints[point].rank;
3444: sremotePoints[sl].index = newLocalPoints[point].index;
3445: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3446: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3447: ++sl;
3448: }
3449: PetscFree2(newLocalPoints,newOwners);
3450: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3451: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3452: }
3453: }
3454: CHKMEMQ;
3455: /* Cleanup */
3456: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3457: ISDestroy(&subvertexIS);
3458: PetscFree(subCells);
3459: return(0);
3460: }
3462: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3463: {
3464: DMLabel label = NULL;
3468: if (labelname) {DMGetLabel(dm, labelname, &label);}
3469: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3470: return(0);
3471: }
3473: /*@C
3474: 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.
3476: Input Parameters:
3477: + dm - The original mesh
3478: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3479: . label - A label name, or NULL
3480: - value - A label value
3482: Output Parameter:
3483: . subdm - The surface mesh
3485: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3487: Level: developer
3489: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3490: @*/
3491: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3492: {
3493: PetscInt dim, cdim, depth;
3499: DMGetDimension(dm, &dim);
3500: DMPlexGetDepth(dm, &depth);
3501: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3502: DMSetType(*subdm, DMPLEX);
3503: DMSetDimension(*subdm, dim-1);
3504: DMGetCoordinateDim(dm, &cdim);
3505: DMSetCoordinateDim(*subdm, cdim);
3506: if (depth == dim) {
3507: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3508: } else {
3509: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3510: }
3511: return(0);
3512: }
3514: /*@
3515: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3517: Input Parameters:
3518: + dm - The original mesh
3519: . cellLabel - The DMLabel marking cells contained in the new mesh
3520: - value - The label value to use
3522: Output Parameter:
3523: . subdm - The new mesh
3525: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3527: Level: developer
3529: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3530: @*/
3531: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3532: {
3533: PetscInt dim;
3539: DMGetDimension(dm, &dim);
3540: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3541: DMSetType(*subdm, DMPLEX);
3542: DMSetDimension(*subdm, dim);
3543: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3544: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3545: return(0);
3546: }
3548: /*@
3549: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3551: Input Parameter:
3552: . dm - The submesh DM
3554: Output Parameter:
3555: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3557: Level: developer
3559: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3560: @*/
3561: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3562: {
3566: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3567: return(0);
3568: }
3570: /*@
3571: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3573: Input Parameters:
3574: + dm - The submesh DM
3575: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3577: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3579: Level: developer
3581: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3582: @*/
3583: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3584: {
3585: DM_Plex *mesh = (DM_Plex *) dm->data;
3586: DMLabel tmp;
3591: tmp = mesh->subpointMap;
3592: mesh->subpointMap = subpointMap;
3593: ++mesh->subpointMap->refct;
3594: DMLabelDestroy(&tmp);
3595: return(0);
3596: }
3598: /*@
3599: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3601: Input Parameter:
3602: . dm - The submesh DM
3604: Output Parameter:
3605: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3607: Note: This IS is guaranteed to be sorted by the construction of the submesh
3609: Level: developer
3611: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3612: @*/
3613: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3614: {
3615: MPI_Comm comm;
3616: DMLabel subpointMap;
3617: IS is;
3618: const PetscInt *opoints;
3619: PetscInt *points, *depths;
3620: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3621: PetscErrorCode ierr;
3626: PetscObjectGetComm((PetscObject)dm,&comm);
3627: *subpointIS = NULL;
3628: DMPlexGetSubpointMap(dm, &subpointMap);
3629: DMPlexGetDepth(dm, &depth);
3630: if (subpointMap && depth >= 0) {
3631: DMPlexGetChart(dm, &pStart, &pEnd);
3632: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3633: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3634: depths[0] = depth;
3635: depths[1] = 0;
3636: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3637: PetscMalloc1(pEnd, &points);
3638: for(d = 0, off = 0; d <= depth; ++d) {
3639: const PetscInt dep = depths[d];
3641: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3642: DMLabelGetStratumSize(subpointMap, dep, &n);
3643: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3644: 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);
3645: } else {
3646: if (!n) {
3647: if (d == 0) {
3648: /* Missing cells */
3649: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3650: } else {
3651: /* Missing faces */
3652: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3653: }
3654: }
3655: }
3656: if (n) {
3657: DMLabelGetStratumIS(subpointMap, dep, &is);
3658: ISGetIndices(is, &opoints);
3659: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3660: ISRestoreIndices(is, &opoints);
3661: ISDestroy(&is);
3662: }
3663: }
3664: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3665: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3666: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3667: }
3668: return(0);
3669: }
3671: /*@
3672: DMPlexGetSubpoint - Return the subpoint corresponding to a point in the original mesh. If the DM
3673: is not a submesh, just return the input point.
3675: Note collective
3677: Input Parameters:
3678: + dm - The submesh DM
3679: - p - The point in the original, from which the submesh was created
3681: Output Parameter:
3682: . subp - The point in the submesh
3684: Level: developer
3686: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS()
3687: @*/
3688: PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp)
3689: {
3690: DMLabel spmap;
3693: *subp = p;
3694: DMPlexGetSubpointMap(dm, &spmap);
3695: if (spmap) {
3696: IS subpointIS;
3697: const PetscInt *subpoints;
3698: PetscInt numSubpoints;
3700: /* TODO Cache the IS, making it look like an index */
3701: DMPlexCreateSubpointIS(dm, &subpointIS);
3702: ISGetLocalSize(subpointIS, &numSubpoints);
3703: ISGetIndices(subpointIS, &subpoints);
3704: PetscFindInt(p, numSubpoints, subpoints, subp);
3705: if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p);
3706: ISRestoreIndices(subpointIS, &subpoints);
3707: ISDestroy(&subpointIS);
3708: }
3709: return(0);
3710: }