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