Actual source code: plexsubmesh.c
petsc-3.9.4 2018-09-11
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;
352: DMPlexGetDepth(dm, &depth);
353: if (depth < 0) return(0);
354: /* Step 1: Expand chart */
355: DMPlexGetChart(dm, &pStart, &pEnd);
356: pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
357: DMPlexSetChart(dmNew, pStart, pEnd);
358: /* Step 2: Set cone and support sizes */
359: for (d = 0; d <= depth; ++d) {
360: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
361: for (p = pStart; p < pEnd; ++p) {
362: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
363: PetscInt size;
365: DMPlexGetConeSize(dm, p, &size);
366: DMPlexSetConeSize(dmNew, newp, size);
367: DMPlexGetSupportSize(dm, p, &size);
368: DMPlexSetSupportSize(dmNew, newp, size);
369: }
370: }
371: return(0);
372: }
374: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
375: {
376: PetscInt *newpoints;
377: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
381: DMPlexGetDepth(dm, &depth);
382: if (depth < 0) return(0);
383: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
384: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
385: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints);
386: /* Step 5: Set cones and supports */
387: DMPlexGetChart(dm, &pStart, &pEnd);
388: for (p = pStart; p < pEnd; ++p) {
389: const PetscInt *points = NULL, *orientations = NULL;
390: PetscInt size,sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
392: DMPlexGetConeSize(dm, p, &size);
393: DMPlexGetCone(dm, p, &points);
394: DMPlexGetConeOrientation(dm, p, &orientations);
395: for (i = 0; i < size; ++i) {
396: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
397: }
398: DMPlexSetCone(dmNew, newp, newpoints);
399: DMPlexSetConeOrientation(dmNew, newp, orientations);
400: DMPlexGetSupportSize(dm, p, &size);
401: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
402: DMPlexGetSupport(dm, p, &points);
403: for (i = 0; i < size; ++i) {
404: newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
405: }
406: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
407: DMPlexSetSupport(dmNew, newp, newpoints);
408: }
409: PetscFree(newpoints);
410: return(0);
411: }
413: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
414: {
415: PetscSection coordSection, newCoordSection;
416: Vec coordinates, newCoordinates;
417: PetscScalar *coords, *newCoords;
418: PetscInt coordSize, sStart, sEnd;
419: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
420: PetscBool hasCells;
424: DMGetCoordinateDim(dm, &dim);
425: DMSetCoordinateDim(dmNew, dim);
426: DMPlexGetDepth(dm, &depth);
427: /* Step 8: Convert coordinates */
428: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
429: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
430: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
431: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
432: DMGetCoordinateSection(dm, &coordSection);
433: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
434: PetscSectionSetNumFields(newCoordSection, 1);
435: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
436: PetscSectionGetChart(coordSection, &sStart, &sEnd);
437: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
438: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
439: if (hasCells) {
440: for (c = cStart; c < cEnd; ++c) {
441: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
443: PetscSectionGetDof(coordSection, c, &dof);
444: PetscSectionSetDof(newCoordSection, cNew, dof);
445: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
446: }
447: }
448: for (v = vStartNew; v < vEndNew; ++v) {
449: PetscSectionSetDof(newCoordSection, v, dim);
450: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
451: }
452: PetscSectionSetUp(newCoordSection);
453: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
454: PetscSectionGetStorageSize(newCoordSection, &coordSize);
455: VecCreate(PETSC_COMM_SELF, &newCoordinates);
456: PetscObjectSetName((PetscObject) newCoordinates, "coordinates");
457: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
458: VecSetBlockSize(newCoordinates, dim);
459: VecSetType(newCoordinates,VECSTANDARD);
460: DMSetCoordinatesLocal(dmNew, newCoordinates);
461: DMGetCoordinatesLocal(dm, &coordinates);
462: VecGetArray(coordinates, &coords);
463: VecGetArray(newCoordinates, &newCoords);
464: if (hasCells) {
465: for (c = cStart; c < cEnd; ++c) {
466: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
468: PetscSectionGetDof(coordSection, c, &dof);
469: PetscSectionGetOffset(coordSection, c, &off);
470: PetscSectionGetOffset(newCoordSection, cNew, &noff);
471: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
472: }
473: }
474: for (v = vStart; v < vEnd; ++v) {
475: PetscInt dof, off, noff, d;
477: PetscSectionGetDof(coordSection, v, &dof);
478: PetscSectionGetOffset(coordSection, v, &off);
479: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
480: for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d];
481: }
482: VecRestoreArray(coordinates, &coords);
483: VecRestoreArray(newCoordinates, &newCoords);
484: VecDestroy(&newCoordinates);
485: PetscSectionDestroy(&newCoordSection);
486: return(0);
487: }
489: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
490: {
491: PetscInt depth = 0;
492: PetscSF sfPoint, sfPointNew;
493: const PetscSFNode *remotePoints;
494: PetscSFNode *gremotePoints;
495: const PetscInt *localPoints;
496: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
497: PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0;
498: PetscErrorCode ierr;
501: DMPlexGetDepth(dm, &depth);
502: /* Step 9: Convert pointSF */
503: DMGetPointSF(dm, &sfPoint);
504: DMGetPointSF(dmNew, &sfPointNew);
505: DMPlexGetChart(dm, &pStart, &pEnd);
506: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
507: totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd;
508: if (numRoots >= 0) {
509: PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
510: for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
511: PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
512: PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
513: PetscMalloc1(numLeaves, &glocalPoints);
514: PetscMalloc1(numLeaves, &gremotePoints);
515: for (l = 0; l < numLeaves; ++l) {
516: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
517: gremotePoints[l].rank = remotePoints[l].rank;
518: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
519: }
520: PetscFree2(newLocation,newRemoteLocation);
521: PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
522: }
523: return(0);
524: }
526: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
527: {
528: PetscSF sfPoint;
529: DMLabel vtkLabel, ghostLabel;
530: const PetscSFNode *leafRemote;
531: const PetscInt *leafLocal;
532: PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f;
533: PetscMPIInt rank;
534: PetscErrorCode ierr;
537: DMPlexGetDepth(dm, &depth);
538: /* Step 10: Convert labels */
539: DMGetNumLabels(dm, &numLabels);
540: for (l = 0; l < numLabels; ++l) {
541: DMLabel label, newlabel;
542: const char *lname;
543: PetscBool isDepth;
544: IS valueIS;
545: const PetscInt *values;
546: PetscInt numValues, val;
548: DMGetLabelName(dm, l, &lname);
549: PetscStrcmp(lname, "depth", &isDepth);
550: if (isDepth) continue;
551: DMCreateLabel(dmNew, lname);
552: DMGetLabel(dm, lname, &label);
553: DMGetLabel(dmNew, lname, &newlabel);
554: DMLabelGetDefaultValue(label,&val);
555: DMLabelSetDefaultValue(newlabel,val);
556: DMLabelGetValueIS(label, &valueIS);
557: ISGetLocalSize(valueIS, &numValues);
558: ISGetIndices(valueIS, &values);
559: for (val = 0; val < numValues; ++val) {
560: IS pointIS;
561: const PetscInt *points;
562: PetscInt numPoints, p;
564: DMLabelGetStratumIS(label, values[val], &pointIS);
565: ISGetLocalSize(pointIS, &numPoints);
566: ISGetIndices(pointIS, &points);
567: for (p = 0; p < numPoints; ++p) {
568: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
570: DMLabelSetValue(newlabel, newpoint, values[val]);
571: }
572: ISRestoreIndices(pointIS, &points);
573: ISDestroy(&pointIS);
574: }
575: ISRestoreIndices(valueIS, &values);
576: ISDestroy(&valueIS);
577: }
578: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
579: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
580: DMGetPointSF(dm, &sfPoint);
581: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
582: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
583: DMCreateLabel(dmNew, "vtk");
584: DMCreateLabel(dmNew, "ghost");
585: DMGetLabel(dmNew, "vtk", &vtkLabel);
586: DMGetLabel(dmNew, "ghost", &ghostLabel);
587: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
588: for (; c < leafLocal[l] && c < cEnd; ++c) {
589: DMLabelSetValue(vtkLabel, c, 1);
590: }
591: if (leafLocal[l] >= cEnd) break;
592: if (leafRemote[l].rank == rank) {
593: DMLabelSetValue(vtkLabel, c, 1);
594: } else {
595: DMLabelSetValue(ghostLabel, c, 2);
596: }
597: }
598: for (; c < cEnd; ++c) {
599: DMLabelSetValue(vtkLabel, c, 1);
600: }
601: if (0) {
602: DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);
603: }
604: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
605: for (f = fStart; f < fEnd; ++f) {
606: PetscInt numCells;
608: DMPlexGetSupportSize(dmNew, f, &numCells);
609: if (numCells < 2) {
610: DMLabelSetValue(ghostLabel, f, 1);
611: } else {
612: const PetscInt *cells = NULL;
613: PetscInt vA, vB;
615: DMPlexGetSupport(dmNew, f, &cells);
616: DMLabelGetValue(vtkLabel, cells[0], &vA);
617: DMLabelGetValue(vtkLabel, cells[1], &vB);
618: if (vA != 1 && vB != 1) {DMLabelSetValue(ghostLabel, f, 1);}
619: }
620: }
621: if (0) {
622: DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);
623: }
624: return(0);
625: }
627: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
628: {
629: DM refTree;
630: PetscSection pSec;
631: PetscInt *parents, *childIDs;
635: DMPlexGetReferenceTree(dm,&refTree);
636: DMPlexSetReferenceTree(dmNew,refTree);
637: DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);
638: if (pSec) {
639: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
640: PetscInt *childIDsShifted;
641: PetscSection pSecShifted;
643: PetscSectionGetChart(pSec,&pStart,&pEnd);
644: DMPlexGetDepth(dm,&depth);
645: pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift);
646: pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift);
647: PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);
648: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);
649: PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);
650: for (p = pStartShifted; p < pEndShifted; p++) {
651: /* start off assuming no children */
652: PetscSectionSetDof(pSecShifted,p,0);
653: }
654: for (p = pStart; p < pEnd; p++) {
655: PetscInt dof;
656: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
658: PetscSectionGetDof(pSec,p,&dof);
659: PetscSectionSetDof(pSecShifted,pNew,dof);
660: }
661: PetscSectionSetUp(pSecShifted);
662: for (p = pStart; p < pEnd; p++) {
663: PetscInt dof;
664: PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift);
666: PetscSectionGetDof(pSec,p,&dof);
667: if (dof) {
668: PetscInt off, offNew;
670: PetscSectionGetOffset(pSec,p,&off);
671: PetscSectionGetOffset(pSecShifted,pNew,&offNew);
672: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift);
673: childIDsShifted[offNew] = childIDs[off];
674: }
675: }
676: DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);
677: PetscFree2(parentsShifted,childIDsShifted);
678: PetscSectionDestroy(&pSecShifted);
679: }
680: return(0);
681: }
683: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
684: {
685: PetscSF sf;
686: IS valueIS;
687: const PetscInt *values, *leaves;
688: PetscInt *depthShift;
689: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
690: PetscBool isper;
691: const PetscReal *maxCell, *L;
692: const DMBoundaryType *bd;
693: PetscErrorCode ierr;
696: DMGetPointSF(dm, &sf);
697: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
698: nleaves = PetscMax(0, nleaves);
699: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
700: /* Count ghost cells */
701: DMLabelGetValueIS(label, &valueIS);
702: ISGetLocalSize(valueIS, &numFS);
703: ISGetIndices(valueIS, &values);
704: Ng = 0;
705: for (fs = 0; fs < numFS; ++fs) {
706: IS faceIS;
707: const PetscInt *faces;
708: PetscInt numFaces, f, numBdFaces = 0;
710: DMLabelGetStratumIS(label, values[fs], &faceIS);
711: ISGetLocalSize(faceIS, &numFaces);
712: ISGetIndices(faceIS, &faces);
713: for (f = 0; f < numFaces; ++f) {
714: PetscInt numChildren;
716: PetscFindInt(faces[f], nleaves, leaves, &loc);
717: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
718: /* non-local and ancestors points don't get to register ghosts */
719: if (loc >= 0 || numChildren) continue;
720: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
721: }
722: Ng += numBdFaces;
723: ISDestroy(&faceIS);
724: }
725: DMPlexGetDepth(dm, &depth);
726: PetscMalloc1(2*(depth+1), &depthShift);
727: for (d = 0; d <= depth; d++) {
728: PetscInt dEnd;
730: DMPlexGetDepthStratum(dm,d,NULL,&dEnd);
731: depthShift[2*d] = dEnd;
732: depthShift[2*d+1] = 0;
733: }
734: if (depth >= 0) depthShift[2*depth+1] = Ng;
735: DMPlexShiftPointSetUp_Internal(depth,depthShift);
736: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
737: /* Step 3: Set cone/support sizes for new points */
738: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
739: DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
740: for (c = cEnd; c < cEnd + Ng; ++c) {
741: DMPlexSetConeSize(gdm, c, 1);
742: }
743: for (fs = 0; fs < numFS; ++fs) {
744: IS faceIS;
745: const PetscInt *faces;
746: PetscInt numFaces, f;
748: DMLabelGetStratumIS(label, values[fs], &faceIS);
749: ISGetLocalSize(faceIS, &numFaces);
750: ISGetIndices(faceIS, &faces);
751: for (f = 0; f < numFaces; ++f) {
752: PetscInt size, numChildren;
754: PetscFindInt(faces[f], nleaves, leaves, &loc);
755: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
756: if (loc >= 0 || numChildren) continue;
757: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
758: DMPlexGetSupportSize(dm, faces[f], &size);
759: if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size);
760: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
761: }
762: ISRestoreIndices(faceIS, &faces);
763: ISDestroy(&faceIS);
764: }
765: /* Step 4: Setup ghosted DM */
766: DMSetUp(gdm);
767: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
768: /* Step 6: Set cones and supports for new points */
769: ghostCell = cEnd;
770: for (fs = 0; fs < numFS; ++fs) {
771: IS faceIS;
772: const PetscInt *faces;
773: PetscInt numFaces, f;
775: DMLabelGetStratumIS(label, values[fs], &faceIS);
776: ISGetLocalSize(faceIS, &numFaces);
777: ISGetIndices(faceIS, &faces);
778: for (f = 0; f < numFaces; ++f) {
779: PetscInt newFace = faces[f] + Ng, numChildren;
781: PetscFindInt(faces[f], nleaves, leaves, &loc);
782: DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);
783: if (loc >= 0 || numChildren) continue;
784: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
785: DMPlexSetCone(gdm, ghostCell, &newFace);
786: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
787: ++ghostCell;
788: }
789: ISRestoreIndices(faceIS, &faces);
790: ISDestroy(&faceIS);
791: }
792: ISRestoreIndices(valueIS, &values);
793: ISDestroy(&valueIS);
794: /* Step 7: Stratify */
795: DMPlexStratify(gdm);
796: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
797: DMPlexShiftSF_Internal(dm, depthShift, gdm);
798: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
799: DMPlexShiftTree_Internal(dm, depthShift, gdm);
800: PetscFree(depthShift);
801: /* Step 7: Periodicity */
802: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
803: DMSetPeriodicity(gdm, isper, maxCell, L, bd);
804: if (numGhostCells) *numGhostCells = Ng;
805: return(0);
806: }
808: /*@C
809: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
811: Collective on dm
813: Input Parameters:
814: + dm - The original DM
815: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
817: Output Parameters:
818: + numGhostCells - The number of ghost cells added to the DM
819: - dmGhosted - The new DM
821: Note: If no label exists of that name, one will be created marking all boundary faces
823: Level: developer
825: .seealso: DMCreate()
826: @*/
827: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
828: {
829: DM gdm;
830: DMLabel label;
831: const char *name = labelName ? labelName : "Face Sets";
832: PetscInt dim;
833: PetscBool flag;
840: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
841: DMSetType(gdm, DMPLEX);
842: DMGetDimension(dm, &dim);
843: DMSetDimension(gdm, dim);
844: DMPlexGetAdjacencyUseCone(dm, &flag);
845: DMPlexSetAdjacencyUseCone(gdm, flag);
846: DMPlexGetAdjacencyUseClosure(dm, &flag);
847: DMPlexSetAdjacencyUseClosure(gdm, flag);
848: DMGetLabel(dm, name, &label);
849: if (!label) {
850: /* Get label for boundary faces */
851: DMCreateLabel(dm, name);
852: DMGetLabel(dm, name, &label);
853: DMPlexMarkBoundaryFaces(dm, 1, label);
854: }
855: DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);
856: DMCopyBoundary(dm, gdm);
857: *dmGhosted = gdm;
858: return(0);
859: }
861: /*
862: We are adding three kinds of points here:
863: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
864: Non-replicated: Points which exist on the fault, but are not replicated
865: Hybrid: Entirely new points, such as cohesive cells
867: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
868: each depth so that the new split/hybrid points can be inserted as a block.
869: */
870: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm)
871: {
872: MPI_Comm comm;
873: IS valueIS;
874: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
875: const PetscInt *values; /* List of depths for which we have replicated points */
876: IS *splitIS;
877: IS *unsplitIS;
878: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
879: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
880: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
881: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
882: const PetscInt **splitPoints; /* Replicated points for each depth */
883: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
884: PetscSection coordSection;
885: Vec coordinates;
886: PetscScalar *coords;
887: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
888: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
889: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
890: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
891: PetscInt *coneNew, *coneONew, *supportNew;
892: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
893: PetscErrorCode ierr;
896: PetscObjectGetComm((PetscObject)dm,&comm);
897: DMGetDimension(dm, &dim);
898: DMPlexGetDepth(dm, &depth);
899: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
900: /* Count split points and add cohesive cells */
901: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
902: PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld);
903: PetscMalloc7(depth+1,&splitIS,depth+1,&unsplitIS,depth+1,&numSplitPoints,depth+1,&numUnsplitPoints,depth+1,&numHybridPoints,depth+1,&splitPoints,depth+1,&unsplitPoints);
904: DMPlexGetHybridBounds(dm, depth >= 0 ? &depthMax[depth] : NULL, depth>1 ? &depthMax[depth-1] : NULL, depth>2 ? &depthMax[1] : NULL, depth >= 0 ? &depthMax[0] : NULL);
905: for (d = 0; d <= depth; ++d) {
906: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
907: depthEnd[d] = pMaxNew[d];
908: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
909: numSplitPoints[d] = 0;
910: numUnsplitPoints[d] = 0;
911: numHybridPoints[d] = 0;
912: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
913: splitPoints[d] = NULL;
914: unsplitPoints[d] = NULL;
915: splitIS[d] = NULL;
916: unsplitIS[d] = NULL;
917: /* we are shifting the existing hybrid points with the stratum behind them, so
918: * the split comes at the end of the normal points, i.e., at depthMax[d] */
919: depthShift[2*d] = depthMax[d];
920: depthShift[2*d+1] = 0;
921: }
922: if (label) {
923: DMLabelGetValueIS(label, &valueIS);
924: ISGetLocalSize(valueIS, &numSP);
925: ISGetIndices(valueIS, &values);
926: }
927: for (sp = 0; sp < numSP; ++sp) {
928: const PetscInt dep = values[sp];
930: if ((dep < 0) || (dep > depth)) continue;
931: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
932: if (splitIS[dep]) {
933: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
934: ISGetIndices(splitIS[dep], &splitPoints[dep]);
935: }
936: DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep]);
937: if (unsplitIS[dep]) {
938: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
939: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
940: }
941: }
942: /* Calculate number of hybrid points */
943: 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 */
944: for (d = 0; d <= depth; ++d) depthShift[2*d+1] = numSplitPoints[d] + numHybridPoints[d];
945: DMPlexShiftPointSetUp_Internal(depth,depthShift);
946: /* the end of the points in this stratum that come before the new points:
947: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
948: * added points */
949: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d],depth,depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
950: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
951: /* Step 3: Set cone/support sizes for new points */
952: for (dep = 0; dep <= depth; ++dep) {
953: for (p = 0; p < numSplitPoints[dep]; ++p) {
954: const PetscInt oldp = splitPoints[dep][p];
955: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
956: const PetscInt splitp = p + pMaxNew[dep];
957: const PetscInt *support;
958: PetscInt coneSize, supportSize, qf, qn, qp, e;
960: DMPlexGetConeSize(dm, oldp, &coneSize);
961: DMPlexSetConeSize(sdm, splitp, coneSize);
962: DMPlexGetSupportSize(dm, oldp, &supportSize);
963: DMPlexSetSupportSize(sdm, splitp, supportSize);
964: if (dep == depth-1) {
965: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
967: /* Add cohesive cells, they are prisms */
968: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
969: } else if (dep == 0) {
970: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
972: DMPlexGetSupport(dm, oldp, &support);
973: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
974: PetscInt val;
976: DMLabelGetValue(label, support[e], &val);
977: if (val == 1) ++qf;
978: if ((val == 1) || (val == (shift + 1))) ++qn;
979: if ((val == 1) || (val == -(shift + 1))) ++qp;
980: }
981: /* Split old vertex: Edges into original vertex and new cohesive edge */
982: DMPlexSetSupportSize(sdm, newp, qn+1);
983: /* Split new vertex: Edges into split vertex and new cohesive edge */
984: DMPlexSetSupportSize(sdm, splitp, qp+1);
985: /* Add hybrid edge */
986: DMPlexSetConeSize(sdm, hybedge, 2);
987: DMPlexSetSupportSize(sdm, hybedge, qf);
988: } else if (dep == dim-2) {
989: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
991: DMPlexGetSupport(dm, oldp, &support);
992: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
993: PetscInt val;
995: DMLabelGetValue(label, support[e], &val);
996: if (val == dim-1) ++qf;
997: if ((val == dim-1) || (val == (shift + dim-1))) ++qn;
998: if ((val == dim-1) || (val == -(shift + dim-1))) ++qp;
999: }
1000: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1001: DMPlexSetSupportSize(sdm, newp, qn+1);
1002: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1003: DMPlexSetSupportSize(sdm, splitp, qp+1);
1004: /* Add hybrid face */
1005: DMPlexSetConeSize(sdm, hybface, 4);
1006: DMPlexSetSupportSize(sdm, hybface, qf);
1007: }
1008: }
1009: }
1010: for (dep = 0; dep <= depth; ++dep) {
1011: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1012: const PetscInt oldp = unsplitPoints[dep][p];
1013: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1014: const PetscInt *support;
1015: PetscInt coneSize, supportSize, qf, e, s;
1017: DMPlexGetConeSize(dm, oldp, &coneSize);
1018: DMPlexGetSupportSize(dm, oldp, &supportSize);
1019: DMPlexGetSupport(dm, oldp, &support);
1020: if (dep == 0) {
1021: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1023: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1024: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1025: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1026: if (e >= 0) ++qf;
1027: }
1028: DMPlexSetSupportSize(sdm, newp, qf+2);
1029: /* Add hybrid edge */
1030: DMPlexSetConeSize(sdm, hybedge, 2);
1031: for (e = 0, qf = 0; e < supportSize; ++e) {
1032: PetscInt val;
1034: DMLabelGetValue(label, support[e], &val);
1035: /* Split and unsplit edges produce hybrid faces */
1036: if (val == 1) ++qf;
1037: if (val == (shift2 + 1)) ++qf;
1038: }
1039: DMPlexSetSupportSize(sdm, hybedge, qf);
1040: } else if (dep == dim-2) {
1041: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1042: PetscInt val;
1044: for (e = 0, qf = 0; e < supportSize; ++e) {
1045: DMLabelGetValue(label, support[e], &val);
1046: if (val == dim-1) qf += 2;
1047: else ++qf;
1048: }
1049: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1050: DMPlexSetSupportSize(sdm, newp, qf+2);
1051: /* Add hybrid face */
1052: for (e = 0, qf = 0; e < supportSize; ++e) {
1053: DMLabelGetValue(label, support[e], &val);
1054: if (val == dim-1) ++qf;
1055: }
1056: DMPlexSetConeSize(sdm, hybface, 4);
1057: DMPlexSetSupportSize(sdm, hybface, qf);
1058: }
1059: }
1060: }
1061: /* Step 4: Setup split DM */
1062: DMSetUp(sdm);
1063: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1064: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1065: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew)*3,&coneNew,PetscMax(maxConeSize, maxConeSizeNew)*3,&coneONew,PetscMax(maxSupportSize, maxSupportSizeNew),&supportNew);
1066: /* Step 6: Set cones and supports for new points */
1067: for (dep = 0; dep <= depth; ++dep) {
1068: for (p = 0; p < numSplitPoints[dep]; ++p) {
1069: const PetscInt oldp = splitPoints[dep][p];
1070: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1071: const PetscInt splitp = p + pMaxNew[dep];
1072: const PetscInt *cone, *support, *ornt;
1073: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1075: DMPlexGetConeSize(dm, oldp, &coneSize);
1076: DMPlexGetCone(dm, oldp, &cone);
1077: DMPlexGetConeOrientation(dm, oldp, &ornt);
1078: DMPlexGetSupportSize(dm, oldp, &supportSize);
1079: DMPlexGetSupport(dm, oldp, &support);
1080: if (dep == depth-1) {
1081: PetscBool hasUnsplit = PETSC_FALSE;
1082: const PetscInt hybcell = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1083: const PetscInt *supportF;
1085: /* Split face: copy in old face to new face to start */
1086: DMPlexGetSupport(sdm, newp, &supportF);
1087: DMPlexSetSupport(sdm, splitp, supportF);
1088: /* Split old face: old vertices/edges in cone so no change */
1089: /* Split new face: new vertices/edges in cone */
1090: for (q = 0; q < coneSize; ++q) {
1091: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1092: if (v < 0) {
1093: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1094: 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);
1095: coneNew[2+q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1096: hasUnsplit = PETSC_TRUE;
1097: } else {
1098: coneNew[2+q] = v + pMaxNew[dep-1];
1099: if (dep > 1) {
1100: const PetscInt *econe;
1101: PetscInt econeSize, r, vs, vu;
1103: DMPlexGetConeSize(dm, cone[q], &econeSize);
1104: DMPlexGetCone(dm, cone[q], &econe);
1105: for (r = 0; r < econeSize; ++r) {
1106: PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs);
1107: PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu);
1108: if (vs >= 0) continue;
1109: 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);
1110: hasUnsplit = PETSC_TRUE;
1111: }
1112: }
1113: }
1114: }
1115: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1116: DMPlexSetConeOrientation(sdm, splitp, ornt);
1117: /* Face support */
1118: for (s = 0; s < supportSize; ++s) {
1119: PetscInt val;
1121: DMLabelGetValue(label, support[s], &val);
1122: if (val < 0) {
1123: /* Split old face: Replace negative side cell with cohesive cell */
1124: DMPlexInsertSupport(sdm, newp, s, hybcell);
1125: } else {
1126: /* Split new face: Replace positive side cell with cohesive cell */
1127: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1128: /* Get orientation for cohesive face */
1129: {
1130: const PetscInt *ncone, *nconeO;
1131: PetscInt nconeSize, nc;
1133: DMPlexGetConeSize(dm, support[s], &nconeSize);
1134: DMPlexGetCone(dm, support[s], &ncone);
1135: DMPlexGetConeOrientation(dm, support[s], &nconeO);
1136: for (nc = 0; nc < nconeSize; ++nc) {
1137: if (ncone[nc] == oldp) {
1138: coneONew[0] = nconeO[nc];
1139: break;
1140: }
1141: }
1142: if (nc >= nconeSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %d in neighboring cell %d", oldp, support[s]);
1143: }
1144: }
1145: }
1146: /* Cohesive cell: Old and new split face, then new cohesive faces */
1147: coneNew[0] = newp; /* Extracted negative side orientation above */
1148: coneNew[1] = splitp;
1149: coneONew[1] = coneONew[0];
1150: for (q = 0; q < coneSize; ++q) {
1151: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1152: if (v < 0) {
1153: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1154: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1155: coneONew[2+q] = 0;
1156: } else {
1157: coneNew[2+q] = v + pMaxNew[dep] + numSplitPoints[dep];
1158: }
1159: coneONew[2+q] = 0;
1160: }
1161: DMPlexSetCone(sdm, hybcell, coneNew);
1162: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1163: /* Label the hybrid cells on the boundary of the split */
1164: if (hasUnsplit) {DMLabelSetValue(label, -hybcell, dim);}
1165: } else if (dep == 0) {
1166: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1168: /* Split old vertex: Edges in old split faces and new cohesive edge */
1169: for (e = 0, qn = 0; e < supportSize; ++e) {
1170: PetscInt val;
1172: DMLabelGetValue(label, support[e], &val);
1173: if ((val == 1) || (val == (shift + 1))) {
1174: supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1175: }
1176: }
1177: supportNew[qn] = hybedge;
1178: DMPlexSetSupport(sdm, newp, supportNew);
1179: /* Split new vertex: Edges in new split faces and new cohesive edge */
1180: for (e = 0, qp = 0; e < supportSize; ++e) {
1181: PetscInt val, edge;
1183: DMLabelGetValue(label, support[e], &val);
1184: if (val == 1) {
1185: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1186: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1187: supportNew[qp++] = edge + pMaxNew[dep+1];
1188: } else if (val == -(shift + 1)) {
1189: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1190: }
1191: }
1192: supportNew[qp] = hybedge;
1193: DMPlexSetSupport(sdm, splitp, supportNew);
1194: /* Hybrid edge: Old and new split vertex */
1195: coneNew[0] = newp;
1196: coneNew[1] = splitp;
1197: DMPlexSetCone(sdm, hybedge, coneNew);
1198: for (e = 0, qf = 0; e < supportSize; ++e) {
1199: PetscInt val, edge;
1201: DMLabelGetValue(label, support[e], &val);
1202: if (val == 1) {
1203: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1204: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1205: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1206: }
1207: }
1208: DMPlexSetSupport(sdm, hybedge, supportNew);
1209: } else if (dep == dim-2) {
1210: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1212: /* Split old edge: old vertices in cone so no change */
1213: /* Split new edge: new vertices in cone */
1214: for (q = 0; q < coneSize; ++q) {
1215: PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v);
1216: if (v < 0) {
1217: PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1218: 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);
1219: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1220: } else {
1221: coneNew[q] = v + pMaxNew[dep-1];
1222: }
1223: }
1224: DMPlexSetCone(sdm, splitp, coneNew);
1225: /* Split old edge: Faces in positive side cells and old split faces */
1226: for (e = 0, q = 0; e < supportSize; ++e) {
1227: PetscInt val;
1229: DMLabelGetValue(label, support[e], &val);
1230: if (val == dim-1) {
1231: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1232: } else if (val == (shift + dim-1)) {
1233: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1234: }
1235: }
1236: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1237: DMPlexSetSupport(sdm, newp, supportNew);
1238: /* Split new edge: Faces in negative side cells and new split faces */
1239: for (e = 0, q = 0; e < supportSize; ++e) {
1240: PetscInt val, face;
1242: DMLabelGetValue(label, support[e], &val);
1243: if (val == dim-1) {
1244: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1245: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1246: supportNew[q++] = face + pMaxNew[dep+1];
1247: } else if (val == -(shift + dim-1)) {
1248: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1249: }
1250: }
1251: supportNew[q++] = p + pMaxNew[dep+1] + numSplitPoints[dep+1];
1252: DMPlexSetSupport(sdm, splitp, supportNew);
1253: /* Hybrid face */
1254: coneNew[0] = newp;
1255: coneNew[1] = splitp;
1256: for (v = 0; v < coneSize; ++v) {
1257: PetscInt vertex;
1258: PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex);
1259: if (vertex < 0) {
1260: PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex);
1261: 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);
1262: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1263: } else {
1264: coneNew[2+v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1265: }
1266: }
1267: DMPlexSetCone(sdm, hybface, coneNew);
1268: for (e = 0, qf = 0; e < supportSize; ++e) {
1269: PetscInt val, face;
1271: DMLabelGetValue(label, support[e], &val);
1272: if (val == dim-1) {
1273: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1274: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]);
1275: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1276: }
1277: }
1278: DMPlexSetSupport(sdm, hybface, supportNew);
1279: }
1280: }
1281: }
1282: for (dep = 0; dep <= depth; ++dep) {
1283: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1284: const PetscInt oldp = unsplitPoints[dep][p];
1285: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1286: const PetscInt *cone, *support, *ornt;
1287: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1289: DMPlexGetConeSize(dm, oldp, &coneSize);
1290: DMPlexGetCone(dm, oldp, &cone);
1291: DMPlexGetConeOrientation(dm, oldp, &ornt);
1292: DMPlexGetSupportSize(dm, oldp, &supportSize);
1293: DMPlexGetSupport(dm, oldp, &support);
1294: if (dep == 0) {
1295: const PetscInt hybedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1297: /* Unsplit vertex */
1298: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1299: for (s = 0, q = 0; s < supportSize; ++s) {
1300: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1301: PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e);
1302: if (e >= 0) {
1303: supportNew[q++] = e + pMaxNew[dep+1];
1304: }
1305: }
1306: supportNew[q++] = hybedge;
1307: supportNew[q++] = hybedge;
1308: if (q != supportSizeNew) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp);
1309: DMPlexSetSupport(sdm, newp, supportNew);
1310: /* Hybrid edge */
1311: coneNew[0] = newp;
1312: coneNew[1] = newp;
1313: DMPlexSetCone(sdm, hybedge, coneNew);
1314: for (e = 0, qf = 0; e < supportSize; ++e) {
1315: PetscInt val, edge;
1317: DMLabelGetValue(label, support[e], &val);
1318: if (val == 1) {
1319: PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge);
1320: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]);
1321: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2];
1322: } else if (val == (shift2 + 1)) {
1323: PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge);
1324: if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a unsplit edge", support[e]);
1325: supportNew[qf++] = edge + pMaxNew[dep+2] + numSplitPoints[dep+2] + numSplitPoints[dep+1];
1326: }
1327: }
1328: DMPlexSetSupport(sdm, hybedge, supportNew);
1329: } else if (dep == dim-2) {
1330: const PetscInt hybface = p + pMaxNew[dep+1] + numSplitPoints[dep+1] + numSplitPoints[dep];
1332: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1333: for (f = 0, qf = 0; f < supportSize; ++f) {
1334: PetscInt val, face;
1336: DMLabelGetValue(label, support[f], &val);
1337: if (val == dim-1) {
1338: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1339: if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[f]);
1340: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1341: supportNew[qf++] = face + pMaxNew[dep+1];
1342: } else {
1343: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1344: }
1345: }
1346: supportNew[qf++] = hybface;
1347: supportNew[qf++] = hybface;
1348: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1349: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d\n", newp, qf, supportSizeNew);
1350: DMPlexSetSupport(sdm, newp, supportNew);
1351: /* Add hybrid face */
1352: coneNew[0] = newp;
1353: coneNew[1] = newp;
1354: PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1355: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[0]);
1356: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1357: PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v);
1358: if (v < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex %d is not an unsplit vertex", cone[1]);
1359: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep-1];
1360: DMPlexSetCone(sdm, hybface, coneNew);
1361: for (f = 0, qf = 0; f < supportSize; ++f) {
1362: PetscInt val, face;
1364: DMLabelGetValue(label, support[f], &val);
1365: if (val == dim-1) {
1366: PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face);
1367: supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2];
1368: }
1369: }
1370: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1371: if (qf != supportSizeNew) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d\n", hybface, qf, supportSizeNew);
1372: DMPlexSetSupport(sdm, hybface, supportNew);
1373: }
1374: }
1375: }
1376: /* Step 6b: Replace split points in negative side cones */
1377: for (sp = 0; sp < numSP; ++sp) {
1378: PetscInt dep = values[sp];
1379: IS pIS;
1380: PetscInt numPoints;
1381: const PetscInt *points;
1383: if (dep >= 0) continue;
1384: DMLabelGetStratumIS(label, dep, &pIS);
1385: if (!pIS) continue;
1386: dep = -dep - shift;
1387: ISGetLocalSize(pIS, &numPoints);
1388: ISGetIndices(pIS, &points);
1389: for (p = 0; p < numPoints; ++p) {
1390: const PetscInt oldp = points[p];
1391: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1392: const PetscInt *cone;
1393: PetscInt coneSize, c;
1394: /* PetscBool replaced = PETSC_FALSE; */
1396: /* Negative edge: replace split vertex */
1397: /* Negative cell: replace split face */
1398: DMPlexGetConeSize(sdm, newp, &coneSize);
1399: DMPlexGetCone(sdm, newp, &cone);
1400: for (c = 0; c < coneSize; ++c) {
1401: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c],depth,depthShift);
1402: PetscInt csplitp, cp, val;
1404: DMLabelGetValue(label, coldp, &val);
1405: if (val == dep-1) {
1406: PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);
1407: if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1);
1408: csplitp = pMaxNew[dep-1] + cp;
1409: DMPlexInsertCone(sdm, newp, c, csplitp);
1410: /* replaced = PETSC_TRUE; */
1411: }
1412: }
1413: /* Cells with only a vertex or edge on the submesh have no replacement */
1414: /* if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
1415: }
1416: ISRestoreIndices(pIS, &points);
1417: ISDestroy(&pIS);
1418: }
1419: /* Step 7: Stratify */
1420: DMPlexStratify(sdm);
1421: /* Step 8: Coordinates */
1422: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1423: DMGetCoordinateSection(sdm, &coordSection);
1424: DMGetCoordinatesLocal(sdm, &coordinates);
1425: VecGetArray(coordinates, &coords);
1426: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1427: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1428: const PetscInt splitp = pMaxNew[0] + v;
1429: PetscInt dof, off, soff, d;
1431: PetscSectionGetDof(coordSection, newp, &dof);
1432: PetscSectionGetOffset(coordSection, newp, &off);
1433: PetscSectionGetOffset(coordSection, splitp, &soff);
1434: for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d];
1435: }
1436: VecRestoreArray(coordinates, &coords);
1437: /* Step 9: SF, if I can figure this out we can split the mesh in parallel */
1438: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1439: /* Step 10: Labels */
1440: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1441: DMGetNumLabels(sdm, &numLabels);
1442: for (dep = 0; dep <= depth; ++dep) {
1443: for (p = 0; p < numSplitPoints[dep]; ++p) {
1444: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1445: const PetscInt splitp = pMaxNew[dep] + p;
1446: PetscInt l;
1448: for (l = 0; l < numLabels; ++l) {
1449: DMLabel mlabel;
1450: const char *lname;
1451: PetscInt val;
1452: PetscBool isDepth;
1454: DMGetLabelName(sdm, l, &lname);
1455: PetscStrcmp(lname, "depth", &isDepth);
1456: if (isDepth) continue;
1457: DMGetLabel(sdm, lname, &mlabel);
1458: DMLabelGetValue(mlabel, newp, &val);
1459: if (val >= 0) {
1460: DMLabelSetValue(mlabel, splitp, val);
1461: }
1462: }
1463: }
1464: }
1465: for (sp = 0; sp < numSP; ++sp) {
1466: const PetscInt dep = values[sp];
1468: if ((dep < 0) || (dep > depth)) continue;
1469: if (splitIS[dep]) {ISRestoreIndices(splitIS[dep], &splitPoints[dep]);}
1470: ISDestroy(&splitIS[dep]);
1471: if (unsplitIS[dep]) {ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);}
1472: ISDestroy(&unsplitIS[dep]);
1473: }
1474: if (label) {
1475: ISRestoreIndices(valueIS, &values);
1476: ISDestroy(&valueIS);
1477: }
1478: for (d = 0; d <= depth; ++d) {
1479: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1480: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1481: }
1482: 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);
1483: PetscFree3(coneNew, coneONew, supportNew);
1484: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1485: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1486: return(0);
1487: }
1489: /*@C
1490: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1492: Collective on dm
1494: Input Parameters:
1495: + dm - The original DM
1496: - label - The label specifying the boundary faces (this could be auto-generated)
1498: Output Parameters:
1499: - dmSplit - The new DM
1501: Level: developer
1503: .seealso: DMCreate(), DMPlexLabelCohesiveComplete()
1504: @*/
1505: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)
1506: {
1507: DM sdm;
1508: PetscInt dim;
1514: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1515: DMSetType(sdm, DMPLEX);
1516: DMGetDimension(dm, &dim);
1517: DMSetDimension(sdm, dim);
1518: switch (dim) {
1519: case 2:
1520: case 3:
1521: DMPlexConstructCohesiveCells_Internal(dm, label, sdm);
1522: break;
1523: default:
1524: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim);
1525: }
1526: *dmSplit = sdm;
1527: return(0);
1528: }
1530: /* Returns the side of the surface for a given cell with a face on the surface */
1531: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1532: {
1533: const PetscInt *cone, *ornt;
1534: PetscInt dim, coneSize, c;
1535: PetscErrorCode ierr;
1538: *pos = PETSC_TRUE;
1539: DMGetDimension(dm, &dim);
1540: DMPlexGetConeSize(dm, cell, &coneSize);
1541: DMPlexGetCone(dm, cell, &cone);
1542: DMPlexGetConeOrientation(dm, cell, &ornt);
1543: for (c = 0; c < coneSize; ++c) {
1544: if (cone[c] == face) {
1545: PetscInt o = ornt[c];
1547: if (subdm) {
1548: const PetscInt *subcone, *subornt;
1549: PetscInt subpoint, subface, subconeSize, sc;
1551: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1552: PetscFindInt(face, numSubpoints, subpoints, &subface);
1553: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1554: DMPlexGetCone(subdm, subpoint, &subcone);
1555: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1556: for (sc = 0; sc < subconeSize; ++sc) {
1557: if (subcone[sc] == subface) {
1558: o = subornt[0];
1559: break;
1560: }
1561: }
1562: 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);
1563: }
1564: if (o >= 0) *pos = PETSC_TRUE;
1565: else *pos = PETSC_FALSE;
1566: break;
1567: }
1568: }
1569: 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);
1570: return(0);
1571: }
1573: /*@
1574: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1575: to complete the surface
1577: Input Parameters:
1578: + dm - The DM
1579: . label - A DMLabel marking the surface
1580: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1581: . flip - Flag to flip the submesh normal and replace points on the other side
1582: - subdm - The subDM associated with the label, or NULL
1584: Output Parameter:
1585: . label - A DMLabel marking all surface points
1587: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1589: Level: developer
1591: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete()
1592: @*/
1593: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm)
1594: {
1595: DMLabel depthLabel;
1596: IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL;
1597: const PetscInt *points, *subpoints;
1598: const PetscInt rev = flip ? -1 : 1;
1599: PetscInt *pMax;
1600: PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val;
1601: PetscErrorCode ierr;
1604: DMPlexGetDepth(dm, &depth);
1605: DMGetDimension(dm, &dim);
1606: pSize = PetscMax(depth, dim) + 1;
1607: PetscMalloc1(pSize,&pMax);
1608: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1609: DMPlexGetDepthLabel(dm, &depthLabel);
1610: DMGetDimension(dm, &dim);
1611: if (subdm) {
1612: DMPlexCreateSubpointIS(subdm, &subpointIS);
1613: if (subpointIS) {
1614: ISGetLocalSize(subpointIS, &numSubpoints);
1615: ISGetIndices(subpointIS, &subpoints);
1616: }
1617: }
1618: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1619: DMLabelGetStratumIS(label, dim-1, &dimIS);
1620: if (!dimIS) {
1621: PetscFree(pMax);
1622: ISDestroy(&subpointIS);
1623: return(0);
1624: }
1625: ISGetLocalSize(dimIS, &numPoints);
1626: ISGetIndices(dimIS, &points);
1627: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
1628: const PetscInt *support;
1629: PetscInt supportSize, s;
1631: DMPlexGetSupportSize(dm, points[p], &supportSize);
1632: if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize);
1633: DMPlexGetSupport(dm, points[p], &support);
1634: for (s = 0; s < supportSize; ++s) {
1635: const PetscInt *cone;
1636: PetscInt coneSize, c;
1637: PetscBool pos;
1639: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
1640: if (pos) {DMLabelSetValue(label, support[s], rev*(shift+dim));}
1641: else {DMLabelSetValue(label, support[s], -rev*(shift+dim));}
1642: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
1643: /* Put faces touching the fault in the label */
1644: DMPlexGetConeSize(dm, support[s], &coneSize);
1645: DMPlexGetCone(dm, support[s], &cone);
1646: for (c = 0; c < coneSize; ++c) {
1647: const PetscInt point = cone[c];
1649: DMLabelGetValue(label, point, &val);
1650: if (val == -1) {
1651: PetscInt *closure = NULL;
1652: PetscInt closureSize, cl;
1654: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1655: for (cl = 0; cl < closureSize*2; cl += 2) {
1656: const PetscInt clp = closure[cl];
1657: PetscInt bval = -1;
1659: DMLabelGetValue(label, clp, &val);
1660: if (blabel) {DMLabelGetValue(blabel, clp, &bval);}
1661: if ((val >= 0) && (val < dim-1) && (bval < 0)) {
1662: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));
1663: break;
1664: }
1665: }
1666: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1667: }
1668: }
1669: }
1670: }
1671: ISRestoreIndices(dimIS, &points);
1672: ISDestroy(&dimIS);
1673: if (subpointIS) {ISRestoreIndices(subpointIS, &subpoints);}
1674: ISDestroy(&subpointIS);
1675: /* Mark boundary points as unsplit */
1676: if (blabel) {
1677: DMLabelGetStratumIS(blabel, 1, &dimIS);
1678: ISGetLocalSize(dimIS, &numPoints);
1679: ISGetIndices(dimIS, &points);
1680: for (p = 0; p < numPoints; ++p) {
1681: const PetscInt point = points[p];
1682: PetscInt val, bval;
1684: DMLabelGetValue(blabel, point, &bval);
1685: if (bval >= 0) {
1686: DMLabelGetValue(label, point, &val);
1687: if ((val < 0) || (val > dim)) {
1688: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
1689: DMLabelClearValue(blabel, point, bval);
1690: }
1691: }
1692: }
1693: for (p = 0; p < numPoints; ++p) {
1694: const PetscInt point = points[p];
1695: PetscInt val, bval;
1697: DMLabelGetValue(blabel, point, &bval);
1698: if (bval >= 0) {
1699: const PetscInt *cone, *support;
1700: PetscInt coneSize, supportSize, s, valA, valB, valE;
1702: /* Mark as unsplit */
1703: DMLabelGetValue(label, point, &val);
1704: 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);
1705: DMLabelClearValue(label, point, val);
1706: DMLabelSetValue(label, point, shift2+val);
1707: /* Check for cross-edge
1708: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
1709: if (val != 0) continue;
1710: DMPlexGetSupport(dm, point, &support);
1711: DMPlexGetSupportSize(dm, point, &supportSize);
1712: for (s = 0; s < supportSize; ++s) {
1713: DMPlexGetCone(dm, support[s], &cone);
1714: DMPlexGetConeSize(dm, support[s], &coneSize);
1715: if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize);
1716: DMLabelGetValue(blabel, cone[0], &valA);
1717: DMLabelGetValue(blabel, cone[1], &valB);
1718: DMLabelGetValue(blabel, support[s], &valE);
1719: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {DMLabelSetValue(blabel, support[s], 2);}
1720: }
1721: }
1722: }
1723: ISRestoreIndices(dimIS, &points);
1724: ISDestroy(&dimIS);
1725: }
1726: /* Search for other cells/faces/edges connected to the fault by a vertex */
1727: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1728: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1729: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1730: cMax = cMax < 0 ? cEnd : cMax;
1731: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1732: DMLabelGetStratumIS(label, 0, &dimIS);
1733: if (blabel) {DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);}
1734: if (dimIS && crossEdgeIS) {
1735: IS vertIS = dimIS;
1737: ISExpand(vertIS, crossEdgeIS, &dimIS);
1738: ISDestroy(&crossEdgeIS);
1739: ISDestroy(&vertIS);
1740: }
1741: if (!dimIS) {
1742: PetscFree(pMax);
1743: return(0);
1744: }
1745: ISGetLocalSize(dimIS, &numPoints);
1746: ISGetIndices(dimIS, &points);
1747: for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */
1748: PetscInt *star = NULL;
1749: PetscInt starSize, s;
1750: PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */
1752: /* All points connected to the fault are inside a cell, so at the top level we will only check cells */
1753: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1754: while (again) {
1755: if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault");
1756: again = 0;
1757: for (s = 0; s < starSize*2; s += 2) {
1758: const PetscInt point = star[s];
1759: const PetscInt *cone;
1760: PetscInt coneSize, c;
1762: if ((point < cStart) || (point >= cMax)) continue;
1763: DMLabelGetValue(label, point, &val);
1764: if (val != -1) continue;
1765: again = again == 1 ? 1 : 2;
1766: DMPlexGetConeSize(dm, point, &coneSize);
1767: DMPlexGetCone(dm, point, &cone);
1768: for (c = 0; c < coneSize; ++c) {
1769: DMLabelGetValue(label, cone[c], &val);
1770: if (val != -1) {
1771: const PetscInt *ccone;
1772: PetscInt cconeSize, cc, side;
1774: 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);
1775: if (val > 0) side = 1;
1776: else side = -1;
1777: DMLabelSetValue(label, point, side*(shift+dim));
1778: /* Mark cell faces which touch the fault */
1779: DMPlexGetConeSize(dm, point, &cconeSize);
1780: DMPlexGetCone(dm, point, &ccone);
1781: for (cc = 0; cc < cconeSize; ++cc) {
1782: PetscInt *closure = NULL;
1783: PetscInt closureSize, cl;
1785: DMLabelGetValue(label, ccone[cc], &val);
1786: if (val != -1) continue;
1787: DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1788: for (cl = 0; cl < closureSize*2; cl += 2) {
1789: const PetscInt clp = closure[cl];
1791: DMLabelGetValue(label, clp, &val);
1792: if (val == -1) continue;
1793: DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));
1794: break;
1795: }
1796: DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);
1797: }
1798: again = 1;
1799: break;
1800: }
1801: }
1802: }
1803: }
1804: /* Classify the rest by cell membership */
1805: for (s = 0; s < starSize*2; s += 2) {
1806: const PetscInt point = star[s];
1808: DMLabelGetValue(label, point, &val);
1809: if (val == -1) {
1810: PetscInt *sstar = NULL;
1811: PetscInt sstarSize, ss;
1812: PetscBool marked = PETSC_FALSE;
1814: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1815: for (ss = 0; ss < sstarSize*2; ss += 2) {
1816: const PetscInt spoint = sstar[ss];
1818: if ((spoint < cStart) || (spoint >= cMax)) continue;
1819: DMLabelGetValue(label, spoint, &val);
1820: if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point);
1821: DMLabelGetValue(depthLabel, point, &dep);
1822: if (val > 0) {
1823: DMLabelSetValue(label, point, shift+dep);
1824: } else {
1825: DMLabelSetValue(label, point, -(shift+dep));
1826: }
1827: marked = PETSC_TRUE;
1828: break;
1829: }
1830: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);
1831: DMLabelGetValue(depthLabel, point, &dep);
1832: if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point);
1833: }
1834: }
1835: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);
1836: }
1837: ISRestoreIndices(dimIS, &points);
1838: ISDestroy(&dimIS);
1839: /* If any faces touching the fault divide cells on either side, split them */
1840: DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);
1841: DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);
1842: ISExpand(facePosIS, faceNegIS, &dimIS);
1843: ISDestroy(&facePosIS);
1844: ISDestroy(&faceNegIS);
1845: ISGetLocalSize(dimIS, &numPoints);
1846: ISGetIndices(dimIS, &points);
1847: for (p = 0; p < numPoints; ++p) {
1848: const PetscInt point = points[p];
1849: const PetscInt *support;
1850: PetscInt supportSize, valA, valB;
1852: DMPlexGetSupportSize(dm, point, &supportSize);
1853: if (supportSize != 2) continue;
1854: DMPlexGetSupport(dm, point, &support);
1855: DMLabelGetValue(label, support[0], &valA);
1856: DMLabelGetValue(label, support[1], &valB);
1857: if ((valA == -1) || (valB == -1)) continue;
1858: if (valA*valB > 0) continue;
1859: /* Split the face */
1860: DMLabelGetValue(label, point, &valA);
1861: DMLabelClearValue(label, point, valA);
1862: DMLabelSetValue(label, point, dim-1);
1863: /* Label its closure:
1864: unmarked: label as unsplit
1865: incident: relabel as split
1866: split: do nothing
1867: */
1868: {
1869: PetscInt *closure = NULL;
1870: PetscInt closureSize, cl;
1872: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1873: for (cl = 0; cl < closureSize*2; cl += 2) {
1874: DMLabelGetValue(label, closure[cl], &valA);
1875: if (valA == -1) { /* Mark as unsplit */
1876: DMLabelGetValue(depthLabel, closure[cl], &dep);
1877: DMLabelSetValue(label, closure[cl], shift2+dep);
1878: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1879: DMLabelGetValue(depthLabel, closure[cl], &dep);
1880: DMLabelClearValue(label, closure[cl], valA);
1881: DMLabelSetValue(label, closure[cl], dep);
1882: }
1883: }
1884: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1885: }
1886: }
1887: ISRestoreIndices(dimIS, &points);
1888: ISDestroy(&dimIS);
1889: PetscFree(pMax);
1890: return(0);
1891: }
1893: /*@
1894: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
1896: Collective on dm
1898: Input Parameters:
1899: + dm - The original DM
1900: - labelName - The label specifying the interface vertices
1902: Output Parameters:
1903: + hybridLabel - The label fully marking the interface
1904: - dmHybrid - The new DM
1906: Level: developer
1908: .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate()
1909: @*/
1910: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid)
1911: {
1912: DM idm;
1913: DMLabel subpointMap, hlabel;
1914: PetscInt dim;
1921: DMGetDimension(dm, &dim);
1922: DMPlexCreateSubmesh(dm, label, 1, &idm);
1923: DMPlexOrient(idm);
1924: DMPlexGetSubpointMap(idm, &subpointMap);
1925: DMLabelDuplicate(subpointMap, &hlabel);
1926: DMLabelClearStratum(hlabel, dim);
1927: DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);
1928: DMDestroy(&idm);
1929: DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);
1930: if (hybridLabel) *hybridLabel = hlabel;
1931: else {DMLabelDestroy(&hlabel);}
1932: return(0);
1933: }
1935: /* Here we need the explicit assumption that:
1937: For any marked cell, the marked vertices constitute a single face
1938: */
1939: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
1940: {
1941: IS subvertexIS = NULL;
1942: const PetscInt *subvertices;
1943: PetscInt *pStart, *pEnd, *pMax, pSize;
1944: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
1945: PetscErrorCode ierr;
1948: *numFaces = 0;
1949: *nFV = 0;
1950: DMPlexGetDepth(dm, &depth);
1951: DMGetDimension(dm, &dim);
1952: pSize = PetscMax(depth, dim) + 1;
1953: PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);
1954: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
1955: for (d = 0; d <= depth; ++d) {
1956: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
1957: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
1958: }
1959: /* Loop over initial vertices and mark all faces in the collective star() */
1960: if (vertexLabel) {DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);}
1961: if (subvertexIS) {
1962: ISGetSize(subvertexIS, &numSubVerticesInitial);
1963: ISGetIndices(subvertexIS, &subvertices);
1964: }
1965: for (v = 0; v < numSubVerticesInitial; ++v) {
1966: const PetscInt vertex = subvertices[v];
1967: PetscInt *star = NULL;
1968: PetscInt starSize, s, numCells = 0, c;
1970: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
1971: for (s = 0; s < starSize*2; s += 2) {
1972: const PetscInt point = star[s];
1973: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
1974: }
1975: for (c = 0; c < numCells; ++c) {
1976: const PetscInt cell = star[c];
1977: PetscInt *closure = NULL;
1978: PetscInt closureSize, cl;
1979: PetscInt cellLoc, numCorners = 0, faceSize = 0;
1981: DMLabelGetValue(subpointMap, cell, &cellLoc);
1982: if (cellLoc == 2) continue;
1983: if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc);
1984: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1985: for (cl = 0; cl < closureSize*2; cl += 2) {
1986: const PetscInt point = closure[cl];
1987: PetscInt vertexLoc;
1989: if ((point >= pStart[0]) && (point < pEnd[0])) {
1990: ++numCorners;
1991: DMLabelGetValue(vertexLabel, point, &vertexLoc);
1992: if (vertexLoc == value) closure[faceSize++] = point;
1993: }
1994: }
1995: if (!(*nFV)) {DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);}
1996: if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
1997: if (faceSize == *nFV) {
1998: const PetscInt *cells = NULL;
1999: PetscInt numCells, nc;
2001: ++(*numFaces);
2002: for (cl = 0; cl < faceSize; ++cl) {
2003: DMLabelSetValue(subpointMap, closure[cl], 0);
2004: }
2005: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2006: for (nc = 0; nc < numCells; ++nc) {
2007: DMLabelSetValue(subpointMap, cells[nc], 2);
2008: }
2009: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2010: }
2011: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2012: }
2013: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2014: }
2015: if (subvertexIS) {
2016: ISRestoreIndices(subvertexIS, &subvertices);
2017: }
2018: ISDestroy(&subvertexIS);
2019: PetscFree3(pStart,pEnd,pMax);
2020: return(0);
2021: }
2023: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm)
2024: {
2025: IS subvertexIS = NULL;
2026: const PetscInt *subvertices;
2027: PetscInt *pStart, *pEnd, *pMax;
2028: PetscInt dim, d, numSubVerticesInitial = 0, v;
2029: PetscErrorCode ierr;
2032: DMGetDimension(dm, &dim);
2033: PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);
2034: DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);
2035: for (d = 0; d <= dim; ++d) {
2036: DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2037: if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]);
2038: }
2039: /* Loop over initial vertices and mark all faces in the collective star() */
2040: if (vertexLabel) {
2041: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2042: if (subvertexIS) {
2043: ISGetSize(subvertexIS, &numSubVerticesInitial);
2044: ISGetIndices(subvertexIS, &subvertices);
2045: }
2046: }
2047: for (v = 0; v < numSubVerticesInitial; ++v) {
2048: const PetscInt vertex = subvertices[v];
2049: PetscInt *star = NULL;
2050: PetscInt starSize, s, numFaces = 0, f;
2052: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2053: for (s = 0; s < starSize*2; s += 2) {
2054: const PetscInt point = star[s];
2055: if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point;
2056: }
2057: for (f = 0; f < numFaces; ++f) {
2058: const PetscInt face = star[f];
2059: PetscInt *closure = NULL;
2060: PetscInt closureSize, c;
2061: PetscInt faceLoc;
2063: DMLabelGetValue(subpointMap, face, &faceLoc);
2064: if (faceLoc == dim-1) continue;
2065: if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc);
2066: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2067: for (c = 0; c < closureSize*2; c += 2) {
2068: const PetscInt point = closure[c];
2069: PetscInt vertexLoc;
2071: if ((point >= pStart[0]) && (point < pEnd[0])) {
2072: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2073: if (vertexLoc != value) break;
2074: }
2075: }
2076: if (c == closureSize*2) {
2077: const PetscInt *support;
2078: PetscInt supportSize, s;
2080: for (c = 0; c < closureSize*2; c += 2) {
2081: const PetscInt point = closure[c];
2083: for (d = 0; d < dim; ++d) {
2084: if ((point >= pStart[d]) && (point < pEnd[d])) {
2085: DMLabelSetValue(subpointMap, point, d);
2086: break;
2087: }
2088: }
2089: }
2090: DMPlexGetSupportSize(dm, face, &supportSize);
2091: DMPlexGetSupport(dm, face, &support);
2092: for (s = 0; s < supportSize; ++s) {
2093: DMLabelSetValue(subpointMap, support[s], dim);
2094: }
2095: }
2096: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2097: }
2098: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2099: }
2100: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subvertices);}
2101: ISDestroy(&subvertexIS);
2102: PetscFree3(pStart,pEnd,pMax);
2103: return(0);
2104: }
2106: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2107: {
2108: DMLabel label = NULL;
2109: const PetscInt *cone;
2110: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2111: PetscErrorCode ierr;
2114: *numFaces = 0;
2115: *nFV = 0;
2116: if (labelname) {DMGetLabel(dm, labelname, &label);}
2117: *subCells = NULL;
2118: DMGetDimension(dm, &dim);
2119: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2120: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2121: if (cMax < 0) return(0);
2122: if (label) {
2123: for (c = cMax; c < cEnd; ++c) {
2124: PetscInt val;
2126: DMLabelGetValue(label, c, &val);
2127: if (val == value) {
2128: ++(*numFaces);
2129: DMPlexGetConeSize(dm, c, &coneSize);
2130: }
2131: }
2132: } else {
2133: *numFaces = cEnd - cMax;
2134: DMPlexGetConeSize(dm, cMax, &coneSize);
2135: }
2136: PetscMalloc1(*numFaces *2, subCells);
2137: if (!(*numFaces)) return(0);
2138: *nFV = hasLagrange ? coneSize/3 : coneSize/2;
2139: for (c = cMax; c < cEnd; ++c) {
2140: const PetscInt *cells;
2141: PetscInt numCells;
2143: if (label) {
2144: PetscInt val;
2146: DMLabelGetValue(label, c, &val);
2147: if (val != value) continue;
2148: }
2149: DMPlexGetCone(dm, c, &cone);
2150: for (p = 0; p < *nFV; ++p) {
2151: DMLabelSetValue(subpointMap, cone[p], 0);
2152: }
2153: /* Negative face */
2154: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2155: /* Not true in parallel
2156: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
2157: for (p = 0; p < numCells; ++p) {
2158: DMLabelSetValue(subpointMap, cells[p], 2);
2159: (*subCells)[subc++] = cells[p];
2160: }
2161: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2162: /* Positive face is not included */
2163: }
2164: return(0);
2165: }
2167: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2168: {
2169: PetscInt *pStart, *pEnd;
2170: PetscInt dim, cMax, cEnd, c, d;
2174: DMGetDimension(dm, &dim);
2175: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2176: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2177: if (cMax < 0) return(0);
2178: PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);
2179: for (d = 0; d <= dim; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
2180: for (c = cMax; c < cEnd; ++c) {
2181: const PetscInt *cone;
2182: PetscInt *closure = NULL;
2183: PetscInt fconeSize, coneSize, closureSize, cl, val;
2185: if (label) {
2186: DMLabelGetValue(label, c, &val);
2187: if (val != value) continue;
2188: }
2189: DMPlexGetConeSize(dm, c, &coneSize);
2190: DMPlexGetCone(dm, c, &cone);
2191: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2192: if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
2193: /* Negative face */
2194: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2195: for (cl = 0; cl < closureSize*2; cl += 2) {
2196: const PetscInt point = closure[cl];
2198: for (d = 0; d <= dim; ++d) {
2199: if ((point >= pStart[d]) && (point < pEnd[d])) {
2200: DMLabelSetValue(subpointMap, point, d);
2201: break;
2202: }
2203: }
2204: }
2205: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2206: /* Cells -- positive face is not included */
2207: for (cl = 0; cl < 1; ++cl) {
2208: const PetscInt *support;
2209: PetscInt supportSize, s;
2211: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2212: /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
2213: DMPlexGetSupport(dm, cone[cl], &support);
2214: for (s = 0; s < supportSize; ++s) {
2215: DMLabelSetValue(subpointMap, support[s], dim);
2216: }
2217: }
2218: }
2219: PetscFree2(pStart, pEnd);
2220: return(0);
2221: }
2223: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2224: {
2225: MPI_Comm comm;
2226: PetscBool posOrient = PETSC_FALSE;
2227: const PetscInt debug = 0;
2228: PetscInt cellDim, faceSize, f;
2232: PetscObjectGetComm((PetscObject)dm,&comm);
2233: DMGetDimension(dm, &cellDim);
2234: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
2236: if (cellDim == 1 && numCorners == 2) {
2237: /* Triangle */
2238: faceSize = numCorners-1;
2239: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2240: } else if (cellDim == 2 && numCorners == 3) {
2241: /* Triangle */
2242: faceSize = numCorners-1;
2243: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2244: } else if (cellDim == 3 && numCorners == 4) {
2245: /* Tetrahedron */
2246: faceSize = numCorners-1;
2247: posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
2248: } else if (cellDim == 1 && numCorners == 3) {
2249: /* Quadratic line */
2250: faceSize = 1;
2251: posOrient = PETSC_TRUE;
2252: } else if (cellDim == 2 && numCorners == 4) {
2253: /* Quads */
2254: faceSize = 2;
2255: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2256: posOrient = PETSC_TRUE;
2257: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2258: posOrient = PETSC_TRUE;
2259: } else {
2260: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2261: posOrient = PETSC_FALSE;
2262: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2263: }
2264: } else if (cellDim == 2 && numCorners == 6) {
2265: /* Quadratic triangle (I hate this) */
2266: /* Edges are determined by the first 2 vertices (corners of edges) */
2267: const PetscInt faceSizeTri = 3;
2268: PetscInt sortedIndices[3], i, iFace;
2269: PetscBool found = PETSC_FALSE;
2270: PetscInt faceVerticesTriSorted[9] = {
2271: 0, 3, 4, /* bottom */
2272: 1, 4, 5, /* right */
2273: 2, 3, 5, /* left */
2274: };
2275: PetscInt faceVerticesTri[9] = {
2276: 0, 3, 4, /* bottom */
2277: 1, 4, 5, /* right */
2278: 2, 5, 3, /* left */
2279: };
2281: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2282: PetscSortInt(faceSizeTri, sortedIndices);
2283: for (iFace = 0; iFace < 3; ++iFace) {
2284: const PetscInt ii = iFace*faceSizeTri;
2285: PetscInt fVertex, cVertex;
2287: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
2288: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
2289: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2290: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2291: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
2292: faceVertices[fVertex] = origVertices[cVertex];
2293: break;
2294: }
2295: }
2296: }
2297: found = PETSC_TRUE;
2298: break;
2299: }
2300: }
2301: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
2302: if (posOriented) *posOriented = PETSC_TRUE;
2303: return(0);
2304: } else if (cellDim == 2 && numCorners == 9) {
2305: /* Quadratic quad (I hate this) */
2306: /* Edges are determined by the first 2 vertices (corners of edges) */
2307: const PetscInt faceSizeQuad = 3;
2308: PetscInt sortedIndices[3], i, iFace;
2309: PetscBool found = PETSC_FALSE;
2310: PetscInt faceVerticesQuadSorted[12] = {
2311: 0, 1, 4, /* bottom */
2312: 1, 2, 5, /* right */
2313: 2, 3, 6, /* top */
2314: 0, 3, 7, /* left */
2315: };
2316: PetscInt faceVerticesQuad[12] = {
2317: 0, 1, 4, /* bottom */
2318: 1, 2, 5, /* right */
2319: 2, 3, 6, /* top */
2320: 3, 0, 7, /* left */
2321: };
2323: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2324: PetscSortInt(faceSizeQuad, sortedIndices);
2325: for (iFace = 0; iFace < 4; ++iFace) {
2326: const PetscInt ii = iFace*faceSizeQuad;
2327: PetscInt fVertex, cVertex;
2329: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
2330: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
2331: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2332: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2333: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
2334: faceVertices[fVertex] = origVertices[cVertex];
2335: break;
2336: }
2337: }
2338: }
2339: found = PETSC_TRUE;
2340: break;
2341: }
2342: }
2343: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
2344: if (posOriented) *posOriented = PETSC_TRUE;
2345: return(0);
2346: } else if (cellDim == 3 && numCorners == 8) {
2347: /* Hexes
2348: A hex is two oriented quads with the normal of the first
2349: pointing up at the second.
2351: 7---6
2352: /| /|
2353: 4---5 |
2354: | 1-|-2
2355: |/ |/
2356: 0---3
2358: Faces are determined by the first 4 vertices (corners of faces) */
2359: const PetscInt faceSizeHex = 4;
2360: PetscInt sortedIndices[4], i, iFace;
2361: PetscBool found = PETSC_FALSE;
2362: PetscInt faceVerticesHexSorted[24] = {
2363: 0, 1, 2, 3, /* bottom */
2364: 4, 5, 6, 7, /* top */
2365: 0, 3, 4, 5, /* front */
2366: 2, 3, 5, 6, /* right */
2367: 1, 2, 6, 7, /* back */
2368: 0, 1, 4, 7, /* left */
2369: };
2370: PetscInt faceVerticesHex[24] = {
2371: 1, 2, 3, 0, /* bottom */
2372: 4, 5, 6, 7, /* top */
2373: 0, 3, 5, 4, /* front */
2374: 3, 2, 6, 5, /* right */
2375: 2, 1, 7, 6, /* back */
2376: 1, 0, 4, 7, /* left */
2377: };
2379: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2380: PetscSortInt(faceSizeHex, sortedIndices);
2381: for (iFace = 0; iFace < 6; ++iFace) {
2382: const PetscInt ii = iFace*faceSizeHex;
2383: PetscInt fVertex, cVertex;
2385: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
2386: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
2387: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
2388: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
2389: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2390: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2391: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
2392: faceVertices[fVertex] = origVertices[cVertex];
2393: break;
2394: }
2395: }
2396: }
2397: found = PETSC_TRUE;
2398: break;
2399: }
2400: }
2401: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2402: if (posOriented) *posOriented = PETSC_TRUE;
2403: return(0);
2404: } else if (cellDim == 3 && numCorners == 10) {
2405: /* Quadratic tet */
2406: /* Faces are determined by the first 3 vertices (corners of faces) */
2407: const PetscInt faceSizeTet = 6;
2408: PetscInt sortedIndices[6], i, iFace;
2409: PetscBool found = PETSC_FALSE;
2410: PetscInt faceVerticesTetSorted[24] = {
2411: 0, 1, 2, 6, 7, 8, /* bottom */
2412: 0, 3, 4, 6, 7, 9, /* front */
2413: 1, 4, 5, 7, 8, 9, /* right */
2414: 2, 3, 5, 6, 8, 9, /* left */
2415: };
2416: PetscInt faceVerticesTet[24] = {
2417: 0, 1, 2, 6, 7, 8, /* bottom */
2418: 0, 4, 3, 6, 7, 9, /* front */
2419: 1, 5, 4, 7, 8, 9, /* right */
2420: 2, 3, 5, 8, 6, 9, /* left */
2421: };
2423: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2424: PetscSortInt(faceSizeTet, sortedIndices);
2425: for (iFace=0; iFace < 4; ++iFace) {
2426: const PetscInt ii = iFace*faceSizeTet;
2427: PetscInt fVertex, cVertex;
2429: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
2430: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
2431: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
2432: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
2433: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2434: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2435: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
2436: faceVertices[fVertex] = origVertices[cVertex];
2437: break;
2438: }
2439: }
2440: }
2441: found = PETSC_TRUE;
2442: break;
2443: }
2444: }
2445: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
2446: if (posOriented) *posOriented = PETSC_TRUE;
2447: return(0);
2448: } else if (cellDim == 3 && numCorners == 27) {
2449: /* Quadratic hexes (I hate this)
2450: A hex is two oriented quads with the normal of the first
2451: pointing up at the second.
2453: 7---6
2454: /| /|
2455: 4---5 |
2456: | 3-|-2
2457: |/ |/
2458: 0---1
2460: Faces are determined by the first 4 vertices (corners of faces) */
2461: const PetscInt faceSizeQuadHex = 9;
2462: PetscInt sortedIndices[9], i, iFace;
2463: PetscBool found = PETSC_FALSE;
2464: PetscInt faceVerticesQuadHexSorted[54] = {
2465: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2466: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2467: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2468: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2469: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2470: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2471: };
2472: PetscInt faceVerticesQuadHex[54] = {
2473: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2474: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2475: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2476: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2477: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2478: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2479: };
2481: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2482: PetscSortInt(faceSizeQuadHex, sortedIndices);
2483: for (iFace = 0; iFace < 6; ++iFace) {
2484: const PetscInt ii = iFace*faceSizeQuadHex;
2485: PetscInt fVertex, cVertex;
2487: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
2488: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
2489: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
2490: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
2491: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2492: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2493: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
2494: faceVertices[fVertex] = origVertices[cVertex];
2495: break;
2496: }
2497: }
2498: }
2499: found = PETSC_TRUE;
2500: break;
2501: }
2502: }
2503: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
2504: if (posOriented) *posOriented = PETSC_TRUE;
2505: return(0);
2506: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2507: if (!posOrient) {
2508: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
2509: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f];
2510: } else {
2511: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
2512: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2513: }
2514: if (posOriented) *posOriented = posOrient;
2515: return(0);
2516: }
2518: /*@
2519: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2520: in faceVertices. The orientation is such that the face normal points out of the cell
2522: Not collective
2524: Input Parameters:
2525: + dm - The original mesh
2526: . cell - The cell mesh point
2527: . faceSize - The number of vertices on the face
2528: . face - The face vertices
2529: . numCorners - The number of vertices on the cell
2530: . indices - Local numbering of face vertices in cell cone
2531: - origVertices - Original face vertices
2533: Output Parameter:
2534: + faceVertices - The face vertices properly oriented
2535: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2537: Level: developer
2539: .seealso: DMPlexGetCone()
2540: @*/
2541: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2542: {
2543: const PetscInt *cone = NULL;
2544: PetscInt coneSize, v, f, v2;
2545: PetscInt oppositeVertex = -1;
2546: PetscErrorCode ierr;
2549: DMPlexGetConeSize(dm, cell, &coneSize);
2550: DMPlexGetCone(dm, cell, &cone);
2551: for (v = 0, v2 = 0; v < coneSize; ++v) {
2552: PetscBool found = PETSC_FALSE;
2554: for (f = 0; f < faceSize; ++f) {
2555: if (face[f] == cone[v]) {
2556: found = PETSC_TRUE; break;
2557: }
2558: }
2559: if (found) {
2560: indices[v2] = v;
2561: origVertices[v2] = cone[v];
2562: ++v2;
2563: } else {
2564: oppositeVertex = v;
2565: }
2566: }
2567: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2568: return(0);
2569: }
2571: /*
2572: DMPlexInsertFace_Internal - Puts a face into the mesh
2574: Not collective
2576: Input Parameters:
2577: + dm - The DMPlex
2578: . numFaceVertex - The number of vertices in the face
2579: . faceVertices - The vertices in the face for dm
2580: . subfaceVertices - The vertices in the face for subdm
2581: . numCorners - The number of vertices in the cell
2582: . cell - A cell in dm containing the face
2583: . subcell - A cell in subdm containing the face
2584: . firstFace - First face in the mesh
2585: - newFacePoint - Next face in the mesh
2587: Output Parameters:
2588: . newFacePoint - Contains next face point number on input, updated on output
2590: Level: developer
2591: */
2592: 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)
2593: {
2594: MPI_Comm comm;
2595: DM_Plex *submesh = (DM_Plex*) subdm->data;
2596: const PetscInt *faces;
2597: PetscInt numFaces, coneSize;
2598: PetscErrorCode ierr;
2601: PetscObjectGetComm((PetscObject)dm,&comm);
2602: DMPlexGetConeSize(subdm, subcell, &coneSize);
2603: if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize);
2604: #if 0
2605: /* Cannot use this because support() has not been constructed yet */
2606: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2607: #else
2608: {
2609: PetscInt f;
2611: numFaces = 0;
2612: DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2613: for (f = firstFace; f < *newFacePoint; ++f) {
2614: PetscInt dof, off, d;
2616: PetscSectionGetDof(submesh->coneSection, f, &dof);
2617: PetscSectionGetOffset(submesh->coneSection, f, &off);
2618: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2619: for (d = 0; d < dof; ++d) {
2620: const PetscInt p = submesh->cones[off+d];
2621: PetscInt v;
2623: for (v = 0; v < numFaceVertices; ++v) {
2624: if (subfaceVertices[v] == p) break;
2625: }
2626: if (v == numFaceVertices) break;
2627: }
2628: if (d == dof) {
2629: numFaces = 1;
2630: ((PetscInt*) faces)[0] = f;
2631: }
2632: }
2633: }
2634: #endif
2635: if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces);
2636: else if (numFaces == 1) {
2637: /* Add the other cell neighbor for this face */
2638: DMPlexSetCone(subdm, subcell, faces);
2639: } else {
2640: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2641: PetscBool posOriented;
2643: DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2644: origVertices = &orientedVertices[numFaceVertices];
2645: indices = &orientedVertices[numFaceVertices*2];
2646: orientedSubVertices = &orientedVertices[numFaceVertices*3];
2647: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2648: /* TODO: I know that routine should return a permutation, not the indices */
2649: for (v = 0; v < numFaceVertices; ++v) {
2650: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2651: for (ov = 0; ov < numFaceVertices; ++ov) {
2652: if (orientedVertices[ov] == vertex) {
2653: orientedSubVertices[ov] = subvertex;
2654: break;
2655: }
2656: }
2657: if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex);
2658: }
2659: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2660: DMPlexSetCone(subdm, subcell, newFacePoint);
2661: DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2662: ++(*newFacePoint);
2663: }
2664: #if 0
2665: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2666: #else
2667: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);
2668: #endif
2669: return(0);
2670: }
2672: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2673: {
2674: MPI_Comm comm;
2675: DMLabel subpointMap;
2676: IS subvertexIS, subcellIS;
2677: const PetscInt *subVertices, *subCells;
2678: PetscInt numSubVertices, firstSubVertex, numSubCells;
2679: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2680: PetscInt vStart, vEnd, c, f;
2681: PetscErrorCode ierr;
2684: PetscObjectGetComm((PetscObject)dm,&comm);
2685: /* Create subpointMap which marks the submesh */
2686: DMLabelCreate("subpoint_map", &subpointMap);
2687: DMPlexSetSubpointMap(subdm, subpointMap);
2688: DMLabelDestroy(&subpointMap);
2689: if (vertexLabel) {DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);}
2690: /* Setup chart */
2691: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
2692: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
2693: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
2694: DMPlexSetVTKCellHeight(subdm, 1);
2695: /* Set cone sizes */
2696: firstSubVertex = numSubCells;
2697: firstSubFace = numSubCells+numSubVertices;
2698: newFacePoint = firstSubFace;
2699: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
2700: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
2701: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
2702: if (subcellIS) {ISGetIndices(subcellIS, &subCells);}
2703: for (c = 0; c < numSubCells; ++c) {
2704: DMPlexSetConeSize(subdm, c, 1);
2705: }
2706: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
2707: DMPlexSetConeSize(subdm, f, nFV);
2708: }
2709: DMSetUp(subdm);
2710: /* Create face cones */
2711: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2712: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2713: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2714: for (c = 0; c < numSubCells; ++c) {
2715: const PetscInt cell = subCells[c];
2716: const PetscInt subcell = c;
2717: PetscInt *closure = NULL;
2718: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
2720: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2721: for (cl = 0; cl < closureSize*2; cl += 2) {
2722: const PetscInt point = closure[cl];
2723: PetscInt subVertex;
2725: if ((point >= vStart) && (point < vEnd)) {
2726: ++numCorners;
2727: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
2728: if (subVertex >= 0) {
2729: closure[faceSize] = point;
2730: subface[faceSize] = firstSubVertex+subVertex;
2731: ++faceSize;
2732: }
2733: }
2734: }
2735: if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize);
2736: if (faceSize == nFV) {
2737: DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
2738: }
2739: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2740: }
2741: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
2742: DMPlexSymmetrize(subdm);
2743: DMPlexStratify(subdm);
2744: /* Build coordinates */
2745: {
2746: PetscSection coordSection, subCoordSection;
2747: Vec coordinates, subCoordinates;
2748: PetscScalar *coords, *subCoords;
2749: PetscInt numComp, coordSize, v;
2750: const char *name;
2752: DMGetCoordinateSection(dm, &coordSection);
2753: DMGetCoordinatesLocal(dm, &coordinates);
2754: DMGetCoordinateSection(subdm, &subCoordSection);
2755: PetscSectionSetNumFields(subCoordSection, 1);
2756: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2757: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2758: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
2759: for (v = 0; v < numSubVertices; ++v) {
2760: const PetscInt vertex = subVertices[v];
2761: const PetscInt subvertex = firstSubVertex+v;
2762: PetscInt dof;
2764: PetscSectionGetDof(coordSection, vertex, &dof);
2765: PetscSectionSetDof(subCoordSection, subvertex, dof);
2766: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2767: }
2768: PetscSectionSetUp(subCoordSection);
2769: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2770: VecCreate(PETSC_COMM_SELF, &subCoordinates);
2771: PetscObjectGetName((PetscObject)coordinates,&name);
2772: PetscObjectSetName((PetscObject)subCoordinates,name);
2773: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
2774: VecSetType(subCoordinates,VECSTANDARD);
2775: if (coordSize) {
2776: VecGetArray(coordinates, &coords);
2777: VecGetArray(subCoordinates, &subCoords);
2778: for (v = 0; v < numSubVertices; ++v) {
2779: const PetscInt vertex = subVertices[v];
2780: const PetscInt subvertex = firstSubVertex+v;
2781: PetscInt dof, off, sdof, soff, d;
2783: PetscSectionGetDof(coordSection, vertex, &dof);
2784: PetscSectionGetOffset(coordSection, vertex, &off);
2785: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
2786: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
2787: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
2788: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
2789: }
2790: VecRestoreArray(coordinates, &coords);
2791: VecRestoreArray(subCoordinates, &subCoords);
2792: }
2793: DMSetCoordinatesLocal(subdm, subCoordinates);
2794: VecDestroy(&subCoordinates);
2795: }
2796: /* Cleanup */
2797: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
2798: ISDestroy(&subvertexIS);
2799: if (subcellIS) {ISRestoreIndices(subcellIS, &subCells);}
2800: ISDestroy(&subcellIS);
2801: return(0);
2802: }
2804: PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
2805: {
2806: PetscInt subPoint;
2809: PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr;
2810: return subPoint < 0 ? subPoint : firstSubPoint+subPoint;
2811: }
2813: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
2814: {
2815: MPI_Comm comm;
2816: DMLabel subpointMap;
2817: IS *subpointIS;
2818: const PetscInt **subpoints;
2819: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
2820: PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v;
2821: PetscMPIInt rank;
2822: PetscErrorCode ierr;
2825: PetscObjectGetComm((PetscObject)dm,&comm);
2826: MPI_Comm_rank(comm, &rank);
2827: /* Create subpointMap which marks the submesh */
2828: DMLabelCreate("subpoint_map", &subpointMap);
2829: DMPlexSetSubpointMap(subdm, subpointMap);
2830: if (cellHeight) {
2831: if (isCohesive) {DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2832: else {DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);}
2833: } else {
2834: DMLabel depth;
2835: IS pointIS;
2836: const PetscInt *points;
2837: PetscInt numPoints;
2839: DMPlexGetDepthLabel(dm, &depth);
2840: DMLabelGetStratumSize(label, value, &numPoints);
2841: DMLabelGetStratumIS(label, value, &pointIS);
2842: ISGetIndices(pointIS, &points);
2843: for (p = 0; p < numPoints; ++p) {
2844: PetscInt *closure = NULL;
2845: PetscInt closureSize, c, pdim;
2847: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2848: for (c = 0; c < closureSize*2; c += 2) {
2849: DMLabelGetValue(depth, closure[c], &pdim);
2850: DMLabelSetValue(subpointMap, closure[c], pdim);
2851: }
2852: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2853: }
2854: ISRestoreIndices(pointIS, &points);
2855: ISDestroy(&pointIS);
2856: }
2857: DMLabelDestroy(&subpointMap);
2858: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
2859: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2860: cMax = (cMax < 0) ? cEnd : cMax;
2861: /* Setup chart */
2862: DMGetDimension(dm, &dim);
2863: PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);
2864: for (d = 0; d <= dim; ++d) {
2865: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
2866: totSubPoints += numSubPoints[d];
2867: }
2868: DMPlexSetChart(subdm, 0, totSubPoints);
2869: DMPlexSetVTKCellHeight(subdm, cellHeight);
2870: /* Set cone sizes */
2871: firstSubPoint[dim] = 0;
2872: firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim];
2873: if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];}
2874: if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];}
2875: for (d = 0; d <= dim; ++d) {
2876: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
2877: if (subpointIS[d]) {ISGetIndices(subpointIS[d], &subpoints[d]);}
2878: }
2879: for (d = 0; d <= dim; ++d) {
2880: for (p = 0; p < numSubPoints[d]; ++p) {
2881: const PetscInt point = subpoints[d][p];
2882: const PetscInt subpoint = firstSubPoint[d] + p;
2883: const PetscInt *cone;
2884: PetscInt coneSize, coneSizeNew, c, val;
2886: DMPlexGetConeSize(dm, point, &coneSize);
2887: DMPlexSetConeSize(subdm, subpoint, coneSize);
2888: if (cellHeight && (d == dim)) {
2889: DMPlexGetCone(dm, point, &cone);
2890: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2891: DMLabelGetValue(subpointMap, cone[c], &val);
2892: if (val >= 0) coneSizeNew++;
2893: }
2894: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
2895: }
2896: }
2897: }
2898: DMSetUp(subdm);
2899: /* Set cones */
2900: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2901: PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);
2902: for (d = 0; d <= dim; ++d) {
2903: for (p = 0; p < numSubPoints[d]; ++p) {
2904: const PetscInt point = subpoints[d][p];
2905: const PetscInt subpoint = firstSubPoint[d] + p;
2906: const PetscInt *cone, *ornt;
2907: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
2909: if (d == dim-1) {
2910: const PetscInt *support, *cone, *ornt;
2911: PetscInt supportSize, coneSize, s, subc;
2913: DMPlexGetSupport(dm, point, &support);
2914: DMPlexGetSupportSize(dm, point, &supportSize);
2915: for (s = 0; s < supportSize; ++s) {
2916: if ((support[s] < cMax) || (support[s] >= cEnd)) continue;
2917: PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);
2918: if (subc >= 0) {
2919: const PetscInt ccell = subpoints[d+1][subc];
2921: DMPlexGetCone(dm, ccell, &cone);
2922: DMPlexGetConeSize(dm, ccell, &coneSize);
2923: DMPlexGetConeOrientation(dm, ccell, &ornt);
2924: for (c = 0; c < coneSize; ++c) {
2925: if (cone[c] == point) {
2926: fornt = ornt[c];
2927: break;
2928: }
2929: }
2930: break;
2931: }
2932: }
2933: }
2934: DMPlexGetConeSize(dm, point, &coneSize);
2935: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
2936: DMPlexGetCone(dm, point, &cone);
2937: DMPlexGetConeOrientation(dm, point, &ornt);
2938: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
2939: PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);
2940: if (subc >= 0) {
2941: coneNew[coneSizeNew] = firstSubPoint[d-1] + subc;
2942: orntNew[coneSizeNew] = ornt[c];
2943: ++coneSizeNew;
2944: }
2945: }
2946: if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize);
2947: if (fornt < 0) {
2948: /* This should be replaced by a call to DMPlexReverseCell() */
2949: #if 0
2950: DMPlexReverseCell(subdm, subpoint);
2951: #else
2952: for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) {
2953: PetscInt faceSize, tmp;
2955: tmp = coneNew[c];
2956: coneNew[c] = coneNew[coneSizeNew-1-c];
2957: coneNew[coneSizeNew-1-c] = tmp;
2958: DMPlexGetConeSize(dm, cone[c], &faceSize);
2959: tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c];
2960: orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c];
2961: orntNew[coneSizeNew-1-c] = tmp;
2962: }
2963: }
2964: DMPlexSetCone(subdm, subpoint, coneNew);
2965: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
2966: #endif
2967: }
2968: }
2969: PetscFree2(coneNew,orntNew);
2970: DMPlexSymmetrize(subdm);
2971: DMPlexStratify(subdm);
2972: /* Build coordinates */
2973: {
2974: PetscSection coordSection, subCoordSection;
2975: Vec coordinates, subCoordinates;
2976: PetscScalar *coords, *subCoords;
2977: PetscInt cdim, numComp, coordSize;
2978: const char *name;
2980: DMGetCoordinateDim(dm, &cdim);
2981: DMGetCoordinateSection(dm, &coordSection);
2982: DMGetCoordinatesLocal(dm, &coordinates);
2983: DMGetCoordinateSection(subdm, &subCoordSection);
2984: PetscSectionSetNumFields(subCoordSection, 1);
2985: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
2986: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
2987: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);
2988: for (v = 0; v < numSubPoints[0]; ++v) {
2989: const PetscInt vertex = subpoints[0][v];
2990: const PetscInt subvertex = firstSubPoint[0]+v;
2991: PetscInt dof;
2993: PetscSectionGetDof(coordSection, vertex, &dof);
2994: PetscSectionSetDof(subCoordSection, subvertex, dof);
2995: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
2996: }
2997: PetscSectionSetUp(subCoordSection);
2998: PetscSectionGetStorageSize(subCoordSection, &coordSize);
2999: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3000: PetscObjectGetName((PetscObject)coordinates,&name);
3001: PetscObjectSetName((PetscObject)subCoordinates,name);
3002: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3003: VecSetBlockSize(subCoordinates, cdim);
3004: VecSetType(subCoordinates,VECSTANDARD);
3005: VecGetArray(coordinates, &coords);
3006: VecGetArray(subCoordinates, &subCoords);
3007: for (v = 0; v < numSubPoints[0]; ++v) {
3008: const PetscInt vertex = subpoints[0][v];
3009: const PetscInt subvertex = firstSubPoint[0]+v;
3010: PetscInt dof, off, sdof, soff, d;
3012: PetscSectionGetDof(coordSection, vertex, &dof);
3013: PetscSectionGetOffset(coordSection, vertex, &off);
3014: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3015: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3016: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3017: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3018: }
3019: VecRestoreArray(coordinates, &coords);
3020: VecRestoreArray(subCoordinates, &subCoords);
3021: DMSetCoordinatesLocal(subdm, subCoordinates);
3022: VecDestroy(&subCoordinates);
3023: }
3024: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3025: {
3026: PetscSF sfPoint, sfPointSub;
3027: IS subpIS;
3028: const PetscSFNode *remotePoints;
3029: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3030: const PetscInt *localPoints, *subpoints;
3031: PetscInt *slocalPoints;
3032: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p;
3033: PetscMPIInt rank;
3035: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3036: DMGetPointSF(dm, &sfPoint);
3037: DMGetPointSF(subdm, &sfPointSub);
3038: DMPlexGetChart(dm, &pStart, &pEnd);
3039: DMPlexGetChart(subdm, NULL, &numSubroots);
3040: DMPlexCreateSubpointIS(subdm, &subpIS);
3041: if (subpIS) {
3042: ISGetIndices(subpIS, &subpoints);
3043: ISGetLocalSize(subpIS, &numSubpoints);
3044: }
3045: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3046: if (numRoots >= 0) {
3047: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3048: for (p = 0; p < pEnd-pStart; ++p) {
3049: newLocalPoints[p].rank = -2;
3050: newLocalPoints[p].index = -2;
3051: }
3052: /* Set subleaves */
3053: for (l = 0; l < numLeaves; ++l) {
3054: const PetscInt point = localPoints[l];
3055: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3057: if (subpoint < 0) continue;
3058: newLocalPoints[point-pStart].rank = rank;
3059: newLocalPoints[point-pStart].index = subpoint;
3060: ++numSubleaves;
3061: }
3062: /* Must put in owned subpoints */
3063: for (p = pStart; p < pEnd; ++p) {
3064: const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints);
3066: if (subpoint < 0) {
3067: newOwners[p-pStart].rank = -3;
3068: newOwners[p-pStart].index = -3;
3069: } else {
3070: newOwners[p-pStart].rank = rank;
3071: newOwners[p-pStart].index = subpoint;
3072: }
3073: }
3074: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3075: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3076: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3077: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3078: PetscMalloc1(numSubleaves, &slocalPoints);
3079: PetscMalloc1(numSubleaves, &sremotePoints);
3080: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3081: const PetscInt point = localPoints[l];
3082: const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints);
3084: if (subpoint < 0) continue;
3085: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3086: slocalPoints[sl] = subpoint;
3087: sremotePoints[sl].rank = newLocalPoints[point].rank;
3088: sremotePoints[sl].index = newLocalPoints[point].index;
3089: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3090: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3091: ++sl;
3092: }
3093: if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves);
3094: PetscFree2(newLocalPoints,newOwners);
3095: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3096: }
3097: if (subpIS) {
3098: ISRestoreIndices(subpIS, &subpoints);
3099: ISDestroy(&subpIS);
3100: }
3101: }
3102: /* Cleanup */
3103: for (d = 0; d <= dim; ++d) {
3104: if (subpointIS[d]) {ISRestoreIndices(subpointIS[d], &subpoints[d]);}
3105: ISDestroy(&subpointIS[d]);
3106: }
3107: PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);
3108: return(0);
3109: }
3111: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3112: {
3116: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);
3117: return(0);
3118: }
3120: /*@
3121: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3123: Input Parameters:
3124: + dm - The original mesh
3125: . vertexLabel - The DMLabel marking vertices contained in the surface
3126: - value - The label value to use
3128: Output Parameter:
3129: . subdm - The surface mesh
3131: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3133: Level: developer
3135: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3136: @*/
3137: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm)
3138: {
3139: PetscInt dim, cdim, depth;
3145: DMGetDimension(dm, &dim);
3146: DMPlexGetDepth(dm, &depth);
3147: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3148: DMSetType(*subdm, DMPLEX);
3149: DMSetDimension(*subdm, dim-1);
3150: DMGetCoordinateDim(dm, &cdim);
3151: DMSetCoordinateDim(*subdm, cdim);
3152: if (depth == dim) {
3153: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);
3154: } else {
3155: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3156: }
3157: return(0);
3158: }
3160: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3161: {
3162: MPI_Comm comm;
3163: DMLabel subpointMap;
3164: IS subvertexIS;
3165: const PetscInt *subVertices;
3166: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3167: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3168: PetscInt cMax, c, f;
3169: PetscErrorCode ierr;
3172: PetscObjectGetComm((PetscObject)dm, &comm);
3173: /* Create subpointMap which marks the submesh */
3174: DMLabelCreate("subpoint_map", &subpointMap);
3175: DMPlexSetSubpointMap(subdm, subpointMap);
3176: DMLabelDestroy(&subpointMap);
3177: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3178: /* Setup chart */
3179: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3180: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3181: DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);
3182: DMPlexSetVTKCellHeight(subdm, 1);
3183: /* Set cone sizes */
3184: firstSubVertex = numSubCells;
3185: firstSubFace = numSubCells+numSubVertices;
3186: newFacePoint = firstSubFace;
3187: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3188: if (subvertexIS) {ISGetIndices(subvertexIS, &subVertices);}
3189: for (c = 0; c < numSubCells; ++c) {
3190: DMPlexSetConeSize(subdm, c, 1);
3191: }
3192: for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) {
3193: DMPlexSetConeSize(subdm, f, nFV);
3194: }
3195: DMSetUp(subdm);
3196: /* Create face cones */
3197: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3198: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
3199: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3200: for (c = 0; c < numSubCells; ++c) {
3201: const PetscInt cell = subCells[c];
3202: const PetscInt subcell = c;
3203: const PetscInt *cone, *cells;
3204: PetscInt numCells, subVertex, p, v;
3206: if (cell < cMax) continue;
3207: DMPlexGetCone(dm, cell, &cone);
3208: for (v = 0; v < nFV; ++v) {
3209: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3210: subface[v] = firstSubVertex+subVertex;
3211: }
3212: DMPlexSetCone(subdm, newFacePoint, subface);
3213: DMPlexSetCone(subdm, subcell, &newFacePoint);
3214: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3215: /* Not true in parallel
3216: if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3217: for (p = 0; p < numCells; ++p) {
3218: PetscInt negsubcell;
3220: if (cells[p] >= cMax) continue;
3221: /* I know this is a crap search */
3222: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3223: if (subCells[negsubcell] == cells[p]) break;
3224: }
3225: if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell);
3226: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3227: }
3228: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3229: ++newFacePoint;
3230: }
3231: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);
3232: DMPlexSymmetrize(subdm);
3233: DMPlexStratify(subdm);
3234: /* Build coordinates */
3235: {
3236: PetscSection coordSection, subCoordSection;
3237: Vec coordinates, subCoordinates;
3238: PetscScalar *coords, *subCoords;
3239: PetscInt cdim, numComp, coordSize, v;
3240: const char *name;
3242: DMGetCoordinateDim(dm, &cdim);
3243: DMGetCoordinateSection(dm, &coordSection);
3244: DMGetCoordinatesLocal(dm, &coordinates);
3245: DMGetCoordinateSection(subdm, &subCoordSection);
3246: PetscSectionSetNumFields(subCoordSection, 1);
3247: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3248: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3249: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);
3250: for (v = 0; v < numSubVertices; ++v) {
3251: const PetscInt vertex = subVertices[v];
3252: const PetscInt subvertex = firstSubVertex+v;
3253: PetscInt dof;
3255: PetscSectionGetDof(coordSection, vertex, &dof);
3256: PetscSectionSetDof(subCoordSection, subvertex, dof);
3257: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3258: }
3259: PetscSectionSetUp(subCoordSection);
3260: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3261: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3262: PetscObjectGetName((PetscObject)coordinates,&name);
3263: PetscObjectSetName((PetscObject)subCoordinates,name);
3264: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3265: VecSetBlockSize(subCoordinates, cdim);
3266: VecSetType(subCoordinates,VECSTANDARD);
3267: VecGetArray(coordinates, &coords);
3268: VecGetArray(subCoordinates, &subCoords);
3269: for (v = 0; v < numSubVertices; ++v) {
3270: const PetscInt vertex = subVertices[v];
3271: const PetscInt subvertex = firstSubVertex+v;
3272: PetscInt dof, off, sdof, soff, d;
3274: PetscSectionGetDof(coordSection, vertex, &dof);
3275: PetscSectionGetOffset(coordSection, vertex, &off);
3276: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3277: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3278: if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof);
3279: for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d];
3280: }
3281: VecRestoreArray(coordinates, &coords);
3282: VecRestoreArray(subCoordinates, &subCoords);
3283: DMSetCoordinatesLocal(subdm, subCoordinates);
3284: VecDestroy(&subCoordinates);
3285: }
3286: /* Build SF */
3287: CHKMEMQ;
3288: {
3289: PetscSF sfPoint, sfPointSub;
3290: const PetscSFNode *remotePoints;
3291: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3292: const PetscInt *localPoints;
3293: PetscInt *slocalPoints;
3294: PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3295: PetscMPIInt rank;
3297: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3298: DMGetPointSF(dm, &sfPoint);
3299: DMGetPointSF(subdm, &sfPointSub);
3300: DMPlexGetChart(dm, &pStart, &pEnd);
3301: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3302: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3303: if (numRoots >= 0) {
3304: /* Only vertices should be shared */
3305: PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);
3306: for (p = 0; p < pEnd-pStart; ++p) {
3307: newLocalPoints[p].rank = -2;
3308: newLocalPoints[p].index = -2;
3309: }
3310: /* Set subleaves */
3311: for (l = 0; l < numLeaves; ++l) {
3312: const PetscInt point = localPoints[l];
3313: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3315: if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point);
3316: if (subPoint < 0) continue;
3317: newLocalPoints[point-pStart].rank = rank;
3318: newLocalPoints[point-pStart].index = subPoint;
3319: ++numSubLeaves;
3320: }
3321: /* Must put in owned subpoints */
3322: for (p = pStart; p < pEnd; ++p) {
3323: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3325: if (subPoint < 0) {
3326: newOwners[p-pStart].rank = -3;
3327: newOwners[p-pStart].index = -3;
3328: } else {
3329: newOwners[p-pStart].rank = rank;
3330: newOwners[p-pStart].index = subPoint;
3331: }
3332: }
3333: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3334: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3335: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3336: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);
3337: PetscMalloc1(numSubLeaves, &slocalPoints);
3338: PetscMalloc1(numSubLeaves, &sremotePoints);
3339: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3340: const PetscInt point = localPoints[l];
3341: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3343: if (subPoint < 0) continue;
3344: if (newLocalPoints[point].rank == rank) {++ll; continue;}
3345: slocalPoints[sl] = subPoint;
3346: sremotePoints[sl].rank = newLocalPoints[point].rank;
3347: sremotePoints[sl].index = newLocalPoints[point].index;
3348: if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point);
3349: if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point);
3350: ++sl;
3351: }
3352: PetscFree2(newLocalPoints,newOwners);
3353: if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves);
3354: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3355: }
3356: }
3357: CHKMEMQ;
3358: /* Cleanup */
3359: if (subvertexIS) {ISRestoreIndices(subvertexIS, &subVertices);}
3360: ISDestroy(&subvertexIS);
3361: PetscFree(subCells);
3362: return(0);
3363: }
3365: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3366: {
3367: DMLabel label = NULL;
3371: if (labelname) {DMGetLabel(dm, labelname, &label);}
3372: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);
3373: return(0);
3374: }
3376: /*@C
3377: 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.
3379: Input Parameters:
3380: + dm - The original mesh
3381: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3382: . label - A label name, or NULL
3383: - value - A label value
3385: Output Parameter:
3386: . subdm - The surface mesh
3388: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3390: Level: developer
3392: .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh()
3393: @*/
3394: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3395: {
3396: PetscInt dim, cdim, depth;
3402: DMGetDimension(dm, &dim);
3403: DMPlexGetDepth(dm, &depth);
3404: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3405: DMSetType(*subdm, DMPLEX);
3406: DMSetDimension(*subdm, dim-1);
3407: DMGetCoordinateDim(dm, &cdim);
3408: DMSetCoordinateDim(*subdm, cdim);
3409: if (depth == dim) {
3410: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3411: } else {
3412: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3413: }
3414: return(0);
3415: }
3417: /*@
3418: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3420: Input Parameters:
3421: + dm - The original mesh
3422: . cellLabel - The DMLabel marking cells contained in the new mesh
3423: - value - The label value to use
3425: Output Parameter:
3426: . subdm - The new mesh
3428: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3430: Level: developer
3432: .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue()
3433: @*/
3434: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3435: {
3436: PetscInt dim;
3442: DMGetDimension(dm, &dim);
3443: DMCreate(PetscObjectComm((PetscObject) dm), subdm);
3444: DMSetType(*subdm, DMPLEX);
3445: DMSetDimension(*subdm, dim);
3446: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3447: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);
3448: return(0);
3449: }
3451: /*@
3452: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3454: Input Parameter:
3455: . dm - The submesh DM
3457: Output Parameter:
3458: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3460: Level: developer
3462: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3463: @*/
3464: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3465: {
3469: *subpointMap = ((DM_Plex*) dm->data)->subpointMap;
3470: return(0);
3471: }
3473: /*@
3474: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3476: Input Parameters:
3477: + dm - The submesh DM
3478: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3480: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3482: Level: developer
3484: .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS()
3485: @*/
3486: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3487: {
3488: DM_Plex *mesh = (DM_Plex *) dm->data;
3489: DMLabel tmp;
3494: tmp = mesh->subpointMap;
3495: mesh->subpointMap = subpointMap;
3496: ++mesh->subpointMap->refct;
3497: DMLabelDestroy(&tmp);
3498: return(0);
3499: }
3501: /*@
3502: DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data
3504: Input Parameter:
3505: . dm - The submesh DM
3507: Output Parameter:
3508: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3510: Note: This IS is guaranteed to be sorted by the construction of the submesh
3512: Level: developer
3514: .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap()
3515: @*/
3516: PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS)
3517: {
3518: MPI_Comm comm;
3519: DMLabel subpointMap;
3520: IS is;
3521: const PetscInt *opoints;
3522: PetscInt *points, *depths;
3523: PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off;
3524: PetscErrorCode ierr;
3529: PetscObjectGetComm((PetscObject)dm,&comm);
3530: *subpointIS = NULL;
3531: DMPlexGetSubpointMap(dm, &subpointMap);
3532: DMPlexGetDepth(dm, &depth);
3533: if (subpointMap && depth >= 0) {
3534: DMPlexGetChart(dm, &pStart, &pEnd);
3535: if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart);
3536: DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);
3537: depths[0] = depth;
3538: depths[1] = 0;
3539: for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;}
3540: PetscMalloc1(pEnd, &points);
3541: for(d = 0, off = 0; d <= depth; ++d) {
3542: const PetscInt dep = depths[d];
3544: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3545: DMLabelGetStratumSize(subpointMap, dep, &n);
3546: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3547: 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);
3548: } else {
3549: if (!n) {
3550: if (d == 0) {
3551: /* Missing cells */
3552: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1;
3553: } else {
3554: /* Missing faces */
3555: for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3556: }
3557: }
3558: }
3559: if (n) {
3560: DMLabelGetStratumIS(subpointMap, dep, &is);
3561: ISGetIndices(is, &opoints);
3562: for(p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3563: ISRestoreIndices(is, &opoints);
3564: ISDestroy(&is);
3565: }
3566: }
3567: DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);
3568: if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd);
3569: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3570: }
3571: return(0);
3572: }