Actual source code: plexsubmesh.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petscsf.h>
5: static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid)
6: {
7: DMPolytopeType ct;
9: DMPlexGetCellType(dm, p, &ct);
10: switch (ct) {
11: case DM_POLYTOPE_POINT_PRISM_TENSOR:
12: case DM_POLYTOPE_SEG_PRISM_TENSOR:
13: case DM_POLYTOPE_TRI_PRISM_TENSOR:
14: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
15: *isHybrid = PETSC_TRUE;
16: break;
17: default:
18: *isHybrid = PETSC_FALSE;
19: break;
20: }
21: return 0;
22: }
24: static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd)
25: {
26: DMLabel ctLabel;
28: if (cStart) *cStart = -1;
29: if (cEnd) *cEnd = -1;
30: DMPlexGetCellTypeLabel(dm, &ctLabel);
31: switch (dim) {
32: case 1:
33: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd);
34: break;
35: case 2:
36: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd);
37: break;
38: case 3:
39: DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd);
40: if (*cStart < 0) DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd);
41: break;
42: default:
43: return 0;
44: }
45: return 0;
46: }
48: static PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label)
49: {
50: PetscInt fStart, fEnd, f;
52: DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
53: for (f = fStart; f < fEnd; ++f) {
54: PetscInt supportSize;
56: DMPlexGetSupportSize(dm, f, &supportSize);
57: if (supportSize == 1) {
58: if (val < 0) {
59: PetscInt *closure = NULL;
60: PetscInt clSize, cl, cval;
62: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
63: for (cl = 0; cl < clSize * 2; cl += 2) {
64: DMLabelGetValue(label, closure[cl], &cval);
65: if (cval < 0) continue;
66: DMLabelSetValue(label, f, cval);
67: break;
68: }
69: if (cl == clSize * 2) DMLabelSetValue(label, f, 1);
70: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure);
71: } else {
72: DMLabelSetValue(label, f, val);
73: }
74: }
75: }
76: return 0;
77: }
79: /*@
80: DMPlexMarkBoundaryFaces - Mark all faces on the boundary
82: Not Collective
84: Input Parameters:
85: + dm - The original DM
86: - val - The marker value, or PETSC_DETERMINE to use some value in the closure (or 1 if none are found)
88: Output Parameter:
89: . label - The DMLabel marking boundary faces with the given value
91: Level: developer
93: .seealso: `DMLabelCreate()`, `DMCreateLabel()`
94: @*/
95: PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label)
96: {
97: DMPlexInterpolatedFlag flg;
100: DMPlexIsInterpolated(dm, &flg);
102: DMPlexMarkBoundaryFaces_Internal(dm, val, 0, label);
103: return 0;
104: }
106: static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells)
107: {
108: IS valueIS;
109: PetscSF sfPoint;
110: const PetscInt *values;
111: PetscInt numValues, v, cStart, cEnd, nroots;
113: DMLabelGetNumValues(label, &numValues);
114: DMLabelGetValueIS(label, &valueIS);
115: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
116: ISGetIndices(valueIS, &values);
117: for (v = 0; v < numValues; ++v) {
118: IS pointIS;
119: const PetscInt *points;
120: PetscInt numPoints, p;
122: DMLabelGetStratumSize(label, values[v], &numPoints);
123: DMLabelGetStratumIS(label, values[v], &pointIS);
124: ISGetIndices(pointIS, &points);
125: for (p = 0; p < numPoints; ++p) {
126: PetscInt q = points[p];
127: PetscInt *closure = NULL;
128: PetscInt closureSize, c;
130: if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */
131: continue;
132: }
133: DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
134: for (c = 0; c < closureSize * 2; c += 2) DMLabelSetValue(label, closure[c], values[v]);
135: DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure);
136: }
137: ISRestoreIndices(pointIS, &points);
138: ISDestroy(&pointIS);
139: }
140: ISRestoreIndices(valueIS, &values);
141: ISDestroy(&valueIS);
142: DMGetPointSF(dm, &sfPoint);
143: PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL);
144: if (nroots >= 0) {
145: DMLabel lblRoots, lblLeaves;
146: IS valueIS, pointIS;
147: const PetscInt *values;
148: PetscInt numValues, v;
150: /* Pull point contributions from remote leaves into local roots */
151: DMLabelGather(label, sfPoint, &lblLeaves);
152: DMLabelGetValueIS(lblLeaves, &valueIS);
153: ISGetLocalSize(valueIS, &numValues);
154: ISGetIndices(valueIS, &values);
155: for (v = 0; v < numValues; ++v) {
156: const PetscInt value = values[v];
158: DMLabelGetStratumIS(lblLeaves, value, &pointIS);
159: DMLabelInsertIS(label, pointIS, value);
160: ISDestroy(&pointIS);
161: }
162: ISRestoreIndices(valueIS, &values);
163: ISDestroy(&valueIS);
164: DMLabelDestroy(&lblLeaves);
165: /* Push point contributions from roots into remote leaves */
166: DMLabelDistribute(label, sfPoint, &lblRoots);
167: DMLabelGetValueIS(lblRoots, &valueIS);
168: ISGetLocalSize(valueIS, &numValues);
169: ISGetIndices(valueIS, &values);
170: for (v = 0; v < numValues; ++v) {
171: const PetscInt value = values[v];
173: DMLabelGetStratumIS(lblRoots, value, &pointIS);
174: DMLabelInsertIS(label, pointIS, value);
175: ISDestroy(&pointIS);
176: }
177: ISRestoreIndices(valueIS, &values);
178: ISDestroy(&valueIS);
179: DMLabelDestroy(&lblRoots);
180: }
181: return 0;
182: }
184: /*@
185: DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface
187: Input Parameters:
188: + dm - The DM
189: - label - A DMLabel marking the surface points
191: Output Parameter:
192: . label - A DMLabel marking all surface points in the transitive closure
194: Level: developer
196: .seealso: `DMPlexLabelCohesiveComplete()`
197: @*/
198: PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label)
199: {
200: DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE);
201: return 0;
202: }
204: /*@
205: DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point
207: Input Parameters:
208: + dm - The DM
209: - label - A DMLabel marking the surface points
211: Output Parameter:
212: . label - A DMLabel incorporating cells
214: Level: developer
216: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
218: .seealso: `DMPlexLabelAddFaceCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
219: @*/
220: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
221: {
222: IS valueIS;
223: const PetscInt *values;
224: PetscInt numValues, v, cStart, cEnd;
226: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
227: DMLabelGetNumValues(label, &numValues);
228: DMLabelGetValueIS(label, &valueIS);
229: ISGetIndices(valueIS, &values);
230: for (v = 0; v < numValues; ++v) {
231: IS pointIS;
232: const PetscInt *points;
233: PetscInt numPoints, p;
235: DMLabelGetStratumSize(label, values[v], &numPoints);
236: DMLabelGetStratumIS(label, values[v], &pointIS);
237: ISGetIndices(pointIS, &points);
238: for (p = 0; p < numPoints; ++p) {
239: PetscInt *closure = NULL;
240: PetscInt closureSize, cl;
242: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
243: for (cl = closureSize - 1; cl > 0; --cl) {
244: const PetscInt cell = closure[cl * 2];
245: if ((cell >= cStart) && (cell < cEnd)) {
246: DMLabelSetValue(label, cell, values[v]);
247: break;
248: }
249: }
250: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure);
251: }
252: ISRestoreIndices(pointIS, &points);
253: ISDestroy(&pointIS);
254: }
255: ISRestoreIndices(valueIS, &values);
256: ISDestroy(&valueIS);
257: return 0;
258: }
260: /*@
261: DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face
263: Input Parameters:
264: + dm - The DM
265: - label - A DMLabel marking the surface points
267: Output Parameter:
268: . label - A DMLabel incorporating cells
270: Level: developer
272: Note: The cells allow FEM boundary conditions to be applied using the cell geometry
274: .seealso: `DMPlexLabelAddCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
275: @*/
276: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
277: {
278: IS valueIS;
279: const PetscInt *values;
280: PetscInt numValues, v, cStart, cEnd, fStart, fEnd;
282: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
283: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
284: DMLabelGetNumValues(label, &numValues);
285: DMLabelGetValueIS(label, &valueIS);
286: ISGetIndices(valueIS, &values);
287: for (v = 0; v < numValues; ++v) {
288: IS pointIS;
289: const PetscInt *points;
290: PetscInt numPoints, p;
292: DMLabelGetStratumSize(label, values[v], &numPoints);
293: DMLabelGetStratumIS(label, values[v], &pointIS);
294: ISGetIndices(pointIS, &points);
295: for (p = 0; p < numPoints; ++p) {
296: const PetscInt face = points[p];
297: PetscInt *closure = NULL;
298: PetscInt closureSize, cl;
300: if ((face < fStart) || (face >= fEnd)) continue;
301: DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
302: for (cl = closureSize - 1; cl > 0; --cl) {
303: const PetscInt cell = closure[cl * 2];
304: if ((cell >= cStart) && (cell < cEnd)) {
305: DMLabelSetValue(label, cell, values[v]);
306: break;
307: }
308: }
309: DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure);
310: }
311: ISRestoreIndices(pointIS, &points);
312: ISDestroy(&pointIS);
313: }
314: ISRestoreIndices(valueIS, &values);
315: ISDestroy(&valueIS);
316: return 0;
317: }
319: /*@
320: DMPlexLabelClearCells - Remove cells from a label
322: Input Parameters:
323: + dm - The DM
324: - label - A DMLabel marking surface points and their adjacent cells
326: Output Parameter:
327: . label - A DMLabel without cells
329: Level: developer
331: Note: This undoes DMPlexLabelAddCells() or DMPlexLabelAddFaceCells()
333: .seealso: `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`, `DMPlexLabelAddCells()`
334: @*/
335: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
336: {
337: IS valueIS;
338: const PetscInt *values;
339: PetscInt numValues, v, cStart, cEnd;
341: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
342: DMLabelGetNumValues(label, &numValues);
343: DMLabelGetValueIS(label, &valueIS);
344: ISGetIndices(valueIS, &values);
345: for (v = 0; v < numValues; ++v) {
346: IS pointIS;
347: const PetscInt *points;
348: PetscInt numPoints, p;
350: DMLabelGetStratumSize(label, values[v], &numPoints);
351: DMLabelGetStratumIS(label, values[v], &pointIS);
352: ISGetIndices(pointIS, &points);
353: for (p = 0; p < numPoints; ++p) {
354: PetscInt point = points[p];
356: if (point >= cStart && point < cEnd) DMLabelClearValue(label, point, values[v]);
357: }
358: ISRestoreIndices(pointIS, &points);
359: ISDestroy(&pointIS);
360: }
361: ISRestoreIndices(valueIS, &values);
362: ISDestroy(&valueIS);
363: return 0;
364: }
366: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
367: * index (skipping first, which is (0,0)) */
368: static inline PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
369: {
370: PetscInt d, off = 0;
372: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
373: for (d = 0; d < depth; d++) {
374: PetscInt firstd = d;
375: PetscInt firstStart = depthShift[2 * d];
376: PetscInt e;
378: for (e = d + 1; e <= depth; e++) {
379: if (depthShift[2 * e] < firstStart) {
380: firstd = e;
381: firstStart = depthShift[2 * d];
382: }
383: }
384: if (firstd != d) {
385: PetscInt swap[2];
387: e = firstd;
388: swap[0] = depthShift[2 * d];
389: swap[1] = depthShift[2 * d + 1];
390: depthShift[2 * d] = depthShift[2 * e];
391: depthShift[2 * d + 1] = depthShift[2 * e + 1];
392: depthShift[2 * e] = swap[0];
393: depthShift[2 * e + 1] = swap[1];
394: }
395: }
396: /* convert (oldstart, added) to (oldstart, newstart) */
397: for (d = 0; d <= depth; d++) {
398: off += depthShift[2 * d + 1];
399: depthShift[2 * d + 1] = depthShift[2 * d] + off;
400: }
401: return 0;
402: }
404: /* depthShift is a list of (old, new) pairs */
405: static inline PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
406: {
407: PetscInt d;
408: PetscInt newOff = 0;
410: for (d = 0; d <= depth; d++) {
411: if (p < depthShift[2 * d]) return p + newOff;
412: else newOff = depthShift[2 * d + 1] - depthShift[2 * d];
413: }
414: return p + newOff;
415: }
417: /* depthShift is a list of (old, new) pairs */
418: static inline PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
419: {
420: PetscInt d;
421: PetscInt newOff = 0;
423: for (d = 0; d <= depth; d++) {
424: if (p < depthShift[2 * d + 1]) return p + newOff;
425: else newOff = depthShift[2 * d] - depthShift[2 * d + 1];
426: }
427: return p + newOff;
428: }
430: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
431: {
432: PetscInt depth = 0, d, pStart, pEnd, p;
433: DMLabel depthLabel;
435: DMPlexGetDepth(dm, &depth);
436: if (depth < 0) return 0;
437: /* Step 1: Expand chart */
438: DMPlexGetChart(dm, &pStart, &pEnd);
439: pEnd = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
440: DMPlexSetChart(dmNew, pStart, pEnd);
441: DMCreateLabel(dmNew, "depth");
442: DMPlexGetDepthLabel(dmNew, &depthLabel);
443: DMCreateLabel(dmNew, "celltype");
444: /* Step 2: Set cone and support sizes */
445: for (d = 0; d <= depth; ++d) {
446: PetscInt pStartNew, pEndNew;
447: IS pIS;
449: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
450: pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
451: pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
452: ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS);
453: DMLabelSetStratumIS(depthLabel, d, pIS);
454: ISDestroy(&pIS);
455: for (p = pStart; p < pEnd; ++p) {
456: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
457: PetscInt size;
458: DMPolytopeType ct;
460: DMPlexGetConeSize(dm, p, &size);
461: DMPlexSetConeSize(dmNew, newp, size);
462: DMPlexGetSupportSize(dm, p, &size);
463: DMPlexSetSupportSize(dmNew, newp, size);
464: DMPlexGetCellType(dm, p, &ct);
465: DMPlexSetCellType(dmNew, newp, ct);
466: }
467: }
468: return 0;
469: }
471: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
472: {
473: PetscInt *newpoints;
474: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
476: DMPlexGetDepth(dm, &depth);
477: if (depth < 0) return 0;
478: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
479: DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew);
480: PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)), &newpoints);
481: /* Step 5: Set cones and supports */
482: DMPlexGetChart(dm, &pStart, &pEnd);
483: for (p = pStart; p < pEnd; ++p) {
484: const PetscInt *points = NULL, *orientations = NULL;
485: PetscInt size, sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
487: DMPlexGetConeSize(dm, p, &size);
488: DMPlexGetCone(dm, p, &points);
489: DMPlexGetConeOrientation(dm, p, &orientations);
490: for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
491: DMPlexSetCone(dmNew, newp, newpoints);
492: DMPlexSetConeOrientation(dmNew, newp, orientations);
493: DMPlexGetSupportSize(dm, p, &size);
494: DMPlexGetSupportSize(dmNew, newp, &sizeNew);
495: DMPlexGetSupport(dm, p, &points);
496: for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
497: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
498: DMPlexSetSupport(dmNew, newp, newpoints);
499: }
500: PetscFree(newpoints);
501: return 0;
502: }
504: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
505: {
506: PetscSection coordSection, newCoordSection;
507: Vec coordinates, newCoordinates;
508: PetscScalar *coords, *newCoords;
509: PetscInt coordSize, sStart, sEnd;
510: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
511: PetscBool hasCells;
513: DMGetCoordinateDim(dm, &dim);
514: DMSetCoordinateDim(dmNew, dim);
515: DMPlexGetDepth(dm, &depth);
516: /* Step 8: Convert coordinates */
517: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
518: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
519: DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);
520: DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);
521: DMGetCoordinateSection(dm, &coordSection);
522: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
523: PetscSectionSetNumFields(newCoordSection, 1);
524: PetscSectionSetFieldComponents(newCoordSection, 0, dim);
525: PetscSectionGetChart(coordSection, &sStart, &sEnd);
526: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
527: PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);
528: if (hasCells) {
529: for (c = cStart; c < cEnd; ++c) {
530: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
532: PetscSectionGetDof(coordSection, c, &dof);
533: PetscSectionSetDof(newCoordSection, cNew, dof);
534: PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);
535: }
536: }
537: for (v = vStartNew; v < vEndNew; ++v) {
538: PetscSectionSetDof(newCoordSection, v, dim);
539: PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
540: }
541: PetscSectionSetUp(newCoordSection);
542: DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);
543: PetscSectionGetStorageSize(newCoordSection, &coordSize);
544: VecCreate(PETSC_COMM_SELF, &newCoordinates);
545: PetscObjectSetName((PetscObject)newCoordinates, "coordinates");
546: VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);
547: VecSetBlockSize(newCoordinates, dim);
548: VecSetType(newCoordinates, VECSTANDARD);
549: DMSetCoordinatesLocal(dmNew, newCoordinates);
550: DMGetCoordinatesLocal(dm, &coordinates);
551: VecGetArray(coordinates, &coords);
552: VecGetArray(newCoordinates, &newCoords);
553: if (hasCells) {
554: for (c = cStart; c < cEnd; ++c) {
555: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
557: PetscSectionGetDof(coordSection, c, &dof);
558: PetscSectionGetOffset(coordSection, c, &off);
559: PetscSectionGetOffset(newCoordSection, cNew, &noff);
560: for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
561: }
562: }
563: for (v = vStart; v < vEnd; ++v) {
564: PetscInt dof, off, noff, d;
566: PetscSectionGetDof(coordSection, v, &dof);
567: PetscSectionGetOffset(coordSection, v, &off);
568: PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);
569: for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
570: }
571: VecRestoreArray(coordinates, &coords);
572: VecRestoreArray(newCoordinates, &newCoords);
573: VecDestroy(&newCoordinates);
574: PetscSectionDestroy(&newCoordSection);
575: return 0;
576: }
578: static PetscErrorCode DMPlexShiftSF_Single(DM dm, PetscInt depthShift[], PetscSF sf, PetscSF sfNew)
579: {
580: const PetscSFNode *remotePoints;
581: PetscSFNode *gremotePoints;
582: const PetscInt *localPoints;
583: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
584: PetscInt numRoots, numLeaves, l, pStart, pEnd, depth = 0, totShift = 0;
586: DMPlexGetDepth(dm, &depth);
587: DMPlexGetChart(dm, &pStart, &pEnd);
588: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
589: totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd;
590: if (numRoots >= 0) {
591: PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation);
592: for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
593: PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
594: PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
595: PetscMalloc1(numLeaves, &glocalPoints);
596: PetscMalloc1(numLeaves, &gremotePoints);
597: for (l = 0; l < numLeaves; ++l) {
598: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
599: gremotePoints[l].rank = remotePoints[l].rank;
600: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
601: }
602: PetscFree2(newLocation, newRemoteLocation);
603: PetscSFSetGraph(sfNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
604: }
605: return 0;
606: }
608: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
609: {
610: PetscSF sfPoint, sfPointNew;
611: PetscBool useNatural;
613: /* Step 9: Convert pointSF */
614: DMGetPointSF(dm, &sfPoint);
615: DMGetPointSF(dmNew, &sfPointNew);
616: DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew);
617: /* Step 9b: Convert naturalSF */
618: DMGetUseNatural(dm, &useNatural);
619: if (useNatural) {
620: PetscSF sfNat, sfNatNew;
622: DMSetUseNatural(dmNew, useNatural);
623: DMGetNaturalSF(dm, &sfNat);
624: DMGetNaturalSF(dmNew, &sfNatNew);
625: DMPlexShiftSF_Single(dm, depthShift, sfNat, sfNatNew);
626: }
627: return 0;
628: }
630: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
631: {
632: PetscInt depth = 0, numLabels, l;
634: DMPlexGetDepth(dm, &depth);
635: /* Step 10: Convert labels */
636: DMGetNumLabels(dm, &numLabels);
637: for (l = 0; l < numLabels; ++l) {
638: DMLabel label, newlabel;
639: const char *lname;
640: PetscBool isDepth, isDim;
641: IS valueIS;
642: const PetscInt *values;
643: PetscInt numValues, val;
645: DMGetLabelName(dm, l, &lname);
646: PetscStrcmp(lname, "depth", &isDepth);
647: if (isDepth) continue;
648: PetscStrcmp(lname, "dim", &isDim);
649: if (isDim) continue;
650: DMCreateLabel(dmNew, lname);
651: DMGetLabel(dm, lname, &label);
652: DMGetLabel(dmNew, lname, &newlabel);
653: DMLabelGetDefaultValue(label, &val);
654: DMLabelSetDefaultValue(newlabel, val);
655: DMLabelGetValueIS(label, &valueIS);
656: ISGetLocalSize(valueIS, &numValues);
657: ISGetIndices(valueIS, &values);
658: for (val = 0; val < numValues; ++val) {
659: IS pointIS;
660: const PetscInt *points;
661: PetscInt numPoints, p;
663: DMLabelGetStratumIS(label, values[val], &pointIS);
664: ISGetLocalSize(pointIS, &numPoints);
665: ISGetIndices(pointIS, &points);
666: for (p = 0; p < numPoints; ++p) {
667: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
669: DMLabelSetValue(newlabel, newpoint, values[val]);
670: }
671: ISRestoreIndices(pointIS, &points);
672: ISDestroy(&pointIS);
673: }
674: ISRestoreIndices(valueIS, &values);
675: ISDestroy(&valueIS);
676: }
677: return 0;
678: }
680: static PetscErrorCode DMPlexCreateVTKLabel_Internal(DM dm, PetscBool createGhostLabel, DM dmNew)
681: {
682: PetscSF sfPoint;
683: DMLabel vtkLabel, ghostLabel = NULL;
684: const PetscSFNode *leafRemote;
685: const PetscInt *leafLocal;
686: PetscInt cellHeight, cStart, cEnd, c, fStart, fEnd, f, numLeaves, l;
687: PetscMPIInt rank;
689: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
690: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
691: DMGetPointSF(dm, &sfPoint);
692: DMPlexGetVTKCellHeight(dmNew, &cellHeight);
693: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
694: PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);
695: DMCreateLabel(dmNew, "vtk");
696: DMGetLabel(dmNew, "vtk", &vtkLabel);
697: if (createGhostLabel) {
698: DMCreateLabel(dmNew, "ghost");
699: DMGetLabel(dmNew, "ghost", &ghostLabel);
700: }
701: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
702: for (; c < leafLocal[l] && c < cEnd; ++c) DMLabelSetValue(vtkLabel, c, 1);
703: if (leafLocal[l] >= cEnd) break;
704: if (leafRemote[l].rank == rank) {
705: DMLabelSetValue(vtkLabel, c, 1);
706: } else if (ghostLabel) DMLabelSetValue(ghostLabel, c, 2);
707: }
708: for (; c < cEnd; ++c) DMLabelSetValue(vtkLabel, c, 1);
709: if (ghostLabel) {
710: DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);
711: for (f = fStart; f < fEnd; ++f) {
712: PetscInt numCells;
714: DMPlexGetSupportSize(dmNew, f, &numCells);
715: if (numCells < 2) {
716: DMLabelSetValue(ghostLabel, f, 1);
717: } else {
718: const PetscInt *cells = NULL;
719: PetscInt vA, vB;
721: DMPlexGetSupport(dmNew, f, &cells);
722: DMLabelGetValue(vtkLabel, cells[0], &vA);
723: DMLabelGetValue(vtkLabel, cells[1], &vB);
724: if (vA != 1 && vB != 1) DMLabelSetValue(ghostLabel, f, 1);
725: }
726: }
727: }
728: return 0;
729: }
731: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
732: {
733: DM refTree;
734: PetscSection pSec;
735: PetscInt *parents, *childIDs;
737: DMPlexGetReferenceTree(dm, &refTree);
738: DMPlexSetReferenceTree(dmNew, refTree);
739: DMPlexGetTree(dm, &pSec, &parents, &childIDs, NULL, NULL);
740: if (pSec) {
741: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
742: PetscInt *childIDsShifted;
743: PetscSection pSecShifted;
745: PetscSectionGetChart(pSec, &pStart, &pEnd);
746: DMPlexGetDepth(dm, &depth);
747: pStartShifted = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
748: pEndShifted = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
749: PetscMalloc2(pEndShifted - pStartShifted, &parentsShifted, pEndShifted - pStartShifted, &childIDsShifted);
750: PetscSectionCreate(PetscObjectComm((PetscObject)dmNew), &pSecShifted);
751: PetscSectionSetChart(pSecShifted, pStartShifted, pEndShifted);
752: for (p = pStartShifted; p < pEndShifted; p++) {
753: /* start off assuming no children */
754: PetscSectionSetDof(pSecShifted, p, 0);
755: }
756: for (p = pStart; p < pEnd; p++) {
757: PetscInt dof;
758: PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);
760: PetscSectionGetDof(pSec, p, &dof);
761: PetscSectionSetDof(pSecShifted, pNew, dof);
762: }
763: PetscSectionSetUp(pSecShifted);
764: for (p = pStart; p < pEnd; p++) {
765: PetscInt dof;
766: PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);
768: PetscSectionGetDof(pSec, p, &dof);
769: if (dof) {
770: PetscInt off, offNew;
772: PetscSectionGetOffset(pSec, p, &off);
773: PetscSectionGetOffset(pSecShifted, pNew, &offNew);
774: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off], depth, depthShift);
775: childIDsShifted[offNew] = childIDs[off];
776: }
777: }
778: DMPlexSetTree(dmNew, pSecShifted, parentsShifted, childIDsShifted);
779: PetscFree2(parentsShifted, childIDsShifted);
780: PetscSectionDestroy(&pSecShifted);
781: }
782: return 0;
783: }
785: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
786: {
787: PetscSF sf;
788: IS valueIS;
789: const PetscInt *values, *leaves;
790: PetscInt *depthShift;
791: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
792: const PetscReal *maxCell, *Lstart, *L;
794: DMGetPointSF(dm, &sf);
795: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
796: nleaves = PetscMax(0, nleaves);
797: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
798: /* Count ghost cells */
799: DMLabelGetValueIS(label, &valueIS);
800: ISGetLocalSize(valueIS, &numFS);
801: ISGetIndices(valueIS, &values);
802: Ng = 0;
803: for (fs = 0; fs < numFS; ++fs) {
804: IS faceIS;
805: const PetscInt *faces;
806: PetscInt numFaces, f, numBdFaces = 0;
808: DMLabelGetStratumIS(label, values[fs], &faceIS);
809: ISGetLocalSize(faceIS, &numFaces);
810: ISGetIndices(faceIS, &faces);
811: for (f = 0; f < numFaces; ++f) {
812: PetscInt numChildren;
814: PetscFindInt(faces[f], nleaves, leaves, &loc);
815: DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
816: /* non-local and ancestors points don't get to register ghosts */
817: if (loc >= 0 || numChildren) continue;
818: if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces;
819: }
820: Ng += numBdFaces;
821: ISRestoreIndices(faceIS, &faces);
822: ISDestroy(&faceIS);
823: }
824: DMPlexGetDepth(dm, &depth);
825: PetscMalloc1(2 * (depth + 1), &depthShift);
826: for (d = 0; d <= depth; d++) {
827: PetscInt dEnd;
829: DMPlexGetDepthStratum(dm, d, NULL, &dEnd);
830: depthShift[2 * d] = dEnd;
831: depthShift[2 * d + 1] = 0;
832: }
833: if (depth >= 0) depthShift[2 * depth + 1] = Ng;
834: DMPlexShiftPointSetUp_Internal(depth, depthShift);
835: DMPlexShiftSizes_Internal(dm, depthShift, gdm);
836: /* Step 3: Set cone/support sizes for new points */
837: DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
838: for (c = cEnd; c < cEnd + Ng; ++c) DMPlexSetConeSize(gdm, c, 1);
839: for (fs = 0; fs < numFS; ++fs) {
840: IS faceIS;
841: const PetscInt *faces;
842: PetscInt numFaces, f;
844: DMLabelGetStratumIS(label, values[fs], &faceIS);
845: ISGetLocalSize(faceIS, &numFaces);
846: ISGetIndices(faceIS, &faces);
847: for (f = 0; f < numFaces; ++f) {
848: PetscInt size, numChildren;
850: PetscFindInt(faces[f], nleaves, leaves, &loc);
851: DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
852: if (loc >= 0 || numChildren) continue;
853: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
854: DMPlexGetSupportSize(dm, faces[f], &size);
856: DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);
857: }
858: ISRestoreIndices(faceIS, &faces);
859: ISDestroy(&faceIS);
860: }
861: /* Step 4: Setup ghosted DM */
862: DMSetUp(gdm);
863: DMPlexShiftPoints_Internal(dm, depthShift, gdm);
864: /* Step 6: Set cones and supports for new points */
865: ghostCell = cEnd;
866: for (fs = 0; fs < numFS; ++fs) {
867: IS faceIS;
868: const PetscInt *faces;
869: PetscInt numFaces, f;
871: DMLabelGetStratumIS(label, values[fs], &faceIS);
872: ISGetLocalSize(faceIS, &numFaces);
873: ISGetIndices(faceIS, &faces);
874: for (f = 0; f < numFaces; ++f) {
875: PetscInt newFace = faces[f] + Ng, numChildren;
877: PetscFindInt(faces[f], nleaves, leaves, &loc);
878: DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL);
879: if (loc >= 0 || numChildren) continue;
880: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
881: DMPlexSetCone(gdm, ghostCell, &newFace);
882: DMPlexInsertSupport(gdm, newFace, 1, ghostCell);
883: ++ghostCell;
884: }
885: ISRestoreIndices(faceIS, &faces);
886: ISDestroy(&faceIS);
887: }
888: ISRestoreIndices(valueIS, &values);
889: ISDestroy(&valueIS);
890: DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);
891: DMPlexShiftSF_Internal(dm, depthShift, gdm);
892: DMPlexShiftLabels_Internal(dm, depthShift, gdm);
893: DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm);
894: DMPlexShiftTree_Internal(dm, depthShift, gdm);
895: PetscFree(depthShift);
896: for (c = cEnd; c < cEnd + Ng; ++c) DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST);
897: /* Step 7: Periodicity */
898: DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
899: DMSetPeriodicity(gdm, maxCell, Lstart, L);
900: if (numGhostCells) *numGhostCells = Ng;
901: return 0;
902: }
904: /*@C
905: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
907: Collective on dm
909: Input Parameters:
910: + dm - The original DM
911: - labelName - The label specifying the boundary faces, or "Face Sets" if this is NULL
913: Output Parameters:
914: + numGhostCells - The number of ghost cells added to the DM
915: - dmGhosted - The new DM
917: Note: If no label exists of that name, one will be created marking all boundary faces
919: Level: developer
921: .seealso: `DMCreate()`
922: @*/
923: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
924: {
925: DM gdm;
926: DMLabel label;
927: const char *name = labelName ? labelName : "Face Sets";
928: PetscInt dim, Ng = 0;
929: PetscBool useCone, useClosure;
934: DMCreate(PetscObjectComm((PetscObject)dm), &gdm);
935: DMSetType(gdm, DMPLEX);
936: DMGetDimension(dm, &dim);
937: DMSetDimension(gdm, dim);
938: DMGetBasicAdjacency(dm, &useCone, &useClosure);
939: DMSetBasicAdjacency(gdm, useCone, useClosure);
940: DMGetLabel(dm, name, &label);
941: if (!label) {
942: /* Get label for boundary faces */
943: DMCreateLabel(dm, name);
944: DMGetLabel(dm, name, &label);
945: DMPlexMarkBoundaryFaces(dm, 1, label);
946: }
947: DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm);
948: DMCopyDisc(dm, gdm);
949: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, gdm);
950: gdm->setfromoptionscalled = dm->setfromoptionscalled;
951: if (numGhostCells) *numGhostCells = Ng;
952: *dmGhosted = gdm;
953: return 0;
954: }
956: static PetscErrorCode DivideCells_Private(DM dm, DMLabel label, DMPlexPointQueue queue)
957: {
958: PetscInt dim, d, shift = 100, *pStart, *pEnd;
960: DMGetDimension(dm, &dim);
961: PetscMalloc2(dim, &pStart, dim, &pEnd);
962: for (d = 0; d < dim; ++d) DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
963: while (!DMPlexPointQueueEmpty(queue)) {
964: PetscInt cell = -1;
965: PetscInt *closure = NULL;
966: PetscInt closureSize, cl, cval;
968: DMPlexPointQueueDequeue(queue, &cell);
969: DMLabelGetValue(label, cell, &cval);
970: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
971: /* Mark points in the cell closure that touch the fault */
972: for (d = 0; d < dim; ++d) {
973: for (cl = 0; cl < closureSize * 2; cl += 2) {
974: const PetscInt clp = closure[cl];
975: PetscInt clval;
977: if ((clp < pStart[d]) || (clp >= pEnd[d])) continue;
978: DMLabelGetValue(label, clp, &clval);
979: if (clval == -1) {
980: const PetscInt *cone;
981: PetscInt coneSize, c;
983: /* If a cone point touches the fault, then this point touches the fault */
984: DMPlexGetCone(dm, clp, &cone);
985: DMPlexGetConeSize(dm, clp, &coneSize);
986: for (c = 0; c < coneSize; ++c) {
987: PetscInt cpval;
989: DMLabelGetValue(label, cone[c], &cpval);
990: if (cpval != -1) {
991: PetscInt dep;
993: DMPlexGetPointDepth(dm, clp, &dep);
994: clval = cval < 0 ? -(shift + dep) : shift + dep;
995: DMLabelSetValue(label, clp, clval);
996: break;
997: }
998: }
999: }
1000: /* Mark neighbor cells through marked faces (these cells must also touch the fault) */
1001: if (d == dim - 1 && clval != -1) {
1002: const PetscInt *support;
1003: PetscInt supportSize, s, nval;
1005: DMPlexGetSupport(dm, clp, &support);
1006: DMPlexGetSupportSize(dm, clp, &supportSize);
1007: for (s = 0; s < supportSize; ++s) {
1008: DMLabelGetValue(label, support[s], &nval);
1009: if (nval == -1) {
1010: DMLabelSetValue(label, support[s], clval < 0 ? clval - 1 : clval + 1);
1011: DMPlexPointQueueEnqueue(queue, support[s]);
1012: }
1013: }
1014: }
1015: }
1016: }
1017: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1018: }
1019: PetscFree2(pStart, pEnd);
1020: return 0;
1021: }
1023: typedef struct {
1024: DM dm;
1025: DMPlexPointQueue queue;
1026: } PointDivision;
1028: static PetscErrorCode divideCell(DMLabel label, PetscInt p, PetscInt val, void *ctx)
1029: {
1030: PointDivision *div = (PointDivision *)ctx;
1031: PetscInt cval = val < 0 ? val - 1 : val + 1;
1032: const PetscInt *support;
1033: PetscInt supportSize, s;
1035: DMPlexGetSupport(div->dm, p, &support);
1036: DMPlexGetSupportSize(div->dm, p, &supportSize);
1037: for (s = 0; s < supportSize; ++s) {
1038: DMLabelSetValue(label, support[s], cval);
1039: DMPlexPointQueueEnqueue(div->queue, support[s]);
1040: }
1041: return 0;
1042: }
1044: /* Mark cells by label propagation */
1045: static PetscErrorCode DMPlexLabelFaultHalo(DM dm, DMLabel faultLabel)
1046: {
1047: DMPlexPointQueue queue = NULL;
1048: PointDivision div;
1049: PetscSF pointSF;
1050: IS pointIS;
1051: const PetscInt *points;
1052: PetscBool empty;
1053: PetscInt dim, shift = 100, n, i;
1055: DMGetDimension(dm, &dim);
1056: DMPlexPointQueueCreate(1024, &queue);
1057: div.dm = dm;
1058: div.queue = queue;
1059: /* Enqueue cells on fault */
1060: DMLabelGetStratumIS(faultLabel, shift + dim, &pointIS);
1061: if (pointIS) {
1062: ISGetLocalSize(pointIS, &n);
1063: ISGetIndices(pointIS, &points);
1064: for (i = 0; i < n; ++i) DMPlexPointQueueEnqueue(queue, points[i]);
1065: ISRestoreIndices(pointIS, &points);
1066: ISDestroy(&pointIS);
1067: }
1068: DMLabelGetStratumIS(faultLabel, -(shift + dim), &pointIS);
1069: if (pointIS) {
1070: ISGetLocalSize(pointIS, &n);
1071: ISGetIndices(pointIS, &points);
1072: for (i = 0; i < n; ++i) DMPlexPointQueueEnqueue(queue, points[i]);
1073: ISRestoreIndices(pointIS, &points);
1074: ISDestroy(&pointIS);
1075: }
1077: DMGetPointSF(dm, &pointSF);
1078: DMLabelPropagateBegin(faultLabel, pointSF);
1079: /* While edge queue is not empty: */
1080: DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
1081: while (!empty) {
1082: DivideCells_Private(dm, faultLabel, queue);
1083: DMLabelPropagatePush(faultLabel, pointSF, divideCell, &div);
1084: DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty);
1085: }
1086: DMLabelPropagateEnd(faultLabel, pointSF);
1087: DMPlexPointQueueDestroy(&queue);
1088: return 0;
1089: }
1091: /*
1092: We are adding three kinds of points here:
1093: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
1094: Non-replicated: Points which exist on the fault, but are not replicated
1095: Ghost: These are shared fault faces which are not owned by this process. These do not produce hybrid cells and do not replicate
1096: Hybrid: Entirely new points, such as cohesive cells
1098: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
1099: each depth so that the new split/hybrid points can be inserted as a block.
1100: */
1101: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
1102: {
1103: MPI_Comm comm;
1104: IS valueIS;
1105: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
1106: const PetscInt *values; /* List of depths for which we have replicated points */
1107: IS *splitIS;
1108: IS *unsplitIS;
1109: IS ghostIS;
1110: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
1111: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
1112: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
1113: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
1114: PetscInt numGhostPoints; /* The number of unowned, shared fault faces */
1115: const PetscInt **splitPoints; /* Replicated points for each depth */
1116: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
1117: const PetscInt *ghostPoints; /* Ghost fault faces */
1118: PetscSection coordSection;
1119: Vec coordinates;
1120: PetscScalar *coords;
1121: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
1122: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
1123: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
1124: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
1125: PetscInt *coneNew, *coneONew, *supportNew;
1126: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
1128: PetscObjectGetComm((PetscObject)dm, &comm);
1129: DMGetDimension(dm, &dim);
1130: DMPlexGetDepth(dm, &depth);
1131: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1132: /* We do not want this label automatically computed, instead we compute it here */
1133: DMCreateLabel(sdm, "celltype");
1134: /* Count split points and add cohesive cells */
1135: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
1136: PetscMalloc5(depth + 1, &depthMax, depth + 1, &depthEnd, 2 * (depth + 1), &depthShift, depth + 1, &pMaxNew, depth + 1, &numHybridPointsOld);
1137: PetscMalloc7(depth + 1, &splitIS, depth + 1, &unsplitIS, depth + 1, &numSplitPoints, depth + 1, &numUnsplitPoints, depth + 1, &numHybridPoints, depth + 1, &splitPoints, depth + 1, &unsplitPoints);
1138: for (d = 0; d <= depth; ++d) {
1139: DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);
1140: DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL);
1141: depthEnd[d] = pMaxNew[d];
1142: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1143: numSplitPoints[d] = 0;
1144: numUnsplitPoints[d] = 0;
1145: numHybridPoints[d] = 0;
1146: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1147: splitPoints[d] = NULL;
1148: unsplitPoints[d] = NULL;
1149: splitIS[d] = NULL;
1150: unsplitIS[d] = NULL;
1151: /* we are shifting the existing hybrid points with the stratum behind them, so
1152: * the split comes at the end of the normal points, i.e., at depthMax[d] */
1153: depthShift[2 * d] = depthMax[d];
1154: depthShift[2 * d + 1] = 0;
1155: }
1156: numGhostPoints = 0;
1157: ghostPoints = NULL;
1158: if (label) {
1159: DMLabelGetValueIS(label, &valueIS);
1160: ISGetLocalSize(valueIS, &numSP);
1161: ISGetIndices(valueIS, &values);
1162: }
1163: for (sp = 0; sp < numSP; ++sp) {
1164: const PetscInt dep = values[sp];
1166: if ((dep < 0) || (dep > depth)) continue;
1167: DMLabelGetStratumIS(label, dep, &splitIS[dep]);
1168: if (splitIS[dep]) {
1169: ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]);
1170: ISGetIndices(splitIS[dep], &splitPoints[dep]);
1171: }
1172: DMLabelGetStratumIS(label, shift2 + dep, &unsplitIS[dep]);
1173: if (unsplitIS[dep]) {
1174: ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]);
1175: ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]);
1176: }
1177: }
1178: DMLabelGetStratumIS(label, shift2 + dim - 1, &ghostIS);
1179: if (ghostIS) {
1180: ISGetLocalSize(ghostIS, &numGhostPoints);
1181: ISGetIndices(ghostIS, &ghostPoints);
1182: }
1183: /* Calculate number of hybrid points */
1184: 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 */
1185: for (d = 0; d <= depth; ++d) depthShift[2 * d + 1] = numSplitPoints[d] + numHybridPoints[d];
1186: DMPlexShiftPointSetUp_Internal(depth, depthShift);
1187: /* the end of the points in this stratum that come before the new points:
1188: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1189: * added points */
1190: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d], depth, depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1191: DMPlexShiftSizes_Internal(dm, depthShift, sdm);
1192: /* Step 3: Set cone/support sizes for new points */
1193: for (dep = 0; dep <= depth; ++dep) {
1194: for (p = 0; p < numSplitPoints[dep]; ++p) {
1195: const PetscInt oldp = splitPoints[dep][p];
1196: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1197: const PetscInt splitp = p + pMaxNew[dep];
1198: const PetscInt *support;
1199: DMPolytopeType ct;
1200: PetscInt coneSize, supportSize, qf, qn, qp, e;
1202: DMPlexGetConeSize(dm, oldp, &coneSize);
1203: DMPlexSetConeSize(sdm, splitp, coneSize);
1204: DMPlexGetSupportSize(dm, oldp, &supportSize);
1205: DMPlexSetSupportSize(sdm, splitp, supportSize);
1206: DMPlexGetCellType(dm, oldp, &ct);
1207: DMPlexSetCellType(sdm, splitp, ct);
1208: if (dep == depth - 1) {
1209: const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1211: /* Add cohesive cells, they are prisms */
1212: DMPlexSetConeSize(sdm, hybcell, 2 + coneSize);
1213: switch (coneSize) {
1214: case 2:
1215: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR);
1216: break;
1217: case 3:
1218: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR);
1219: break;
1220: case 4:
1221: DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR);
1222: break;
1223: }
1224: /* Shared fault faces with only one support cell now have two with the cohesive cell */
1225: /* TODO Check thaat oldp has rootdegree == 1 */
1226: if (supportSize == 1) {
1227: const PetscInt *support;
1228: PetscInt val;
1230: DMPlexGetSupport(dm, oldp, &support);
1231: DMLabelGetValue(label, support[0], &val);
1232: if (val < 0) DMPlexSetSupportSize(sdm, splitp, 2);
1233: else DMPlexSetSupportSize(sdm, newp, 2);
1234: }
1235: } else if (dep == 0) {
1236: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1238: DMPlexGetSupport(dm, oldp, &support);
1239: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1240: PetscInt val;
1242: DMLabelGetValue(label, support[e], &val);
1243: if (val == 1) ++qf;
1244: if ((val == 1) || (val == (shift + 1))) ++qn;
1245: if ((val == 1) || (val == -(shift + 1))) ++qp;
1246: }
1247: /* Split old vertex: Edges into original vertex and new cohesive edge */
1248: DMPlexSetSupportSize(sdm, newp, qn + 1);
1249: /* Split new vertex: Edges into split vertex and new cohesive edge */
1250: DMPlexSetSupportSize(sdm, splitp, qp + 1);
1251: /* Add hybrid edge */
1252: DMPlexSetConeSize(sdm, hybedge, 2);
1253: DMPlexSetSupportSize(sdm, hybedge, qf);
1254: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1255: } else if (dep == dim - 2) {
1256: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1258: DMPlexGetSupport(dm, oldp, &support);
1259: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1260: PetscInt val;
1262: DMLabelGetValue(label, support[e], &val);
1263: if (val == dim - 1) ++qf;
1264: if ((val == dim - 1) || (val == (shift + dim - 1))) ++qn;
1265: if ((val == dim - 1) || (val == -(shift + dim - 1))) ++qp;
1266: }
1267: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1268: DMPlexSetSupportSize(sdm, newp, qn + 1);
1269: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1270: DMPlexSetSupportSize(sdm, splitp, qp + 1);
1271: /* Add hybrid face */
1272: DMPlexSetConeSize(sdm, hybface, 4);
1273: DMPlexSetSupportSize(sdm, hybface, qf);
1274: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1275: }
1276: }
1277: }
1278: for (dep = 0; dep <= depth; ++dep) {
1279: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1280: const PetscInt oldp = unsplitPoints[dep][p];
1281: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1282: const PetscInt *support;
1283: PetscInt coneSize, supportSize, qf, e, s;
1285: DMPlexGetConeSize(dm, oldp, &coneSize);
1286: DMPlexGetSupportSize(dm, oldp, &supportSize);
1287: DMPlexGetSupport(dm, oldp, &support);
1288: if (dep == 0) {
1289: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1291: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1292: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1293: PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e);
1294: if (e >= 0) ++qf;
1295: }
1296: DMPlexSetSupportSize(sdm, newp, qf + 2);
1297: /* Add hybrid edge */
1298: DMPlexSetConeSize(sdm, hybedge, 2);
1299: for (e = 0, qf = 0; e < supportSize; ++e) {
1300: PetscInt val;
1302: DMLabelGetValue(label, support[e], &val);
1303: /* Split and unsplit edges produce hybrid faces */
1304: if (val == 1) ++qf;
1305: if (val == (shift2 + 1)) ++qf;
1306: }
1307: DMPlexSetSupportSize(sdm, hybedge, qf);
1308: DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR);
1309: } else if (dep == dim - 2) {
1310: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1311: PetscInt val;
1313: for (e = 0, qf = 0; e < supportSize; ++e) {
1314: DMLabelGetValue(label, support[e], &val);
1315: if (val == dim - 1) qf += 2;
1316: else ++qf;
1317: }
1318: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1319: DMPlexSetSupportSize(sdm, newp, qf + 2);
1320: /* Add hybrid face */
1321: for (e = 0, qf = 0; e < supportSize; ++e) {
1322: DMLabelGetValue(label, support[e], &val);
1323: if (val == dim - 1) ++qf;
1324: }
1325: DMPlexSetConeSize(sdm, hybface, 4);
1326: DMPlexSetSupportSize(sdm, hybface, qf);
1327: DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR);
1328: }
1329: }
1330: }
1331: /* Step 4: Setup split DM */
1332: DMSetUp(sdm);
1333: DMPlexShiftPoints_Internal(dm, depthShift, sdm);
1334: DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew);
1335: PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneNew, PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneONew, PetscMax(maxSupportSize, maxSupportSizeNew), &supportNew);
1336: /* Step 6: Set cones and supports for new points */
1337: for (dep = 0; dep <= depth; ++dep) {
1338: for (p = 0; p < numSplitPoints[dep]; ++p) {
1339: const PetscInt oldp = splitPoints[dep][p];
1340: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1341: const PetscInt splitp = p + pMaxNew[dep];
1342: const PetscInt *cone, *support, *ornt;
1343: DMPolytopeType ct;
1344: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1346: DMPlexGetCellType(dm, oldp, &ct);
1347: DMPlexGetConeSize(dm, oldp, &coneSize);
1348: DMPlexGetCone(dm, oldp, &cone);
1349: DMPlexGetConeOrientation(dm, oldp, &ornt);
1350: DMPlexGetSupportSize(dm, oldp, &supportSize);
1351: DMPlexGetSupport(dm, oldp, &support);
1352: if (dep == depth - 1) {
1353: PetscBool hasUnsplit = PETSC_FALSE;
1354: const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1355: const PetscInt *supportF;
1357: coneONew[0] = coneONew[1] = -1000;
1358: /* Split face: copy in old face to new face to start */
1359: DMPlexGetSupport(sdm, newp, &supportF);
1360: DMPlexSetSupport(sdm, splitp, supportF);
1361: /* Split old face: old vertices/edges in cone so no change */
1362: /* Split new face: new vertices/edges in cone */
1363: for (q = 0; q < coneSize; ++q) {
1364: PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1365: if (v < 0) {
1366: PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1368: coneNew[2 + q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1369: hasUnsplit = PETSC_TRUE;
1370: } else {
1371: coneNew[2 + q] = v + pMaxNew[dep - 1];
1372: if (dep > 1) {
1373: const PetscInt *econe;
1374: PetscInt econeSize, r, vs, vu;
1376: DMPlexGetConeSize(dm, cone[q], &econeSize);
1377: DMPlexGetCone(dm, cone[q], &econe);
1378: for (r = 0; r < econeSize; ++r) {
1379: PetscFindInt(econe[r], numSplitPoints[dep - 2], splitPoints[dep - 2], &vs);
1380: PetscFindInt(econe[r], numUnsplitPoints[dep - 2], unsplitPoints[dep - 2], &vu);
1381: if (vs >= 0) continue;
1383: hasUnsplit = PETSC_TRUE;
1384: }
1385: }
1386: }
1387: }
1388: DMPlexSetCone(sdm, splitp, &coneNew[2]);
1389: DMPlexSetConeOrientation(sdm, splitp, ornt);
1390: /* Face support */
1391: PetscInt vals[2];
1393: DMLabelGetValue(label, support[0], &vals[0]);
1394: if (supportSize > 1) DMLabelGetValue(label, support[1], &vals[1]);
1395: else vals[1] = -vals[0];
1398: for (s = 0; s < 2; ++s) {
1399: if (s >= supportSize) {
1400: if (vals[s] < 0) {
1401: /* Ghost old face: Replace negative side cell with cohesive cell */
1402: DMPlexInsertSupport(sdm, newp, s, hybcell);
1403: } else {
1404: /* Ghost new face: Replace positive side cell with cohesive cell */
1405: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1406: }
1407: } else {
1408: if (vals[s] < 0) {
1409: /* Split old face: Replace negative side cell with cohesive cell */
1410: DMPlexInsertSupport(sdm, newp, s, hybcell);
1411: } else {
1412: /* Split new face: Replace positive side cell with cohesive cell */
1413: DMPlexInsertSupport(sdm, splitp, s, hybcell);
1414: }
1415: }
1416: }
1417: /* Get orientation for cohesive face using the positive side cell */
1418: {
1419: const PetscInt *ncone, *nconeO;
1420: PetscInt nconeSize, nc, ocell;
1421: PetscBool flip = PETSC_FALSE;
1423: if (supportSize > 1) {
1424: ocell = vals[0] < 0 ? support[1] : support[0];
1425: } else {
1426: ocell = support[0];
1427: flip = vals[0] < 0 ? PETSC_TRUE : PETSC_FALSE;
1428: }
1429: DMPlexGetConeSize(dm, ocell, &nconeSize);
1430: DMPlexGetCone(dm, ocell, &ncone);
1431: DMPlexGetConeOrientation(dm, ocell, &nconeO);
1432: for (nc = 0; nc < nconeSize; ++nc) {
1433: if (ncone[nc] == oldp) {
1434: coneONew[0] = flip ? -(nconeO[nc] + 1) : nconeO[nc];
1435: break;
1436: }
1437: }
1439: }
1440: /* Cohesive cell: Old and new split face, then new cohesive faces */
1441: {
1442: const PetscInt No = DMPolytopeTypeGetNumArrangments(ct) / 2;
1444: }
1445: const PetscInt *arr = DMPolytopeTypeGetArrangment(ct, coneONew[0]);
1447: coneNew[0] = newp; /* Extracted negative side orientation above */
1448: coneNew[1] = splitp;
1449: coneONew[1] = coneONew[0];
1450: for (q = 0; q < coneSize; ++q) {
1451: /* Hybrid faces must follow order from oriented end face */
1452: const PetscInt qa = arr[q * 2 + 0];
1453: const PetscInt qo = arr[q * 2 + 1];
1454: DMPolytopeType ft = dep == 2 ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1456: PetscFindInt(cone[qa], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1457: if (v < 0) {
1458: PetscFindInt(cone[qa], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1459: coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1460: } else {
1461: coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep];
1462: }
1463: coneONew[2 + q] = DMPolytopeTypeComposeOrientation(ft, qo, ornt[qa]);
1464: }
1465: DMPlexSetCone(sdm, hybcell, coneNew);
1466: DMPlexSetConeOrientation(sdm, hybcell, coneONew);
1467: /* Label the hybrid cells on the boundary of the split */
1468: if (hasUnsplit) DMLabelSetValue(label, -hybcell, dim);
1469: } else if (dep == 0) {
1470: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1472: /* Split old vertex: Edges in old split faces and new cohesive edge */
1473: for (e = 0, qn = 0; e < supportSize; ++e) {
1474: PetscInt val;
1476: DMLabelGetValue(label, support[e], &val);
1477: if ((val == 1) || (val == (shift + 1))) supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1478: }
1479: supportNew[qn] = hybedge;
1480: DMPlexSetSupport(sdm, newp, supportNew);
1481: /* Split new vertex: Edges in new split faces and new cohesive edge */
1482: for (e = 0, qp = 0; e < supportSize; ++e) {
1483: PetscInt val, edge;
1485: DMLabelGetValue(label, support[e], &val);
1486: if (val == 1) {
1487: PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1489: supportNew[qp++] = edge + pMaxNew[dep + 1];
1490: } else if (val == -(shift + 1)) {
1491: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1492: }
1493: }
1494: supportNew[qp] = hybedge;
1495: DMPlexSetSupport(sdm, splitp, supportNew);
1496: /* Hybrid edge: Old and new split vertex */
1497: coneNew[0] = newp;
1498: coneNew[1] = splitp;
1499: DMPlexSetCone(sdm, hybedge, coneNew);
1500: for (e = 0, qf = 0; e < supportSize; ++e) {
1501: PetscInt val, edge;
1503: DMLabelGetValue(label, support[e], &val);
1504: if (val == 1) {
1505: PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1507: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1508: }
1509: }
1510: DMPlexSetSupport(sdm, hybedge, supportNew);
1511: } else if (dep == dim - 2) {
1512: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1514: /* Split old edge: old vertices in cone so no change */
1515: /* Split new edge: new vertices in cone */
1516: for (q = 0; q < coneSize; ++q) {
1517: PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v);
1518: if (v < 0) {
1519: PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1521: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1522: } else {
1523: coneNew[q] = v + pMaxNew[dep - 1];
1524: }
1525: }
1526: DMPlexSetCone(sdm, splitp, coneNew);
1527: /* Split old edge: Faces in positive side cells and old split faces */
1528: for (e = 0, q = 0; e < supportSize; ++e) {
1529: PetscInt val;
1531: DMLabelGetValue(label, support[e], &val);
1532: if (val == dim - 1) {
1533: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1534: } else if (val == (shift + dim - 1)) {
1535: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1536: }
1537: }
1538: supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1539: DMPlexSetSupport(sdm, newp, supportNew);
1540: /* Split new edge: Faces in negative side cells and new split faces */
1541: for (e = 0, q = 0; e < supportSize; ++e) {
1542: PetscInt val, face;
1544: DMLabelGetValue(label, support[e], &val);
1545: if (val == dim - 1) {
1546: PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1548: supportNew[q++] = face + pMaxNew[dep + 1];
1549: } else if (val == -(shift + dim - 1)) {
1550: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1551: }
1552: }
1553: supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1554: DMPlexSetSupport(sdm, splitp, supportNew);
1555: /* Hybrid face */
1556: coneNew[0] = newp;
1557: coneNew[1] = splitp;
1558: for (v = 0; v < coneSize; ++v) {
1559: PetscInt vertex;
1560: PetscFindInt(cone[v], numSplitPoints[dep - 1], splitPoints[dep - 1], &vertex);
1561: if (vertex < 0) {
1562: PetscFindInt(cone[v], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &vertex);
1564: coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1565: } else {
1566: coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
1567: }
1568: }
1569: DMPlexSetCone(sdm, hybface, coneNew);
1570: for (e = 0, qf = 0; e < supportSize; ++e) {
1571: PetscInt val, face;
1573: DMLabelGetValue(label, support[e], &val);
1574: if (val == dim - 1) {
1575: PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1577: supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1578: }
1579: }
1580: DMPlexSetSupport(sdm, hybface, supportNew);
1581: }
1582: }
1583: }
1584: for (dep = 0; dep <= depth; ++dep) {
1585: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1586: const PetscInt oldp = unsplitPoints[dep][p];
1587: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1588: const PetscInt *cone, *support;
1589: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
1591: DMPlexGetConeSize(dm, oldp, &coneSize);
1592: DMPlexGetCone(dm, oldp, &cone);
1593: DMPlexGetSupportSize(dm, oldp, &supportSize);
1594: DMPlexGetSupport(dm, oldp, &support);
1595: if (dep == 0) {
1596: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1598: /* Unsplit vertex */
1599: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1600: for (s = 0, q = 0; s < supportSize; ++s) {
1601: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
1602: PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e);
1603: if (e >= 0) supportNew[q++] = e + pMaxNew[dep + 1];
1604: }
1605: supportNew[q++] = hybedge;
1606: supportNew[q++] = hybedge;
1608: DMPlexSetSupport(sdm, newp, supportNew);
1609: /* Hybrid edge */
1610: coneNew[0] = newp;
1611: coneNew[1] = newp;
1612: DMPlexSetCone(sdm, hybedge, coneNew);
1613: for (e = 0, qf = 0; e < supportSize; ++e) {
1614: PetscInt val, edge;
1616: DMLabelGetValue(label, support[e], &val);
1617: if (val == 1) {
1618: PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge);
1620: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1621: } else if (val == (shift2 + 1)) {
1622: PetscFindInt(support[e], numUnsplitPoints[dep + 1], unsplitPoints[dep + 1], &edge);
1624: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2] + numSplitPoints[dep + 1];
1625: }
1626: }
1627: DMPlexSetSupport(sdm, hybedge, supportNew);
1628: } else if (dep == dim - 2) {
1629: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1631: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
1632: for (f = 0, qf = 0; f < supportSize; ++f) {
1633: PetscInt val, face;
1635: DMLabelGetValue(label, support[f], &val);
1636: if (val == dim - 1) {
1637: PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1639: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1640: supportNew[qf++] = face + pMaxNew[dep + 1];
1641: } else {
1642: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
1643: }
1644: }
1645: supportNew[qf++] = hybface;
1646: supportNew[qf++] = hybface;
1647: DMPlexGetSupportSize(sdm, newp, &supportSizeNew);
1649: DMPlexSetSupport(sdm, newp, supportNew);
1650: /* Add hybrid face */
1651: coneNew[0] = newp;
1652: coneNew[1] = newp;
1653: PetscFindInt(cone[0], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1655: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1656: PetscFindInt(cone[1], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v);
1658: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1659: DMPlexSetCone(sdm, hybface, coneNew);
1660: for (f = 0, qf = 0; f < supportSize; ++f) {
1661: PetscInt val, face;
1663: DMLabelGetValue(label, support[f], &val);
1664: if (val == dim - 1) {
1665: PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face);
1666: supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1667: }
1668: }
1669: DMPlexGetSupportSize(sdm, hybface, &supportSizeNew);
1671: DMPlexSetSupport(sdm, hybface, supportNew);
1672: }
1673: }
1674: }
1675: /* Step 6b: Replace split points in negative side cones */
1676: for (sp = 0; sp < numSP; ++sp) {
1677: PetscInt dep = values[sp];
1678: IS pIS;
1679: PetscInt numPoints;
1680: const PetscInt *points;
1682: if (dep >= 0) continue;
1683: DMLabelGetStratumIS(label, dep, &pIS);
1684: if (!pIS) continue;
1685: dep = -dep - shift;
1686: ISGetLocalSize(pIS, &numPoints);
1687: ISGetIndices(pIS, &points);
1688: for (p = 0; p < numPoints; ++p) {
1689: const PetscInt oldp = points[p];
1690: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
1691: const PetscInt *cone;
1692: PetscInt coneSize, c;
1693: /* PetscBool replaced = PETSC_FALSE; */
1695: /* Negative edge: replace split vertex */
1696: /* Negative cell: replace split face */
1697: DMPlexGetConeSize(sdm, newp, &coneSize);
1698: DMPlexGetCone(sdm, newp, &cone);
1699: for (c = 0; c < coneSize; ++c) {
1700: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c], depth, depthShift);
1701: PetscInt csplitp, cp, val;
1703: DMLabelGetValue(label, coldp, &val);
1704: if (val == dep - 1) {
1705: PetscFindInt(coldp, numSplitPoints[dep - 1], splitPoints[dep - 1], &cp);
1707: csplitp = pMaxNew[dep - 1] + cp;
1708: DMPlexInsertCone(sdm, newp, c, csplitp);
1709: /* replaced = PETSC_TRUE; */
1710: }
1711: }
1712: /* Cells with only a vertex or edge on the submesh have no replacement */
1714: }
1715: ISRestoreIndices(pIS, &points);
1716: ISDestroy(&pIS);
1717: }
1718: /* Step 7: Coordinates */
1719: DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);
1720: DMGetCoordinateSection(sdm, &coordSection);
1721: DMGetCoordinatesLocal(sdm, &coordinates);
1722: VecGetArray(coordinates, &coords);
1723: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
1724: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
1725: const PetscInt splitp = pMaxNew[0] + v;
1726: PetscInt dof, off, soff, d;
1728: PetscSectionGetDof(coordSection, newp, &dof);
1729: PetscSectionGetOffset(coordSection, newp, &off);
1730: PetscSectionGetOffset(coordSection, splitp, &soff);
1731: for (d = 0; d < dof; ++d) coords[soff + d] = coords[off + d];
1732: }
1733: VecRestoreArray(coordinates, &coords);
1734: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
1735: DMPlexShiftSF_Internal(dm, depthShift, sdm);
1736: /* TODO We need to associate the ghost points with the correct replica */
1737: /* Step 9: Labels */
1738: DMPlexShiftLabels_Internal(dm, depthShift, sdm);
1739: DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm);
1740: DMGetNumLabels(sdm, &numLabels);
1741: for (dep = 0; dep <= depth; ++dep) {
1742: for (p = 0; p < numSplitPoints[dep]; ++p) {
1743: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
1744: const PetscInt splitp = pMaxNew[dep] + p;
1745: PetscInt l;
1747: if (splitLabel) {
1748: const PetscInt val = 100 + dep;
1750: DMLabelSetValue(splitLabel, newp, val);
1751: DMLabelSetValue(splitLabel, splitp, -val);
1752: }
1753: for (l = 0; l < numLabels; ++l) {
1754: DMLabel mlabel;
1755: const char *lname;
1756: PetscInt val;
1757: PetscBool isDepth;
1759: DMGetLabelName(sdm, l, &lname);
1760: PetscStrcmp(lname, "depth", &isDepth);
1761: if (isDepth) continue;
1762: DMGetLabel(sdm, lname, &mlabel);
1763: DMLabelGetValue(mlabel, newp, &val);
1764: if (val >= 0) DMLabelSetValue(mlabel, splitp, val);
1765: }
1766: }
1767: }
1768: for (sp = 0; sp < numSP; ++sp) {
1769: const PetscInt dep = values[sp];
1771: if ((dep < 0) || (dep > depth)) continue;
1772: if (splitIS[dep]) ISRestoreIndices(splitIS[dep], &splitPoints[dep]);
1773: ISDestroy(&splitIS[dep]);
1774: if (unsplitIS[dep]) ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);
1775: ISDestroy(&unsplitIS[dep]);
1776: }
1777: if (ghostIS) ISRestoreIndices(ghostIS, &ghostPoints);
1778: ISDestroy(&ghostIS);
1779: if (label) {
1780: ISRestoreIndices(valueIS, &values);
1781: ISDestroy(&valueIS);
1782: }
1783: for (d = 0; d <= depth; ++d) {
1784: DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);
1785: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
1786: }
1787: PetscFree3(coneNew, coneONew, supportNew);
1788: PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);
1789: PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);
1790: return 0;
1791: }
1793: /*@C
1794: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
1796: Collective on dm
1798: Input Parameters:
1799: + dm - The original DM
1800: - label - The label specifying the boundary faces (this could be auto-generated)
1802: Output Parameters:
1803: + splitLabel - The label containing the split points, or NULL if no output is desired
1804: - dmSplit - The new DM
1806: Level: developer
1808: .seealso: `DMCreate()`, `DMPlexLabelCohesiveComplete()`
1809: @*/
1810: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
1811: {
1812: DM sdm;
1813: PetscInt dim;
1817: DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
1818: DMSetType(sdm, DMPLEX);
1819: DMGetDimension(dm, &dim);
1820: DMSetDimension(sdm, dim);
1821: switch (dim) {
1822: case 2:
1823: case 3:
1824: DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm);
1825: break;
1826: default:
1827: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %" PetscInt_FMT, dim);
1828: }
1829: *dmSplit = sdm;
1830: return 0;
1831: }
1833: /* Returns the side of the surface for a given cell with a face on the surface */
1834: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
1835: {
1836: const PetscInt *cone, *ornt;
1837: PetscInt dim, coneSize, c;
1839: *pos = PETSC_TRUE;
1840: DMGetDimension(dm, &dim);
1841: DMPlexGetConeSize(dm, cell, &coneSize);
1842: DMPlexGetCone(dm, cell, &cone);
1843: DMPlexGetConeOrientation(dm, cell, &ornt);
1844: for (c = 0; c < coneSize; ++c) {
1845: if (cone[c] == face) {
1846: PetscInt o = ornt[c];
1848: if (subdm) {
1849: const PetscInt *subcone, *subornt;
1850: PetscInt subpoint, subface, subconeSize, sc;
1852: PetscFindInt(cell, numSubpoints, subpoints, &subpoint);
1853: PetscFindInt(face, numSubpoints, subpoints, &subface);
1854: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
1855: DMPlexGetCone(subdm, subpoint, &subcone);
1856: DMPlexGetConeOrientation(subdm, subpoint, &subornt);
1857: for (sc = 0; sc < subconeSize; ++sc) {
1858: if (subcone[sc] == subface) {
1859: o = subornt[0];
1860: break;
1861: }
1862: }
1864: }
1865: if (o >= 0) *pos = PETSC_TRUE;
1866: else *pos = PETSC_FALSE;
1867: break;
1868: }
1869: }
1871: return 0;
1872: }
1874: static PetscErrorCode CheckFaultEdge_Private(DM dm, DMLabel label)
1875: {
1876: IS facePosIS, faceNegIS, dimIS;
1877: const PetscInt *points;
1878: PetscInt dim, numPoints, p, shift = 100, shift2 = 200;
1880: DMGetDimension(dm, &dim);
1881: /* If any faces touching the fault divide cells on either side, split them */
1882: DMLabelGetStratumIS(label, shift + dim - 1, &facePosIS);
1883: DMLabelGetStratumIS(label, -(shift + dim - 1), &faceNegIS);
1884: if (!facePosIS || !faceNegIS) {
1885: ISDestroy(&facePosIS);
1886: ISDestroy(&faceNegIS);
1887: return 0;
1888: }
1889: ISExpand(facePosIS, faceNegIS, &dimIS);
1890: ISDestroy(&facePosIS);
1891: ISDestroy(&faceNegIS);
1892: ISGetLocalSize(dimIS, &numPoints);
1893: ISGetIndices(dimIS, &points);
1894: for (p = 0; p < numPoints; ++p) {
1895: const PetscInt point = points[p];
1896: const PetscInt *support;
1897: PetscInt supportSize, valA, valB;
1899: DMPlexGetSupportSize(dm, point, &supportSize);
1900: if (supportSize != 2) continue;
1901: DMPlexGetSupport(dm, point, &support);
1902: DMLabelGetValue(label, support[0], &valA);
1903: DMLabelGetValue(label, support[1], &valB);
1904: if ((valA == -1) || (valB == -1)) continue;
1905: if (valA * valB > 0) continue;
1906: /* Check that this face is not incident on only unsplit faces, meaning has at least one split face */
1907: {
1908: PetscInt *closure = NULL;
1909: PetscBool split = PETSC_FALSE;
1910: PetscInt closureSize, cl;
1912: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1913: for (cl = 0; cl < closureSize * 2; cl += 2) {
1914: DMLabelGetValue(label, closure[cl], &valA);
1915: if ((valA >= 0) && (valA <= dim)) {
1916: split = PETSC_TRUE;
1917: break;
1918: }
1919: }
1920: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1921: if (!split) continue;
1922: }
1923: /* Split the face */
1924: DMLabelGetValue(label, point, &valA);
1925: DMLabelClearValue(label, point, valA);
1926: DMLabelSetValue(label, point, dim - 1);
1927: /* Label its closure:
1928: unmarked: label as unsplit
1929: incident: relabel as split
1930: split: do nothing
1931: */
1932: {
1933: PetscInt *closure = NULL;
1934: PetscInt closureSize, cl, dep;
1936: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1937: for (cl = 0; cl < closureSize * 2; cl += 2) {
1938: DMLabelGetValue(label, closure[cl], &valA);
1939: if (valA == -1) { /* Mark as unsplit */
1940: DMPlexGetPointDepth(dm, closure[cl], &dep);
1941: DMLabelSetValue(label, closure[cl], shift2 + dep);
1942: } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) {
1943: DMPlexGetPointDepth(dm, closure[cl], &dep);
1944: DMLabelClearValue(label, closure[cl], valA);
1945: DMLabelSetValue(label, closure[cl], dep);
1946: }
1947: }
1948: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
1949: }
1950: }
1951: ISRestoreIndices(dimIS, &points);
1952: ISDestroy(&dimIS);
1953: return 0;
1954: }
1956: /*@
1957: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
1958: to complete the surface
1960: Input Parameters:
1961: + dm - The DM
1962: . label - A DMLabel marking the surface
1963: . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically
1964: . bvalue - Value of DMLabel marking the vertices on the boundary
1965: . flip - Flag to flip the submesh normal and replace points on the other side
1966: - subdm - The subDM associated with the label, or NULL
1968: Output Parameter:
1969: . label - A DMLabel marking all surface points
1971: Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
1973: Level: developer
1975: .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexLabelComplete()`
1976: @*/
1977: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscInt bvalue, PetscBool flip, DM subdm)
1978: {
1979: DMLabel depthLabel;
1980: IS dimIS, subpointIS = NULL;
1981: const PetscInt *points, *subpoints;
1982: const PetscInt rev = flip ? -1 : 1;
1983: PetscInt shift = 100, shift2 = 200, shift3 = 300, dim, depth, numPoints, numSubpoints, p, val;
1985: DMPlexGetDepth(dm, &depth);
1986: DMGetDimension(dm, &dim);
1987: DMPlexGetDepthLabel(dm, &depthLabel);
1988: if (subdm) {
1989: DMPlexGetSubpointIS(subdm, &subpointIS);
1990: if (subpointIS) {
1991: ISGetLocalSize(subpointIS, &numSubpoints);
1992: ISGetIndices(subpointIS, &subpoints);
1993: }
1994: }
1995: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
1996: DMLabelGetStratumIS(label, dim - 1, &dimIS);
1997: if (!dimIS) goto divide;
1998: ISGetLocalSize(dimIS, &numPoints);
1999: ISGetIndices(dimIS, &points);
2000: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
2001: const PetscInt *support;
2002: PetscInt supportSize, s;
2004: DMPlexGetSupportSize(dm, points[p], &supportSize);
2005: #if 0
2006: if (supportSize != 2) {
2007: const PetscInt *lp;
2008: PetscInt Nlp, pind;
2010: /* Check that for a cell with a single support face, that face is in the SF */
2011: /* THis check only works for the remote side. We would need root side information */
2012: PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL);
2013: PetscFindInt(points[p], Nlp, lp, &pind);
2015: }
2016: #endif
2017: DMPlexGetSupport(dm, points[p], &support);
2018: for (s = 0; s < supportSize; ++s) {
2019: const PetscInt *cone;
2020: PetscInt coneSize, c;
2021: PetscBool pos;
2023: GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);
2024: if (pos) DMLabelSetValue(label, support[s], rev * (shift + dim));
2025: else DMLabelSetValue(label, support[s], -rev * (shift + dim));
2026: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
2027: /* Put faces touching the fault in the label */
2028: DMPlexGetConeSize(dm, support[s], &coneSize);
2029: DMPlexGetCone(dm, support[s], &cone);
2030: for (c = 0; c < coneSize; ++c) {
2031: const PetscInt point = cone[c];
2033: DMLabelGetValue(label, point, &val);
2034: if (val == -1) {
2035: PetscInt *closure = NULL;
2036: PetscInt closureSize, cl;
2038: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2039: for (cl = 0; cl < closureSize * 2; cl += 2) {
2040: const PetscInt clp = closure[cl];
2041: PetscInt bval = -1;
2043: DMLabelGetValue(label, clp, &val);
2044: if (blabel) DMLabelGetValue(blabel, clp, &bval);
2045: if ((val >= 0) && (val < dim - 1) && (bval < 0)) {
2046: DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift + dim - 1 : -(shift + dim - 1));
2047: break;
2048: }
2049: }
2050: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2051: }
2052: }
2053: }
2054: }
2055: ISRestoreIndices(dimIS, &points);
2056: ISDestroy(&dimIS);
2057: /* Mark boundary points as unsplit */
2058: if (blabel) {
2059: IS bdIS;
2061: DMLabelGetStratumIS(blabel, bvalue, &bdIS);
2062: ISGetLocalSize(bdIS, &numPoints);
2063: ISGetIndices(bdIS, &points);
2064: for (p = 0; p < numPoints; ++p) {
2065: const PetscInt point = points[p];
2066: PetscInt val, bval;
2068: DMLabelGetValue(blabel, point, &bval);
2069: if (bval >= 0) {
2070: DMLabelGetValue(label, point, &val);
2071: if ((val < 0) || (val > dim)) {
2072: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
2073: DMLabelClearValue(blabel, point, bval);
2074: }
2075: }
2076: }
2077: for (p = 0; p < numPoints; ++p) {
2078: const PetscInt point = points[p];
2079: PetscInt val, bval;
2081: DMLabelGetValue(blabel, point, &bval);
2082: if (bval >= 0) {
2083: const PetscInt *cone, *support;
2084: PetscInt coneSize, supportSize, s, valA, valB, valE;
2086: /* Mark as unsplit */
2087: DMLabelGetValue(label, point, &val);
2089: DMLabelClearValue(label, point, val);
2090: DMLabelSetValue(label, point, shift2 + val);
2091: /* Check for cross-edge
2092: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
2093: if (val != 0) continue;
2094: DMPlexGetSupport(dm, point, &support);
2095: DMPlexGetSupportSize(dm, point, &supportSize);
2096: for (s = 0; s < supportSize; ++s) {
2097: DMPlexGetCone(dm, support[s], &cone);
2098: DMPlexGetConeSize(dm, support[s], &coneSize);
2100: DMLabelGetValue(blabel, cone[0], &valA);
2101: DMLabelGetValue(blabel, cone[1], &valB);
2102: DMLabelGetValue(blabel, support[s], &valE);
2103: if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) DMLabelSetValue(blabel, support[s], 2);
2104: }
2105: }
2106: }
2107: ISRestoreIndices(bdIS, &points);
2108: ISDestroy(&bdIS);
2109: }
2110: /* Mark ghost fault cells */
2111: {
2112: PetscSF sf;
2113: const PetscInt *leaves;
2114: PetscInt Nl, l;
2116: DMGetPointSF(dm, &sf);
2117: PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL);
2118: DMLabelGetStratumIS(label, dim - 1, &dimIS);
2119: if (!dimIS) goto divide;
2120: ISGetLocalSize(dimIS, &numPoints);
2121: ISGetIndices(dimIS, &points);
2122: if (Nl > 0) {
2123: for (p = 0; p < numPoints; ++p) {
2124: const PetscInt point = points[p];
2125: PetscInt val;
2127: PetscFindInt(point, Nl, leaves, &l);
2128: if (l >= 0) {
2129: PetscInt *closure = NULL;
2130: PetscInt closureSize, cl;
2132: DMLabelGetValue(label, point, &val);
2134: DMLabelClearValue(label, point, val);
2135: DMLabelSetValue(label, point, shift3 + val);
2136: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2137: for (cl = 2; cl < closureSize * 2; cl += 2) {
2138: const PetscInt clp = closure[cl];
2140: DMLabelGetValue(label, clp, &val);
2142: DMLabelClearValue(label, clp, val);
2143: DMLabelSetValue(label, clp, shift3 + val);
2144: }
2145: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2146: }
2147: }
2148: }
2149: ISRestoreIndices(dimIS, &points);
2150: ISDestroy(&dimIS);
2151: }
2152: divide:
2153: if (subpointIS) ISRestoreIndices(subpointIS, &subpoints);
2154: DMPlexLabelFaultHalo(dm, label);
2155: CheckFaultEdge_Private(dm, label);
2156: return 0;
2157: }
2159: /* Check that no cell have all vertices on the fault */
2160: PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2161: {
2162: IS subpointIS;
2163: const PetscInt *dmpoints;
2164: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
2166: if (!label) return 0;
2167: DMLabelGetDefaultValue(label, &defaultValue);
2168: DMPlexGetSubpointIS(subdm, &subpointIS);
2169: if (!subpointIS) return 0;
2170: DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);
2171: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2172: ISGetIndices(subpointIS, &dmpoints);
2173: for (c = cStart; c < cEnd; ++c) {
2174: PetscBool invalidCell = PETSC_TRUE;
2175: PetscInt *closure = NULL;
2176: PetscInt closureSize, cl;
2178: DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2179: for (cl = 0; cl < closureSize * 2; cl += 2) {
2180: PetscInt value = 0;
2182: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2183: DMLabelGetValue(label, closure[cl], &value);
2184: if (value == defaultValue) {
2185: invalidCell = PETSC_FALSE;
2186: break;
2187: }
2188: }
2189: DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);
2190: if (invalidCell) {
2191: ISRestoreIndices(subpointIS, &dmpoints);
2192: ISDestroy(&subpointIS);
2193: DMDestroy(&subdm);
2194: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %" PetscInt_FMT " has all of its vertices on the submesh.", dmpoints[c]);
2195: }
2196: }
2197: ISRestoreIndices(subpointIS, &dmpoints);
2198: return 0;
2199: }
2201: /*@
2202: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
2204: Collective on dm
2206: Input Parameters:
2207: + dm - The original DM
2208: . label - The label specifying the interface vertices
2209: . bdlabel - The optional label specifying the interface boundary vertices
2210: - bdvalue - Value of optional label specifying the interface boundary vertices
2212: Output Parameters:
2213: + hybridLabel - The label fully marking the interface, or NULL if no output is desired
2214: . splitLabel - The label containing the split points, or NULL if no output is desired
2215: . dmInterface - The new interface DM, or NULL
2216: - dmHybrid - The new DM with cohesive cells
2218: Note: The hybridLabel indicates what parts of the original mesh impinged on the division surface. For points
2219: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2220: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2221: one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the
2222: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
2224: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2225: mesh. The label value is $\pm 100+dim$ for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2226: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
2228: The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using
2229: DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh.
2231: Level: developer
2233: .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexLabelCohesiveComplete()`, `DMPlexGetSubpointMap()`, `DMCreate()`
2234: @*/
2235: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, PetscInt bdvalue, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2236: {
2237: DM idm;
2238: DMLabel subpointMap, hlabel, slabel = NULL;
2239: PetscInt dim;
2248: DMGetDimension(dm, &dim);
2249: DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);
2250: DMPlexCheckValidSubmesh_Private(dm, label, idm);
2251: DMPlexOrient(idm);
2252: DMPlexGetSubpointMap(idm, &subpointMap);
2253: DMLabelDuplicate(subpointMap, &hlabel);
2254: DMLabelClearStratum(hlabel, dim);
2255: if (splitLabel) {
2256: const char *name;
2257: char sname[PETSC_MAX_PATH_LEN];
2259: PetscObjectGetName((PetscObject)hlabel, &name);
2260: PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);
2261: PetscStrcat(sname, " split");
2262: DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);
2263: }
2264: DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, bdvalue, PETSC_FALSE, idm);
2265: if (dmInterface) {
2266: *dmInterface = idm;
2267: } else DMDestroy(&idm);
2268: DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);
2269: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmHybrid);
2270: if (hybridLabel) *hybridLabel = hlabel;
2271: else DMLabelDestroy(&hlabel);
2272: if (splitLabel) *splitLabel = slabel;
2273: {
2274: DM cdm;
2275: DMLabel ctLabel;
2277: /* We need to somehow share the celltype label with the coordinate dm */
2278: DMGetCoordinateDM(*dmHybrid, &cdm);
2279: DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel);
2280: DMSetLabel(cdm, ctLabel);
2281: }
2282: return 0;
2283: }
2285: /* Here we need the explicit assumption that:
2287: For any marked cell, the marked vertices constitute a single face
2288: */
2289: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2290: {
2291: IS subvertexIS = NULL;
2292: const PetscInt *subvertices;
2293: PetscInt *pStart, *pEnd, pSize;
2294: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2296: *numFaces = 0;
2297: *nFV = 0;
2298: DMPlexGetDepth(dm, &depth);
2299: DMGetDimension(dm, &dim);
2300: pSize = PetscMax(depth, dim) + 1;
2301: PetscMalloc2(pSize, &pStart, pSize, &pEnd);
2302: for (d = 0; d <= depth; ++d) DMPlexGetSimplexOrBoxCells(dm, depth - d, &pStart[d], &pEnd[d]);
2303: /* Loop over initial vertices and mark all faces in the collective star() */
2304: if (vertexLabel) DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2305: if (subvertexIS) {
2306: ISGetSize(subvertexIS, &numSubVerticesInitial);
2307: ISGetIndices(subvertexIS, &subvertices);
2308: }
2309: for (v = 0; v < numSubVerticesInitial; ++v) {
2310: const PetscInt vertex = subvertices[v];
2311: PetscInt *star = NULL;
2312: PetscInt starSize, s, numCells = 0, c;
2314: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2315: for (s = 0; s < starSize * 2; s += 2) {
2316: const PetscInt point = star[s];
2317: if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point;
2318: }
2319: for (c = 0; c < numCells; ++c) {
2320: const PetscInt cell = star[c];
2321: PetscInt *closure = NULL;
2322: PetscInt closureSize, cl;
2323: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2325: DMLabelGetValue(subpointMap, cell, &cellLoc);
2326: if (cellLoc == 2) continue;
2328: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2329: for (cl = 0; cl < closureSize * 2; cl += 2) {
2330: const PetscInt point = closure[cl];
2331: PetscInt vertexLoc;
2333: if ((point >= pStart[0]) && (point < pEnd[0])) {
2334: ++numCorners;
2335: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2336: if (vertexLoc == value) closure[faceSize++] = point;
2337: }
2338: }
2339: if (!(*nFV)) DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);
2341: if (faceSize == *nFV) {
2342: const PetscInt *cells = NULL;
2343: PetscInt numCells, nc;
2345: ++(*numFaces);
2346: for (cl = 0; cl < faceSize; ++cl) DMLabelSetValue(subpointMap, closure[cl], 0);
2347: DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);
2348: for (nc = 0; nc < numCells; ++nc) DMLabelSetValue(subpointMap, cells[nc], 2);
2349: DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);
2350: }
2351: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
2352: }
2353: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2354: }
2355: if (subvertexIS) ISRestoreIndices(subvertexIS, &subvertices);
2356: ISDestroy(&subvertexIS);
2357: PetscFree2(pStart, pEnd);
2358: return 0;
2359: }
2361: static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm)
2362: {
2363: IS subvertexIS = NULL;
2364: const PetscInt *subvertices;
2365: PetscInt *pStart, *pEnd;
2366: PetscInt dim, d, numSubVerticesInitial = 0, v;
2368: DMGetDimension(dm, &dim);
2369: PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd);
2370: for (d = 0; d <= dim; ++d) DMPlexGetSimplexOrBoxCells(dm, dim - d, &pStart[d], &pEnd[d]);
2371: /* Loop over initial vertices and mark all faces in the collective star() */
2372: if (vertexLabel) {
2373: DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);
2374: if (subvertexIS) {
2375: ISGetSize(subvertexIS, &numSubVerticesInitial);
2376: ISGetIndices(subvertexIS, &subvertices);
2377: }
2378: }
2379: for (v = 0; v < numSubVerticesInitial; ++v) {
2380: const PetscInt vertex = subvertices[v];
2381: PetscInt *star = NULL;
2382: PetscInt starSize, s, numFaces = 0, f;
2384: DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2385: for (s = 0; s < starSize * 2; s += 2) {
2386: const PetscInt point = star[s];
2387: PetscInt faceLoc;
2389: if ((point >= pStart[dim - 1]) && (point < pEnd[dim - 1])) {
2390: if (markedFaces) {
2391: DMLabelGetValue(vertexLabel, point, &faceLoc);
2392: if (faceLoc < 0) continue;
2393: }
2394: star[numFaces++] = point;
2395: }
2396: }
2397: for (f = 0; f < numFaces; ++f) {
2398: const PetscInt face = star[f];
2399: PetscInt *closure = NULL;
2400: PetscInt closureSize, c;
2401: PetscInt faceLoc;
2403: DMLabelGetValue(subpointMap, face, &faceLoc);
2404: if (faceLoc == dim - 1) continue;
2406: DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2407: for (c = 0; c < closureSize * 2; c += 2) {
2408: const PetscInt point = closure[c];
2409: PetscInt vertexLoc;
2411: if ((point >= pStart[0]) && (point < pEnd[0])) {
2412: DMLabelGetValue(vertexLabel, point, &vertexLoc);
2413: if (vertexLoc != value) break;
2414: }
2415: }
2416: if (c == closureSize * 2) {
2417: const PetscInt *support;
2418: PetscInt supportSize, s;
2420: for (c = 0; c < closureSize * 2; c += 2) {
2421: const PetscInt point = closure[c];
2423: for (d = 0; d < dim; ++d) {
2424: if ((point >= pStart[d]) && (point < pEnd[d])) {
2425: DMLabelSetValue(subpointMap, point, d);
2426: break;
2427: }
2428: }
2429: }
2430: DMPlexGetSupportSize(dm, face, &supportSize);
2431: DMPlexGetSupport(dm, face, &support);
2432: for (s = 0; s < supportSize; ++s) DMLabelSetValue(subpointMap, support[s], dim);
2433: }
2434: DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);
2435: }
2436: DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);
2437: }
2438: if (subvertexIS) ISRestoreIndices(subvertexIS, &subvertices);
2439: ISDestroy(&subvertexIS);
2440: PetscFree2(pStart, pEnd);
2441: return 0;
2442: }
2444: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
2445: {
2446: DMLabel label = NULL;
2447: const PetscInt *cone;
2448: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
2450: *numFaces = 0;
2451: *nFV = 0;
2452: if (labelname) DMGetLabel(dm, labelname, &label);
2453: *subCells = NULL;
2454: DMGetDimension(dm, &dim);
2455: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2456: if (cMax < 0) return 0;
2457: if (label) {
2458: for (c = cMax; c < cEnd; ++c) {
2459: PetscInt val;
2461: DMLabelGetValue(label, c, &val);
2462: if (val == value) {
2463: ++(*numFaces);
2464: DMPlexGetConeSize(dm, c, &coneSize);
2465: }
2466: }
2467: } else {
2468: *numFaces = cEnd - cMax;
2469: DMPlexGetConeSize(dm, cMax, &coneSize);
2470: }
2471: PetscMalloc1(*numFaces * 2, subCells);
2472: if (!(*numFaces)) return 0;
2473: *nFV = hasLagrange ? coneSize / 3 : coneSize / 2;
2474: for (c = cMax; c < cEnd; ++c) {
2475: const PetscInt *cells;
2476: PetscInt numCells;
2478: if (label) {
2479: PetscInt val;
2481: DMLabelGetValue(label, c, &val);
2482: if (val != value) continue;
2483: }
2484: DMPlexGetCone(dm, c, &cone);
2485: for (p = 0; p < *nFV; ++p) DMLabelSetValue(subpointMap, cone[p], 0);
2486: /* Negative face */
2487: DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);
2488: /* Not true in parallel
2490: for (p = 0; p < numCells; ++p) {
2491: DMLabelSetValue(subpointMap, cells[p], 2);
2492: (*subCells)[subc++] = cells[p];
2493: }
2494: DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);
2495: /* Positive face is not included */
2496: }
2497: return 0;
2498: }
2500: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
2501: {
2502: PetscInt *pStart, *pEnd;
2503: PetscInt dim, cMax, cEnd, c, d;
2505: DMGetDimension(dm, &dim);
2506: DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);
2507: if (cMax < 0) return 0;
2508: PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd);
2509: for (d = 0; d <= dim; ++d) DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);
2510: for (c = cMax; c < cEnd; ++c) {
2511: const PetscInt *cone;
2512: PetscInt *closure = NULL;
2513: PetscInt fconeSize, coneSize, closureSize, cl, val;
2515: if (label) {
2516: DMLabelGetValue(label, c, &val);
2517: if (val != value) continue;
2518: }
2519: DMPlexGetConeSize(dm, c, &coneSize);
2520: DMPlexGetCone(dm, c, &cone);
2521: DMPlexGetConeSize(dm, cone[0], &fconeSize);
2523: /* Negative face */
2524: DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2525: for (cl = 0; cl < closureSize * 2; cl += 2) {
2526: const PetscInt point = closure[cl];
2528: for (d = 0; d <= dim; ++d) {
2529: if ((point >= pStart[d]) && (point < pEnd[d])) {
2530: DMLabelSetValue(subpointMap, point, d);
2531: break;
2532: }
2533: }
2534: }
2535: DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);
2536: /* Cells -- positive face is not included */
2537: for (cl = 0; cl < 1; ++cl) {
2538: const PetscInt *support;
2539: PetscInt supportSize, s;
2541: DMPlexGetSupportSize(dm, cone[cl], &supportSize);
2543: DMPlexGetSupport(dm, cone[cl], &support);
2544: for (s = 0; s < supportSize; ++s) DMLabelSetValue(subpointMap, support[s], dim);
2545: }
2546: }
2547: PetscFree2(pStart, pEnd);
2548: return 0;
2549: }
2551: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2552: {
2553: MPI_Comm comm;
2554: PetscBool posOrient = PETSC_FALSE;
2555: const PetscInt debug = 0;
2556: PetscInt cellDim, faceSize, f;
2558: PetscObjectGetComm((PetscObject)dm, &comm);
2559: DMGetDimension(dm, &cellDim);
2560: if (debug) PetscPrintf(comm, "cellDim: %" PetscInt_FMT " numCorners: %" PetscInt_FMT "\n", cellDim, numCorners);
2562: if (cellDim == 1 && numCorners == 2) {
2563: /* Triangle */
2564: faceSize = numCorners - 1;
2565: posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2566: } else if (cellDim == 2 && numCorners == 3) {
2567: /* Triangle */
2568: faceSize = numCorners - 1;
2569: posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2570: } else if (cellDim == 3 && numCorners == 4) {
2571: /* Tetrahedron */
2572: faceSize = numCorners - 1;
2573: posOrient = (oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
2574: } else if (cellDim == 1 && numCorners == 3) {
2575: /* Quadratic line */
2576: faceSize = 1;
2577: posOrient = PETSC_TRUE;
2578: } else if (cellDim == 2 && numCorners == 4) {
2579: /* Quads */
2580: faceSize = 2;
2581: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
2582: posOrient = PETSC_TRUE;
2583: } else if ((indices[0] == 3) && (indices[1] == 0)) {
2584: posOrient = PETSC_TRUE;
2585: } else {
2586: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
2587: posOrient = PETSC_FALSE;
2588: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
2589: }
2590: } else if (cellDim == 2 && numCorners == 6) {
2591: /* Quadratic triangle (I hate this) */
2592: /* Edges are determined by the first 2 vertices (corners of edges) */
2593: const PetscInt faceSizeTri = 3;
2594: PetscInt sortedIndices[3], i, iFace;
2595: PetscBool found = PETSC_FALSE;
2596: PetscInt faceVerticesTriSorted[9] = {
2597: 0, 3, 4, /* bottom */
2598: 1, 4, 5, /* right */
2599: 2, 3, 5, /* left */
2600: };
2601: PetscInt faceVerticesTri[9] = {
2602: 0, 3, 4, /* bottom */
2603: 1, 4, 5, /* right */
2604: 2, 5, 3, /* left */
2605: };
2607: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
2608: PetscSortInt(faceSizeTri, sortedIndices);
2609: for (iFace = 0; iFace < 3; ++iFace) {
2610: const PetscInt ii = iFace * faceSizeTri;
2611: PetscInt fVertex, cVertex;
2613: if ((sortedIndices[0] == faceVerticesTriSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTriSorted[ii + 1])) {
2614: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
2615: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
2616: if (indices[cVertex] == faceVerticesTri[ii + fVertex]) {
2617: faceVertices[fVertex] = origVertices[cVertex];
2618: break;
2619: }
2620: }
2621: }
2622: found = PETSC_TRUE;
2623: break;
2624: }
2625: }
2627: if (posOriented) *posOriented = PETSC_TRUE;
2628: return 0;
2629: } else if (cellDim == 2 && numCorners == 9) {
2630: /* Quadratic quad (I hate this) */
2631: /* Edges are determined by the first 2 vertices (corners of edges) */
2632: const PetscInt faceSizeQuad = 3;
2633: PetscInt sortedIndices[3], i, iFace;
2634: PetscBool found = PETSC_FALSE;
2635: PetscInt faceVerticesQuadSorted[12] = {
2636: 0, 1, 4, /* bottom */
2637: 1, 2, 5, /* right */
2638: 2, 3, 6, /* top */
2639: 0, 3, 7, /* left */
2640: };
2641: PetscInt faceVerticesQuad[12] = {
2642: 0, 1, 4, /* bottom */
2643: 1, 2, 5, /* right */
2644: 2, 3, 6, /* top */
2645: 3, 0, 7, /* left */
2646: };
2648: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
2649: PetscSortInt(faceSizeQuad, sortedIndices);
2650: for (iFace = 0; iFace < 4; ++iFace) {
2651: const PetscInt ii = iFace * faceSizeQuad;
2652: PetscInt fVertex, cVertex;
2654: if ((sortedIndices[0] == faceVerticesQuadSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadSorted[ii + 1])) {
2655: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
2656: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
2657: if (indices[cVertex] == faceVerticesQuad[ii + fVertex]) {
2658: faceVertices[fVertex] = origVertices[cVertex];
2659: break;
2660: }
2661: }
2662: }
2663: found = PETSC_TRUE;
2664: break;
2665: }
2666: }
2668: if (posOriented) *posOriented = PETSC_TRUE;
2669: return 0;
2670: } else if (cellDim == 3 && numCorners == 8) {
2671: /* Hexes
2672: A hex is two oriented quads with the normal of the first
2673: pointing up at the second.
2675: 7---6
2676: /| /|
2677: 4---5 |
2678: | 1-|-2
2679: |/ |/
2680: 0---3
2682: Faces are determined by the first 4 vertices (corners of faces) */
2683: const PetscInt faceSizeHex = 4;
2684: PetscInt sortedIndices[4], i, iFace;
2685: PetscBool found = PETSC_FALSE;
2686: PetscInt faceVerticesHexSorted[24] = {
2687: 0, 1, 2, 3, /* bottom */
2688: 4, 5, 6, 7, /* top */
2689: 0, 3, 4, 5, /* front */
2690: 2, 3, 5, 6, /* right */
2691: 1, 2, 6, 7, /* back */
2692: 0, 1, 4, 7, /* left */
2693: };
2694: PetscInt faceVerticesHex[24] = {
2695: 1, 2, 3, 0, /* bottom */
2696: 4, 5, 6, 7, /* top */
2697: 0, 3, 5, 4, /* front */
2698: 3, 2, 6, 5, /* right */
2699: 2, 1, 7, 6, /* back */
2700: 1, 0, 4, 7, /* left */
2701: };
2703: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
2704: PetscSortInt(faceSizeHex, sortedIndices);
2705: for (iFace = 0; iFace < 6; ++iFace) {
2706: const PetscInt ii = iFace * faceSizeHex;
2707: PetscInt fVertex, cVertex;
2709: if ((sortedIndices[0] == faceVerticesHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesHexSorted[ii + 3])) {
2710: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
2711: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
2712: if (indices[cVertex] == faceVerticesHex[ii + fVertex]) {
2713: faceVertices[fVertex] = origVertices[cVertex];
2714: break;
2715: }
2716: }
2717: }
2718: found = PETSC_TRUE;
2719: break;
2720: }
2721: }
2723: if (posOriented) *posOriented = PETSC_TRUE;
2724: return 0;
2725: } else if (cellDim == 3 && numCorners == 10) {
2726: /* Quadratic tet */
2727: /* Faces are determined by the first 3 vertices (corners of faces) */
2728: const PetscInt faceSizeTet = 6;
2729: PetscInt sortedIndices[6], i, iFace;
2730: PetscBool found = PETSC_FALSE;
2731: PetscInt faceVerticesTetSorted[24] = {
2732: 0, 1, 2, 6, 7, 8, /* bottom */
2733: 0, 3, 4, 6, 7, 9, /* front */
2734: 1, 4, 5, 7, 8, 9, /* right */
2735: 2, 3, 5, 6, 8, 9, /* left */
2736: };
2737: PetscInt faceVerticesTet[24] = {
2738: 0, 1, 2, 6, 7, 8, /* bottom */
2739: 0, 4, 3, 6, 7, 9, /* front */
2740: 1, 5, 4, 7, 8, 9, /* right */
2741: 2, 3, 5, 8, 6, 9, /* left */
2742: };
2744: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
2745: PetscSortInt(faceSizeTet, sortedIndices);
2746: for (iFace = 0; iFace < 4; ++iFace) {
2747: const PetscInt ii = iFace * faceSizeTet;
2748: PetscInt fVertex, cVertex;
2750: if ((sortedIndices[0] == faceVerticesTetSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTetSorted[ii + 1]) && (sortedIndices[2] == faceVerticesTetSorted[ii + 2]) && (sortedIndices[3] == faceVerticesTetSorted[ii + 3])) {
2751: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
2752: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
2753: if (indices[cVertex] == faceVerticesTet[ii + fVertex]) {
2754: faceVertices[fVertex] = origVertices[cVertex];
2755: break;
2756: }
2757: }
2758: }
2759: found = PETSC_TRUE;
2760: break;
2761: }
2762: }
2764: if (posOriented) *posOriented = PETSC_TRUE;
2765: return 0;
2766: } else if (cellDim == 3 && numCorners == 27) {
2767: /* Quadratic hexes (I hate this)
2768: A hex is two oriented quads with the normal of the first
2769: pointing up at the second.
2771: 7---6
2772: /| /|
2773: 4---5 |
2774: | 3-|-2
2775: |/ |/
2776: 0---1
2778: Faces are determined by the first 4 vertices (corners of faces) */
2779: const PetscInt faceSizeQuadHex = 9;
2780: PetscInt sortedIndices[9], i, iFace;
2781: PetscBool found = PETSC_FALSE;
2782: PetscInt faceVerticesQuadHexSorted[54] = {
2783: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
2784: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2785: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
2786: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
2787: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
2788: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
2789: };
2790: PetscInt faceVerticesQuadHex[54] = {
2791: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
2792: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
2793: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
2794: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
2795: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
2796: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
2797: };
2799: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
2800: PetscSortInt(faceSizeQuadHex, sortedIndices);
2801: for (iFace = 0; iFace < 6; ++iFace) {
2802: const PetscInt ii = iFace * faceSizeQuadHex;
2803: PetscInt fVertex, cVertex;
2805: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesQuadHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesQuadHexSorted[ii + 3])) {
2806: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
2807: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
2808: if (indices[cVertex] == faceVerticesQuadHex[ii + fVertex]) {
2809: faceVertices[fVertex] = origVertices[cVertex];
2810: break;
2811: }
2812: }
2813: }
2814: found = PETSC_TRUE;
2815: break;
2816: }
2817: }
2819: if (posOriented) *posOriented = PETSC_TRUE;
2820: return 0;
2821: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
2822: if (!posOrient) {
2823: if (debug) PetscPrintf(comm, " Reversing initial face orientation\n");
2824: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize - 1 - f];
2825: } else {
2826: if (debug) PetscPrintf(comm, " Keeping initial face orientation\n");
2827: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
2828: }
2829: if (posOriented) *posOriented = posOrient;
2830: return 0;
2831: }
2833: /*@
2834: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
2835: in faceVertices. The orientation is such that the face normal points out of the cell
2837: Not collective
2839: Input Parameters:
2840: + dm - The original mesh
2841: . cell - The cell mesh point
2842: . faceSize - The number of vertices on the face
2843: . face - The face vertices
2844: . numCorners - The number of vertices on the cell
2845: . indices - Local numbering of face vertices in cell cone
2846: - origVertices - Original face vertices
2848: Output Parameters:
2849: + faceVertices - The face vertices properly oriented
2850: - posOriented - PETSC_TRUE if the face was oriented with outward normal
2852: Level: developer
2854: .seealso: `DMPlexGetCone()`
2855: @*/
2856: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
2857: {
2858: const PetscInt *cone = NULL;
2859: PetscInt coneSize, v, f, v2;
2860: PetscInt oppositeVertex = -1;
2862: DMPlexGetConeSize(dm, cell, &coneSize);
2863: DMPlexGetCone(dm, cell, &cone);
2864: for (v = 0, v2 = 0; v < coneSize; ++v) {
2865: PetscBool found = PETSC_FALSE;
2867: for (f = 0; f < faceSize; ++f) {
2868: if (face[f] == cone[v]) {
2869: found = PETSC_TRUE;
2870: break;
2871: }
2872: }
2873: if (found) {
2874: indices[v2] = v;
2875: origVertices[v2] = cone[v];
2876: ++v2;
2877: } else {
2878: oppositeVertex = v;
2879: }
2880: }
2881: DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
2882: return 0;
2883: }
2885: /*
2886: DMPlexInsertFace_Internal - Puts a face into the mesh
2888: Not collective
2890: Input Parameters:
2891: + dm - The DMPlex
2892: . numFaceVertex - The number of vertices in the face
2893: . faceVertices - The vertices in the face for dm
2894: . subfaceVertices - The vertices in the face for subdm
2895: . numCorners - The number of vertices in the cell
2896: . cell - A cell in dm containing the face
2897: . subcell - A cell in subdm containing the face
2898: . firstFace - First face in the mesh
2899: - newFacePoint - Next face in the mesh
2901: Output Parameters:
2902: . newFacePoint - Contains next face point number on input, updated on output
2904: Level: developer
2905: */
2906: 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)
2907: {
2908: MPI_Comm comm;
2909: DM_Plex *submesh = (DM_Plex *)subdm->data;
2910: const PetscInt *faces;
2911: PetscInt numFaces, coneSize;
2913: PetscObjectGetComm((PetscObject)dm, &comm);
2914: DMPlexGetConeSize(subdm, subcell, &coneSize);
2916: #if 0
2917: /* Cannot use this because support() has not been constructed yet */
2918: DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2919: #else
2920: {
2921: PetscInt f;
2923: numFaces = 0;
2924: DMGetWorkArray(subdm, 1, MPIU_INT, (void **)&faces);
2925: for (f = firstFace; f < *newFacePoint; ++f) {
2926: PetscInt dof, off, d;
2928: PetscSectionGetDof(submesh->coneSection, f, &dof);
2929: PetscSectionGetOffset(submesh->coneSection, f, &off);
2930: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
2931: for (d = 0; d < dof; ++d) {
2932: const PetscInt p = submesh->cones[off + d];
2933: PetscInt v;
2935: for (v = 0; v < numFaceVertices; ++v) {
2936: if (subfaceVertices[v] == p) break;
2937: }
2938: if (v == numFaceVertices) break;
2939: }
2940: if (d == dof) {
2941: numFaces = 1;
2942: ((PetscInt *)faces)[0] = f;
2943: }
2944: }
2945: }
2946: #endif
2948: if (numFaces == 1) {
2949: /* Add the other cell neighbor for this face */
2950: DMPlexSetCone(subdm, subcell, faces);
2951: } else {
2952: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
2953: PetscBool posOriented;
2955: DMGetWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2956: origVertices = &orientedVertices[numFaceVertices];
2957: indices = &orientedVertices[numFaceVertices * 2];
2958: orientedSubVertices = &orientedVertices[numFaceVertices * 3];
2959: DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);
2960: /* TODO: I know that routine should return a permutation, not the indices */
2961: for (v = 0; v < numFaceVertices; ++v) {
2962: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
2963: for (ov = 0; ov < numFaceVertices; ++ov) {
2964: if (orientedVertices[ov] == vertex) {
2965: orientedSubVertices[ov] = subvertex;
2966: break;
2967: }
2968: }
2970: }
2971: DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);
2972: DMPlexSetCone(subdm, subcell, newFacePoint);
2973: DMRestoreWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);
2974: ++(*newFacePoint);
2975: }
2976: #if 0
2977: DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);
2978: #else
2979: DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **)&faces);
2980: #endif
2981: return 0;
2982: }
2984: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
2985: {
2986: MPI_Comm comm;
2987: DMLabel subpointMap;
2988: IS subvertexIS, subcellIS;
2989: const PetscInt *subVertices, *subCells;
2990: PetscInt numSubVertices, firstSubVertex, numSubCells;
2991: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
2992: PetscInt vStart, vEnd, c, f;
2994: PetscObjectGetComm((PetscObject)dm, &comm);
2995: /* Create subpointMap which marks the submesh */
2996: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
2997: DMPlexSetSubpointMap(subdm, subpointMap);
2998: DMLabelDestroy(&subpointMap);
2999: if (vertexLabel) DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);
3000: /* Setup chart */
3001: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3002: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3003: DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices);
3004: DMPlexSetVTKCellHeight(subdm, 1);
3005: /* Set cone sizes */
3006: firstSubVertex = numSubCells;
3007: firstSubFace = numSubCells + numSubVertices;
3008: newFacePoint = firstSubFace;
3009: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3010: if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
3011: DMLabelGetStratumIS(subpointMap, 2, &subcellIS);
3012: if (subcellIS) ISGetIndices(subcellIS, &subCells);
3013: for (c = 0; c < numSubCells; ++c) DMPlexSetConeSize(subdm, c, 1);
3014: for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) DMPlexSetConeSize(subdm, f, nFV);
3015: DMSetUp(subdm);
3016: /* Create face cones */
3017: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3018: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3019: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3020: for (c = 0; c < numSubCells; ++c) {
3021: const PetscInt cell = subCells[c];
3022: const PetscInt subcell = c;
3023: PetscInt *closure = NULL;
3024: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
3026: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
3027: for (cl = 0; cl < closureSize * 2; cl += 2) {
3028: const PetscInt point = closure[cl];
3029: PetscInt subVertex;
3031: if ((point >= vStart) && (point < vEnd)) {
3032: ++numCorners;
3033: PetscFindInt(point, numSubVertices, subVertices, &subVertex);
3034: if (subVertex >= 0) {
3035: closure[faceSize] = point;
3036: subface[faceSize] = firstSubVertex + subVertex;
3037: ++faceSize;
3038: }
3039: }
3040: }
3042: if (faceSize == nFV) DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);
3043: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
3044: }
3045: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3046: DMPlexSymmetrize(subdm);
3047: DMPlexStratify(subdm);
3048: /* Build coordinates */
3049: {
3050: PetscSection coordSection, subCoordSection;
3051: Vec coordinates, subCoordinates;
3052: PetscScalar *coords, *subCoords;
3053: PetscInt numComp, coordSize, v;
3054: const char *name;
3056: DMGetCoordinateSection(dm, &coordSection);
3057: DMGetCoordinatesLocal(dm, &coordinates);
3058: DMGetCoordinateSection(subdm, &subCoordSection);
3059: PetscSectionSetNumFields(subCoordSection, 1);
3060: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3061: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3062: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices);
3063: for (v = 0; v < numSubVertices; ++v) {
3064: const PetscInt vertex = subVertices[v];
3065: const PetscInt subvertex = firstSubVertex + v;
3066: PetscInt dof;
3068: PetscSectionGetDof(coordSection, vertex, &dof);
3069: PetscSectionSetDof(subCoordSection, subvertex, dof);
3070: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3071: }
3072: PetscSectionSetUp(subCoordSection);
3073: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3074: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3075: PetscObjectGetName((PetscObject)coordinates, &name);
3076: PetscObjectSetName((PetscObject)subCoordinates, name);
3077: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3078: VecSetType(subCoordinates, VECSTANDARD);
3079: if (coordSize) {
3080: VecGetArray(coordinates, &coords);
3081: VecGetArray(subCoordinates, &subCoords);
3082: for (v = 0; v < numSubVertices; ++v) {
3083: const PetscInt vertex = subVertices[v];
3084: const PetscInt subvertex = firstSubVertex + v;
3085: PetscInt dof, off, sdof, soff, d;
3087: PetscSectionGetDof(coordSection, vertex, &dof);
3088: PetscSectionGetOffset(coordSection, vertex, &off);
3089: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3090: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3092: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3093: }
3094: VecRestoreArray(coordinates, &coords);
3095: VecRestoreArray(subCoordinates, &subCoords);
3096: }
3097: DMSetCoordinatesLocal(subdm, subCoordinates);
3098: VecDestroy(&subCoordinates);
3099: }
3100: /* Cleanup */
3101: if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
3102: ISDestroy(&subvertexIS);
3103: if (subcellIS) ISRestoreIndices(subcellIS, &subCells);
3104: ISDestroy(&subcellIS);
3105: return 0;
3106: }
3108: /* TODO: Fix this to properly propagate up error conditions it may find */
3109: static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3110: {
3111: PetscInt subPoint;
3114: PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3115: if (ierr) return -1;
3116: return subPoint < 0 ? subPoint : firstSubPoint + subPoint;
3117: }
3119: /* TODO: Fix this to properly propagate up error conditions it may find */
3120: static inline PetscInt DMPlexFilterPointPerm_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[], const PetscInt subIndices[])
3121: {
3122: PetscInt subPoint;
3125: PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3126: if (ierr) return -1;
3127: return subPoint < 0 ? subPoint : firstSubPoint + (subIndices ? subIndices[subPoint] : subPoint);
3128: }
3130: static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm)
3131: {
3132: DMLabel depthLabel;
3133: PetscInt Nl, l, d;
3135: // Reset depth label for fast lookup
3136: DMPlexGetDepthLabel(dm, &depthLabel);
3137: DMLabelMakeAllInvalid_Internal(depthLabel);
3138: DMGetNumLabels(dm, &Nl);
3139: for (l = 0; l < Nl; ++l) {
3140: DMLabel label, newlabel;
3141: const char *lname;
3142: PetscBool isDepth, isDim, isCelltype, isVTK;
3143: IS valueIS;
3144: const PetscInt *values;
3145: PetscInt Nv, v;
3147: DMGetLabelName(dm, l, &lname);
3148: PetscStrcmp(lname, "depth", &isDepth);
3149: PetscStrcmp(lname, "dim", &isDim);
3150: PetscStrcmp(lname, "celltype", &isCelltype);
3151: PetscStrcmp(lname, "vtk", &isVTK);
3152: if (isDepth || isDim || isCelltype || isVTK) continue;
3153: DMCreateLabel(subdm, lname);
3154: DMGetLabel(dm, lname, &label);
3155: DMGetLabel(subdm, lname, &newlabel);
3156: DMLabelGetDefaultValue(label, &v);
3157: DMLabelSetDefaultValue(newlabel, v);
3158: DMLabelGetValueIS(label, &valueIS);
3159: ISGetLocalSize(valueIS, &Nv);
3160: ISGetIndices(valueIS, &values);
3161: for (v = 0; v < Nv; ++v) {
3162: IS pointIS;
3163: const PetscInt *points;
3164: PetscInt Np, p;
3166: DMLabelGetStratumIS(label, values[v], &pointIS);
3167: ISGetLocalSize(pointIS, &Np);
3168: ISGetIndices(pointIS, &points);
3169: for (p = 0; p < Np; ++p) {
3170: const PetscInt point = points[p];
3171: PetscInt subp;
3173: DMPlexGetPointDepth(dm, point, &d);
3174: subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]);
3175: if (subp >= 0) DMLabelSetValue(newlabel, subp, values[v]);
3176: }
3177: ISRestoreIndices(pointIS, &points);
3178: ISDestroy(&pointIS);
3179: }
3180: ISRestoreIndices(valueIS, &values);
3181: ISDestroy(&valueIS);
3182: }
3183: return 0;
3184: }
3186: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm)
3187: {
3188: MPI_Comm comm;
3189: DMLabel subpointMap;
3190: IS *subpointIS;
3191: const PetscInt **subpoints;
3192: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3193: PetscInt totSubPoints = 0, maxConeSize, dim, sdim, cdim, p, d, v;
3194: PetscMPIInt rank;
3196: PetscObjectGetComm((PetscObject)dm, &comm);
3197: MPI_Comm_rank(comm, &rank);
3198: /* Create subpointMap which marks the submesh */
3199: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3200: DMPlexSetSubpointMap(subdm, subpointMap);
3201: if (cellHeight) {
3202: if (isCohesive) DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);
3203: else DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);
3204: } else {
3205: DMLabel depth;
3206: IS pointIS;
3207: const PetscInt *points;
3208: PetscInt numPoints = 0;
3210: DMPlexGetDepthLabel(dm, &depth);
3211: DMLabelGetStratumIS(label, value, &pointIS);
3212: if (pointIS) {
3213: ISGetIndices(pointIS, &points);
3214: ISGetLocalSize(pointIS, &numPoints);
3215: }
3216: for (p = 0; p < numPoints; ++p) {
3217: PetscInt *closure = NULL;
3218: PetscInt closureSize, c, pdim;
3220: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3221: for (c = 0; c < closureSize * 2; c += 2) {
3222: DMLabelGetValue(depth, closure[c], &pdim);
3223: DMLabelSetValue(subpointMap, closure[c], pdim);
3224: }
3225: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
3226: }
3227: if (pointIS) ISRestoreIndices(pointIS, &points);
3228: ISDestroy(&pointIS);
3229: }
3230: /* Setup chart */
3231: DMGetDimension(dm, &dim);
3232: DMGetCoordinateDim(dm, &cdim);
3233: PetscMalloc4(dim + 1, &numSubPoints, dim + 1, &firstSubPoint, dim + 1, &subpointIS, dim + 1, &subpoints);
3234: for (d = 0; d <= dim; ++d) {
3235: DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);
3236: totSubPoints += numSubPoints[d];
3237: }
3238: // Determine submesh dimension
3239: DMGetDimension(subdm, &sdim);
3240: if (sdim > 0) {
3241: // Calling function knows what dimension to use, and we include neighboring cells as well
3242: sdim = dim;
3243: } else {
3244: // We reset the subdimension based on what is being selected
3245: PetscInt lsdim;
3246: for (lsdim = dim; lsdim >= 0; --lsdim)
3247: if (numSubPoints[lsdim]) break;
3248: MPI_Allreduce(&lsdim, &sdim, 1, MPIU_INT, MPIU_MAX, comm);
3249: DMSetDimension(subdm, sdim);
3250: DMSetCoordinateDim(subdm, cdim);
3251: }
3252: DMPlexSetChart(subdm, 0, totSubPoints);
3253: DMPlexSetVTKCellHeight(subdm, cellHeight);
3254: /* Set cone sizes */
3255: firstSubPoint[sdim] = 0;
3256: firstSubPoint[0] = firstSubPoint[sdim] + numSubPoints[sdim];
3257: if (sdim > 1) firstSubPoint[sdim - 1] = firstSubPoint[0] + numSubPoints[0];
3258: if (sdim > 2) firstSubPoint[sdim - 2] = firstSubPoint[sdim - 1] + numSubPoints[sdim - 1];
3259: for (d = 0; d <= sdim; ++d) {
3260: DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);
3261: if (subpointIS[d]) ISGetIndices(subpointIS[d], &subpoints[d]);
3262: }
3263: /* We do not want this label automatically computed, instead we compute it here */
3264: DMCreateLabel(subdm, "celltype");
3265: for (d = 0; d <= sdim; ++d) {
3266: for (p = 0; p < numSubPoints[d]; ++p) {
3267: const PetscInt point = subpoints[d][p];
3268: const PetscInt subpoint = firstSubPoint[d] + p;
3269: const PetscInt *cone;
3270: PetscInt coneSize;
3272: DMPlexGetConeSize(dm, point, &coneSize);
3273: if (cellHeight && (d == sdim)) {
3274: PetscInt coneSizeNew, c, val;
3276: DMPlexGetCone(dm, point, &cone);
3277: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3278: DMLabelGetValue(subpointMap, cone[c], &val);
3279: if (val >= 0) coneSizeNew++;
3280: }
3281: DMPlexSetConeSize(subdm, subpoint, coneSizeNew);
3282: DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);
3283: } else {
3284: DMPolytopeType ct;
3286: DMPlexSetConeSize(subdm, subpoint, coneSize);
3287: DMPlexGetCellType(dm, point, &ct);
3288: DMPlexSetCellType(subdm, subpoint, ct);
3289: }
3290: }
3291: }
3292: DMLabelDestroy(&subpointMap);
3293: DMSetUp(subdm);
3294: /* Set cones */
3295: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3296: PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew);
3297: for (d = 0; d <= sdim; ++d) {
3298: for (p = 0; p < numSubPoints[d]; ++p) {
3299: const PetscInt point = subpoints[d][p];
3300: const PetscInt subpoint = firstSubPoint[d] + p;
3301: const PetscInt *cone, *ornt;
3302: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3304: if (d == sdim - 1) {
3305: const PetscInt *support, *cone, *ornt;
3306: PetscInt supportSize, coneSize, s, subc;
3308: DMPlexGetSupport(dm, point, &support);
3309: DMPlexGetSupportSize(dm, point, &supportSize);
3310: for (s = 0; s < supportSize; ++s) {
3311: PetscBool isHybrid = PETSC_FALSE;
3313: DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);
3314: if (!isHybrid) continue;
3315: PetscFindInt(support[s], numSubPoints[d + 1], subpoints[d + 1], &subc);
3316: if (subc >= 0) {
3317: const PetscInt ccell = subpoints[d + 1][subc];
3319: DMPlexGetCone(dm, ccell, &cone);
3320: DMPlexGetConeSize(dm, ccell, &coneSize);
3321: DMPlexGetConeOrientation(dm, ccell, &ornt);
3322: for (c = 0; c < coneSize; ++c) {
3323: if (cone[c] == point) {
3324: fornt = ornt[c];
3325: break;
3326: }
3327: }
3328: break;
3329: }
3330: }
3331: }
3332: DMPlexGetConeSize(dm, point, &coneSize);
3333: DMPlexGetConeSize(subdm, subpoint, &subconeSize);
3334: DMPlexGetCone(dm, point, &cone);
3335: DMPlexGetConeOrientation(dm, point, &ornt);
3336: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3337: PetscFindInt(cone[c], numSubPoints[d - 1], subpoints[d - 1], &subc);
3338: if (subc >= 0) {
3339: coneNew[coneSizeNew] = firstSubPoint[d - 1] + subc;
3340: orntNew[coneSizeNew] = ornt[c];
3341: ++coneSizeNew;
3342: }
3343: }
3345: DMPlexSetCone(subdm, subpoint, coneNew);
3346: DMPlexSetConeOrientation(subdm, subpoint, orntNew);
3347: if (fornt < 0) DMPlexOrientPoint(subdm, subpoint, fornt);
3348: }
3349: }
3350: PetscFree2(coneNew, orntNew);
3351: DMPlexSymmetrize(subdm);
3352: DMPlexStratify(subdm);
3353: /* Build coordinates */
3354: {
3355: PetscSection coordSection, subCoordSection;
3356: Vec coordinates, subCoordinates;
3357: PetscScalar *coords, *subCoords;
3358: PetscInt cdim, numComp, coordSize;
3359: const char *name;
3361: DMGetCoordinateDim(dm, &cdim);
3362: DMGetCoordinateSection(dm, &coordSection);
3363: DMGetCoordinatesLocal(dm, &coordinates);
3364: DMGetCoordinateSection(subdm, &subCoordSection);
3365: PetscSectionSetNumFields(subCoordSection, 1);
3366: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3367: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3368: PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0] + numSubPoints[0]);
3369: for (v = 0; v < numSubPoints[0]; ++v) {
3370: const PetscInt vertex = subpoints[0][v];
3371: const PetscInt subvertex = firstSubPoint[0] + v;
3372: PetscInt dof;
3374: PetscSectionGetDof(coordSection, vertex, &dof);
3375: PetscSectionSetDof(subCoordSection, subvertex, dof);
3376: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3377: }
3378: PetscSectionSetUp(subCoordSection);
3379: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3380: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3381: PetscObjectGetName((PetscObject)coordinates, &name);
3382: PetscObjectSetName((PetscObject)subCoordinates, name);
3383: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3384: VecSetBlockSize(subCoordinates, cdim);
3385: VecSetType(subCoordinates, VECSTANDARD);
3386: VecGetArray(coordinates, &coords);
3387: VecGetArray(subCoordinates, &subCoords);
3388: for (v = 0; v < numSubPoints[0]; ++v) {
3389: const PetscInt vertex = subpoints[0][v];
3390: const PetscInt subvertex = firstSubPoint[0] + v;
3391: PetscInt dof, off, sdof, soff, d;
3393: PetscSectionGetDof(coordSection, vertex, &dof);
3394: PetscSectionGetOffset(coordSection, vertex, &off);
3395: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3396: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3398: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3399: }
3400: VecRestoreArray(coordinates, &coords);
3401: VecRestoreArray(subCoordinates, &subCoords);
3402: DMSetCoordinatesLocal(subdm, subCoordinates);
3403: VecDestroy(&subCoordinates);
3404: }
3405: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
3406: {
3407: PetscSF sfPoint, sfPointSub;
3408: IS subpIS;
3409: const PetscSFNode *remotePoints;
3410: PetscSFNode *sremotePoints = NULL, *newLocalPoints = NULL, *newOwners = NULL;
3411: const PetscInt *localPoints, *subpoints, *rootdegree;
3412: PetscInt *slocalPoints = NULL, *sortedPoints = NULL, *sortedIndices = NULL;
3413: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl = 0, ll = 0, pStart, pEnd, p;
3414: PetscMPIInt rank, size;
3416: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
3417: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
3418: DMGetPointSF(dm, &sfPoint);
3419: DMGetPointSF(subdm, &sfPointSub);
3420: DMPlexGetChart(dm, &pStart, &pEnd);
3421: DMPlexGetChart(subdm, NULL, &numSubroots);
3422: DMPlexGetSubpointIS(subdm, &subpIS);
3423: if (subpIS) {
3424: PetscBool sorted = PETSC_TRUE;
3426: ISGetIndices(subpIS, &subpoints);
3427: ISGetLocalSize(subpIS, &numSubpoints);
3428: for (p = 1; p < numSubpoints; ++p) sorted = sorted && (subpoints[p] >= subpoints[p - 1]) ? PETSC_TRUE : PETSC_FALSE;
3429: if (!sorted) {
3430: PetscMalloc2(numSubpoints, &sortedPoints, numSubpoints, &sortedIndices);
3431: for (p = 0; p < numSubpoints; ++p) sortedIndices[p] = p;
3432: PetscArraycpy(sortedPoints, subpoints, numSubpoints);
3433: PetscSortIntWithArray(numSubpoints, sortedPoints, sortedIndices);
3434: }
3435: }
3436: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3437: if (numRoots >= 0) {
3438: PetscSFComputeDegreeBegin(sfPoint, &rootdegree);
3439: PetscSFComputeDegreeEnd(sfPoint, &rootdegree);
3440: PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners);
3441: for (p = 0; p < pEnd - pStart; ++p) {
3442: newLocalPoints[p].rank = -2;
3443: newLocalPoints[p].index = -2;
3444: }
3445: /* Set subleaves */
3446: for (l = 0; l < numLeaves; ++l) {
3447: const PetscInt point = localPoints[l];
3448: const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);
3450: if (subpoint < 0) continue;
3451: newLocalPoints[point - pStart].rank = rank;
3452: newLocalPoints[point - pStart].index = subpoint;
3453: ++numSubleaves;
3454: }
3455: /* Must put in owned subpoints */
3456: for (p = pStart; p < pEnd; ++p) {
3457: newOwners[p - pStart].rank = -3;
3458: newOwners[p - pStart].index = -3;
3459: }
3460: for (p = 0; p < numSubpoints; ++p) {
3461: /* Hold on to currently owned points */
3462: if (rootdegree[subpoints[p] - pStart]) newOwners[subpoints[p] - pStart].rank = rank + size;
3463: else newOwners[subpoints[p] - pStart].rank = rank;
3464: newOwners[subpoints[p] - pStart].index = p;
3465: }
3466: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3467: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3468: for (p = pStart; p < pEnd; ++p)
3469: if (newOwners[p - pStart].rank >= size) newOwners[p - pStart].rank -= size;
3470: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3471: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3472: PetscMalloc1(numSubleaves, &slocalPoints);
3473: PetscMalloc1(numSubleaves, &sremotePoints);
3474: for (l = 0; l < numLeaves; ++l) {
3475: const PetscInt point = localPoints[l];
3476: const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);
3478: if (subpoint < 0) continue;
3479: if (newLocalPoints[point].rank == rank) {
3480: ++ll;
3481: continue;
3482: }
3483: slocalPoints[sl] = subpoint;
3484: sremotePoints[sl].rank = newLocalPoints[point].rank;
3485: sremotePoints[sl].index = newLocalPoints[point].index;
3488: ++sl;
3489: }
3491: PetscFree2(newLocalPoints, newOwners);
3492: PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3493: }
3494: if (subpIS) ISRestoreIndices(subpIS, &subpoints);
3495: PetscFree2(sortedPoints, sortedIndices);
3496: }
3497: /* Filter labels */
3498: DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm);
3499: /* Cleanup */
3500: for (d = 0; d <= sdim; ++d) {
3501: if (subpointIS[d]) ISRestoreIndices(subpointIS[d], &subpoints[d]);
3502: ISDestroy(&subpointIS[d]);
3503: }
3504: PetscFree4(numSubPoints, firstSubPoint, subpointIS, subpoints);
3505: return 0;
3506: }
3508: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
3509: {
3510: DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);
3511: return 0;
3512: }
3514: /*@
3515: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
3517: Input Parameters:
3518: + dm - The original mesh
3519: . vertexLabel - The DMLabel marking points contained in the surface
3520: . value - The label value to use
3521: - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked
3523: Output Parameter:
3524: . subdm - The surface mesh
3526: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3528: Level: developer
3530: .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`
3531: @*/
3532: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
3533: {
3534: DMPlexInterpolatedFlag interpolated;
3535: PetscInt dim, cdim;
3539: DMGetDimension(dm, &dim);
3540: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3541: DMSetType(*subdm, DMPLEX);
3542: DMSetDimension(*subdm, dim - 1);
3543: DMGetCoordinateDim(dm, &cdim);
3544: DMSetCoordinateDim(*subdm, cdim);
3545: DMPlexIsInterpolated(dm, &interpolated);
3547: if (interpolated) {
3548: DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);
3549: } else {
3550: DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);
3551: }
3552: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3553: return 0;
3554: }
3556: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
3557: {
3558: MPI_Comm comm;
3559: DMLabel subpointMap;
3560: IS subvertexIS;
3561: const PetscInt *subVertices;
3562: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
3563: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
3564: PetscInt c, f;
3566: PetscObjectGetComm((PetscObject)dm, &comm);
3567: /* Create subpointMap which marks the submesh */
3568: DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);
3569: DMPlexSetSubpointMap(subdm, subpointMap);
3570: DMLabelDestroy(&subpointMap);
3571: DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);
3572: /* Setup chart */
3573: DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);
3574: DMLabelGetStratumSize(subpointMap, 2, &numSubCells);
3575: DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices);
3576: DMPlexSetVTKCellHeight(subdm, 1);
3577: /* Set cone sizes */
3578: firstSubVertex = numSubCells;
3579: firstSubFace = numSubCells + numSubVertices;
3580: newFacePoint = firstSubFace;
3581: DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);
3582: if (subvertexIS) ISGetIndices(subvertexIS, &subVertices);
3583: for (c = 0; c < numSubCells; ++c) DMPlexSetConeSize(subdm, c, 1);
3584: for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) DMPlexSetConeSize(subdm, f, nFV);
3585: DMSetUp(subdm);
3586: /* Create face cones */
3587: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
3588: DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3589: for (c = 0; c < numSubCells; ++c) {
3590: const PetscInt cell = subCells[c];
3591: const PetscInt subcell = c;
3592: const PetscInt *cone, *cells;
3593: PetscBool isHybrid = PETSC_FALSE;
3594: PetscInt numCells, subVertex, p, v;
3596: DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);
3597: if (!isHybrid) continue;
3598: DMPlexGetCone(dm, cell, &cone);
3599: for (v = 0; v < nFV; ++v) {
3600: PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);
3601: subface[v] = firstSubVertex + subVertex;
3602: }
3603: DMPlexSetCone(subdm, newFacePoint, subface);
3604: DMPlexSetCone(subdm, subcell, &newFacePoint);
3605: DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);
3606: /* Not true in parallel
3608: for (p = 0; p < numCells; ++p) {
3609: PetscInt negsubcell;
3610: PetscBool isHybrid = PETSC_FALSE;
3612: DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);
3613: if (isHybrid) continue;
3614: /* I know this is a crap search */
3615: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
3616: if (subCells[negsubcell] == cells[p]) break;
3617: }
3619: DMPlexSetCone(subdm, negsubcell, &newFacePoint);
3620: }
3621: DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);
3622: ++newFacePoint;
3623: }
3624: DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface);
3625: DMPlexSymmetrize(subdm);
3626: DMPlexStratify(subdm);
3627: /* Build coordinates */
3628: {
3629: PetscSection coordSection, subCoordSection;
3630: Vec coordinates, subCoordinates;
3631: PetscScalar *coords, *subCoords;
3632: PetscInt cdim, numComp, coordSize, v;
3633: const char *name;
3635: DMGetCoordinateDim(dm, &cdim);
3636: DMGetCoordinateSection(dm, &coordSection);
3637: DMGetCoordinatesLocal(dm, &coordinates);
3638: DMGetCoordinateSection(subdm, &subCoordSection);
3639: PetscSectionSetNumFields(subCoordSection, 1);
3640: PetscSectionGetFieldComponents(coordSection, 0, &numComp);
3641: PetscSectionSetFieldComponents(subCoordSection, 0, numComp);
3642: PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices);
3643: for (v = 0; v < numSubVertices; ++v) {
3644: const PetscInt vertex = subVertices[v];
3645: const PetscInt subvertex = firstSubVertex + v;
3646: PetscInt dof;
3648: PetscSectionGetDof(coordSection, vertex, &dof);
3649: PetscSectionSetDof(subCoordSection, subvertex, dof);
3650: PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);
3651: }
3652: PetscSectionSetUp(subCoordSection);
3653: PetscSectionGetStorageSize(subCoordSection, &coordSize);
3654: VecCreate(PETSC_COMM_SELF, &subCoordinates);
3655: PetscObjectGetName((PetscObject)coordinates, &name);
3656: PetscObjectSetName((PetscObject)subCoordinates, name);
3657: VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);
3658: VecSetBlockSize(subCoordinates, cdim);
3659: VecSetType(subCoordinates, VECSTANDARD);
3660: VecGetArray(coordinates, &coords);
3661: VecGetArray(subCoordinates, &subCoords);
3662: for (v = 0; v < numSubVertices; ++v) {
3663: const PetscInt vertex = subVertices[v];
3664: const PetscInt subvertex = firstSubVertex + v;
3665: PetscInt dof, off, sdof, soff, d;
3667: PetscSectionGetDof(coordSection, vertex, &dof);
3668: PetscSectionGetOffset(coordSection, vertex, &off);
3669: PetscSectionGetDof(subCoordSection, subvertex, &sdof);
3670: PetscSectionGetOffset(subCoordSection, subvertex, &soff);
3672: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3673: }
3674: VecRestoreArray(coordinates, &coords);
3675: VecRestoreArray(subCoordinates, &subCoords);
3676: DMSetCoordinatesLocal(subdm, subCoordinates);
3677: VecDestroy(&subCoordinates);
3678: }
3679: /* Build SF */
3680: CHKMEMQ;
3681: {
3682: PetscSF sfPoint, sfPointSub;
3683: const PetscSFNode *remotePoints;
3684: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
3685: const PetscInt *localPoints;
3686: PetscInt *slocalPoints;
3687: PetscInt numRoots, numLeaves, numSubRoots = numSubCells + numSubFaces + numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
3688: PetscMPIInt rank;
3690: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
3691: DMGetPointSF(dm, &sfPoint);
3692: DMGetPointSF(subdm, &sfPointSub);
3693: DMPlexGetChart(dm, &pStart, &pEnd);
3694: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3695: PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
3696: if (numRoots >= 0) {
3697: /* Only vertices should be shared */
3698: PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners);
3699: for (p = 0; p < pEnd - pStart; ++p) {
3700: newLocalPoints[p].rank = -2;
3701: newLocalPoints[p].index = -2;
3702: }
3703: /* Set subleaves */
3704: for (l = 0; l < numLeaves; ++l) {
3705: const PetscInt point = localPoints[l];
3706: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3709: if (subPoint < 0) continue;
3710: newLocalPoints[point - pStart].rank = rank;
3711: newLocalPoints[point - pStart].index = subPoint;
3712: ++numSubLeaves;
3713: }
3714: /* Must put in owned subpoints */
3715: for (p = pStart; p < pEnd; ++p) {
3716: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
3718: if (subPoint < 0) {
3719: newOwners[p - pStart].rank = -3;
3720: newOwners[p - pStart].index = -3;
3721: } else {
3722: newOwners[p - pStart].rank = rank;
3723: newOwners[p - pStart].index = subPoint;
3724: }
3725: }
3726: PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3727: PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);
3728: PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3729: PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints, MPI_REPLACE);
3730: PetscMalloc1(numSubLeaves, &slocalPoints);
3731: PetscMalloc1(numSubLeaves, &sremotePoints);
3732: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
3733: const PetscInt point = localPoints[l];
3734: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
3736: if (subPoint < 0) continue;
3737: if (newLocalPoints[point].rank == rank) {
3738: ++ll;
3739: continue;
3740: }
3741: slocalPoints[sl] = subPoint;
3742: sremotePoints[sl].rank = newLocalPoints[point].rank;
3743: sremotePoints[sl].index = newLocalPoints[point].index;
3746: ++sl;
3747: }
3748: PetscFree2(newLocalPoints, newOwners);
3750: PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);
3751: }
3752: }
3753: CHKMEMQ;
3754: /* Cleanup */
3755: if (subvertexIS) ISRestoreIndices(subvertexIS, &subVertices);
3756: ISDestroy(&subvertexIS);
3757: PetscFree(subCells);
3758: return 0;
3759: }
3761: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
3762: {
3763: DMLabel label = NULL;
3765: if (labelname) DMGetLabel(dm, labelname, &label);
3766: DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);
3767: return 0;
3768: }
3770: /*@C
3771: DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label can be given to restrict the cells.
3773: Input Parameters:
3774: + dm - The original mesh
3775: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
3776: . label - A label name, or NULL
3777: - value - A label value
3779: Output Parameter:
3780: . subdm - The surface mesh
3782: Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap().
3784: Level: developer
3786: .seealso: `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()`
3787: @*/
3788: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
3789: {
3790: PetscInt dim, cdim, depth;
3794: DMGetDimension(dm, &dim);
3795: DMPlexGetDepth(dm, &depth);
3796: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3797: DMSetType(*subdm, DMPLEX);
3798: DMSetDimension(*subdm, dim - 1);
3799: DMGetCoordinateDim(dm, &cdim);
3800: DMSetCoordinateDim(*subdm, cdim);
3801: if (depth == dim) {
3802: DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);
3803: } else {
3804: DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);
3805: }
3806: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3807: return 0;
3808: }
3810: /*@
3811: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
3813: Input Parameters:
3814: + dm - The original mesh
3815: . cellLabel - The DMLabel marking cells contained in the new mesh
3816: - value - The label value to use
3818: Output Parameter:
3819: . subdm - The new mesh
3821: Notes:
3822: This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`.
3824: Level: developer
3826: .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`, `DMPlexCreateSubmesh()`
3827: @*/
3828: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm)
3829: {
3830: PetscInt dim, overlap;
3834: DMGetDimension(dm, &dim);
3835: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
3836: DMSetType(*subdm, DMPLEX);
3837: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
3838: DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);
3839: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm);
3840: // It is possible to obtain a surface mesh where some faces are in SF
3841: // We should either mark the mesh as having an overlap, or delete these from the SF
3842: DMPlexGetOverlap(dm, &overlap);
3843: if (!overlap) {
3844: PetscSF sf;
3845: const PetscInt *leaves;
3846: PetscInt cStart, cEnd, Nl;
3847: PetscBool hasSubcell = PETSC_FALSE, ghasSubcell;
3849: DMPlexGetHeightStratum(*subdm, 0, &cStart, &cEnd);
3850: DMGetPointSF(*subdm, &sf);
3851: PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL);
3852: for (PetscInt l = 0; l < Nl; ++l) {
3853: const PetscInt point = leaves ? leaves[l] : l;
3855: if (point >= cStart && point < cEnd) {
3856: hasSubcell = PETSC_TRUE;
3857: break;
3858: }
3859: }
3860: MPIU_Allreduce(&hasSubcell, &ghasSubcell, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
3861: if (ghasSubcell) DMPlexSetOverlap(*subdm, NULL, 1);
3862: }
3863: return 0;
3864: }
3866: /*@
3867: DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values
3869: Input Parameter:
3870: . dm - The submesh DM
3872: Output Parameter:
3873: . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3875: Level: developer
3877: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
3878: @*/
3879: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
3880: {
3883: *subpointMap = ((DM_Plex *)dm->data)->subpointMap;
3884: return 0;
3885: }
3887: /*@
3888: DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values
3890: Input Parameters:
3891: + dm - The submesh DM
3892: - subpointMap - The DMLabel of all the points from the original mesh in this submesh
3894: Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh()
3896: Level: developer
3898: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
3899: @*/
3900: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
3901: {
3902: DM_Plex *mesh = (DM_Plex *)dm->data;
3903: DMLabel tmp;
3906: tmp = mesh->subpointMap;
3907: mesh->subpointMap = subpointMap;
3908: PetscObjectReference((PetscObject)mesh->subpointMap);
3909: DMLabelDestroy(&tmp);
3910: return 0;
3911: }
3913: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
3914: {
3915: DMLabel spmap;
3916: PetscInt depth, d;
3918: DMPlexGetSubpointMap(dm, &spmap);
3919: DMPlexGetDepth(dm, &depth);
3920: if (spmap && depth >= 0) {
3921: DM_Plex *mesh = (DM_Plex *)dm->data;
3922: PetscInt *points, *depths;
3923: PetscInt pStart, pEnd, p, off;
3925: DMPlexGetChart(dm, &pStart, &pEnd);
3927: PetscMalloc1(pEnd, &points);
3928: DMGetWorkArray(dm, depth + 1, MPIU_INT, &depths);
3929: depths[0] = depth;
3930: depths[1] = 0;
3931: for (d = 2; d <= depth; ++d) depths[d] = depth + 1 - d;
3932: for (d = 0, off = 0; d <= depth; ++d) {
3933: const PetscInt dep = depths[d];
3934: PetscInt depStart, depEnd, n;
3936: DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);
3937: DMLabelGetStratumSize(spmap, dep, &n);
3938: if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */
3940: } else {
3941: if (!n) {
3942: if (d == 0) {
3943: /* Missing cells */
3944: for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = -1;
3945: } else {
3946: /* Missing faces */
3947: for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = PETSC_MAX_INT;
3948: }
3949: }
3950: }
3951: if (n) {
3952: IS is;
3953: const PetscInt *opoints;
3955: DMLabelGetStratumIS(spmap, dep, &is);
3956: ISGetIndices(is, &opoints);
3957: for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
3958: ISRestoreIndices(is, &opoints);
3959: ISDestroy(&is);
3960: }
3961: }
3962: DMRestoreWorkArray(dm, depth + 1, MPIU_INT, &depths);
3964: ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);
3965: PetscObjectStateGet((PetscObject)spmap, &mesh->subpointState);
3966: }
3967: return 0;
3968: }
3970: /*@
3971: DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data
3973: Input Parameter:
3974: . dm - The submesh DM
3976: Output Parameter:
3977: . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh
3979: Note: This IS is guaranteed to be sorted by the construction of the submesh
3981: Level: developer
3983: .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()`
3984: @*/
3985: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
3986: {
3987: DM_Plex *mesh = (DM_Plex *)dm->data;
3988: DMLabel spmap;
3989: PetscObjectState state;
3993: DMPlexGetSubpointMap(dm, &spmap);
3994: PetscObjectStateGet((PetscObject)spmap, &state);
3995: if (state != mesh->subpointState || !mesh->subpointIS) DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);
3996: *subpointIS = mesh->subpointIS;
3997: return 0;
3998: }
4000: /*@
4001: DMGetEnclosureRelation - Get the relationship between dmA and dmB
4003: Input Parameters:
4004: + dmA - The first DM
4005: - dmB - The second DM
4007: Output Parameter:
4008: . rel - The relation of dmA to dmB
4010: Level: intermediate
4012: .seealso: `DMGetEnclosurePoint()`
4013: @*/
4014: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
4015: {
4016: DM plexA, plexB, sdm;
4017: DMLabel spmap;
4018: PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;
4021: *rel = DM_ENC_NONE;
4022: if (!dmA || !dmB) return 0;
4025: if (dmA == dmB) {
4026: *rel = DM_ENC_EQUALITY;
4027: return 0;
4028: }
4029: DMConvert(dmA, DMPLEX, &plexA);
4030: DMConvert(dmB, DMPLEX, &plexB);
4031: DMPlexGetChart(plexA, &pStartA, &pEndA);
4032: DMPlexGetChart(plexB, &pStartB, &pEndB);
4033: /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
4034: - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
4035: if ((pStartA == pStartB) && (pEndA == pEndB)) {
4036: *rel = DM_ENC_EQUALITY;
4037: goto end;
4038: }
4039: NpA = pEndA - pStartA;
4040: NpB = pEndB - pStartB;
4041: if (NpA == NpB) goto end;
4042: sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
4043: DMPlexGetSubpointMap(sdm, &spmap);
4044: if (!spmap) goto end;
4045: /* TODO Check the space mapped to by subpointMap is same size as dm */
4046: if (NpA > NpB) {
4047: *rel = DM_ENC_SUPERMESH;
4048: } else {
4049: *rel = DM_ENC_SUBMESH;
4050: }
4051: end:
4052: DMDestroy(&plexA);
4053: DMDestroy(&plexB);
4054: return 0;
4055: }
4057: /*@
4058: DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB
4060: Input Parameters:
4061: + dmA - The first DM
4062: . dmB - The second DM
4063: . etype - The type of enclosure relation that dmA has to dmB
4064: - pB - A point of dmB
4066: Output Parameter:
4067: . pA - The corresponding point of dmA
4069: Level: intermediate
4071: .seealso: `DMGetEnclosureRelation()`
4072: @*/
4073: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
4074: {
4075: DM sdm;
4076: IS subpointIS;
4077: const PetscInt *subpoints;
4078: PetscInt numSubpoints;
4080: /* TODO Cache the IS, making it look like an index */
4081: switch (etype) {
4082: case DM_ENC_SUPERMESH:
4083: sdm = dmB;
4084: DMPlexGetSubpointIS(sdm, &subpointIS);
4085: ISGetIndices(subpointIS, &subpoints);
4086: *pA = subpoints[pB];
4087: ISRestoreIndices(subpointIS, &subpoints);
4088: break;
4089: case DM_ENC_SUBMESH:
4090: sdm = dmA;
4091: DMPlexGetSubpointIS(sdm, &subpointIS);
4092: ISGetLocalSize(subpointIS, &numSubpoints);
4093: ISGetIndices(subpointIS, &subpoints);
4094: PetscFindInt(pB, numSubpoints, subpoints, pA);
4095: if (*pA < 0) {
4096: DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");
4097: DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");
4098: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB);
4099: }
4100: ISRestoreIndices(subpointIS, &subpoints);
4101: break;
4102: case DM_ENC_EQUALITY:
4103: case DM_ENC_NONE:
4104: *pA = pB;
4105: break;
4106: case DM_ENC_UNKNOWN: {
4107: DMEnclosureType enc;
4109: DMGetEnclosureRelation(dmA, dmB, &enc);
4110: DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);
4111: } break;
4112: default:
4113: SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int)etype);
4114: }
4115: return 0;
4116: }