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