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] = 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, PetscBool useCone)
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, useCone, &closureSize, &closure));
234: for (c = 0; c < closureSize * 2; c += 2) PetscCall(DMLabelSetValue(label, closure[c], values[v]));
235: PetscCall(DMPlexRestoreTransitiveClosure(dm, q, useCone, &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, we add their transitive closure
287: Input Parameters:
288: + dm - The `DM`
289: - label - A `DMLabel` marking the points
291: Output Parameter:
292: . label - A `DMLabel` marking all 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, PETSC_TRUE));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: DMPlexLabelCompleteStar - Starting with a label marking points, we add their star
308: Input Parameters:
309: + dm - The `DM`
310: - label - A `DMLabel` marking the points
312: Output Parameter:
313: . label - A `DMLabel` marking all points in the star
315: Level: developer
317: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelComplete()`
318: @*/
319: PetscErrorCode DMPlexLabelCompleteStar(DM dm, DMLabel label)
320: {
321: PetscFunctionBegin;
322: PetscCall(DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE, PETSC_FALSE));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscInt label_defval_private = -1;
327: static PetscInt label_errval_private = -2;
328: static void MPIAPI label_value_check(void *a, void *b, int *len, MPI_Datatype *datatype)
329: {
330: const int N = *len;
332: if (*datatype == MPIU_INT) {
333: PetscInt *A = (PetscInt *)a;
334: PetscInt *B = (PetscInt *)b;
336: for (int i = 0; i < N; i++) {
337: // Propagate errors
338: if (A[i] == label_errval_private || B[i] == label_errval_private) {
339: B[i] = label_errval_private;
340: continue;
341: }
342: // Default values do not propagate
343: if (A[i] == label_defval_private) continue;
344: // Override default values
345: if (B[i] == label_defval_private) {
346: B[i] = A[i];
347: continue;
348: }
349: B[i] = A[i] != B[i] ? label_errval_private : B[i];
350: }
351: }
352: }
354: /*@
355: DMPlexCheckLabel - Check that points matched by the pointSF have the same value in the label
357: Input Parameters:
358: + dm - The `DM`
359: . reduceop - The MPI reduction operation to use for comparison, or `MPI_OP_NULL` for the default
360: - label - A `DMLabel` marking the points
362: Level: advanced
364: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelCohesiveComplete()`
365: @*/
366: PetscErrorCode DMPlexCheckLabel(DM dm, MPI_Op reduceop, DMLabel label)
367: {
368: PetscSF sf;
369: IS valueIS;
370: MPI_Op lreduceop = reduceop;
371: const PetscInt *leaves, *values, *degree;
372: const PetscSFNode *remotes;
373: PetscInt *rvalues, *lvalues;
374: PetscInt Nr, Nl, Nv;
375: PetscBool mismatch = PETSC_FALSE, gmismatch;
376: MPI_Comm comm;
378: PetscFunctionBegin;
379: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
380: PetscCall(DMLabelGetDefaultValue(label, &label_defval_private));
381: label_errval_private = label_defval_private - 1;
383: PetscCall(DMLabelGetValueIS(label, &valueIS));
384: PetscCall(ISGetLocalSize(valueIS, &Nv));
385: PetscCall(ISGetIndices(valueIS, &values));
386: for (PetscInt v = 0; v < Nv; ++v) {
387: PetscCheck(values[v] != label_errval_private, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label value %" PetscInt_FMT " matches the choice for the error value", values[v]);
388: }
390: PetscCall(DMGetPointSF(dm, &sf));
391: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &leaves, &remotes));
392: PetscCall(PetscSFComputeDegreeBegin(sf, °ree));
393: PetscCall(PetscSFComputeDegreeEnd(sf, °ree));
394: PetscCall(PetscMalloc2(Nr, &rvalues, Nr, &lvalues));
395: PetscCall(PetscSFView(sf, NULL));
397: for (PetscInt l = 0; l < Nl; ++l) PetscCall(DMLabelGetValue(label, leaves[l], &lvalues[leaves[l]]));
398: for (PetscInt r = 0; r < Nr; ++r) rvalues[r] = label_defval_private;
399: if (reduceop == MPI_OP_NULL) PetscCallMPI(MPI_Op_create(label_value_check, PETSC_TRUE, &lreduceop));
400: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, lvalues, rvalues, lreduceop));
401: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, lvalues, rvalues, lreduceop));
402: if (reduceop == MPI_OP_NULL) PetscCallMPI(MPI_Op_free(&lreduceop));
404: PetscCall(ISGetIndices(valueIS, &values));
405: for (PetscInt i = 0; i < Nv; ++i) {
406: const PetscInt val = values[i];
407: IS stratumIS;
408: const PetscInt *points;
409: PetscInt Ns;
411: PetscCall(DMLabelGetStratumIS(label, val, &stratumIS));
412: PetscCall(ISGetLocalSize(stratumIS, &Ns));
413: PetscCall(ISGetIndices(stratumIS, &points));
414: for (PetscInt s = 0; s < Ns; ++s) {
415: const PetscInt point = points[s];
416: PetscInt val;
418: // Check only shared points
419: if (degree[point]) {
420: PetscCall(DMLabelGetValue(label, point, &val));
421: if (val != rvalues[point]) {
422: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Label value mismatch (%" PetscInt_FMT " != %" PetscInt_FMT ") at point %" PetscInt_FMT "\n", PetscGlobalRank, val, rvalues[point], point));
423: mismatch = PETSC_TRUE;
424: }
425: }
426: }
427: PetscCall(ISRestoreIndices(stratumIS, &points));
428: PetscCall(ISDestroy(&stratumIS));
429: }
430: PetscCall(ISRestoreIndices(valueIS, &values));
431: PetscCall(ISDestroy(&valueIS));
432: PetscCall(PetscFree2(rvalues, lvalues));
433: PetscCallMPI(MPIU_Allreduce(&mismatch, &gmismatch, 1, MPI_C_BOOL, MPI_LOR, comm));
434: PetscCheck(!gmismatch, comm, PETSC_ERR_ARG_WRONG, "Label value mismatch detected");
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: DMPlexReconcileLabel - Force points matched by the pointSF to have the same value in the label
441: Input Parameters:
442: + dm - The `DM`
443: . reduceop - The MPI reduction operation to use for comparison, or `MPI_OP_NULL` to use the root value
444: . newValue - Value to additionally give to newly added points, or `PETSC_DETERMINE` use only the normal value
445: - label - A `DMLabel` marking the points
447: Level: advanced
449: Note:
450: Newly added values can come from isolated embedded surface points on processes without a surface face.
452: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelCohesiveComplete()`
453: @*/
454: PetscErrorCode DMPlexReconcileLabel(DM dm, MPI_Op reduceop, PetscInt newValue, DMLabel label)
455: {
456: const PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
457: const PetscInt shift3 = 300;
458: PetscSF sf;
459: IS valueIS;
460: const PetscInt *leaves, *values, *degree;
461: const PetscSFNode *remotes;
462: PetscInt *rvalues, *lvalues;
463: PetscInt Nr, Nl, Nv;
464: MPI_Comm comm;
466: PetscFunctionBegin;
467: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
468: PetscCall(DMGetPointSF(dm, &sf));
469: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &leaves, &remotes));
470: if (Nr < 0) PetscFunctionReturn(PETSC_SUCCESS);
471: PetscCall(DMViewFromOptions(dm, NULL, "-reconcile_view"));
472: PetscCall(PetscSFViewFromOptions(sf, NULL, "-reconcile_view"));
473: PetscCall(DMLabelViewFromOptions(label, NULL, "-reconcile_view"));
474: PetscCall(DMLabelGetValueIS(label, &valueIS));
475: PetscCall(ISGetLocalSize(valueIS, &Nv));
476: PetscCall(ISGetIndices(valueIS, &values));
478: PetscCall(PetscSFComputeDegreeBegin(sf, °ree));
479: PetscCall(PetscSFComputeDegreeEnd(sf, °ree));
480: PetscCall(PetscMalloc2(Nr, &rvalues, Nr, &lvalues));
481: for (PetscInt i = 0; i < Nr; ++i) rvalues[i] = -1;
482: // First set root values of shared points
483: PetscCall(ISGetIndices(valueIS, &values));
484: for (PetscInt i = 0; i < Nv; ++i) {
485: const PetscInt val = values[i];
486: IS stratumIS;
487: const PetscInt *points;
488: PetscInt Ns;
490: PetscCall(DMLabelGetStratumIS(label, val, &stratumIS));
491: if (!stratumIS) continue;
492: PetscCall(ISGetLocalSize(stratumIS, &Ns));
493: PetscCall(ISGetIndices(stratumIS, &points));
494: // Set shared points
495: for (PetscInt s = 0; s < Ns; ++s) {
496: if (degree[points[s]]) PetscCall(DMLabelGetValue(label, points[s], &rvalues[points[s]]));
497: }
498: PetscCall(ISRestoreIndices(stratumIS, &points));
499: PetscCall(ISDestroy(&stratumIS));
500: }
501: // Reduce in leaf values
502: if (reduceop != MPI_OP_NULL) {
503: for (PetscInt l = 0; l < Nl; ++l) PetscCall(DMLabelGetValue(label, leaves[l], &lvalues[leaves[l]]));
504: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, lvalues, rvalues, reduceop));
505: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, lvalues, rvalues, reduceop));
506: // Update root values
507: for (PetscInt i = 0; i < Nv; ++i) {
508: const PetscInt val = values[i];
509: IS stratumIS;
510: const PetscInt *points;
511: PetscInt Ns;
513: PetscCall(DMLabelGetStratumIS(label, val, &stratumIS));
514: if (!stratumIS) continue;
515: PetscCall(ISGetLocalSize(stratumIS, &Ns));
516: PetscCall(ISGetIndices(stratumIS, &points));
517: // Set shared points
518: for (PetscInt s = 0; s < Ns; ++s) {
519: const PetscInt point = points[s];
521: if (degree[point]) {
522: PetscInt val;
524: PetscCall(DMLabelGetValue(label, point, &val));
525: // Do not discard or push ghost designation
526: if (val != rvalues[point] && val < shift3 && rvalues[point] < shift3) {
527: PetscCall(DMLabelClearValue(label, point, val));
528: PetscCall(DMLabelSetValue(label, point, rvalues[point]));
529: }
530: }
531: }
532: PetscCall(ISRestoreIndices(stratumIS, &points));
533: PetscCall(ISDestroy(&stratumIS));
534: }
535: }
536: // Check for root values which were not labeled, but now are
537: for (PetscInt point = 0; point < Nr; ++point) {
538: if (degree[point]) {
539: PetscInt val;
541: PetscCall(DMLabelGetValue(label, point, &val));
542: if (val == -1 && val != rvalues[point]) {
543: if (newValue != PETSC_DETERMINE) PetscCall(DMLabelSetValue(label, point, newValue));
544: PetscCall(DMLabelSetValue(label, point, rvalues[point]));
545: }
546: }
547: }
548: // Broadcast root values to leaves
549: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, rvalues, lvalues, MPI_REPLACE));
550: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, rvalues, lvalues, MPI_REPLACE));
551: // Update leaf values
552: for (PetscInt l = 0; l < Nl; ++l) {
553: const PetscInt point = leaves[l];
554: PetscInt val;
556: PetscCall(DMLabelGetValue(label, point, &val));
557: // Do not discard or push ghost designation
558: if (val != lvalues[point] && val < shift3 && lvalues[point] < shift3) {
559: PetscCall(DMLabelClearValue(label, point, val));
560: PetscCall(DMLabelSetValue(label, point, lvalues[point]));
561: }
562: }
564: if (debug) {
565: PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
566: PetscViewer selfviewer;
567: const char *name;
568: PetscMPIInt rank;
570: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
571: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
572: PetscCall(PetscObjectGetName((PetscObject)label, &name));
573: PetscCallMPI(MPI_Comm_rank(comm, &rank));
574: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: Reconciled label %s\n", rank, name));
575: PetscCall(DMLabelView(label, selfviewer));
576: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
577: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
578: }
579: PetscCall(ISRestoreIndices(valueIS, &values));
580: PetscCall(ISDestroy(&valueIS));
581: PetscCall(PetscFree2(rvalues, lvalues));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: /*@
586: DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point
588: Input Parameters:
589: + dm - The `DM`
590: - label - A `DMLabel` marking the surface points
592: Output Parameter:
593: . label - A `DMLabel` incorporating cells
595: Level: developer
597: Note:
598: The cells allow FEM boundary conditions to be applied using the cell geometry
600: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelAddFaceCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
601: @*/
602: PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label)
603: {
604: IS valueIS;
605: const PetscInt *values;
606: PetscInt numValues, v, csStart, csEnd, chStart, chEnd;
608: PetscFunctionBegin;
609: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &csStart, &csEnd));
610: PetscCall(DMPlexGetHeightStratum(dm, 0, &chStart, &chEnd));
611: PetscCall(DMLabelGetNumValues(label, &numValues));
612: PetscCall(DMLabelGetValueIS(label, &valueIS));
613: PetscCall(ISGetIndices(valueIS, &values));
614: for (v = 0; v < numValues; ++v) {
615: IS pointIS;
616: const PetscInt *points;
617: PetscInt numPoints, p;
619: PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints));
620: PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS));
621: PetscCall(ISGetIndices(pointIS, &points));
622: for (p = 0; p < numPoints; ++p) {
623: const PetscInt point = points[p];
624: PetscInt *closure = NULL;
625: PetscInt closureSize, cl, h, cStart, cEnd;
626: DMPolytopeType ct;
628: // If the point is a hybrid, allow hybrid cells
629: PetscCall(DMPlexGetCellType(dm, point, &ct));
630: PetscCall(DMPlexGetPointHeight(dm, point, &h));
631: if (DMPolytopeTypeIsHybrid(ct)) {
632: cStart = chStart;
633: cEnd = chEnd;
634: } else {
635: cStart = csStart;
636: cEnd = csEnd;
637: }
639: PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure));
640: for (cl = closureSize - 1; cl > 0; --cl) {
641: const PetscInt cell = closure[cl * 2];
642: if (cell >= cStart && cell < cEnd) {
643: PetscCall(DMLabelSetValue(label, cell, values[v]));
644: break;
645: }
646: }
647: PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure));
648: }
649: PetscCall(ISRestoreIndices(pointIS, &points));
650: PetscCall(ISDestroy(&pointIS));
651: }
652: PetscCall(ISRestoreIndices(valueIS, &values));
653: PetscCall(ISDestroy(&valueIS));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*@
658: DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face
660: Input Parameters:
661: + dm - The `DM`
662: - label - A `DMLabel` marking the surface points
664: Output Parameter:
665: . label - A `DMLabel` incorporating cells
667: Level: developer
669: Note:
670: The cells allow FEM boundary conditions to be applied using the cell geometry
672: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelAddCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`
673: @*/
674: PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label)
675: {
676: IS valueIS;
677: const PetscInt *values;
678: PetscInt numValues, v, cStart, cEnd, fStart, fEnd;
680: PetscFunctionBegin;
681: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
682: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
683: PetscCall(DMLabelGetNumValues(label, &numValues));
684: PetscCall(DMLabelGetValueIS(label, &valueIS));
685: PetscCall(ISGetIndices(valueIS, &values));
686: for (v = 0; v < numValues; ++v) {
687: IS pointIS;
688: const PetscInt *points;
689: PetscInt numPoints, p;
691: PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints));
692: PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS));
693: PetscCall(ISGetIndices(pointIS, &points));
694: for (p = 0; p < numPoints; ++p) {
695: const PetscInt face = points[p];
696: PetscInt *closure = NULL;
697: PetscInt closureSize, cl;
699: if ((face < fStart) || (face >= fEnd)) continue;
700: PetscCall(DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure));
701: for (cl = closureSize - 1; cl > 0; --cl) {
702: const PetscInt cell = closure[cl * 2];
703: if (cell >= cStart && cell < cEnd) {
704: PetscCall(DMLabelSetValue(label, cell, values[v]));
705: break;
706: }
707: }
708: PetscCall(DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure));
709: }
710: PetscCall(ISRestoreIndices(pointIS, &points));
711: PetscCall(ISDestroy(&pointIS));
712: }
713: PetscCall(ISRestoreIndices(valueIS, &values));
714: PetscCall(ISDestroy(&valueIS));
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: /*@
719: DMPlexLabelClearCells - Remove cells from a label
721: Input Parameters:
722: + dm - The `DM`
723: - label - A `DMLabel` marking surface points and their adjacent cells
725: Output Parameter:
726: . label - A `DMLabel` without cells
728: Level: developer
730: Note:
731: This undoes `DMPlexLabelAddCells()` or `DMPlexLabelAddFaceCells()`
733: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`, `DMPlexLabelAddCells()`
734: @*/
735: PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label)
736: {
737: IS valueIS;
738: const PetscInt *values;
739: PetscInt numValues, v, cStart, cEnd;
741: PetscFunctionBegin;
742: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
743: PetscCall(DMLabelGetNumValues(label, &numValues));
744: PetscCall(DMLabelGetValueIS(label, &valueIS));
745: PetscCall(ISGetIndices(valueIS, &values));
746: for (v = 0; v < numValues; ++v) {
747: IS pointIS;
748: const PetscInt *points;
749: PetscInt numPoints, p;
751: PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints));
752: PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS));
753: PetscCall(ISGetIndices(pointIS, &points));
754: for (p = 0; p < numPoints; ++p) {
755: PetscInt point = points[p];
757: if (point >= cStart && point < cEnd) PetscCall(DMLabelClearValue(label, point, values[v]));
758: }
759: PetscCall(ISRestoreIndices(pointIS, &points));
760: PetscCall(ISDestroy(&pointIS));
761: }
762: PetscCall(ISRestoreIndices(valueIS, &values));
763: PetscCall(ISDestroy(&valueIS));
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending
768: * index (skipping first, which is (0,0)) */
769: static inline PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[])
770: {
771: PetscInt d, off = 0;
773: PetscFunctionBegin;
774: /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */
775: for (d = 0; d < depth; d++) {
776: PetscInt firstd = d;
777: PetscInt firstStart = depthShift[2 * d];
778: PetscInt e;
780: for (e = d + 1; e <= depth; e++) {
781: if (depthShift[2 * e] < firstStart) {
782: firstd = e;
783: firstStart = depthShift[2 * d];
784: }
785: }
786: if (firstd != d) {
787: PetscInt swap[2];
789: e = firstd;
790: swap[0] = depthShift[2 * d];
791: swap[1] = depthShift[2 * d + 1];
792: depthShift[2 * d] = depthShift[2 * e];
793: depthShift[2 * d + 1] = depthShift[2 * e + 1];
794: depthShift[2 * e] = swap[0];
795: depthShift[2 * e + 1] = swap[1];
796: }
797: }
798: /* convert (oldstart, added) to (oldstart, newstart) */
799: for (d = 0; d <= depth; d++) {
800: off += depthShift[2 * d + 1];
801: depthShift[2 * d + 1] = depthShift[2 * d] + off;
802: }
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: /* depthShift is a list of (old, new) pairs */
807: static inline PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
808: {
809: PetscInt d;
810: PetscInt newOff = 0;
812: for (d = 0; d <= depth; d++) {
813: if (p < depthShift[2 * d]) return p + newOff;
814: else newOff = depthShift[2 * d + 1] - depthShift[2 * d];
815: }
816: return p + newOff;
817: }
819: /* depthShift is a list of (old, new) pairs */
820: static inline PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[])
821: {
822: PetscInt d;
823: PetscInt newOff = 0;
825: for (d = 0; d <= depth; d++) {
826: if (p < depthShift[2 * d + 1]) return p + newOff;
827: else newOff = depthShift[2 * d] - depthShift[2 * d + 1];
828: }
829: return p + newOff;
830: }
832: static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew)
833: {
834: PetscInt depth = 0, d, pStart, pEnd, p;
835: DMLabel depthLabel;
837: PetscFunctionBegin;
838: PetscCall(DMPlexGetDepth(dm, &depth));
839: if (depth < 0) PetscFunctionReturn(PETSC_SUCCESS);
840: /* Step 1: Expand chart */
841: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
842: pEnd = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
843: PetscCall(DMPlexSetChart(dmNew, pStart, pEnd));
844: PetscCall(DMCreateLabel(dmNew, "depth"));
845: PetscCall(DMPlexGetDepthLabel(dmNew, &depthLabel));
846: PetscCall(DMCreateLabel(dmNew, "celltype"));
847: /* Step 2: Set cone and support sizes */
848: for (d = 0; d <= depth; ++d) {
849: PetscInt pStartNew, pEndNew;
850: IS pIS;
852: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
853: pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
854: pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
855: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS));
856: PetscCall(DMLabelSetStratumIS(depthLabel, d, pIS));
857: PetscCall(ISDestroy(&pIS));
858: for (p = pStart; p < pEnd; ++p) {
859: PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
860: PetscInt size;
861: DMPolytopeType ct;
863: PetscCall(DMPlexGetConeSize(dm, p, &size));
864: PetscCall(DMPlexSetConeSize(dmNew, newp, size));
865: PetscCall(DMPlexGetSupportSize(dm, p, &size));
866: PetscCall(DMPlexSetSupportSize(dmNew, newp, size));
867: PetscCall(DMPlexGetCellType(dm, p, &ct));
868: PetscCall(DMPlexSetCellType(dmNew, newp, ct));
869: }
870: }
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew)
875: {
876: PetscInt *newpoints;
877: PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p;
879: PetscFunctionBegin;
880: PetscCall(DMPlexGetDepth(dm, &depth));
881: if (depth < 0) PetscFunctionReturn(PETSC_SUCCESS);
882: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
883: PetscCall(DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew));
884: PetscCall(PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)), &newpoints));
885: /* Step 5: Set cones and supports */
886: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
887: for (p = pStart; p < pEnd; ++p) {
888: const PetscInt *points = NULL, *orientations = NULL;
889: PetscInt size, sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift);
891: PetscCall(DMPlexGetConeSize(dm, p, &size));
892: PetscCall(DMPlexGetCone(dm, p, &points));
893: PetscCall(DMPlexGetConeOrientation(dm, p, &orientations));
894: for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
895: PetscCall(DMPlexSetCone(dmNew, newp, newpoints));
896: PetscCall(DMPlexSetConeOrientation(dmNew, newp, orientations));
897: PetscCall(DMPlexGetSupportSize(dm, p, &size));
898: PetscCall(DMPlexGetSupportSize(dmNew, newp, &sizeNew));
899: PetscCall(DMPlexGetSupport(dm, p, &points));
900: for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift);
901: for (i = size; i < sizeNew; ++i) newpoints[i] = 0;
902: PetscCall(DMPlexSetSupport(dmNew, newp, newpoints));
903: }
904: PetscCall(PetscFree(newpoints));
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew)
909: {
910: PetscSection coordSection, newCoordSection;
911: Vec coordinates, newCoordinates;
912: PetscScalar *coords, *newCoords;
913: PetscInt coordSize, sStart, sEnd;
914: PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v;
915: PetscBool hasCells;
917: PetscFunctionBegin;
918: PetscCall(DMGetCoordinateDim(dm, &dim));
919: PetscCall(DMSetCoordinateDim(dmNew, dim));
920: PetscCall(DMPlexGetDepth(dm, &depth));
921: /* Step 8: Convert coordinates */
922: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
923: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
924: PetscCall(DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew));
925: PetscCall(DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew));
926: PetscCall(DMGetCoordinateSection(dm, &coordSection));
927: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection));
928: PetscCall(PetscSectionSetNumFields(newCoordSection, 1));
929: PetscCall(PetscSectionSetFieldComponents(newCoordSection, 0, dim));
930: PetscCall(PetscSectionGetChart(coordSection, &sStart, &sEnd));
931: hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE;
932: PetscCall(PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew));
933: if (hasCells) {
934: for (c = cStart; c < cEnd; ++c) {
935: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof;
937: PetscCall(PetscSectionGetDof(coordSection, c, &dof));
938: PetscCall(PetscSectionSetDof(newCoordSection, cNew, dof));
939: PetscCall(PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof));
940: }
941: }
942: for (v = vStartNew; v < vEndNew; ++v) {
943: PetscCall(PetscSectionSetDof(newCoordSection, v, dim));
944: PetscCall(PetscSectionSetFieldDof(newCoordSection, v, 0, dim));
945: }
946: PetscCall(PetscSectionSetUp(newCoordSection));
947: PetscCall(DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection));
948: PetscCall(PetscSectionGetStorageSize(newCoordSection, &coordSize));
949: PetscCall(VecCreate(PETSC_COMM_SELF, &newCoordinates));
950: PetscCall(PetscObjectSetName((PetscObject)newCoordinates, "coordinates"));
951: PetscCall(VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE));
952: PetscCall(VecSetBlockSize(newCoordinates, dim));
953: PetscCall(VecSetType(newCoordinates, VECSTANDARD));
954: PetscCall(DMSetCoordinatesLocal(dmNew, newCoordinates));
955: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
956: PetscCall(VecGetArray(coordinates, &coords));
957: PetscCall(VecGetArray(newCoordinates, &newCoords));
958: if (hasCells) {
959: for (c = cStart; c < cEnd; ++c) {
960: PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d;
962: PetscCall(PetscSectionGetDof(coordSection, c, &dof));
963: PetscCall(PetscSectionGetOffset(coordSection, c, &off));
964: PetscCall(PetscSectionGetOffset(newCoordSection, cNew, &noff));
965: for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
966: }
967: }
968: for (v = vStart; v < vEnd; ++v) {
969: PetscInt dof, off, noff, d;
971: PetscCall(PetscSectionGetDof(coordSection, v, &dof));
972: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
973: PetscCall(PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff));
974: for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d];
975: }
976: PetscCall(VecRestoreArray(coordinates, &coords));
977: PetscCall(VecRestoreArray(newCoordinates, &newCoords));
978: PetscCall(VecDestroy(&newCoordinates));
979: PetscCall(PetscSectionDestroy(&newCoordSection));
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: static PetscErrorCode DMPlexShiftSF_Single(DM dm, PetscInt depthShift[], PetscSF sf, PetscSF sfNew)
984: {
985: const PetscSFNode *remotePoints;
986: PetscSFNode *gremotePoints;
987: const PetscInt *localPoints;
988: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
989: PetscInt numRoots, numLeaves, l, pStart, pEnd, depth = 0, totShift = 0;
991: PetscFunctionBegin;
992: PetscCall(DMPlexGetDepth(dm, &depth));
993: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
994: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints));
995: totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd;
996: if (numRoots >= 0) {
997: PetscCall(PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation));
998: for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift);
999: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
1000: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
1001: PetscCall(PetscMalloc1(numLeaves, &glocalPoints));
1002: PetscCall(PetscMalloc1(numLeaves, &gremotePoints));
1003: for (l = 0; l < numLeaves; ++l) {
1004: glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift);
1005: gremotePoints[l].rank = remotePoints[l].rank;
1006: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
1007: }
1008: PetscCall(PetscFree2(newLocation, newRemoteLocation));
1009: PetscCall(PetscSFSetGraph(sfNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER));
1010: }
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew)
1015: {
1016: PetscSF sfPoint, sfPointNew;
1017: PetscBool useNatural;
1019: PetscFunctionBegin;
1020: /* Step 9: Convert pointSF */
1021: PetscCall(DMGetPointSF(dm, &sfPoint));
1022: PetscCall(DMGetPointSF(dmNew, &sfPointNew));
1023: PetscCall(DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew));
1024: /* Step 9b: Convert naturalSF */
1025: PetscCall(DMGetUseNatural(dm, &useNatural));
1026: if (useNatural) {
1027: PetscSF sfNat, sfNatNew;
1029: PetscCall(DMSetUseNatural(dmNew, useNatural));
1030: PetscCall(DMGetNaturalSF(dm, &sfNat));
1031: PetscCall(DMGetNaturalSF(dmNew, &sfNatNew));
1032: PetscCall(DMPlexShiftSF_Single(dm, depthShift, sfNat, sfNatNew));
1033: }
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew)
1038: {
1039: PetscInt depth = 0, numLabels, l;
1041: PetscFunctionBegin;
1042: PetscCall(DMPlexGetDepth(dm, &depth));
1043: /* Step 10: Convert labels */
1044: PetscCall(DMGetNumLabels(dm, &numLabels));
1045: for (l = 0; l < numLabels; ++l) {
1046: DMLabel label, newlabel;
1047: const char *lname;
1048: PetscBool isDepth, isDim;
1049: IS valueIS;
1050: const PetscInt *values;
1051: PetscInt numValues, val;
1053: PetscCall(DMGetLabelName(dm, l, &lname));
1054: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
1055: if (isDepth) continue;
1056: PetscCall(PetscStrcmp(lname, "dim", &isDim));
1057: if (isDim) continue;
1058: PetscCall(DMCreateLabel(dmNew, lname));
1059: PetscCall(DMGetLabel(dm, lname, &label));
1060: PetscCall(DMGetLabel(dmNew, lname, &newlabel));
1061: PetscCall(DMLabelGetDefaultValue(label, &val));
1062: PetscCall(DMLabelSetDefaultValue(newlabel, val));
1063: PetscCall(DMLabelGetValueIS(label, &valueIS));
1064: PetscCall(ISGetLocalSize(valueIS, &numValues));
1065: PetscCall(ISGetIndices(valueIS, &values));
1066: for (val = 0; val < numValues; ++val) {
1067: IS pointIS;
1068: const PetscInt *points;
1069: PetscInt numPoints, p;
1071: PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS));
1072: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1073: PetscCall(ISGetIndices(pointIS, &points));
1074: for (p = 0; p < numPoints; ++p) {
1075: const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift);
1077: PetscCall(DMLabelSetValue(newlabel, newpoint, values[val]));
1078: }
1079: PetscCall(ISRestoreIndices(pointIS, &points));
1080: PetscCall(ISDestroy(&pointIS));
1081: }
1082: PetscCall(ISRestoreIndices(valueIS, &values));
1083: PetscCall(ISDestroy(&valueIS));
1084: }
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: static PetscErrorCode DMPlexCreateVTKLabel_Internal(DM dm, PetscBool createGhostLabel, DM dmNew)
1089: {
1090: PetscSF sfPoint;
1091: DMLabel vtkLabel, ghostLabel = NULL;
1092: const PetscSFNode *leafRemote;
1093: const PetscInt *leafLocal;
1094: PetscInt cellHeight, cStart, cEnd, c, fStart, fEnd, f, numLeaves, l;
1095: PetscMPIInt rank;
1097: PetscFunctionBegin;
1098: /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */
1099: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1100: PetscCall(DMGetPointSF(dm, &sfPoint));
1101: PetscCall(DMPlexGetVTKCellHeight(dmNew, &cellHeight));
1102: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
1103: PetscCall(PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote));
1104: PetscCall(DMCreateLabel(dmNew, "vtk"));
1105: PetscCall(DMGetLabel(dmNew, "vtk", &vtkLabel));
1106: if (createGhostLabel) {
1107: PetscCall(DMCreateLabel(dmNew, "ghost"));
1108: PetscCall(DMGetLabel(dmNew, "ghost", &ghostLabel));
1109: }
1110: for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) {
1111: for (; c < leafLocal[l] && c < cEnd; ++c) PetscCall(DMLabelSetValue(vtkLabel, c, 1));
1112: if (leafLocal[l] >= cEnd) break;
1113: if (leafRemote[l].rank == rank) PetscCall(DMLabelSetValue(vtkLabel, c, 1));
1114: else if (ghostLabel) PetscCall(DMLabelSetValue(ghostLabel, c, 2));
1115: }
1116: for (; c < cEnd; ++c) PetscCall(DMLabelSetValue(vtkLabel, c, 1));
1117: if (ghostLabel) {
1118: PetscCall(DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd));
1119: for (f = fStart; f < fEnd; ++f) {
1120: PetscInt numCells;
1122: PetscCall(DMPlexGetSupportSize(dmNew, f, &numCells));
1123: if (numCells < 2) PetscCall(DMLabelSetValue(ghostLabel, f, 1));
1124: else {
1125: const PetscInt *cells = NULL;
1126: PetscInt vA, vB;
1128: PetscCall(DMPlexGetSupport(dmNew, f, &cells));
1129: PetscCall(DMLabelGetValue(vtkLabel, cells[0], &vA));
1130: PetscCall(DMLabelGetValue(vtkLabel, cells[1], &vB));
1131: if (vA != 1 && vB != 1) PetscCall(DMLabelSetValue(ghostLabel, f, 1));
1132: }
1133: }
1134: }
1135: PetscFunctionReturn(PETSC_SUCCESS);
1136: }
1138: static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew)
1139: {
1140: DM refTree;
1141: PetscSection pSec;
1142: PetscInt *parents, *childIDs;
1144: PetscFunctionBegin;
1145: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
1146: PetscCall(DMPlexSetReferenceTree(dmNew, refTree));
1147: PetscCall(DMPlexGetTree(dm, &pSec, &parents, &childIDs, NULL, NULL));
1148: if (pSec) {
1149: PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth;
1150: PetscInt *childIDsShifted;
1151: PetscSection pSecShifted;
1153: PetscCall(PetscSectionGetChart(pSec, &pStart, &pEnd));
1154: PetscCall(DMPlexGetDepth(dm, &depth));
1155: pStartShifted = DMPlexShiftPoint_Internal(pStart, depth, depthShift);
1156: pEndShifted = DMPlexShiftPoint_Internal(pEnd, depth, depthShift);
1157: PetscCall(PetscMalloc2(pEndShifted - pStartShifted, &parentsShifted, pEndShifted - pStartShifted, &childIDsShifted));
1158: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmNew), &pSecShifted));
1159: PetscCall(PetscSectionSetChart(pSecShifted, pStartShifted, pEndShifted));
1160: for (p = pStartShifted; p < pEndShifted; p++) {
1161: /* start off assuming no children */
1162: PetscCall(PetscSectionSetDof(pSecShifted, p, 0));
1163: }
1164: for (p = pStart; p < pEnd; p++) {
1165: PetscInt dof;
1166: PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);
1168: PetscCall(PetscSectionGetDof(pSec, p, &dof));
1169: PetscCall(PetscSectionSetDof(pSecShifted, pNew, dof));
1170: }
1171: PetscCall(PetscSectionSetUp(pSecShifted));
1172: for (p = pStart; p < pEnd; p++) {
1173: PetscInt dof;
1174: PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift);
1176: PetscCall(PetscSectionGetDof(pSec, p, &dof));
1177: if (dof) {
1178: PetscInt off, offNew;
1180: PetscCall(PetscSectionGetOffset(pSec, p, &off));
1181: PetscCall(PetscSectionGetOffset(pSecShifted, pNew, &offNew));
1182: parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off], depth, depthShift);
1183: childIDsShifted[offNew] = childIDs[off];
1184: }
1185: }
1186: PetscCall(DMPlexSetTree(dmNew, pSecShifted, parentsShifted, childIDsShifted));
1187: PetscCall(PetscFree2(parentsShifted, childIDsShifted));
1188: PetscCall(PetscSectionDestroy(&pSecShifted));
1189: }
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm)
1194: {
1195: PetscSF sf;
1196: IS valueIS;
1197: const PetscInt *values, *leaves;
1198: PetscInt *depthShift;
1199: PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c;
1200: const PetscReal *maxCell, *Lstart, *L;
1202: PetscFunctionBegin;
1203: PetscCall(DMGetPointSF(dm, &sf));
1204: PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
1205: nleaves = PetscMax(0, nleaves);
1206: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
1207: /* Count ghost cells */
1208: PetscCall(DMLabelGetValueIS(label, &valueIS));
1209: PetscCall(ISGetLocalSize(valueIS, &numFS));
1210: PetscCall(ISGetIndices(valueIS, &values));
1211: Ng = 0;
1212: for (fs = 0; fs < numFS; ++fs) {
1213: IS faceIS;
1214: const PetscInt *faces;
1215: PetscInt numFaces, f, numBdFaces = 0;
1217: PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS));
1218: PetscCall(ISGetLocalSize(faceIS, &numFaces));
1219: PetscCall(ISGetIndices(faceIS, &faces));
1220: for (f = 0; f < numFaces; ++f) {
1221: PetscInt numChildren;
1223: PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc));
1224: PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL));
1225: /* non-local and ancestors points don't get to register ghosts */
1226: if (loc >= 0 || numChildren) continue;
1227: if (faces[f] >= fStart && faces[f] < fEnd) ++numBdFaces;
1228: }
1229: Ng += numBdFaces;
1230: PetscCall(ISRestoreIndices(faceIS, &faces));
1231: PetscCall(ISDestroy(&faceIS));
1232: }
1233: PetscCall(DMPlexGetDepth(dm, &depth));
1234: PetscCall(PetscMalloc1(2 * (depth + 1), &depthShift));
1235: for (d = 0; d <= depth; d++) {
1236: PetscInt dEnd;
1238: PetscCall(DMPlexGetDepthStratum(dm, d, NULL, &dEnd));
1239: depthShift[2 * d] = dEnd;
1240: depthShift[2 * d + 1] = 0;
1241: }
1242: if (depth >= 0) depthShift[2 * depth + 1] = Ng;
1243: PetscCall(DMPlexShiftPointSetUp_Internal(depth, depthShift));
1244: PetscCall(DMPlexShiftSizes_Internal(dm, depthShift, gdm));
1245: /* Step 3: Set cone/support sizes for new points */
1246: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
1247: for (c = cEnd; c < cEnd + Ng; ++c) PetscCall(DMPlexSetConeSize(gdm, c, 1));
1248: for (fs = 0; fs < numFS; ++fs) {
1249: IS faceIS;
1250: const PetscInt *faces;
1251: PetscInt numFaces, f;
1253: PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS));
1254: PetscCall(ISGetLocalSize(faceIS, &numFaces));
1255: PetscCall(ISGetIndices(faceIS, &faces));
1256: for (f = 0; f < numFaces; ++f) {
1257: PetscInt size, numChildren;
1259: PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc));
1260: PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL));
1261: if (loc >= 0 || numChildren) continue;
1262: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
1263: PetscCall(DMPlexGetSupportSize(dm, faces[f], &size));
1264: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %" PetscInt_FMT " with %" PetscInt_FMT " support cells", faces[f], size);
1265: PetscCall(DMPlexSetSupportSize(gdm, faces[f] + Ng, 2));
1266: }
1267: PetscCall(ISRestoreIndices(faceIS, &faces));
1268: PetscCall(ISDestroy(&faceIS));
1269: }
1270: /* Step 4: Setup ghosted DM */
1271: PetscCall(DMSetUp(gdm));
1272: PetscCall(DMPlexShiftPoints_Internal(dm, depthShift, gdm));
1273: /* Step 6: Set cones and supports for new points */
1274: ghostCell = cEnd;
1275: for (fs = 0; fs < numFS; ++fs) {
1276: IS faceIS;
1277: const PetscInt *faces;
1278: PetscInt numFaces, f;
1280: PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS));
1281: PetscCall(ISGetLocalSize(faceIS, &numFaces));
1282: PetscCall(ISGetIndices(faceIS, &faces));
1283: for (f = 0; f < numFaces; ++f) {
1284: PetscInt newFace = faces[f] + Ng, numChildren;
1286: PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc));
1287: PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL));
1288: if (loc >= 0 || numChildren) continue;
1289: if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue;
1290: PetscCall(DMPlexSetCone(gdm, ghostCell, &newFace));
1291: PetscCall(DMPlexInsertSupport(gdm, newFace, 1, ghostCell));
1292: ++ghostCell;
1293: }
1294: PetscCall(ISRestoreIndices(faceIS, &faces));
1295: PetscCall(ISDestroy(&faceIS));
1296: }
1297: PetscCall(ISRestoreIndices(valueIS, &values));
1298: PetscCall(ISDestroy(&valueIS));
1299: PetscCall(DMPlexShiftCoordinates_Internal(dm, depthShift, gdm));
1300: PetscCall(DMPlexShiftSF_Internal(dm, depthShift, gdm));
1301: PetscCall(DMPlexShiftLabels_Internal(dm, depthShift, gdm));
1302: PetscCall(DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm));
1303: PetscCall(DMPlexShiftTree_Internal(dm, depthShift, gdm));
1304: PetscCall(PetscFree(depthShift));
1305: for (c = cEnd; c < cEnd + Ng; ++c) PetscCall(DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST));
1306: /* Step 7: Periodicity */
1307: PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L));
1308: PetscCall(DMSetPeriodicity(gdm, maxCell, Lstart, L));
1309: if (numGhostCells) *numGhostCells = Ng;
1310: PetscFunctionReturn(PETSC_SUCCESS);
1311: }
1313: /*@
1314: DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face
1316: Collective
1318: Input Parameters:
1319: + dm - The original `DM`
1320: - labelName - The label specifying the boundary faces, or "Face Sets" if this is `NULL`
1322: Output Parameters:
1323: + numGhostCells - The number of ghost cells added to the `DM`
1324: - dmGhosted - The new `DM`
1326: Level: developer
1328: Note:
1329: If no label exists of that name, one will be created marking all boundary faces
1331: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1332: @*/
1333: PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted)
1334: {
1335: DM gdm;
1336: DMLabel label;
1337: const char *name = labelName ? labelName : "Face Sets";
1338: PetscInt dim, Ng = 0;
1339: PetscBool useCone, useClosure;
1341: PetscFunctionBegin;
1343: if (numGhostCells) PetscAssertPointer(numGhostCells, 3);
1344: PetscAssertPointer(dmGhosted, 4);
1345: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &gdm));
1346: PetscCall(DMSetType(gdm, DMPLEX));
1347: PetscCall(DMGetDimension(dm, &dim));
1348: PetscCall(DMSetDimension(gdm, dim));
1349: PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
1350: PetscCall(DMSetBasicAdjacency(gdm, useCone, useClosure));
1351: PetscCall(DMGetLabel(dm, name, &label));
1352: if (!label) {
1353: /* Get label for boundary faces */
1354: PetscCall(DMCreateLabel(dm, name));
1355: PetscCall(DMGetLabel(dm, name, &label));
1356: PetscCall(DMPlexMarkBoundaryFaces(dm, 1, label));
1357: }
1358: PetscCall(DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm));
1359: PetscCall(DMCopyDisc(dm, gdm));
1360: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, gdm));
1361: gdm->setfromoptionscalled = dm->setfromoptionscalled;
1362: if (numGhostCells) *numGhostCells = Ng;
1363: *dmGhosted = gdm;
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: static PetscErrorCode DivideCells_Private(DM dm, DMLabel label, DMPlexPointQueue queue)
1368: {
1369: const PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
1370: const PetscInt shift = 100;
1371: PetscInt dim, d, *pStart, *pEnd;
1372: MPI_Comm comm;
1373: PetscMPIInt rank;
1375: PetscFunctionBegin;
1376: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1377: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1378: PetscCall(DMGetDimension(dm, &dim));
1379: PetscCall(PetscMalloc2(dim, &pStart, dim, &pEnd));
1380: for (d = 0; d < dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]));
1381: while (!DMPlexPointQueueEmpty(queue)) {
1382: PetscInt cell = -1;
1383: PetscInt *closure = NULL;
1384: PetscInt closureSize, cl, cval, depth;
1386: PetscCall(DMPlexPointQueueDequeue(queue, &cell));
1387: PetscCall(DMLabelGetValue(label, cell, &cval));
1388: PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
1389: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]Dequeued fault cell %" PetscInt_FMT " depth %" PetscInt_FMT "\n", rank, cell, depth));
1390: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1391: /* Mark points in the cell closure that touch the fault */
1392: for (d = 0; d < dim; ++d) {
1393: for (cl = 0; cl < closureSize * 2; cl += 2) {
1394: const PetscInt clp = closure[cl];
1395: PetscInt clval;
1397: if ((clp < pStart[d]) || (clp >= pEnd[d])) continue;
1398: PetscCall(DMLabelGetValue(label, clp, &clval));
1399: if (clval == -1) {
1400: const PetscInt *cone;
1401: PetscInt coneSize, c;
1403: /* If a cone point touches the fault, then this point touches the fault */
1404: PetscCall(DMPlexGetCone(dm, clp, &cone));
1405: PetscCall(DMPlexGetConeSize(dm, clp, &coneSize));
1406: for (c = 0; c < coneSize; ++c) {
1407: PetscInt cpval;
1409: PetscCall(DMLabelGetValue(label, cone[c], &cpval));
1410: if (cpval != -1) {
1411: PetscInt dep;
1413: PetscCall(DMPlexGetPointDepth(dm, clp, &dep));
1414: clval = cval < 0 ? -(shift + dep) : shift + dep;
1415: PetscCall(DMLabelSetValue(label, clp, clval));
1416: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d]Labeled point %" PetscInt_FMT " with value %" PetscInt_FMT "\n", rank, clp, clval));
1417: break;
1418: }
1419: }
1420: }
1421: /* Mark neighbor cells through marked faces (these cells must also touch the fault) */
1422: if (d == dim - 1 && clval != -1) {
1423: const PetscInt *support;
1424: PetscInt supportSize, s, nval;
1426: PetscCall(DMPlexGetSupport(dm, clp, &support));
1427: PetscCall(DMPlexGetSupportSize(dm, clp, &supportSize));
1428: for (s = 0; s < supportSize; ++s) {
1429: PetscCall(DMLabelGetValue(label, support[s], &nval));
1430: if (nval == -1) {
1431: PetscCall(DMLabelSetValue(label, support[s], clval < 0 ? clval - 1 : clval + 1));
1432: PetscCall(DMPlexPointQueueEnqueue(queue, support[s]));
1433: if (debug) {
1434: PetscCall(PetscSynchronizedPrintf(comm, "[%d]Neighbor: Labeled point %" PetscInt_FMT " with value %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, support[s], clval < 0 ? clval - 1 : clval + 1, clval));
1435: PetscCall(PetscSynchronizedPrintf(comm, "[%d]Enqueued point %" PetscInt_FMT "\n", rank, support[s]));
1436: }
1437: }
1438: }
1439: }
1440: }
1441: }
1442: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
1443: }
1444: PetscCall(PetscFree2(pStart, pEnd));
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: typedef struct {
1449: DM dm;
1450: DMPlexPointQueue queue;
1451: } PointDivision;
1453: static PetscErrorCode divideCell(DMLabel label, PetscInt p, PetscInt val, PetscCtx ctx)
1454: {
1455: PointDivision *div = (PointDivision *)ctx;
1456: PetscInt cval = val < 0 ? val - 1 : val + 1;
1457: const PetscInt debug = ((DM_Plex *)div->dm->data)->printCohesive;
1458: const PetscInt *support;
1459: PetscInt supportSize, s;
1460: PetscMPIInt rank;
1462: PetscFunctionBegin;
1463: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)div->dm), &rank));
1464: if (cval < 100) {
1465: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)div->dm), "[%d]DivideCell: Ignored fault point %" PetscInt_FMT " with value %" PetscInt_FMT "\n", rank, p, cval));
1466: PetscFunctionReturn(PETSC_SUCCESS);
1467: }
1468: PetscCall(DMPlexGetSupport(div->dm, p, &support));
1469: PetscCall(DMPlexGetSupportSize(div->dm, p, &supportSize));
1470: for (s = 0; s < supportSize; ++s) {
1471: PetscCall(DMLabelSetValue(label, support[s], cval));
1472: PetscCall(DMPlexPointQueueEnqueue(div->queue, support[s]));
1473: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)div->dm), "[%d]DivideCell: Enqueued point %" PetscInt_FMT " from parent %" PetscInt_FMT " with value %" PetscInt_FMT "\n", rank, support[s], p, cval));
1474: }
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: /* Mark cells by label propagation */
1479: static PetscErrorCode DMPlexLabelFaultHalo(DM dm, DMLabel faultLabel)
1480: {
1481: const PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
1482: const PetscInt shift = 100;
1483: DMPlexPointQueue queue = NULL;
1484: PointDivision div;
1485: PetscSF pointSF;
1486: IS pointIS;
1487: const PetscInt *points;
1488: PetscBool empty;
1489: PetscInt dim, n;
1490: PetscMPIInt rank;
1492: PetscFunctionBegin;
1493: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1494: PetscCall(DMGetDimension(dm, &dim));
1495: PetscCall(DMPlexPointQueueCreate(1024, &queue));
1496: div.dm = dm;
1497: div.queue = queue;
1498: /* Enqueue cells on fault */
1499: PetscCall(DMLabelGetStratumIS(faultLabel, shift + dim, &pointIS));
1500: if (pointIS) {
1501: PetscCall(ISGetLocalSize(pointIS, &n));
1502: PetscCall(ISGetIndices(pointIS, &points));
1503: for (PetscInt i = 0; i < n; ++i) {
1504: PetscCall(DMPlexPointQueueEnqueue(queue, points[i]));
1505: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]Init: Enqueued point %" PetscInt_FMT " with value %" PetscInt_FMT "\n", rank, points[i], shift + dim));
1506: }
1507: PetscCall(ISRestoreIndices(pointIS, &points));
1508: PetscCall(ISDestroy(&pointIS));
1509: }
1510: PetscCall(DMLabelGetStratumIS(faultLabel, -(shift + dim), &pointIS));
1511: if (pointIS) {
1512: PetscCall(ISGetLocalSize(pointIS, &n));
1513: PetscCall(ISGetIndices(pointIS, &points));
1514: for (PetscInt i = 0; i < n; ++i) {
1515: PetscCall(DMPlexPointQueueEnqueue(queue, points[i]));
1516: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]Init: Enqueued point %" PetscInt_FMT " with value %" PetscInt_FMT "\n", rank, points[i], -(shift + dim)));
1517: }
1518: PetscCall(ISRestoreIndices(pointIS, &points));
1519: PetscCall(ISDestroy(&pointIS));
1520: }
1521: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1523: PetscCall(DMGetPointSF(dm, &pointSF));
1524: PetscCall(DMLabelPropagateBegin(faultLabel, pointSF));
1525: /* While edge queue is not empty: */
1526: PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
1527: while (!empty) {
1528: PetscCall(DivideCells_Private(dm, faultLabel, queue));
1529: PetscCall(DMLabelPropagatePush(faultLabel, pointSF, divideCell, &div));
1530: PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty));
1531: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1532: }
1533: PetscCall(DMLabelPropagateEnd(faultLabel, pointSF));
1534: PetscCall(DMPlexPointQueueDestroy(&queue));
1535: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
1536: PetscFunctionReturn(PETSC_SUCCESS);
1537: }
1539: /*
1540: We are adding three kinds of points here:
1541: Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault
1542: Non-replicated: Points which exist on the fault, but are not replicated
1543: Ghost: These are shared fault faces which are not owned by this process. These do not produce hybrid cells and do not replicate
1544: Hybrid: Entirely new points, such as cohesive cells
1546: When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at
1547: each depth so that the new split/hybrid points can be inserted as a block.
1548: */
1549: static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm)
1550: {
1551: MPI_Comm comm;
1552: IS valueIS;
1553: PetscInt numSP = 0; /* The number of depths for which we have replicated points */
1554: const PetscInt *values; /* List of depths for which we have replicated points */
1555: IS *splitIS;
1556: IS *unsplitIS;
1557: IS ghostIS;
1558: PetscInt *numSplitPoints; /* The number of replicated points at each depth */
1559: PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */
1560: PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */
1561: PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */
1562: PetscInt numGhostPoints; /* The number of unowned, shared fault faces */
1563: const PetscInt **splitPoints; /* Replicated points for each depth */
1564: const PetscInt **unsplitPoints; /* Non-replicated points for each depth */
1565: const PetscInt *ghostPoints; /* Ghost fault faces */
1566: PetscSection coordSection;
1567: Vec coordinates;
1568: PetscScalar *coords;
1569: PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */
1570: PetscInt *depthEnd; /* The point limit at each depth in the original mesh */
1571: PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */
1572: PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */
1573: PetscInt *coneNew, *coneONew, *supportNew;
1574: PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v;
1576: PetscFunctionBegin;
1577: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1578: PetscCall(DMGetDimension(dm, &dim));
1579: PetscCall(DMPlexGetDepth(dm, &depth));
1580: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1581: /* We do not want this label automatically computed, instead we compute it here */
1582: PetscCall(DMCreateLabel(sdm, "celltype"));
1583: /* Count split points and add cohesive cells */
1584: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
1585: PetscCall(PetscMalloc5(depth + 1, &depthMax, depth + 1, &depthEnd, 2 * (depth + 1), &depthShift, depth + 1, &pMaxNew, depth + 1, &numHybridPointsOld));
1586: PetscCall(PetscMalloc7(depth + 1, &splitIS, depth + 1, &unsplitIS, depth + 1, &numSplitPoints, depth + 1, &numUnsplitPoints, depth + 1, &numHybridPoints, depth + 1, &splitPoints, depth + 1, &unsplitPoints));
1587: for (d = 0; d <= depth; ++d) {
1588: PetscCall(DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]));
1589: PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL));
1590: depthEnd[d] = pMaxNew[d];
1591: depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d];
1592: numSplitPoints[d] = 0;
1593: numUnsplitPoints[d] = 0;
1594: numHybridPoints[d] = 0;
1595: numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d];
1596: splitPoints[d] = NULL;
1597: unsplitPoints[d] = NULL;
1598: splitIS[d] = NULL;
1599: unsplitIS[d] = NULL;
1600: /* we are shifting the existing hybrid points with the stratum behind them, so
1601: * the split comes at the end of the normal points, i.e., at depthMax[d] */
1602: depthShift[2 * d] = depthMax[d];
1603: depthShift[2 * d + 1] = 0;
1604: }
1605: numGhostPoints = 0;
1606: ghostPoints = NULL;
1607: if (label) {
1608: PetscCall(DMLabelGetValueIS(label, &valueIS));
1609: PetscCall(ISGetLocalSize(valueIS, &numSP));
1610: PetscCall(ISGetIndices(valueIS, &values));
1611: }
1612: for (sp = 0; sp < numSP; ++sp) {
1613: const PetscInt dep = values[sp];
1615: if ((dep < 0) || (dep > depth)) continue;
1616: PetscCall(DMLabelGetStratumIS(label, dep, &splitIS[dep]));
1617: if (splitIS[dep]) {
1618: PetscCall(ISGetLocalSize(splitIS[dep], &numSplitPoints[dep]));
1619: PetscCall(ISGetIndices(splitIS[dep], &splitPoints[dep]));
1620: }
1621: PetscCall(DMLabelGetStratumIS(label, shift2 + dep, &unsplitIS[dep]));
1622: if (unsplitIS[dep]) {
1623: PetscCall(ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep]));
1624: PetscCall(ISGetIndices(unsplitIS[dep], &unsplitPoints[dep]));
1625: }
1626: }
1627: PetscCall(DMLabelGetStratumIS(label, shift2 + dim - 1, &ghostIS));
1628: if (ghostIS) {
1629: PetscCall(ISGetLocalSize(ghostIS, &numGhostPoints));
1630: PetscCall(ISGetIndices(ghostIS, &ghostPoints));
1631: }
1632: /* Calculate number of hybrid points */
1633: 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 */
1634: for (d = 0; d <= depth; ++d) depthShift[2 * d + 1] = numSplitPoints[d] + numHybridPoints[d];
1635: PetscCall(DMPlexShiftPointSetUp_Internal(depth, depthShift));
1636: /* the end of the points in this stratum that come before the new points:
1637: * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly
1638: * added points */
1639: for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d], depth, depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]);
1640: PetscCall(DMPlexShiftSizes_Internal(dm, depthShift, sdm));
1641: /* Step 3: Set cone/support sizes for new points */
1642: for (dep = 0; dep <= depth; ++dep) {
1643: for (p = 0; p < numSplitPoints[dep]; ++p) {
1644: const PetscInt oldp = splitPoints[dep][p];
1645: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1646: const PetscInt splitp = p + pMaxNew[dep];
1647: const PetscInt *support;
1648: DMPolytopeType ct;
1649: PetscInt coneSize, supportSize, qf, qn, qp, e;
1651: PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize));
1652: PetscCall(DMPlexSetConeSize(sdm, splitp, coneSize));
1653: PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize));
1654: PetscCall(DMPlexSetSupportSize(sdm, splitp, supportSize));
1655: PetscCall(DMPlexGetCellType(dm, oldp, &ct));
1656: PetscCall(DMPlexSetCellType(sdm, splitp, ct));
1657: if (dep == depth - 1) {
1658: const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1660: /* Add cohesive cells, they are prisms */
1661: PetscCall(DMPlexSetConeSize(sdm, hybcell, 2 + coneSize));
1662: switch (coneSize) {
1663: case 2:
1664: PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR));
1665: break;
1666: case 3:
1667: PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR));
1668: break;
1669: case 4:
1670: PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR));
1671: break;
1672: }
1673: /* Shared fault faces with only one support cell now have two with the cohesive cell */
1674: /* TODO Check thaat oldp has rootdegree == 1 */
1675: if (supportSize == 1) {
1676: const PetscInt *support;
1677: PetscInt val;
1679: PetscCall(DMPlexGetSupport(dm, oldp, &support));
1680: PetscCall(DMLabelGetValue(label, support[0], &val));
1681: if (val < 0) PetscCall(DMPlexSetSupportSize(sdm, splitp, 2));
1682: else PetscCall(DMPlexSetSupportSize(sdm, newp, 2));
1683: }
1684: } else if (dep == 0) {
1685: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1687: PetscCall(DMPlexGetSupport(dm, oldp, &support));
1688: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1689: PetscInt val;
1691: PetscCall(DMLabelGetValue(label, support[e], &val));
1692: if (val == 1) ++qf;
1693: if ((val == 1) || (val == (shift + 1))) ++qn;
1694: if ((val == 1) || (val == -(shift + 1))) ++qp;
1695: }
1696: /* Split old vertex: Edges into original vertex and new cohesive edge */
1697: PetscCall(DMPlexSetSupportSize(sdm, newp, qn + 1));
1698: /* Split new vertex: Edges into split vertex and new cohesive edge */
1699: PetscCall(DMPlexSetSupportSize(sdm, splitp, qp + 1));
1700: /* Add hybrid edge */
1701: PetscCall(DMPlexSetConeSize(sdm, hybedge, 2));
1702: PetscCall(DMPlexSetSupportSize(sdm, hybedge, qf));
1703: PetscCall(DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR));
1704: } else if (dep == dim - 2) {
1705: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1707: PetscCall(DMPlexGetSupport(dm, oldp, &support));
1708: for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) {
1709: PetscInt val;
1711: PetscCall(DMLabelGetValue(label, support[e], &val));
1712: if (val == dim - 1) ++qf;
1713: if ((val == dim - 1) || (val == (shift + dim - 1))) ++qn;
1714: if ((val == dim - 1) || (val == -(shift + dim - 1))) ++qp;
1715: }
1716: /* Split old edge: Faces into original edge and cohesive face (positive side?) */
1717: PetscCall(DMPlexSetSupportSize(sdm, newp, qn + 1));
1718: /* Split new edge: Faces into split edge and cohesive face (negative side?) */
1719: PetscCall(DMPlexSetSupportSize(sdm, splitp, qp + 1));
1720: /* Add hybrid face */
1721: PetscCall(DMPlexSetConeSize(sdm, hybface, 4));
1722: PetscCall(DMPlexSetSupportSize(sdm, hybface, qf));
1723: PetscCall(DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR));
1724: }
1725: }
1726: }
1727: for (dep = 0; dep <= depth; ++dep) {
1728: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
1729: const PetscInt oldp = unsplitPoints[dep][p];
1730: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1731: const PetscInt *support;
1732: PetscInt coneSize, supportSize, qf, e, s;
1734: PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize));
1735: PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize));
1736: PetscCall(DMPlexGetSupport(dm, oldp, &support));
1737: if (dep == 0) {
1738: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1740: /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */
1741: for (s = 0, qf = 0; s < supportSize; ++s, ++qf) {
1742: PetscCall(PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e));
1743: if (e >= 0) ++qf;
1744: }
1745: PetscCall(DMPlexSetSupportSize(sdm, newp, qf + 2));
1746: /* Add hybrid edge */
1747: PetscCall(DMPlexSetConeSize(sdm, hybedge, 2));
1748: for (e = 0, qf = 0; e < supportSize; ++e) {
1749: PetscInt val;
1751: PetscCall(DMLabelGetValue(label, support[e], &val));
1752: /* Split and unsplit edges produce hybrid faces */
1753: if (val == 1) ++qf;
1754: if (val == (shift2 + 1)) ++qf;
1755: }
1756: PetscCall(DMPlexSetSupportSize(sdm, hybedge, qf));
1757: PetscCall(DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR));
1758: } else if (dep == dim - 2) {
1759: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
1760: PetscInt val;
1762: for (e = 0, qf = 0; e < supportSize; ++e) {
1763: PetscCall(DMLabelGetValue(label, support[e], &val));
1764: if (val == dim - 1) qf += 2;
1765: else ++qf;
1766: }
1767: /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */
1768: PetscCall(DMPlexSetSupportSize(sdm, newp, qf + 2));
1769: /* Add hybrid face */
1770: for (e = 0, qf = 0; e < supportSize; ++e) {
1771: PetscCall(DMLabelGetValue(label, support[e], &val));
1772: if (val == dim - 1) ++qf;
1773: }
1774: PetscCall(DMPlexSetConeSize(sdm, hybface, 4));
1775: PetscCall(DMPlexSetSupportSize(sdm, hybface, qf));
1776: PetscCall(DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR));
1777: }
1778: }
1779: }
1780: /* Step 4: Setup split DM */
1781: PetscCall(DMSetUp(sdm));
1782: PetscCall(DMPlexShiftPoints_Internal(dm, depthShift, sdm));
1783: PetscCall(DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew));
1784: PetscCall(PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneNew, PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneONew, PetscMax(maxSupportSize, maxSupportSizeNew), &supportNew));
1785: /* Step 6: Set cones and supports for new points */
1786: for (dep = 0; dep <= depth; ++dep) {
1787: for (p = 0; p < numSplitPoints[dep]; ++p) {
1788: const PetscInt oldp = splitPoints[dep][p];
1789: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
1790: const PetscInt splitp = p + pMaxNew[dep];
1791: const PetscInt *cone, *support, *ornt;
1792: DMPolytopeType ct;
1793: PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s;
1795: PetscCall(DMPlexGetCellType(dm, oldp, &ct));
1796: PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize));
1797: PetscCall(DMPlexGetCone(dm, oldp, &cone));
1798: PetscCall(DMPlexGetConeOrientation(dm, oldp, &ornt));
1799: PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize));
1800: PetscCall(DMPlexGetSupport(dm, oldp, &support));
1801: if (dep == depth - 1) {
1802: PetscBool hasUnsplit = PETSC_FALSE;
1803: const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1804: const PetscInt *supportF;
1806: coneONew[0] = coneONew[1] = -1000;
1807: /* Split face: copy in old face to new face to start */
1808: PetscCall(DMPlexGetSupport(sdm, newp, &supportF));
1809: PetscCall(DMPlexSetSupport(sdm, splitp, supportF));
1810: /* Split old face: old vertices/edges in cone so no change */
1811: /* Split new face: new vertices/edges in cone */
1812: for (q = 0; q < coneSize; ++q) {
1813: PetscCall(PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v));
1814: if (v < 0) {
1815: PetscCall(PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v));
1816: 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);
1817: coneNew[2 + q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1818: hasUnsplit = PETSC_TRUE;
1819: } else {
1820: coneNew[2 + q] = v + pMaxNew[dep - 1];
1821: if (dep > 1) {
1822: const PetscInt *econe;
1823: PetscInt econeSize, r, vs, vu;
1825: PetscCall(DMPlexGetConeSize(dm, cone[q], &econeSize));
1826: PetscCall(DMPlexGetCone(dm, cone[q], &econe));
1827: for (r = 0; r < econeSize; ++r) {
1828: PetscCall(PetscFindInt(econe[r], numSplitPoints[dep - 2], splitPoints[dep - 2], &vs));
1829: PetscCall(PetscFindInt(econe[r], numUnsplitPoints[dep - 2], unsplitPoints[dep - 2], &vu));
1830: if (vs >= 0) continue;
1831: 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);
1832: hasUnsplit = PETSC_TRUE;
1833: }
1834: }
1835: }
1836: }
1837: PetscCall(DMPlexSetCone(sdm, splitp, &coneNew[2]));
1838: PetscCall(DMPlexSetConeOrientation(sdm, splitp, ornt));
1839: /* Face support */
1840: PetscInt vals[2];
1842: PetscCall(DMLabelGetValue(label, support[0], &vals[0]));
1843: if (supportSize > 1) PetscCall(DMLabelGetValue(label, support[1], &vals[1]));
1844: else vals[1] = -vals[0];
1845: PetscCheck(vals[0] * vals[1] < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid support labels %" PetscInt_FMT " %" PetscInt_FMT, vals[0], vals[1]);
1847: for (s = 0; s < 2; ++s) {
1848: if (s >= supportSize) {
1849: if (vals[s] < 0) {
1850: /* Ghost old face: Replace negative side cell with cohesive cell */
1851: PetscCall(DMPlexInsertSupport(sdm, newp, s, hybcell));
1852: } else {
1853: /* Ghost new face: Replace positive side cell with cohesive cell */
1854: PetscCall(DMPlexInsertSupport(sdm, splitp, s, hybcell));
1855: }
1856: } else {
1857: if (vals[s] < 0) {
1858: /* Split old face: Replace negative side cell with cohesive cell */
1859: PetscCall(DMPlexInsertSupport(sdm, newp, s, hybcell));
1860: } else {
1861: /* Split new face: Replace positive side cell with cohesive cell */
1862: PetscCall(DMPlexInsertSupport(sdm, splitp, s, hybcell));
1863: }
1864: }
1865: }
1866: /* Get orientation for cohesive face using the positive side cell */
1867: {
1868: const PetscInt *ncone, *nconeO;
1869: PetscInt nconeSize, nc, ocell;
1870: PetscBool flip = PETSC_FALSE;
1872: if (supportSize > 1) {
1873: ocell = vals[0] < 0 ? support[1] : support[0];
1874: } else {
1875: ocell = support[0];
1876: flip = vals[0] < 0 ? PETSC_TRUE : PETSC_FALSE;
1877: }
1878: PetscCall(DMPlexGetConeSize(dm, ocell, &nconeSize));
1879: PetscCall(DMPlexGetCone(dm, ocell, &ncone));
1880: PetscCall(DMPlexGetConeOrientation(dm, ocell, &nconeO));
1881: for (nc = 0; nc < nconeSize; ++nc) {
1882: if (ncone[nc] == oldp) {
1883: coneONew[0] = flip ? -(nconeO[nc] + 1) : nconeO[nc];
1884: break;
1885: }
1886: }
1887: PetscCheck(nc < nconeSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %" PetscInt_FMT " in neighboring cell %" PetscInt_FMT, oldp, ocell);
1888: }
1889: /* Cohesive cell: Old and new split face, then new cohesive faces */
1890: {
1891: const PetscInt No = DMPolytopeTypeGetNumArrangements(ct) / 2;
1892: PetscCheck(coneONew[0] >= -No && coneONew[0] < No, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid %s orientation %" PetscInt_FMT, DMPolytopeTypes[ct], coneONew[0]);
1893: }
1894: const PetscInt *arr = DMPolytopeTypeGetArrangement(ct, coneONew[0]);
1896: coneNew[0] = newp; /* Extracted negative side orientation above */
1897: coneNew[1] = splitp;
1898: coneONew[1] = coneONew[0];
1899: for (q = 0; q < coneSize; ++q) {
1900: /* Hybrid faces must follow order from oriented end face */
1901: const PetscInt qa = arr[q * 2 + 0];
1902: const PetscInt qo = arr[q * 2 + 1];
1903: DMPolytopeType ft = dep == 2 ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
1905: PetscCall(PetscFindInt(cone[qa], numSplitPoints[dep - 1], splitPoints[dep - 1], &v));
1906: if (v < 0) {
1907: PetscCall(PetscFindInt(cone[qa], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v));
1908: coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
1909: } else {
1910: coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep];
1911: }
1912: coneONew[2 + q] = DMPolytopeTypeComposeOrientation(ft, qo, ornt[qa]);
1913: }
1914: PetscCall(DMPlexSetCone(sdm, hybcell, coneNew));
1915: PetscCall(DMPlexSetConeOrientation(sdm, hybcell, coneONew));
1916: /* Label the hybrid cells on the boundary of the split */
1917: if (hasUnsplit) PetscCall(DMLabelSetValue(label, -hybcell, dim));
1918: } else if (dep == 0) {
1919: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1921: /* Split old vertex: Edges in old split faces and new cohesive edge */
1922: for (e = 0, qn = 0; e < supportSize; ++e) {
1923: PetscInt val;
1925: PetscCall(DMLabelGetValue(label, support[e], &val));
1926: if ((val == 1) || (val == (shift + 1))) supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1927: }
1928: supportNew[qn] = hybedge;
1929: PetscCall(DMPlexSetSupport(sdm, newp, supportNew));
1930: /* Split new vertex: Edges in new split faces and new cohesive edge */
1931: for (e = 0, qp = 0; e < supportSize; ++e) {
1932: PetscInt val, edge;
1934: PetscCall(DMLabelGetValue(label, support[e], &val));
1935: if (val == 1) {
1936: PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge));
1937: PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]);
1938: supportNew[qp++] = edge + pMaxNew[dep + 1];
1939: } else if (val == -(shift + 1)) {
1940: supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1941: }
1942: }
1943: supportNew[qp] = hybedge;
1944: PetscCall(DMPlexSetSupport(sdm, splitp, supportNew));
1945: /* Hybrid edge: Old and new split vertex */
1946: coneNew[0] = newp;
1947: coneNew[1] = splitp;
1948: PetscCall(DMPlexSetCone(sdm, hybedge, coneNew));
1949: for (e = 0, qf = 0; e < supportSize; ++e) {
1950: PetscInt val, edge;
1952: PetscCall(DMLabelGetValue(label, support[e], &val));
1953: if (val == 1) {
1954: PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge));
1955: PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]);
1956: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
1957: }
1958: }
1959: PetscCall(DMPlexSetSupport(sdm, hybedge, supportNew));
1960: } else if (dep == dim - 2) {
1961: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1963: /* Split old edge: old vertices in cone so no change */
1964: /* Split new edge: new vertices in cone */
1965: for (q = 0; q < coneSize; ++q) {
1966: PetscCall(PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v));
1967: if (v < 0) {
1968: PetscCall(PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v));
1969: 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);
1970: coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/;
1971: } else {
1972: coneNew[q] = v + pMaxNew[dep - 1];
1973: }
1974: }
1975: PetscCall(DMPlexSetCone(sdm, splitp, coneNew));
1976: /* Split old edge: Faces in positive side cells and old split faces */
1977: for (e = 0, q = 0; e < supportSize; ++e) {
1978: PetscInt val;
1980: PetscCall(DMLabelGetValue(label, support[e], &val));
1981: if (val == dim - 1) {
1982: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1983: } else if (val == (shift + dim - 1)) {
1984: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
1985: }
1986: }
1987: supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
1988: PetscCall(DMPlexSetSupport(sdm, newp, supportNew));
1989: /* Split new edge: Faces in negative side cells and new split faces */
1990: for (e = 0, q = 0; e < supportSize; ++e) {
1991: PetscInt val, face;
1993: PetscCall(DMLabelGetValue(label, support[e], &val));
1994: if (val == dim - 1) {
1995: PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face));
1996: PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[e]);
1997: supportNew[q++] = face + pMaxNew[dep + 1];
1998: } else if (val == -(shift + dim - 1)) {
1999: supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/;
2000: }
2001: }
2002: supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1];
2003: PetscCall(DMPlexSetSupport(sdm, splitp, supportNew));
2004: /* Hybrid face */
2005: coneNew[0] = newp;
2006: coneNew[1] = splitp;
2007: for (v = 0; v < coneSize; ++v) {
2008: PetscInt vertex;
2009: PetscCall(PetscFindInt(cone[v], numSplitPoints[dep - 1], splitPoints[dep - 1], &vertex));
2010: if (vertex < 0) {
2011: PetscCall(PetscFindInt(cone[v], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &vertex));
2012: 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);
2013: coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
2014: } else {
2015: coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep];
2016: }
2017: }
2018: PetscCall(DMPlexSetCone(sdm, hybface, coneNew));
2019: for (e = 0, qf = 0; e < supportSize; ++e) {
2020: PetscInt val, face;
2022: PetscCall(DMLabelGetValue(label, support[e], &val));
2023: if (val == dim - 1) {
2024: PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face));
2025: PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[e]);
2026: supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
2027: }
2028: }
2029: PetscCall(DMPlexSetSupport(sdm, hybface, supportNew));
2030: }
2031: }
2032: }
2033: for (dep = 0; dep <= depth; ++dep) {
2034: for (p = 0; p < numUnsplitPoints[dep]; ++p) {
2035: const PetscInt oldp = unsplitPoints[dep][p];
2036: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/;
2037: const PetscInt *cone, *support;
2038: PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s;
2040: PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize));
2041: PetscCall(DMPlexGetCone(dm, oldp, &cone));
2042: PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize));
2043: PetscCall(DMPlexGetSupport(dm, oldp, &support));
2044: if (dep == 0) {
2045: const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
2047: /* Unsplit vertex */
2048: PetscCall(DMPlexGetSupportSize(sdm, newp, &supportSizeNew));
2049: for (s = 0, q = 0; s < supportSize; ++s) {
2050: supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/;
2051: PetscCall(PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e));
2052: if (e >= 0) supportNew[q++] = e + pMaxNew[dep + 1];
2053: }
2054: supportNew[q++] = hybedge;
2055: supportNew[q++] = hybedge;
2056: PetscCheck(q == supportSizeNew, comm, PETSC_ERR_ARG_WRONG, "Support size %" PetscInt_FMT " != %" PetscInt_FMT " for vertex %" PetscInt_FMT, q, supportSizeNew, newp);
2057: PetscCall(DMPlexSetSupport(sdm, newp, supportNew));
2058: /* Hybrid edge */
2059: coneNew[0] = newp;
2060: coneNew[1] = newp;
2061: PetscCall(DMPlexSetCone(sdm, hybedge, coneNew));
2062: for (e = 0, qf = 0; e < supportSize; ++e) {
2063: PetscInt val, edge;
2065: PetscCall(DMLabelGetValue(label, support[e], &val));
2066: if (val == 1) {
2067: PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge));
2068: PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]);
2069: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
2070: } else if (val == (shift2 + 1)) {
2071: PetscCall(PetscFindInt(support[e], numUnsplitPoints[dep + 1], unsplitPoints[dep + 1], &edge));
2072: PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a unsplit edge", support[e]);
2073: supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2] + numSplitPoints[dep + 1];
2074: }
2075: }
2076: PetscCall(DMPlexSetSupport(sdm, hybedge, supportNew));
2077: } else if (dep == dim - 2) {
2078: const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep];
2080: /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */
2081: for (f = 0, qf = 0; f < supportSize; ++f) {
2082: PetscInt val, face;
2084: PetscCall(DMLabelGetValue(label, support[f], &val));
2085: if (val == dim - 1) {
2086: PetscCall(PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face));
2087: PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[f]);
2088: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
2089: supportNew[qf++] = face + pMaxNew[dep + 1];
2090: } else {
2091: supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/;
2092: }
2093: }
2094: supportNew[qf++] = hybface;
2095: supportNew[qf++] = hybface;
2096: PetscCall(DMPlexGetSupportSize(sdm, newp, &supportSizeNew));
2097: PetscCheck(qf == supportSizeNew, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %" PetscInt_FMT " is %" PetscInt_FMT " != %" PetscInt_FMT, newp, qf, supportSizeNew);
2098: PetscCall(DMPlexSetSupport(sdm, newp, supportNew));
2099: /* Add hybrid face */
2100: coneNew[0] = newp;
2101: coneNew[1] = newp;
2102: PetscCall(PetscFindInt(cone[0], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v));
2103: PetscCheck(v >= 0, comm, PETSC_ERR_ARG_WRONG, "Vertex %" PetscInt_FMT " is not an unsplit vertex", cone[0]);
2104: coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
2105: PetscCall(PetscFindInt(cone[1], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v));
2106: PetscCheck(v >= 0, comm, PETSC_ERR_ARG_WRONG, "Vertex %" PetscInt_FMT " is not an unsplit vertex", cone[1]);
2107: coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1];
2108: PetscCall(DMPlexSetCone(sdm, hybface, coneNew));
2109: for (f = 0, qf = 0; f < supportSize; ++f) {
2110: PetscInt val, face;
2112: PetscCall(DMLabelGetValue(label, support[f], &val));
2113: if (val == dim - 1) {
2114: PetscCall(PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face));
2115: supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2];
2116: }
2117: }
2118: PetscCall(DMPlexGetSupportSize(sdm, hybface, &supportSizeNew));
2119: PetscCheck(qf == supportSizeNew, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %" PetscInt_FMT " is %" PetscInt_FMT " != %" PetscInt_FMT, hybface, qf, supportSizeNew);
2120: PetscCall(DMPlexSetSupport(sdm, hybface, supportNew));
2121: }
2122: }
2123: }
2124: /* Step 6b: Replace split points in negative side cones */
2125: for (sp = 0; sp < numSP; ++sp) {
2126: PetscInt dep = values[sp];
2127: IS pIS;
2128: PetscInt numPoints;
2129: const PetscInt *points;
2131: if (dep >= 0) continue;
2132: PetscCall(DMLabelGetStratumIS(label, dep, &pIS));
2133: if (!pIS) continue;
2134: dep = -dep - shift;
2135: PetscCall(ISGetLocalSize(pIS, &numPoints));
2136: PetscCall(ISGetIndices(pIS, &points));
2137: for (p = 0; p < numPoints; ++p) {
2138: const PetscInt oldp = points[p];
2139: const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/;
2140: const PetscInt *cone;
2141: PetscInt coneSize, c;
2142: /* PetscBool replaced = PETSC_FALSE; */
2144: /* Negative edge: replace split vertex */
2145: /* Negative cell: replace split face */
2146: PetscCall(DMPlexGetConeSize(sdm, newp, &coneSize));
2147: PetscCall(DMPlexGetCone(sdm, newp, &cone));
2148: for (c = 0; c < coneSize; ++c) {
2149: const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c], depth, depthShift);
2150: PetscInt csplitp, cp, val;
2152: PetscCall(DMLabelGetValue(label, coldp, &val));
2153: if (val == dep - 1) {
2154: PetscCall(PetscFindInt(coldp, numSplitPoints[dep - 1], splitPoints[dep - 1], &cp));
2155: PetscCheck(cp >= 0, comm, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " is not a split point of dimension %" PetscInt_FMT, oldp, dep - 1);
2156: csplitp = pMaxNew[dep - 1] + cp;
2157: PetscCall(DMPlexInsertCone(sdm, newp, c, csplitp));
2158: /* replaced = PETSC_TRUE; */
2159: }
2160: }
2161: /* Cells with only a vertex or edge on the submesh have no replacement */
2162: /* PetscCheck(replaced,comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */
2163: }
2164: PetscCall(ISRestoreIndices(pIS, &points));
2165: PetscCall(ISDestroy(&pIS));
2166: }
2167: PetscCall(DMPlexReorderCohesiveSupports(sdm));
2168: /* Step 7: Coordinates */
2169: PetscCall(DMPlexShiftCoordinates_Internal(dm, depthShift, sdm));
2170: PetscCall(DMGetCoordinateSection(sdm, &coordSection));
2171: PetscCall(DMGetCoordinatesLocal(sdm, &coordinates));
2172: PetscCall(VecGetArray(coordinates, &coords));
2173: for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) {
2174: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/;
2175: const PetscInt splitp = pMaxNew[0] + v;
2176: PetscInt dof, off, soff, d;
2178: PetscCall(PetscSectionGetDof(coordSection, newp, &dof));
2179: PetscCall(PetscSectionGetOffset(coordSection, newp, &off));
2180: PetscCall(PetscSectionGetOffset(coordSection, splitp, &soff));
2181: for (d = 0; d < dof; ++d) coords[soff + d] = coords[off + d];
2182: }
2183: PetscCall(VecRestoreArray(coordinates, &coords));
2184: /* Step 8: SF, if I can figure this out we can split the mesh in parallel */
2185: PetscCall(DMPlexShiftSF_Internal(dm, depthShift, sdm));
2186: /* TODO We need to associate the ghost points with the correct replica */
2187: /* Step 9: Labels */
2188: PetscCall(DMPlexShiftLabels_Internal(dm, depthShift, sdm));
2189: PetscCall(DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm));
2190: PetscCall(DMGetNumLabels(sdm, &numLabels));
2191: for (dep = 0; dep <= depth; ++dep) {
2192: for (p = 0; p < numSplitPoints[dep]; ++p) {
2193: const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/;
2194: const PetscInt splitp = pMaxNew[dep] + p;
2195: PetscInt l;
2197: if (splitLabel) {
2198: const PetscInt val = 100 + dep;
2200: PetscCall(DMLabelSetValue(splitLabel, newp, val));
2201: PetscCall(DMLabelSetValue(splitLabel, splitp, -val));
2202: }
2203: for (l = 0; l < numLabels; ++l) {
2204: DMLabel mlabel;
2205: const char *lname;
2206: PetscInt val;
2207: PetscBool isDepth;
2209: PetscCall(DMGetLabelName(sdm, l, &lname));
2210: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
2211: if (isDepth) continue;
2212: PetscCall(DMGetLabel(sdm, lname, &mlabel));
2213: PetscCall(DMLabelGetValue(mlabel, newp, &val));
2214: if (val >= 0) PetscCall(DMLabelSetValue(mlabel, splitp, val));
2215: }
2216: }
2217: }
2218: for (sp = 0; sp < numSP; ++sp) {
2219: const PetscInt dep = values[sp];
2221: if ((dep < 0) || (dep > depth)) continue;
2222: if (splitIS[dep]) PetscCall(ISRestoreIndices(splitIS[dep], &splitPoints[dep]));
2223: PetscCall(ISDestroy(&splitIS[dep]));
2224: if (unsplitIS[dep]) PetscCall(ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]));
2225: PetscCall(ISDestroy(&unsplitIS[dep]));
2226: }
2227: if (ghostIS) PetscCall(ISRestoreIndices(ghostIS, &ghostPoints));
2228: PetscCall(ISDestroy(&ghostIS));
2229: if (label) {
2230: PetscCall(ISRestoreIndices(valueIS, &values));
2231: PetscCall(ISDestroy(&valueIS));
2232: }
2233: for (d = 0; d <= depth; ++d) {
2234: PetscCall(DMPlexGetDepthStratum(sdm, d, NULL, &pEnd));
2235: pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d];
2236: }
2237: PetscCall(PetscFree3(coneNew, coneONew, supportNew));
2238: PetscCall(PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld));
2239: PetscCall(PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints));
2240: PetscFunctionReturn(PETSC_SUCCESS);
2241: }
2243: /*@C
2244: DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface
2246: Collective
2248: Input Parameters:
2249: + dm - The original `DM`
2250: - label - The `DMLabel` specifying the boundary faces (this could be auto-generated)
2252: Output Parameters:
2253: + splitLabel - The `DMLabel` containing the split points, or `NULL` if no output is desired
2254: - dmSplit - The new `DM`
2256: Level: developer
2258: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexLabelCohesiveComplete()`
2259: @*/
2260: PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit)
2261: {
2262: DM sdm;
2263: PetscInt dim;
2265: PetscFunctionBegin;
2267: PetscAssertPointer(dmSplit, 4);
2268: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm));
2269: PetscCall(DMSetType(sdm, DMPLEX));
2270: PetscCall(DMGetDimension(dm, &dim));
2271: PetscCall(DMSetDimension(sdm, dim));
2272: switch (dim) {
2273: case 2:
2274: case 3:
2275: PetscCall(DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm));
2276: break;
2277: default:
2278: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %" PetscInt_FMT, dim);
2279: }
2280: *dmSplit = sdm;
2281: PetscFunctionReturn(PETSC_SUCCESS);
2282: }
2284: static PetscErrorCode ViewCohesiveLabel_Static(MPI_Comm comm, DMLabel label)
2285: {
2286: PetscViewer viewer;
2288: PetscFunctionBegin;
2289: PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_(comm), PETSC_COMM_SELF, &viewer));
2290: PetscCall(DMLabelView(label, viewer));
2291: PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_(comm), PETSC_COMM_SELF, &viewer));
2292: PetscFunctionReturn(PETSC_SUCCESS);
2293: }
2295: /* Returns the side of the surface for a given cell with a face on the surface */
2296: static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos)
2297: {
2298: const PetscInt *cone, *ornt;
2299: PetscInt dim, coneSize, c;
2301: PetscFunctionBegin;
2302: *pos = PETSC_TRUE;
2303: PetscCall(DMGetDimension(dm, &dim));
2304: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
2305: PetscCall(DMPlexGetCone(dm, cell, &cone));
2306: PetscCall(DMPlexGetConeOrientation(dm, cell, &ornt));
2307: for (c = 0; c < coneSize; ++c) {
2308: if (cone[c] == face) {
2309: PetscInt o = ornt[c];
2311: if (subdm) {
2312: const PetscInt *subcone, *subornt;
2313: PetscInt subpoint, subface, subconeSize, sc;
2315: PetscCall(PetscFindInt(cell, numSubpoints, subpoints, &subpoint));
2316: PetscCall(PetscFindInt(face, numSubpoints, subpoints, &subface));
2317: PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize));
2318: PetscCall(DMPlexGetCone(subdm, subpoint, &subcone));
2319: PetscCall(DMPlexGetConeOrientation(subdm, subpoint, &subornt));
2320: for (sc = 0; sc < subconeSize; ++sc) {
2321: if (subcone[sc] == subface) {
2322: o = subornt[0];
2323: break;
2324: }
2325: }
2326: 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);
2327: }
2328: if (o >= 0) *pos = PETSC_TRUE;
2329: else *pos = PETSC_FALSE;
2330: break;
2331: }
2332: }
2333: 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);
2334: PetscFunctionReturn(PETSC_SUCCESS);
2335: }
2337: static PetscErrorCode CheckFaultEdge_Private(DM dm, DMLabel label, PetscBool split)
2338: {
2339: const PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
2340: const PetscInt shift = 100;
2341: const PetscInt shift2 = 200;
2342: IS facePosIS, faceNegIS, dimIS;
2343: const PetscInt *points;
2344: PetscInt *closure = NULL, *inclosure = NULL;
2345: PetscInt dim, numPoints;
2347: PetscFunctionBegin;
2348: PetscCall(DMGetDimension(dm, &dim));
2349: // If any faces touching the fault divide cells on either side,
2350: // either split them, or unsplit the connection
2351: PetscCall(DMLabelGetStratumIS(label, shift + dim - 1, &facePosIS));
2352: PetscCall(DMLabelGetStratumIS(label, -(shift + dim - 1), &faceNegIS));
2353: if (!facePosIS || !faceNegIS) {
2354: PetscCall(ISDestroy(&facePosIS));
2355: PetscCall(ISDestroy(&faceNegIS));
2356: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
2357: PetscFunctionReturn(PETSC_SUCCESS);
2358: }
2359: PetscCall(ISExpand(facePosIS, faceNegIS, &dimIS));
2360: PetscCall(ISDestroy(&facePosIS));
2361: PetscCall(ISDestroy(&faceNegIS));
2362: PetscCall(ISGetLocalSize(dimIS, &numPoints));
2363: PetscCall(ISGetIndices(dimIS, &points));
2364: for (PetscInt p = 0; p < numPoints; ++p) {
2365: const PetscInt point = points[p];
2366: const PetscInt *support;
2367: PetscInt supportSize, valA, valB;
2369: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
2370: if (supportSize != 2) continue;
2371: PetscCall(DMPlexGetSupport(dm, point, &support));
2372: PetscCall(DMLabelGetValue(label, support[0], &valA));
2373: PetscCall(DMLabelGetValue(label, support[1], &valB));
2374: if ((valA == -1) || (valB == -1)) continue;
2375: if (valA * valB > 0) continue;
2376: // Check that this face is not incident on only unsplit faces,
2377: // meaning has at least one split face
2378: {
2379: PetscBool split = PETSC_FALSE;
2380: PetscInt Ncl, val;
2382: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &Ncl, &closure));
2383: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2384: PetscCall(DMLabelGetValue(label, closure[cl], &val));
2385: if (val >= 0 && val <= dim) {
2386: split = PETSC_TRUE;
2387: break;
2388: }
2389: }
2390: if (!split) continue;
2391: }
2392: if (debug)
2393: PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d]Point %" PetscInt_FMT " is impinging (%" PetscInt_FMT ":%" PetscInt_FMT ", %" PetscInt_FMT ":%" PetscInt_FMT ")\n", PetscGlobalRank, point, support[0], valA, support[1], valB));
2394: if (split) {
2395: // Split the face
2396: PetscCall(DMLabelGetValue(label, point, &valA));
2397: PetscCall(DMLabelClearValue(label, point, valA));
2398: PetscCall(DMLabelSetValue(label, point, dim - 1));
2399: /* Label its closure:
2400: unmarked: label as unsplit
2401: incident: relabel as split
2402: split: do nothing */
2403: {
2404: PetscInt closureSize, cl, dep;
2406: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure));
2407: for (cl = 0; cl < closureSize * 2; cl += 2) {
2408: PetscCall(DMLabelGetValue(label, closure[cl], &valA));
2409: if (valA == -1) { /* Mark as unsplit */
2410: PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep));
2411: PetscCall(DMLabelSetValue(label, closure[cl], shift2 + dep));
2412: } else if ((valA >= shift && valA < shift2) || (valA <= -shift && valA > -shift2)) {
2413: PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep));
2414: PetscCall(DMLabelClearValue(label, closure[cl], valA));
2415: PetscCall(DMLabelSetValue(label, closure[cl], dep));
2416: }
2417: }
2418: }
2419: } else {
2420: // Unsplit the incident faces and their closures
2421: PetscInt Ncl, dep, val;
2423: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &Ncl, &closure));
2424: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2425: PetscCall(DMLabelGetValue(label, closure[cl], &val));
2426: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Point %" PetscInt_FMT ":%" PetscInt_FMT "\n", PetscGlobalRank, closure[cl], val));
2427: if (val >= 0 && val <= dim) {
2428: PetscInt Nincl, inval, indep;
2430: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Point %" PetscInt_FMT " is being unsplit\n", PetscGlobalRank, closure[cl]));
2431: PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep));
2432: PetscCall(DMLabelClearValue(label, closure[cl], val));
2433: PetscCall(DMLabelSetValue(label, closure[cl], shift2 + dep));
2435: PetscCall(DMPlexGetTransitiveClosure(dm, closure[cl], PETSC_TRUE, &Nincl, &inclosure));
2436: for (PetscInt incl = 0; incl < Nincl * 2; incl += 2) {
2437: PetscCall(DMLabelGetValue(label, inclosure[cl], &inval));
2438: if (inval >= 0 && inval <= dim) {
2439: if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Point %" PetscInt_FMT " is being unsplit\n", PetscGlobalRank, inclosure[incl]));
2440: PetscCall(DMPlexGetPointDepth(dm, inclosure[incl], &indep));
2441: PetscCall(DMLabelClearValue(label, inclosure[incl], inval));
2442: PetscCall(DMLabelSetValue(label, inclosure[incl], shift2 + indep));
2443: }
2444: }
2445: }
2446: }
2447: }
2448: }
2449: PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &inclosure));
2450: PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure));
2451: PetscCall(ISRestoreIndices(dimIS, &points));
2452: PetscCall(ISDestroy(&dimIS));
2453: if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL));
2454: PetscFunctionReturn(PETSC_SUCCESS);
2455: }
2457: static void MPIAPI DMPlexLabelCohesiveValueReduce_Private(void *a, void *b, int *len, MPI_Datatype *datatype)
2458: {
2459: const int N = *len;
2461: if (*datatype == MPIU_INT) {
2462: PetscInt *A = (PetscInt *)a;
2463: PetscInt *B = (PetscInt *)b;
2465: for (int i = 0; i < N; i++) {
2466: // Default values do not propagate
2467: if (A[i] == -1) continue;
2468: // Override default values
2469: if (B[i] == -1) {
2470: B[i] = A[i];
2471: continue;
2472: }
2473: // Unsplit points are preferred
2474: if (A[i] == B[i] + 200) {
2475: B[i] = A[i];
2476: continue;
2477: }
2478: if (B[i] == A[i] + 200) continue;
2479: // Ignore ghost designation
2480: if (B[i] >= 300) {
2481: const PetscInt a = A[i] - 300;
2482: const PetscInt b = B[i] < 0 ? (B[i] < -100 ? -B[i] - 100 : -B[i]) : (B[i] >= 200 ? B[i] - 200 : (B[i] > 100 ? B[i] - 100 : B[i]));
2484: if (a == b) {
2485: B[i] = A[i];
2486: continue;
2487: }
2488: }
2489: if (A[i] >= 300) {
2490: const PetscInt a = A[i] - 300;
2491: const PetscInt b = B[i] < 0 ? (B[i] < -100 ? -B[i] - 100 : -B[i]) : (B[i] >= 200 ? B[i] - 200 : (B[i] > 100 ? B[i] - 100 : B[i]));
2493: if (a == b) continue;
2494: }
2495: // Process can disagree about the surface side on the boundary
2496: if (A[i] == -B[i]) continue;
2497: // Error
2498: if (A[i] != B[i]) {
2499: printf("ERROR A[%d] = %d, B[%d] = %d\n", i, (int)A[i], i, (int)B[i]);
2500: MPI_Abort(MPI_COMM_WORLD, 1);
2501: }
2502: }
2503: }
2504: }
2506: /*@
2507: DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces
2508: to complete the surface
2510: Input Parameters:
2511: + dm - The `DM`
2512: . label - A `DMLabel` marking the surface
2513: . blabel - A `DMLabel` marking the vertices on the boundary which will not be duplicated, or `NULL` to find them automatically
2514: . bvalue - Value of `DMLabel` marking the vertices on the boundary
2515: . flip - Flag to flip the submesh normal and replace points on the other side
2516: . split - Split faces impinging on the surface, rather than clamping the surface boundary
2517: - subdm - The `DM` associated with the label, or `NULL`
2519: Output Parameter:
2520: . label - A `DMLabel` marking all surface points
2522: Level: developer
2524: Note:
2525: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation.
2527: Points are marked with their dimension, combined with a shift based on the type of interation with the surface. For points on the surface itself, the shift is zero. Mesh points impinging on the surface have a shoft of 100, and then are negated for points on the negative side of the fault. Points on the surface boundary, called unsplit, are shifted by 200. Cells on the surface that are not owned by this process are shifted by 300.
2529: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexConstructCohesiveCells()`, `DMPlexLabelComplete()`
2530: @*/
2531: PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscInt bvalue, PetscBool flip, PetscBool split, DM subdm)
2532: {
2533: const PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
2534: const PetscInt shift = 100;
2535: const PetscInt shift2 = 200;
2536: const PetscInt shift3 = split ? 300 : 0;
2537: DMLabel depthLabel;
2538: IS dimIS, subpointIS = NULL;
2539: const PetscInt *points, *subpoints;
2540: const PetscInt rev = flip ? -1 : 1;
2541: PetscInt dim, depth, numPoints, numSubpoints, p, val;
2542: MPI_Comm comm;
2544: PetscFunctionBegin;
2545: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2546: if (debug) PetscCall(ViewCohesiveLabel_Static(comm, label));
2547: PetscCall(DMPlexGetDepth(dm, &depth));
2548: PetscCall(DMGetDimension(dm, &dim));
2549: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
2550: if (subdm) {
2551: PetscCall(DMPlexGetSubpointIS(subdm, &subpointIS));
2552: if (subpointIS) {
2553: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
2554: PetscCall(ISGetIndices(subpointIS, &subpoints));
2555: }
2556: }
2557: /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */
2558: PetscCall(DMLabelGetStratumIS(label, dim - 1, &dimIS));
2559: if (!dimIS) goto divide;
2560: PetscCall(ISGetLocalSize(dimIS, &numPoints));
2561: PetscCall(ISGetIndices(dimIS, &points));
2562: for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */
2563: const PetscInt *support;
2564: PetscInt supportSize, s;
2566: PetscCall(DMPlexGetSupportSize(dm, points[p], &supportSize));
2567: #if 0
2568: if (supportSize != 2) {
2569: const PetscInt *lp;
2570: PetscInt Nlp, pind;
2572: /* Check that for a cell with a single support face, that face is in the SF */
2573: /* THis check only works for the remote side. We would need root side information */
2574: PetscCall(PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL));
2575: PetscCall(PetscFindInt(points[p], Nlp, lp, &pind));
2576: 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);
2577: }
2578: #endif
2579: PetscCall(DMPlexGetSupport(dm, points[p], &support));
2580: for (s = 0; s < supportSize; ++s) {
2581: const PetscInt *cone;
2582: PetscInt coneSize, c;
2583: PetscBool pos;
2585: PetscCall(GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos));
2586: if (pos) PetscCall(DMLabelSetValue(label, support[s], rev * (shift + dim)));
2587: else PetscCall(DMLabelSetValue(label, support[s], -rev * (shift + dim)));
2588: if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Cell %" PetscInt_FMT " on %s side\n", PetscGlobalRank, support[s], pos ? "positive" : "negative"));
2589: if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE;
2590: /* Put faces touching the fault in the label */
2591: PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize));
2592: PetscCall(DMPlexGetCone(dm, support[s], &cone));
2593: for (c = 0; c < coneSize; ++c) {
2594: const PetscInt point = cone[c];
2596: PetscCall(DMLabelGetValue(label, point, &val));
2597: if (val == -1) {
2598: PetscInt *closure = NULL;
2599: PetscInt closureSize, cl;
2601: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure));
2602: for (cl = 0; cl < closureSize * 2; cl += 2) {
2603: const PetscInt clp = closure[cl];
2604: PetscInt bval = -1;
2606: PetscCall(DMLabelGetValue(label, clp, &val));
2607: if (blabel) PetscCall(DMLabelGetValue(blabel, clp, &bval));
2608: if (val >= 0 && val < dim - 1 && bval < 0) {
2609: PetscCall(DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift + dim - 1 : -(shift + dim - 1)));
2610: break;
2611: }
2612: }
2613: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure));
2614: }
2615: }
2616: }
2617: }
2618: PetscCall(ISRestoreIndices(dimIS, &points));
2619: PetscCall(ISDestroy(&dimIS));
2620: /* Mark boundary points as unsplit */
2621: if (blabel) {
2622: IS bdIS;
2624: PetscCall(DMLabelGetStratumIS(blabel, bvalue, &bdIS));
2625: PetscCall(ISGetLocalSize(bdIS, &numPoints));
2626: PetscCall(ISGetIndices(bdIS, &points));
2627: for (p = 0; p < numPoints; ++p) {
2628: const PetscInt point = points[p];
2629: PetscInt val, bval;
2631: PetscCall(DMLabelGetValue(blabel, point, &bval));
2632: if (bval >= 0) {
2633: PetscCall(DMLabelGetValue(label, point, &val));
2634: if ((val < 0) || (val > dim)) {
2635: /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */
2636: PetscCall(DMLabelClearValue(blabel, point, bval));
2637: }
2638: }
2639: }
2640: for (p = 0; p < numPoints; ++p) {
2641: const PetscInt point = points[p];
2642: PetscInt val, bval;
2644: PetscCall(DMLabelGetValue(blabel, point, &bval));
2645: if (bval >= 0) {
2646: const PetscInt *cone, *support;
2647: PetscInt coneSize, supportSize, s, valA, valB, valE;
2649: /* Mark as unsplit */
2650: PetscCall(DMLabelGetValue(label, point, &val));
2651: 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);
2652: PetscCall(DMLabelClearValue(label, point, val));
2653: PetscCall(DMLabelSetValue(label, point, shift2 + val));
2654: /* Check for cross-edge
2655: A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */
2656: if (val != 0) continue;
2657: PetscCall(DMPlexGetSupport(dm, point, &support));
2658: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
2659: for (s = 0; s < supportSize; ++s) {
2660: PetscCall(DMPlexGetCone(dm, support[s], &cone));
2661: PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize));
2662: PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %" PetscInt_FMT " has %" PetscInt_FMT " vertices != 2", support[s], coneSize);
2663: PetscCall(DMLabelGetValue(blabel, cone[0], &valA));
2664: PetscCall(DMLabelGetValue(blabel, cone[1], &valB));
2665: PetscCall(DMLabelGetValue(blabel, support[s], &valE));
2666: if (valE < 0 && valA >= 0 && valB >= 0 && cone[0] != cone[1]) PetscCall(DMLabelSetValue(blabel, support[s], 2));
2667: }
2668: }
2669: }
2670: PetscCall(ISRestoreIndices(bdIS, &points));
2671: PetscCall(ISDestroy(&bdIS));
2672: }
2673: /* Mark ghost fault cells */
2674: {
2675: PetscSF sf;
2676: const PetscInt *leaves;
2677: PetscInt Nl, l;
2679: PetscCall(DMGetPointSF(dm, &sf));
2680: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
2681: PetscCall(DMLabelGetStratumIS(label, dim - 1, &dimIS));
2682: if (!dimIS) goto divide;
2683: PetscCall(ISGetLocalSize(dimIS, &numPoints));
2684: PetscCall(ISGetIndices(dimIS, &points));
2685: if (Nl > 0) {
2686: for (p = 0; p < numPoints; ++p) {
2687: const PetscInt point = points[p];
2688: PetscInt val;
2690: PetscCall(PetscFindInt(point, Nl, leaves, &l));
2691: if (l >= 0) {
2692: PetscInt *closure = NULL;
2693: PetscInt closureSize, cl;
2695: PetscCall(DMLabelGetValue(label, point, &val));
2696: 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);
2697: PetscCall(DMLabelClearValue(label, point, val));
2698: PetscCall(DMLabelSetValue(label, point, shift3 + val));
2699: PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure));
2700: for (cl = 2; cl < closureSize * 2; cl += 2) {
2701: const PetscInt clp = closure[cl];
2703: PetscCall(DMLabelGetValue(label, clp, &val));
2704: 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);
2705: PetscCall(DMLabelClearValue(label, clp, val));
2706: PetscCall(DMLabelSetValue(label, clp, shift3 + val));
2707: }
2708: PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure));
2709: }
2710: }
2711: }
2712: PetscCall(ISRestoreIndices(dimIS, &points));
2713: PetscCall(ISDestroy(&dimIS));
2714: }
2715: divide:
2716: if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL));
2717: if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &subpoints));
2718: if (debug) PetscCall(ViewCohesiveLabel_Static(comm, label));
2719: PetscCall(DMPlexLabelFaultHalo(dm, label));
2720: PetscCall(CheckFaultEdge_Private(dm, label, split));
2721: {
2722: MPI_Op reduceop;
2723: IS newpointIS;
2724: const PetscInt newValue = 1000;
2725: const PetscInt *newpoints;
2726: PetscInt n, defVal;
2728: PetscCallMPI(MPI_Op_create(DMPlexLabelCohesiveValueReduce_Private, PETSC_TRUE, &reduceop));
2729: PetscCall(DMPlexReconcileLabel(dm, reduceop, newValue, label));
2730: PetscCallMPI(MPI_Op_free(&reduceop));
2731: // New points may have connected points
2732: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2733: PetscCall(DMLabelGetStratumIS(label, newValue, &newpointIS));
2734: if (newpointIS) {
2735: PetscCall(ISGetLocalSize(newpointIS, &n));
2736: PetscCall(ISGetIndices(newpointIS, &newpoints));
2737: for (PetscInt i = 0; i < n; ++i) {
2738: PetscBool found = PETSC_FALSE;
2739: const PetscInt point = newpoints[i];
2740: const PetscInt *cone;
2741: PetscInt cS, val;
2743: PetscCall(DMLabelClearValue(label, point, newValue));
2744: PetscCall(DMLabelGetValue(label, point, &val));
2745: // Find point on fault
2746: PetscCall(DMPlexGetCone(dm, point, &cone));
2747: PetscCall(DMPlexGetConeSize(dm, point, &cS));
2748: for (PetscInt c = 0; c < cS; ++c) {
2749: PetscInt cval;
2751: PetscCall(DMLabelGetValue(label, cone[c], &cval));
2752: if (0 <= cval && cval < 100) {
2753: const PetscInt *supp;
2754: PetscInt sS;
2756: found = PETSC_TRUE;
2757: // Add support of fault point to label
2758: PetscCall(DMPlexGetSupport(dm, cone[c], &supp));
2759: PetscCall(DMPlexGetSupportSize(dm, cone[c], &sS));
2760: for (PetscInt s = 0; s < sS; ++s) {
2761: PetscInt sval, pdepth;
2763: PetscCall(DMLabelGetValue(label, supp[s], &sval));
2764: if (sval == defVal) {
2765: PetscCall(DMPlexGetPointDepth(dm, supp[s], &pdepth));
2766: PetscCall(DMLabelSetValue(label, supp[s], val >= 0 ? shift + pdepth : -shift - pdepth));
2767: }
2768: }
2769: }
2770: }
2771: PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "New label point %" PetscInt_FMT " must be connected to fault", point);
2772: }
2773: PetscCall(ISRestoreIndices(newpointIS, &newpoints));
2774: PetscCall(ISDestroy(&newpointIS));
2775: }
2776: }
2777: if (debug) PetscCall(ViewCohesiveLabel_Static(comm, label));
2778: PetscFunctionReturn(PETSC_SUCCESS);
2779: }
2781: /* Check that no cell have all vertices on the fault */
2782: static PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm)
2783: {
2784: IS subpointIS;
2785: const PetscInt *dmpoints;
2786: PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd;
2788: PetscFunctionBegin;
2789: if (!label) PetscFunctionReturn(PETSC_SUCCESS);
2790: PetscCall(DMLabelGetDefaultValue(label, &defaultValue));
2791: PetscCall(DMPlexGetSubpointIS(subdm, &subpointIS));
2792: if (!subpointIS) PetscFunctionReturn(PETSC_SUCCESS);
2793: PetscCall(DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd));
2794: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
2795: PetscCall(ISGetIndices(subpointIS, &dmpoints));
2796: for (c = cStart; c < cEnd; ++c) {
2797: PetscBool invalidCell = PETSC_TRUE;
2798: PetscInt *closure = NULL;
2799: PetscInt closureSize, cl;
2801: PetscCall(DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure));
2802: for (cl = 0; cl < closureSize * 2; cl += 2) {
2803: PetscInt value = 0;
2805: if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue;
2806: PetscCall(DMLabelGetValue(label, closure[cl], &value));
2807: if (value == defaultValue) {
2808: invalidCell = PETSC_FALSE;
2809: break;
2810: }
2811: }
2812: PetscCall(DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure));
2813: if (invalidCell) {
2814: PetscCall(ISRestoreIndices(subpointIS, &dmpoints));
2815: PetscCall(ISDestroy(&subpointIS));
2816: PetscCall(DMDestroy(&subdm));
2817: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %" PetscInt_FMT " has all of its vertices on the submesh.", dmpoints[c]);
2818: }
2819: }
2820: PetscCall(ISRestoreIndices(subpointIS, &dmpoints));
2821: PetscFunctionReturn(PETSC_SUCCESS);
2822: }
2824: /*@
2825: DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface
2827: Collective
2829: Input Parameters:
2830: + dm - The original `DM`
2831: . label - The label specifying the interface vertices
2832: . bdlabel - The optional label specifying the interface boundary vertices
2833: - bdvalue - Value of optional label specifying the interface boundary vertices
2835: Output Parameters:
2836: + hybridLabel - The label fully marking the interface, or `NULL` if no output is desired
2837: . splitLabel - The label containing the split points, or `NULL` if no output is desired
2838: . dmInterface - The new interface `DM`, or `NULL`
2839: - dmHybrid - The new `DM` with cohesive cells
2841: Level: developer
2843: Note:
2844: The hybridLabel indicates what parts of the original mesh impinged on the division surface. For points
2845: directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be
2846: 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with
2847: 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
2848: surface also hits vertex 3, it would be 9 (-101) in hybridLabel.
2850: The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original
2851: mesh. The label value is $\pm 100+dim$ for each point. For example, if two edges 10 and 14 in the hybrid resulting from
2852: splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel.
2854: The dmInterface is a `DM` built from the original division surface. It has a label which can be retrieved using
2855: `DMPlexGetSubpointMap()` which maps each point back to the point in the surface of the original mesh.
2857: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexConstructCohesiveCells()`, `DMPlexLabelCohesiveComplete()`, `DMPlexGetSubpointMap()`, `DMCreate()`
2858: @*/
2859: PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, PetscInt bdvalue, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid)
2860: {
2861: DM idm;
2862: DMLabel subpointMap, hlabel, slabel = NULL;
2863: PetscInt dim;
2865: PetscFunctionBegin;
2867: if (label) PetscAssertPointer(label, 2);
2868: if (bdlabel) PetscAssertPointer(bdlabel, 3);
2869: if (hybridLabel) PetscAssertPointer(hybridLabel, 5);
2870: if (splitLabel) PetscAssertPointer(splitLabel, 6);
2871: if (dmInterface) PetscAssertPointer(dmInterface, 7);
2872: PetscAssertPointer(dmHybrid, 8);
2873: PetscCall(DMGetDimension(dm, &dim));
2874: PetscCall(DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm));
2875: PetscCall(DMPlexCheckValidSubmesh_Private(dm, label, idm));
2876: PetscCall(DMPlexOrient(idm));
2877: PetscCall(DMPlexGetSubpointMap(idm, &subpointMap));
2878: PetscCall(DMLabelDuplicate(subpointMap, &hlabel));
2879: PetscCall(DMLabelClearStratum(hlabel, dim));
2880: if (splitLabel) {
2881: const char *name;
2882: char sname[PETSC_MAX_PATH_LEN];
2884: PetscCall(PetscObjectGetName((PetscObject)hlabel, &name));
2885: PetscCall(PetscStrncpy(sname, name, sizeof(sname)));
2886: PetscCall(PetscStrlcat(sname, " split", sizeof(sname)));
2887: PetscCall(DMLabelCreate(PETSC_COMM_SELF, sname, &slabel));
2888: }
2889: PetscCall(DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, bdvalue, PETSC_FALSE, PETSC_TRUE, idm));
2890: if (dmInterface) {
2891: *dmInterface = idm;
2892: } else PetscCall(DMDestroy(&idm));
2893: PetscCall(DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid));
2894: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmHybrid));
2895: if (hybridLabel) *hybridLabel = hlabel;
2896: else PetscCall(DMLabelDestroy(&hlabel));
2897: if (splitLabel) *splitLabel = slabel;
2898: {
2899: DM cdm;
2900: DMLabel ctLabel;
2902: /* We need to somehow share the celltype label with the coordinate dm */
2903: PetscCall(DMGetCoordinateDM(*dmHybrid, &cdm));
2904: PetscCall(DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel));
2905: PetscCall(DMSetLabel(cdm, ctLabel));
2906: }
2907: PetscFunctionReturn(PETSC_SUCCESS);
2908: }
2910: /* Here we need the explicit assumption that:
2912: For any marked cell, the marked vertices constitute a single face
2913: */
2914: static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm)
2915: {
2916: IS subvertexIS = NULL;
2917: const PetscInt *subvertices;
2918: PetscInt *pStart, *pEnd, pSize;
2919: PetscInt depth, dim, d, numSubVerticesInitial = 0, v;
2921: PetscFunctionBegin;
2922: *numFaces = 0;
2923: *nFV = 0;
2924: PetscCall(DMPlexGetDepth(dm, &depth));
2925: PetscCall(DMGetDimension(dm, &dim));
2926: pSize = PetscMax(depth, dim) + 1;
2927: PetscCall(PetscMalloc2(pSize, &pStart, pSize, &pEnd));
2928: for (d = 0; d <= depth; ++d) PetscCall(DMPlexGetSimplexOrBoxCells(dm, depth - d, &pStart[d], &pEnd[d]));
2929: /* Loop over initial vertices and mark all faces in the collective star() */
2930: if (vertexLabel) PetscCall(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS));
2931: if (subvertexIS) {
2932: PetscCall(ISGetSize(subvertexIS, &numSubVerticesInitial));
2933: PetscCall(ISGetIndices(subvertexIS, &subvertices));
2934: }
2935: for (v = 0; v < numSubVerticesInitial; ++v) {
2936: const PetscInt vertex = subvertices[v];
2937: PetscInt *star = NULL;
2938: PetscInt starSize, s, numCells = 0, c;
2940: PetscCall(DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star));
2941: for (s = 0; s < starSize * 2; s += 2) {
2942: const PetscInt point = star[s];
2943: if (point >= pStart[depth] && point < pEnd[depth]) star[numCells++] = point;
2944: }
2945: for (c = 0; c < numCells; ++c) {
2946: const PetscInt cell = star[c];
2947: PetscInt *closure = NULL;
2948: PetscInt closureSize, cl;
2949: PetscInt cellLoc, numCorners = 0, faceSize = 0;
2951: PetscCall(DMLabelGetValue(subpointMap, cell, &cellLoc));
2952: if (cellLoc == 2) continue;
2953: PetscCheck(cellLoc < 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " has dimension %" PetscInt_FMT " in the surface label", cell, cellLoc);
2954: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
2955: for (cl = 0; cl < closureSize * 2; cl += 2) {
2956: const PetscInt point = closure[cl];
2957: PetscInt vertexLoc;
2959: if (point >= pStart[0] && point < pEnd[0]) {
2960: ++numCorners;
2961: PetscCall(DMLabelGetValue(vertexLabel, point, &vertexLoc));
2962: if (vertexLoc == value) closure[faceSize++] = point;
2963: }
2964: }
2965: if (!*nFV) PetscCall(DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV));
2966: PetscCheck(faceSize <= *nFV, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %" PetscInt_FMT " of an element on the surface", faceSize);
2967: if (faceSize == *nFV) {
2968: const PetscInt *cells = NULL;
2969: PetscInt numCells, nc;
2971: ++(*numFaces);
2972: for (cl = 0; cl < faceSize; ++cl) PetscCall(DMLabelSetValue(subpointMap, closure[cl], 0));
2973: PetscCall(DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells));
2974: for (nc = 0; nc < numCells; ++nc) PetscCall(DMLabelSetValue(subpointMap, cells[nc], 2));
2975: PetscCall(DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells));
2976: }
2977: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
2978: }
2979: PetscCall(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star));
2980: }
2981: if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subvertices));
2982: PetscCall(ISDestroy(&subvertexIS));
2983: PetscCall(PetscFree2(pStart, pEnd));
2984: PetscFunctionReturn(PETSC_SUCCESS);
2985: }
2987: PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, PetscBool addCells, DMLabel subpointMap, DM subdm)
2988: {
2989: IS subvertexIS = NULL;
2990: const PetscInt *subvertices;
2991: PetscInt *pStart, *pEnd;
2992: PetscInt dim, d, numSubVerticesInitial = 0, v;
2994: PetscFunctionBegin;
2995: PetscCall(DMGetDimension(dm, &dim));
2996: PetscCall(PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd));
2997: for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetSimplexOrBoxCells(dm, dim - d, &pStart[d], &pEnd[d]));
2998: /* Loop over initial vertices and mark all faces in the collective star() */
2999: if (vertexLabel) {
3000: PetscCall(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS));
3001: if (subvertexIS) {
3002: PetscCall(ISGetSize(subvertexIS, &numSubVerticesInitial));
3003: PetscCall(ISGetIndices(subvertexIS, &subvertices));
3004: }
3005: }
3006: for (v = 0; v < numSubVerticesInitial; ++v) {
3007: const PetscInt vertex = subvertices[v];
3008: PetscInt *star = NULL;
3009: PetscInt starSize, s, numFaces = 0, f;
3011: PetscCall(DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star));
3012: for (s = 0; s < starSize * 2; s += 2) {
3013: const PetscInt point = star[s];
3014: PetscInt faceLoc;
3016: if (point >= pStart[dim - 1] && point < pEnd[dim - 1]) {
3017: if (markedFaces) {
3018: PetscCall(DMLabelGetValue(vertexLabel, point, &faceLoc));
3019: if (faceLoc < 0) continue;
3020: }
3021: star[numFaces++] = point;
3022: }
3023: }
3024: for (f = 0; f < numFaces; ++f) {
3025: const PetscInt face = star[f];
3026: PetscInt *closure = NULL;
3027: PetscInt closureSize, c;
3028: PetscInt faceLoc;
3030: PetscCall(DMLabelGetValue(subpointMap, face, &faceLoc));
3031: if (faceLoc == dim - 1) continue;
3032: PetscCheck(faceLoc < 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %" PetscInt_FMT " has dimension %" PetscInt_FMT " in the surface label", face, faceLoc);
3033: PetscCall(DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure));
3034: for (c = 0; c < closureSize * 2; c += 2) {
3035: const PetscInt point = closure[c];
3036: PetscInt vertexLoc;
3038: if (point >= pStart[0] && point < pEnd[0]) {
3039: PetscCall(DMLabelGetValue(vertexLabel, point, &vertexLoc));
3040: if (vertexLoc != value) break;
3041: }
3042: }
3043: if (c == closureSize * 2) {
3044: const PetscInt *support;
3045: PetscInt supportSize, s;
3047: for (c = 0; c < closureSize * 2; c += 2) {
3048: const PetscInt point = closure[c];
3050: for (d = 0; d < dim; ++d) {
3051: if (point >= pStart[d] && point < pEnd[d]) {
3052: PetscCall(DMLabelSetValue(subpointMap, point, d));
3053: break;
3054: }
3055: }
3056: }
3057: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
3058: PetscCall(DMPlexGetSupport(dm, face, &support));
3059: if (addCells)
3060: for (s = 0; s < supportSize; ++s) PetscCall(DMLabelSetValue(subpointMap, support[s], dim));
3061: }
3062: PetscCall(DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure));
3063: }
3064: PetscCall(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star));
3065: }
3066: if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subvertices));
3067: PetscCall(ISDestroy(&subvertexIS));
3068: PetscCall(PetscFree2(pStart, pEnd));
3069: PetscFunctionReturn(PETSC_SUCCESS);
3070: }
3072: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm)
3073: {
3074: DMLabel label = NULL;
3075: const PetscInt *cone;
3076: PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1;
3078: PetscFunctionBegin;
3079: *numFaces = 0;
3080: *nFV = 0;
3081: if (labelname) PetscCall(DMGetLabel(dm, labelname, &label));
3082: *subCells = NULL;
3083: PetscCall(DMGetDimension(dm, &dim));
3084: PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd));
3085: if (cMax < 0) PetscFunctionReturn(PETSC_SUCCESS);
3086: if (label) {
3087: for (c = cMax; c < cEnd; ++c) {
3088: PetscInt val;
3090: PetscCall(DMLabelGetValue(label, c, &val));
3091: if (val == value) {
3092: ++(*numFaces);
3093: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
3094: }
3095: }
3096: } else {
3097: *numFaces = cEnd - cMax;
3098: PetscCall(DMPlexGetConeSize(dm, cMax, &coneSize));
3099: }
3100: PetscCall(PetscMalloc1(*numFaces * 2, subCells));
3101: if (!*numFaces) PetscFunctionReturn(PETSC_SUCCESS);
3102: *nFV = hasLagrange ? coneSize / 3 : coneSize / 2;
3103: for (c = cMax; c < cEnd; ++c) {
3104: const PetscInt *cells;
3105: PetscInt numCells;
3107: if (label) {
3108: PetscInt val;
3110: PetscCall(DMLabelGetValue(label, c, &val));
3111: if (val != value) continue;
3112: }
3113: PetscCall(DMPlexGetCone(dm, c, &cone));
3114: for (p = 0; p < *nFV; ++p) PetscCall(DMLabelSetValue(subpointMap, cone[p], 0));
3115: /* Negative face */
3116: PetscCall(DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells));
3117: /* Not true in parallel
3118: PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
3119: for (p = 0; p < numCells; ++p) {
3120: PetscCall(DMLabelSetValue(subpointMap, cells[p], 2));
3121: (*subCells)[subc++] = cells[p];
3122: }
3123: PetscCall(DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells));
3124: /* Positive face is not included */
3125: }
3126: PetscFunctionReturn(PETSC_SUCCESS);
3127: }
3129: static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm)
3130: {
3131: PetscInt *pStart, *pEnd;
3132: PetscInt dim, cMax, cEnd, c, d;
3134: PetscFunctionBegin;
3135: PetscCall(DMGetDimension(dm, &dim));
3136: PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd));
3137: if (cMax < 0) PetscFunctionReturn(PETSC_SUCCESS);
3138: PetscCall(PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd));
3139: for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]));
3140: for (c = cMax; c < cEnd; ++c) {
3141: const PetscInt *cone;
3142: PetscInt *closure = NULL;
3143: PetscInt fconeSize, coneSize, closureSize, cl, val;
3145: if (label) {
3146: PetscCall(DMLabelGetValue(label, c, &val));
3147: if (val != value) continue;
3148: }
3149: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
3150: PetscCall(DMPlexGetCone(dm, c, &cone));
3151: PetscCall(DMPlexGetConeSize(dm, cone[0], &fconeSize));
3152: PetscCheck(coneSize == (fconeSize ? fconeSize : 1) + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells");
3153: /* Negative face */
3154: PetscCall(DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure));
3155: for (cl = 0; cl < closureSize * 2; cl += 2) {
3156: const PetscInt point = closure[cl];
3158: for (d = 0; d <= dim; ++d) {
3159: if (point >= pStart[d] && point < pEnd[d]) {
3160: PetscCall(DMLabelSetValue(subpointMap, point, d));
3161: break;
3162: }
3163: }
3164: }
3165: PetscCall(DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure));
3166: /* Cells -- positive face is not included */
3167: for (cl = 0; cl < 1; ++cl) {
3168: const PetscInt *support;
3169: PetscInt supportSize, s;
3171: PetscCall(DMPlexGetSupportSize(dm, cone[cl], &supportSize));
3172: /* PetscCheck(supportSize == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */
3173: PetscCall(DMPlexGetSupport(dm, cone[cl], &support));
3174: for (s = 0; s < supportSize; ++s) PetscCall(DMLabelSetValue(subpointMap, support[s], dim));
3175: }
3176: }
3177: PetscCall(PetscFree2(pStart, pEnd));
3178: PetscFunctionReturn(PETSC_SUCCESS);
3179: }
3181: static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
3182: {
3183: MPI_Comm comm;
3184: PetscBool posOrient = PETSC_FALSE;
3185: const PetscInt debug = 0;
3186: PetscInt cellDim, faceSize, f;
3188: PetscFunctionBegin;
3189: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3190: PetscCall(DMGetDimension(dm, &cellDim));
3191: if (debug) PetscCall(PetscPrintf(comm, "cellDim: %" PetscInt_FMT " numCorners: %" PetscInt_FMT "\n", cellDim, numCorners));
3193: if (cellDim == 1 && numCorners == 2) {
3194: /* Triangle */
3195: faceSize = numCorners - 1;
3196: posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
3197: } else if (cellDim == 2 && numCorners == 3) {
3198: /* Triangle */
3199: faceSize = numCorners - 1;
3200: posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
3201: } else if (cellDim == 3 && numCorners == 4) {
3202: /* Tetrahedron */
3203: faceSize = numCorners - 1;
3204: posOrient = (oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE;
3205: } else if (cellDim == 1 && numCorners == 3) {
3206: /* Quadratic line */
3207: faceSize = 1;
3208: posOrient = PETSC_TRUE;
3209: } else if (cellDim == 2 && numCorners == 4) {
3210: /* Quads */
3211: faceSize = 2;
3212: if (indices[1] > indices[0] && indices[1] - indices[0] == 1) {
3213: posOrient = PETSC_TRUE;
3214: } else if (indices[0] == 3 && indices[1] == 0) {
3215: posOrient = PETSC_TRUE;
3216: } else {
3217: PetscCheck((indices[0] > indices[1] && indices[0] - indices[1] == 1) || (indices[0] == 0 && indices[1] == 3), comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
3218: posOrient = PETSC_FALSE;
3219: }
3220: } else if (cellDim == 2 && numCorners == 6) {
3221: /* Quadratic triangle (I hate this) */
3222: /* Edges are determined by the first 2 vertices (corners of edges) */
3223: const PetscInt faceSizeTri = 3;
3224: PetscInt sortedIndices[3], i, iFace;
3225: PetscBool found = PETSC_FALSE;
3226: PetscInt faceVerticesTriSorted[9] = {
3227: 0, 3, 4, /* bottom */
3228: 1, 4, 5, /* right */
3229: 2, 3, 5, /* left */
3230: };
3231: PetscInt faceVerticesTri[9] = {
3232: 0, 3, 4, /* bottom */
3233: 1, 4, 5, /* right */
3234: 2, 5, 3, /* left */
3235: };
3237: for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
3238: PetscCall(PetscSortInt(faceSizeTri, sortedIndices));
3239: for (iFace = 0; iFace < 3; ++iFace) {
3240: const PetscInt ii = iFace * faceSizeTri;
3241: PetscInt fVertex, cVertex;
3243: if (sortedIndices[0] == faceVerticesTriSorted[ii + 0] && sortedIndices[1] == faceVerticesTriSorted[ii + 1]) {
3244: for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
3245: for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
3246: if (indices[cVertex] == faceVerticesTri[ii + fVertex]) {
3247: faceVertices[fVertex] = origVertices[cVertex];
3248: break;
3249: }
3250: }
3251: }
3252: found = PETSC_TRUE;
3253: break;
3254: }
3255: }
3256: PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
3257: if (posOriented) *posOriented = PETSC_TRUE;
3258: PetscFunctionReturn(PETSC_SUCCESS);
3259: } else if (cellDim == 2 && numCorners == 9) {
3260: /* Quadratic quad (I hate this) */
3261: /* Edges are determined by the first 2 vertices (corners of edges) */
3262: const PetscInt faceSizeQuad = 3;
3263: PetscInt sortedIndices[3], i, iFace;
3264: PetscBool found = PETSC_FALSE;
3265: PetscInt faceVerticesQuadSorted[12] = {
3266: 0, 1, 4, /* bottom */
3267: 1, 2, 5, /* right */
3268: 2, 3, 6, /* top */
3269: 0, 3, 7, /* left */
3270: };
3271: PetscInt faceVerticesQuad[12] = {
3272: 0, 1, 4, /* bottom */
3273: 1, 2, 5, /* right */
3274: 2, 3, 6, /* top */
3275: 3, 0, 7, /* left */
3276: };
3278: for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
3279: PetscCall(PetscSortInt(faceSizeQuad, sortedIndices));
3280: for (iFace = 0; iFace < 4; ++iFace) {
3281: const PetscInt ii = iFace * faceSizeQuad;
3282: PetscInt fVertex, cVertex;
3284: if (sortedIndices[0] == faceVerticesQuadSorted[ii + 0] && sortedIndices[1] == faceVerticesQuadSorted[ii + 1]) {
3285: for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
3286: for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
3287: if (indices[cVertex] == faceVerticesQuad[ii + fVertex]) {
3288: faceVertices[fVertex] = origVertices[cVertex];
3289: break;
3290: }
3291: }
3292: }
3293: found = PETSC_TRUE;
3294: break;
3295: }
3296: }
3297: PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
3298: if (posOriented) *posOriented = PETSC_TRUE;
3299: PetscFunctionReturn(PETSC_SUCCESS);
3300: } else if (cellDim == 3 && numCorners == 8) {
3301: /* Hexes
3302: A hex is two oriented quads with the normal of the first
3303: pointing up at the second.
3305: 7---6
3306: /| /|
3307: 4---5 |
3308: | 1-|-2
3309: |/ |/
3310: 0---3
3312: Faces are determined by the first 4 vertices (corners of faces) */
3313: const PetscInt faceSizeHex = 4;
3314: PetscInt sortedIndices[4], i, iFace;
3315: PetscBool found = PETSC_FALSE;
3316: PetscInt faceVerticesHexSorted[24] = {
3317: 0, 1, 2, 3, /* bottom */
3318: 4, 5, 6, 7, /* top */
3319: 0, 3, 4, 5, /* front */
3320: 2, 3, 5, 6, /* right */
3321: 1, 2, 6, 7, /* back */
3322: 0, 1, 4, 7, /* left */
3323: };
3324: PetscInt faceVerticesHex[24] = {
3325: 1, 2, 3, 0, /* bottom */
3326: 4, 5, 6, 7, /* top */
3327: 0, 3, 5, 4, /* front */
3328: 3, 2, 6, 5, /* right */
3329: 2, 1, 7, 6, /* back */
3330: 1, 0, 4, 7, /* left */
3331: };
3333: for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
3334: PetscCall(PetscSortInt(faceSizeHex, sortedIndices));
3335: for (iFace = 0; iFace < 6; ++iFace) {
3336: const PetscInt ii = iFace * faceSizeHex;
3337: PetscInt fVertex, cVertex;
3339: if (sortedIndices[0] == faceVerticesHexSorted[ii + 0] && sortedIndices[1] == faceVerticesHexSorted[ii + 1] && sortedIndices[2] == faceVerticesHexSorted[ii + 2] && sortedIndices[3] == faceVerticesHexSorted[ii + 3]) {
3340: for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
3341: for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
3342: if (indices[cVertex] == faceVerticesHex[ii + fVertex]) {
3343: faceVertices[fVertex] = origVertices[cVertex];
3344: break;
3345: }
3346: }
3347: }
3348: found = PETSC_TRUE;
3349: break;
3350: }
3351: }
3352: PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
3353: if (posOriented) *posOriented = PETSC_TRUE;
3354: PetscFunctionReturn(PETSC_SUCCESS);
3355: } else if (cellDim == 3 && numCorners == 10) {
3356: /* Quadratic tet */
3357: /* Faces are determined by the first 3 vertices (corners of faces) */
3358: const PetscInt faceSizeTet = 6;
3359: PetscInt sortedIndices[6], i, iFace;
3360: PetscBool found = PETSC_FALSE;
3361: PetscInt faceVerticesTetSorted[24] = {
3362: 0, 1, 2, 6, 7, 8, /* bottom */
3363: 0, 3, 4, 6, 7, 9, /* front */
3364: 1, 4, 5, 7, 8, 9, /* right */
3365: 2, 3, 5, 6, 8, 9, /* left */
3366: };
3367: PetscInt faceVerticesTet[24] = {
3368: 0, 1, 2, 6, 7, 8, /* bottom */
3369: 0, 4, 3, 6, 7, 9, /* front */
3370: 1, 5, 4, 7, 8, 9, /* right */
3371: 2, 3, 5, 8, 6, 9, /* left */
3372: };
3374: for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
3375: PetscCall(PetscSortInt(faceSizeTet, sortedIndices));
3376: for (iFace = 0; iFace < 4; ++iFace) {
3377: const PetscInt ii = iFace * faceSizeTet;
3378: PetscInt fVertex, cVertex;
3380: if (sortedIndices[0] == faceVerticesTetSorted[ii + 0] && sortedIndices[1] == faceVerticesTetSorted[ii + 1] && sortedIndices[2] == faceVerticesTetSorted[ii + 2] && sortedIndices[3] == faceVerticesTetSorted[ii + 3]) {
3381: for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
3382: for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
3383: if (indices[cVertex] == faceVerticesTet[ii + fVertex]) {
3384: faceVertices[fVertex] = origVertices[cVertex];
3385: break;
3386: }
3387: }
3388: }
3389: found = PETSC_TRUE;
3390: break;
3391: }
3392: }
3393: PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
3394: if (posOriented) *posOriented = PETSC_TRUE;
3395: PetscFunctionReturn(PETSC_SUCCESS);
3396: } else if (cellDim == 3 && numCorners == 27) {
3397: /* Quadratic hexes (I hate this)
3398: A hex is two oriented quads with the normal of the first
3399: pointing up at the second.
3401: 7---6
3402: /| /|
3403: 4---5 |
3404: | 3-|-2
3405: |/ |/
3406: 0---1
3408: Faces are determined by the first 4 vertices (corners of faces) */
3409: const PetscInt faceSizeQuadHex = 9;
3410: PetscInt sortedIndices[9], i, iFace;
3411: PetscBool found = PETSC_FALSE;
3412: PetscInt faceVerticesQuadHexSorted[54] = {
3413: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
3414: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
3415: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
3416: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
3417: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
3418: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
3419: };
3420: PetscInt faceVerticesQuadHex[54] = {
3421: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
3422: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
3423: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
3424: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
3425: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
3426: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
3427: };
3429: for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
3430: PetscCall(PetscSortInt(faceSizeQuadHex, sortedIndices));
3431: for (iFace = 0; iFace < 6; ++iFace) {
3432: const PetscInt ii = iFace * faceSizeQuadHex;
3433: PetscInt fVertex, cVertex;
3435: if (sortedIndices[0] == faceVerticesQuadHexSorted[ii + 0] && sortedIndices[1] == faceVerticesQuadHexSorted[ii + 1] && sortedIndices[2] == faceVerticesQuadHexSorted[ii + 2] && sortedIndices[3] == faceVerticesQuadHexSorted[ii + 3]) {
3436: for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
3437: for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
3438: if (indices[cVertex] == faceVerticesQuadHex[ii + fVertex]) {
3439: faceVertices[fVertex] = origVertices[cVertex];
3440: break;
3441: }
3442: }
3443: }
3444: found = PETSC_TRUE;
3445: break;
3446: }
3447: }
3448: PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
3449: if (posOriented) *posOriented = PETSC_TRUE;
3450: PetscFunctionReturn(PETSC_SUCCESS);
3451: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
3452: if (!posOrient) {
3453: if (debug) PetscCall(PetscPrintf(comm, " Reversing initial face orientation\n"));
3454: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize - 1 - f];
3455: } else {
3456: if (debug) PetscCall(PetscPrintf(comm, " Keeping initial face orientation\n"));
3457: for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f];
3458: }
3459: if (posOriented) *posOriented = posOrient;
3460: PetscFunctionReturn(PETSC_SUCCESS);
3461: }
3463: /*@
3464: DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices,
3465: in faceVertices. The orientation is such that the face normal points out of the cell
3467: Not Collective
3469: Input Parameters:
3470: + dm - The original mesh
3471: . cell - The cell mesh point
3472: . faceSize - The number of vertices on the face
3473: . face - The face vertices
3474: . numCorners - The number of vertices on the cell
3475: . indices - Local numbering of face vertices in cell cone
3476: - origVertices - Original face vertices
3478: Output Parameters:
3479: + faceVertices - The face vertices properly oriented
3480: - posOriented - `PETSC_TRUE` if the face was oriented with outward normal
3482: Level: developer
3484: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetCone()`
3485: @*/
3486: PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
3487: {
3488: const PetscInt *cone = NULL;
3489: PetscInt coneSize, v, f, v2;
3490: PetscInt oppositeVertex = -1;
3492: PetscFunctionBegin;
3493: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3494: PetscCall(DMPlexGetCone(dm, cell, &cone));
3495: for (v = 0, v2 = 0; v < coneSize; ++v) {
3496: PetscBool found = PETSC_FALSE;
3498: for (f = 0; f < faceSize; ++f) {
3499: if (face[f] == cone[v]) {
3500: found = PETSC_TRUE;
3501: break;
3502: }
3503: }
3504: if (found) {
3505: indices[v2] = v;
3506: origVertices[v2] = cone[v];
3507: ++v2;
3508: } else {
3509: oppositeVertex = v;
3510: }
3511: }
3512: PetscCall(DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented));
3513: PetscFunctionReturn(PETSC_SUCCESS);
3514: }
3516: /*
3517: DMPlexInsertFace_Internal - Puts a face into the mesh
3519: Not Collective
3521: Input Parameters:
3522: + dm - The `DMPLEX`
3523: . numFaceVertex - The number of vertices in the face
3524: . faceVertices - The vertices in the face for `dm`
3525: . subfaceVertices - The vertices in the face for subdm
3526: . numCorners - The number of vertices in the `cell`
3527: . cell - A cell in `dm` containing the face
3528: . subcell - A cell in subdm containing the face
3529: . firstFace - First face in the mesh
3530: - newFacePoint - Next face in the mesh
3532: Output Parameter:
3533: . newFacePoint - Contains next face point number on input, updated on output
3535: Level: developer
3536: */
3537: 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)
3538: {
3539: MPI_Comm comm;
3540: DM_Plex *submesh = (DM_Plex *)subdm->data;
3541: const PetscInt *faces;
3542: PetscInt numFaces, coneSize;
3544: PetscFunctionBegin;
3545: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3546: PetscCall(DMPlexGetConeSize(subdm, subcell, &coneSize));
3547: PetscCheck(coneSize == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %" PetscInt_FMT " is %" PetscInt_FMT " != 1", cell, coneSize);
3548: #if 0
3549: /* Cannot use this because support() has not been constructed yet */
3550: PetscCall(DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces));
3551: #else
3552: {
3553: PetscInt f;
3555: numFaces = 0;
3556: PetscCall(DMGetWorkArray(subdm, 1, MPIU_INT, (void **)&faces));
3557: for (f = firstFace; f < *newFacePoint; ++f) {
3558: PetscInt dof, off, d;
3560: PetscCall(PetscSectionGetDof(submesh->coneSection, f, &dof));
3561: PetscCall(PetscSectionGetOffset(submesh->coneSection, f, &off));
3562: /* Yes, I know this is quadratic, but I expect the sizes to be <5 */
3563: for (d = 0; d < dof; ++d) {
3564: const PetscInt p = submesh->cones[off + d];
3565: PetscInt v;
3567: for (v = 0; v < numFaceVertices; ++v) {
3568: if (subfaceVertices[v] == p) break;
3569: }
3570: if (v == numFaceVertices) break;
3571: }
3572: if (d == dof) {
3573: numFaces = 1;
3574: ((PetscInt *)faces)[0] = f;
3575: }
3576: }
3577: }
3578: #endif
3579: PetscCheck(numFaces <= 1, comm, PETSC_ERR_ARG_WRONG, "Vertex set had %" PetscInt_FMT " faces, not one", numFaces);
3580: if (numFaces == 1) {
3581: /* Add the other cell neighbor for this face */
3582: PetscCall(DMPlexSetCone(subdm, subcell, faces));
3583: } else {
3584: PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov;
3585: PetscBool posOriented;
3587: PetscCall(DMGetWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices));
3588: origVertices = &orientedVertices[numFaceVertices];
3589: indices = &orientedVertices[numFaceVertices * 2];
3590: orientedSubVertices = &orientedVertices[numFaceVertices * 3];
3591: PetscCall(DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented));
3592: /* TODO: I know that routine should return a permutation, not the indices */
3593: for (v = 0; v < numFaceVertices; ++v) {
3594: const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v];
3595: for (ov = 0; ov < numFaceVertices; ++ov) {
3596: if (orientedVertices[ov] == vertex) {
3597: orientedSubVertices[ov] = subvertex;
3598: break;
3599: }
3600: }
3601: PetscCheck(ov != numFaceVertices, comm, PETSC_ERR_PLIB, "Could not find face vertex %" PetscInt_FMT " in orientated set", vertex);
3602: }
3603: PetscCall(DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices));
3604: PetscCall(DMPlexSetCone(subdm, subcell, newFacePoint));
3605: PetscCall(DMRestoreWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices));
3606: ++(*newFacePoint);
3607: }
3608: #if 0
3609: PetscCall(DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces));
3610: #else
3611: PetscCall(DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **)&faces));
3612: #endif
3613: PetscFunctionReturn(PETSC_SUCCESS);
3614: }
3616: static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm)
3617: {
3618: MPI_Comm comm;
3619: DMLabel subpointMap;
3620: IS subvertexIS, subcellIS;
3621: const PetscInt *subVertices, *subCells;
3622: PetscInt numSubVertices, firstSubVertex, numSubCells;
3623: PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0;
3624: PetscInt vStart, vEnd, c, f;
3626: PetscFunctionBegin;
3627: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3628: /* Create subpointMap which marks the submesh */
3629: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap));
3630: PetscCall(DMPlexSetSubpointMap(subdm, subpointMap));
3631: if (vertexLabel) PetscCall(DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm));
3632: /* Setup chart */
3633: PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices));
3634: PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells));
3635: PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices));
3636: PetscCall(DMPlexSetVTKCellHeight(subdm, 1));
3637: /* Set cone sizes */
3638: firstSubVertex = numSubCells;
3639: firstSubFace = numSubCells + numSubVertices;
3640: newFacePoint = firstSubFace;
3641: PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS));
3642: if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices));
3643: PetscCall(DMLabelGetStratumIS(subpointMap, 2, &subcellIS));
3644: if (subcellIS) PetscCall(ISGetIndices(subcellIS, &subCells));
3645: for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1));
3646: for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV));
3647: PetscCall(DMSetUp(subdm));
3648: PetscCall(DMLabelDestroy(&subpointMap));
3649: /* Create face cones */
3650: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
3651: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
3652: PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface));
3653: for (c = 0; c < numSubCells; ++c) {
3654: const PetscInt cell = subCells[c];
3655: const PetscInt subcell = c;
3656: PetscInt *closure = NULL;
3657: PetscInt closureSize, cl, numCorners = 0, faceSize = 0;
3659: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
3660: for (cl = 0; cl < closureSize * 2; cl += 2) {
3661: const PetscInt point = closure[cl];
3662: PetscInt subVertex;
3664: if (point >= vStart && point < vEnd) {
3665: ++numCorners;
3666: PetscCall(PetscFindInt(point, numSubVertices, subVertices, &subVertex));
3667: if (subVertex >= 0) {
3668: closure[faceSize] = point;
3669: subface[faceSize] = firstSubVertex + subVertex;
3670: ++faceSize;
3671: }
3672: }
3673: }
3674: PetscCheck(faceSize <= nFV, comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %" PetscInt_FMT " of an element on the surface", faceSize);
3675: if (faceSize == nFV) PetscCall(DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint));
3676: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
3677: }
3678: PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface));
3679: PetscCall(DMPlexSymmetrize(subdm));
3680: PetscCall(DMPlexStratify(subdm));
3681: /* Build coordinates */
3682: {
3683: PetscSection coordSection, subCoordSection;
3684: Vec coordinates, subCoordinates;
3685: PetscScalar *coords, *subCoords;
3686: PetscInt numComp, coordSize, v;
3687: const char *name;
3689: PetscCall(DMGetCoordinateSection(dm, &coordSection));
3690: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
3691: PetscCall(DMGetCoordinateSection(subdm, &subCoordSection));
3692: PetscCall(PetscSectionSetNumFields(subCoordSection, 1));
3693: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp));
3694: PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp));
3695: PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices));
3696: for (v = 0; v < numSubVertices; ++v) {
3697: const PetscInt vertex = subVertices[v];
3698: const PetscInt subvertex = firstSubVertex + v;
3699: PetscInt dof;
3701: PetscCall(PetscSectionGetDof(coordSection, vertex, &dof));
3702: PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof));
3703: PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof));
3704: }
3705: PetscCall(PetscSectionSetUp(subCoordSection));
3706: PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize));
3707: PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates));
3708: PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
3709: PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name));
3710: PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE));
3711: PetscCall(VecSetType(subCoordinates, VECSTANDARD));
3712: if (coordSize) {
3713: PetscCall(VecGetArray(coordinates, &coords));
3714: PetscCall(VecGetArray(subCoordinates, &subCoords));
3715: for (v = 0; v < numSubVertices; ++v) {
3716: const PetscInt vertex = subVertices[v];
3717: const PetscInt subvertex = firstSubVertex + v;
3718: PetscInt dof, off, sdof, soff, d;
3720: PetscCall(PetscSectionGetDof(coordSection, vertex, &dof));
3721: PetscCall(PetscSectionGetOffset(coordSection, vertex, &off));
3722: PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof));
3723: PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff));
3724: 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);
3725: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
3726: }
3727: PetscCall(VecRestoreArray(coordinates, &coords));
3728: PetscCall(VecRestoreArray(subCoordinates, &subCoords));
3729: }
3730: PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates));
3731: PetscCall(VecDestroy(&subCoordinates));
3732: }
3733: /* Cleanup */
3734: if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices));
3735: PetscCall(ISDestroy(&subvertexIS));
3736: if (subcellIS) PetscCall(ISRestoreIndices(subcellIS, &subCells));
3737: PetscCall(ISDestroy(&subcellIS));
3738: PetscFunctionReturn(PETSC_SUCCESS);
3739: }
3741: /* TODO: Fix this to properly propagate up error conditions it may find */
3742: static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[])
3743: {
3744: PetscInt subPoint;
3745: PetscErrorCode ierr;
3747: ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3748: if (ierr) return -1;
3749: return subPoint < 0 ? subPoint : firstSubPoint + subPoint;
3750: }
3752: /* TODO: Fix this to properly propagate up error conditions it may find */
3753: static inline PetscInt DMPlexFilterPointPerm_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[], const PetscInt subIndices[])
3754: {
3755: PetscInt subPoint;
3756: PetscErrorCode ierr;
3758: ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint);
3759: if (ierr) return -1;
3760: return subPoint < 0 ? subPoint : firstSubPoint + (subIndices ? subIndices[subPoint] : subPoint);
3761: }
3763: static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm)
3764: {
3765: DMLabel depthLabel;
3766: PetscInt Nl, l, d;
3768: PetscFunctionBegin;
3769: // Reset depth label for fast lookup
3770: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
3771: PetscCall(DMLabelMakeAllInvalid_Internal(depthLabel));
3772: PetscCall(DMGetNumLabels(dm, &Nl));
3773: for (l = 0; l < Nl; ++l) {
3774: DMLabel label, newlabel;
3775: const char *lname;
3776: PetscBool isDepth, isDim, isCelltype, isVTK;
3777: IS valueIS;
3778: const PetscInt *values;
3779: PetscInt Nv, v;
3781: PetscCall(DMGetLabelName(dm, l, &lname));
3782: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
3783: PetscCall(PetscStrcmp(lname, "dim", &isDim));
3784: PetscCall(PetscStrcmp(lname, "celltype", &isCelltype));
3785: PetscCall(PetscStrcmp(lname, "vtk", &isVTK));
3786: if (isDepth || isDim || isCelltype || isVTK) continue;
3787: PetscCall(DMCreateLabel(subdm, lname));
3788: PetscCall(DMGetLabel(dm, lname, &label));
3789: PetscCall(DMGetLabel(subdm, lname, &newlabel));
3790: PetscCall(DMLabelGetDefaultValue(label, &v));
3791: PetscCall(DMLabelSetDefaultValue(newlabel, v));
3792: PetscCall(DMLabelGetValueIS(label, &valueIS));
3793: PetscCall(ISGetLocalSize(valueIS, &Nv));
3794: PetscCall(ISGetIndices(valueIS, &values));
3795: for (v = 0; v < Nv; ++v) {
3796: IS pointIS;
3797: const PetscInt *points;
3798: PetscInt Np, p;
3800: PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS));
3801: PetscCall(ISGetLocalSize(pointIS, &Np));
3802: PetscCall(ISGetIndices(pointIS, &points));
3803: for (p = 0; p < Np; ++p) {
3804: const PetscInt point = points[p];
3805: PetscInt subp, subdepth;
3807: PetscCall(DMPlexGetPointDepth(dm, point, &d));
3808: PetscCall(DMPlexGetDepth(subdm, &subdepth));
3809: if (d > subdepth) continue;
3810: subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]);
3811: if (subp >= 0) PetscCall(DMLabelSetValue(newlabel, subp, values[v]));
3812: }
3813: PetscCall(ISRestoreIndices(pointIS, &points));
3814: PetscCall(ISDestroy(&pointIS));
3815: }
3816: PetscCall(ISRestoreIndices(valueIS, &values));
3817: PetscCall(ISDestroy(&valueIS));
3818: }
3819: PetscFunctionReturn(PETSC_SUCCESS);
3820: }
3822: static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, PetscSF *ownershipTransferSF, DM subdm)
3823: {
3824: MPI_Comm comm;
3825: DMLabel subpointMap;
3826: IS *subpointIS;
3827: const PetscInt **subpoints;
3828: PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew;
3829: PetscInt totSubPoints = 0, maxConeSize, dim, sdim, cdim, p, d, coordinate_type;
3830: PetscMPIInt rank;
3832: PetscFunctionBegin;
3833: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
3834: PetscCallMPI(MPI_Comm_rank(comm, &rank));
3835: /* Create subpointMap which marks the submesh */
3836: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap));
3837: PetscCall(DMPlexSetSubpointMap(subdm, subpointMap));
3838: if (cellHeight) {
3839: if (isCohesive) PetscCall(DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm));
3840: else PetscCall(DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, PETSC_TRUE, subpointMap, subdm));
3841: } else {
3842: DMLabel depth;
3843: IS pointIS;
3844: const PetscInt *points;
3845: PetscInt numPoints = 0, pStart, pEnd;
3846: PetscBool *isGhost = NULL;
3848: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
3849: PetscCall(DMPlexGetDepthLabel(dm, &depth));
3850: if (label) PetscCall(DMLabelGetStratumIS(label, value, &pointIS));
3851: else {
3852: PetscInt cStart, cEnd;
3854: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3855: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cStart, cStart, 1, &pointIS));
3856: }
3857: if (pointIS) {
3858: PetscCall(ISGetIndices(pointIS, &points));
3859: PetscCall(ISGetLocalSize(pointIS, &numPoints));
3860: }
3861: if (ignoreLabelHalo) {
3862: PetscSF pointSF;
3863: PetscInt nleaves;
3864: const PetscInt *ilocal = NULL;
3866: PetscCall(DMGetPointSF(dm, &pointSF));
3867: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, NULL));
3868: PetscCall(PetscMalloc1(pEnd - pStart, &isGhost));
3869: for (p = 0; p < pEnd - pStart; ++p) isGhost[p] = PETSC_FALSE;
3870: for (p = 0; p < nleaves; ++p) isGhost[(ilocal ? ilocal[p] : p) - pStart] = PETSC_TRUE;
3871: }
3872: for (p = 0; p < numPoints; ++p) {
3873: PetscInt *closure = NULL;
3874: PetscInt closureSize, c, pdim;
3876: if (isGhost && isGhost[points[p] - pStart]) continue;
3877: PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure));
3878: for (c = 0; c < closureSize * 2; c += 2) {
3879: PetscCall(DMLabelGetValue(depth, closure[c], &pdim));
3880: PetscCall(DMLabelSetValue(subpointMap, closure[c], pdim));
3881: }
3882: PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure));
3883: }
3884: PetscCall(PetscFree(isGhost));
3885: if (pointIS) PetscCall(ISRestoreIndices(pointIS, &points));
3886: PetscCall(ISDestroy(&pointIS));
3887: }
3888: /* Setup chart */
3889: PetscCall(DMGetDimension(dm, &dim));
3890: PetscCall(DMGetCoordinateDim(dm, &cdim));
3891: PetscCall(PetscMalloc4(dim + 1, &numSubPoints, dim + 1, &firstSubPoint, dim + 1, &subpointIS, dim + 1, &subpoints));
3892: for (d = 0; d <= dim; ++d) {
3893: PetscCall(DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]));
3894: totSubPoints += numSubPoints[d];
3895: }
3896: // Determine submesh dimension
3897: PetscCall(DMGetDimension(subdm, &sdim));
3898: if (sdim > 0) {
3899: // Calling function knows what dimension to use, and we include neighboring cells as well
3900: sdim = dim;
3901: } else {
3902: // We reset the subdimension based on what is being selected
3903: PetscInt lsdim;
3904: for (lsdim = dim; lsdim >= 0; --lsdim)
3905: if (numSubPoints[lsdim]) break;
3906: PetscCallMPI(MPIU_Allreduce(&lsdim, &sdim, 1, MPIU_INT, MPI_MAX, comm));
3907: PetscCall(DMSetDimension(subdm, sdim));
3908: PetscCall(DMSetCoordinateDim(subdm, cdim));
3909: }
3910: PetscCall(DMPlexSetChart(subdm, 0, totSubPoints));
3911: PetscCall(DMPlexSetVTKCellHeight(subdm, cellHeight));
3912: /* Set cone sizes */
3913: firstSubPoint[sdim] = 0;
3914: firstSubPoint[0] = firstSubPoint[sdim] + numSubPoints[sdim];
3915: if (sdim > 1) firstSubPoint[sdim - 1] = firstSubPoint[0] + numSubPoints[0];
3916: if (sdim > 2) firstSubPoint[sdim - 2] = firstSubPoint[sdim - 1] + numSubPoints[sdim - 1];
3917: for (d = 0; d <= sdim; ++d) {
3918: PetscCall(DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]));
3919: if (subpointIS[d]) PetscCall(ISGetIndices(subpointIS[d], &subpoints[d]));
3920: }
3921: /* We do not want this label automatically computed, instead we compute it here */
3922: PetscCall(DMCreateLabel(subdm, "celltype"));
3923: for (d = 0; d <= sdim; ++d) {
3924: for (p = 0; p < numSubPoints[d]; ++p) {
3925: const PetscInt point = subpoints[d][p];
3926: const PetscInt subpoint = firstSubPoint[d] + p;
3927: const PetscInt *cone;
3928: PetscInt coneSize;
3930: PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
3931: if (cellHeight && (d == sdim)) {
3932: PetscInt coneSizeNew, c, val;
3934: PetscCall(DMPlexGetCone(dm, point, &cone));
3935: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3936: PetscCall(DMLabelGetValue(subpointMap, cone[c], &val));
3937: if (val >= 0) coneSizeNew++;
3938: }
3939: PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSizeNew));
3940: PetscCall(DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST));
3941: } else {
3942: DMPolytopeType ct;
3944: PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSize));
3945: PetscCall(DMPlexGetCellType(dm, point, &ct));
3946: PetscCall(DMPlexSetCellType(subdm, subpoint, ct));
3947: }
3948: }
3949: }
3950: PetscCall(DMLabelDestroy(&subpointMap));
3951: PetscCall(DMSetUp(subdm));
3952: /* Set cones */
3953: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
3954: PetscCall(PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew));
3955: for (d = 0; d <= sdim; ++d) {
3956: for (p = 0; p < numSubPoints[d]; ++p) {
3957: const PetscInt point = subpoints[d][p];
3958: const PetscInt subpoint = firstSubPoint[d] + p;
3959: const PetscInt *cone, *ornt;
3960: PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0;
3962: if (d == sdim - 1) {
3963: const PetscInt *support, *cone, *ornt;
3964: PetscInt supportSize, coneSize, s, subc;
3966: PetscCall(DMPlexGetSupport(dm, point, &support));
3967: PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
3968: for (s = 0; s < supportSize; ++s) {
3969: PetscBool isHybrid = PETSC_FALSE;
3971: PetscCall(DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid));
3972: if (!isHybrid) continue;
3973: PetscCall(PetscFindInt(support[s], numSubPoints[d + 1], subpoints[d + 1], &subc));
3974: if (subc >= 0) {
3975: const PetscInt ccell = subpoints[d + 1][subc];
3977: PetscCall(DMPlexGetCone(dm, ccell, &cone));
3978: PetscCall(DMPlexGetConeSize(dm, ccell, &coneSize));
3979: PetscCall(DMPlexGetConeOrientation(dm, ccell, &ornt));
3980: for (c = 0; c < coneSize; ++c) {
3981: if (cone[c] == point) {
3982: fornt = ornt[c];
3983: break;
3984: }
3985: }
3986: break;
3987: }
3988: }
3989: }
3990: PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
3991: PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize));
3992: PetscCall(DMPlexGetCone(dm, point, &cone));
3993: PetscCall(DMPlexGetConeOrientation(dm, point, &ornt));
3994: for (c = 0, coneSizeNew = 0; c < coneSize; ++c) {
3995: PetscCall(PetscFindInt(cone[c], numSubPoints[d - 1], subpoints[d - 1], &subc));
3996: if (subc >= 0) {
3997: coneNew[coneSizeNew] = firstSubPoint[d - 1] + subc;
3998: orntNew[coneSizeNew] = ornt[c];
3999: ++coneSizeNew;
4000: }
4001: }
4002: PetscCheck(coneSizeNew == subconeSize, comm, PETSC_ERR_PLIB, "Number of cone points located %" PetscInt_FMT " does not match subcone size %" PetscInt_FMT, coneSizeNew, subconeSize);
4003: PetscCall(DMPlexSetCone(subdm, subpoint, coneNew));
4004: PetscCall(DMPlexSetConeOrientation(subdm, subpoint, orntNew));
4005: if (fornt < 0) PetscCall(DMPlexOrientPoint(subdm, subpoint, fornt));
4006: }
4007: }
4008: PetscCall(PetscFree2(coneNew, orntNew));
4009: PetscCall(DMPlexSymmetrize(subdm));
4010: PetscCall(DMPlexStratify(subdm));
4011: /* Build coordinates */
4012: for (coordinate_type = 0; coordinate_type < 2; ++coordinate_type) {
4013: DM coordDM = NULL, subCoordDM = NULL;
4014: PetscSection coordSection = NULL, subCoordSection = NULL;
4015: Vec coordinates = NULL, subCoordinates = NULL;
4016: PetscScalar *coords = NULL, *subCoords = NULL;
4017: PetscInt bs, numComp, coordSize, firstP, lastP, firstSubP = totSubPoints, lastSubP = -1, numFields;
4018: const char *name;
4019: PetscBool localized = (PetscBool)coordinate_type;
4021: if (!dm->coordinates[coordinate_type].dm) continue;
4022: if (!localized) {
4023: PetscCall(DMGetCoordinateDM(dm, &coordDM));
4024: PetscCall(DMGetCoordinateDM(subdm, &subCoordDM));
4025: } else {
4026: {
4027: PetscInt localizationHeight;
4028: PetscBool sparseLocalize;
4030: PetscCall(DMGetSparseLocalize(dm, &sparseLocalize));
4031: PetscCall(DMSetSparseLocalize(subdm, sparseLocalize));
4032: PetscCall(DMGetCoordinateDM(dm, &coordDM));
4033: PetscCall(DMGetCoordinateDM(subdm, &subCoordDM));
4034: PetscCall(DMPlexGetMaxProjectionHeight(coordDM, &localizationHeight));
4035: PetscCall(DMPlexSetMaxProjectionHeight(subCoordDM, localizationHeight - (dim - sdim) - cellHeight));
4036: PetscUseTypeMethod(subdm, createcellcoordinatedm, &subCoordDM);
4037: PetscCall(DMSetCellCoordinateDM(subdm, subCoordDM));
4038: PetscCall(DMDestroy(&subCoordDM));
4039: }
4040: PetscCall(DMGetCellCoordinateDM(dm, &coordDM));
4041: PetscCall(DMGetCellCoordinateDM(subdm, &subCoordDM));
4042: }
4043: PetscCall(DMGetNumFields(coordDM, &numFields));
4044: if (numFields > 0) {
4045: PetscFE fe = NULL;
4046: PetscSpace P = NULL;
4047: PetscClassId id;
4048: PetscInt degree;
4050: PetscCall(DMGetField(coordDM, 0, NULL, (PetscObject *)&fe));
4051: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
4052: if (id == PETSCFE_CLASSID) {
4053: if (sdim == dim && cellHeight == 0) {
4054: /* TODO: Handle Field labels correctly */
4055: PetscCall(DMSetField(subCoordDM, 0, NULL, (PetscObject)fe));
4056: PetscCall(DMCreateDS(subCoordDM));
4057: } else {
4058: /* TODO: Reconstruct the lower-dimensional FE more robustly */
4059: PetscCall(PetscFEGetBasisSpace(fe, &P));
4060: PetscCall(PetscSpaceGetDegree(P, °ree, NULL));
4061: PetscCall(DMPlexCreateCoordinateSpace(subdm, degree, localized, PETSC_FALSE));
4062: }
4063: }
4064: }
4065: if (!localized) {
4066: PetscCall(DMGetCoordinateSection(dm, &coordSection));
4067: PetscCall(DMGetCoordinateSection(subdm, &subCoordSection));
4068: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
4069: } else {
4070: PetscCall(DMGetCellCoordinateSection(dm, &coordSection));
4071: PetscCall(DMGetCellCoordinateSection(subdm, &subCoordSection));
4072: PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates));
4073: }
4074: PetscCall(PetscSectionSetNumFields(subCoordSection, 1));
4075: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp));
4076: PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp));
4077: PetscCall(PetscSectionGetChart(coordSection, &firstP, &lastP));
4078: for (d = 0; d <= sdim; ++d) {
4079: for (p = 0; p < numSubPoints[d]; ++p) {
4080: const PetscInt point = subpoints[d][p];
4081: const PetscInt subpoint = firstSubPoint[d] + p;
4083: if (point >= firstP && point < lastP) {
4084: firstSubP = PetscMin(firstSubP, subpoint);
4085: lastSubP = PetscMax(lastSubP, subpoint);
4086: }
4087: }
4088: }
4089: lastSubP += 1;
4090: if (firstSubP == totSubPoints) {
4091: /* Zero if there is no coordinate point. */
4092: firstSubP = 0;
4093: lastSubP = 0;
4094: }
4095: PetscCall(PetscSectionSetChart(subCoordSection, firstSubP, lastSubP));
4096: for (d = 0; d <= sdim; ++d) {
4097: for (p = 0; p < numSubPoints[d]; ++p) {
4098: const PetscInt point = subpoints[d][p];
4099: const PetscInt subpoint = firstSubPoint[d] + p;
4100: PetscInt dof;
4102: if (point >= firstP && point < lastP) {
4103: PetscCall(PetscSectionGetDof(coordSection, point, &dof));
4104: PetscCall(PetscSectionSetDof(subCoordSection, subpoint, dof));
4105: PetscCall(PetscSectionSetFieldDof(subCoordSection, subpoint, 0, dof));
4106: }
4107: }
4108: }
4109: PetscCall(PetscSectionSetUp(subCoordSection));
4110: PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize));
4111: PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates));
4112: PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
4113: PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name));
4114: PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE));
4115: PetscCall(VecGetBlockSize(coordinates, &bs));
4116: PetscCall(VecSetBlockSize(subCoordinates, bs));
4117: PetscCall(VecSetType(subCoordinates, VECSTANDARD));
4118: PetscCall(VecGetArray(coordinates, &coords));
4119: PetscCall(VecGetArray(subCoordinates, &subCoords));
4120: for (d = 0; d <= sdim; ++d) {
4121: for (p = 0; p < numSubPoints[d]; ++p) {
4122: const PetscInt point = subpoints[d][p];
4123: const PetscInt subpoint = firstSubPoint[d] + p;
4124: PetscInt dof, off, sdof, soff, d;
4126: if (point >= firstP && point < lastP) {
4127: PetscCall(PetscSectionGetDof(coordSection, point, &dof));
4128: PetscCall(PetscSectionGetOffset(coordSection, point, &off));
4129: PetscCall(PetscSectionGetDof(subCoordSection, subpoint, &sdof));
4130: PetscCall(PetscSectionGetOffset(subCoordSection, subpoint, &soff));
4131: PetscCheck(dof == sdof, comm, PETSC_ERR_PLIB, "Coordinate dimension %" PetscInt_FMT " on subpoint %" PetscInt_FMT ", point %" PetscInt_FMT " should be %" PetscInt_FMT, sdof, subpoint, point, dof);
4132: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
4133: }
4134: }
4135: }
4136: PetscCall(VecRestoreArray(coordinates, &coords));
4137: PetscCall(VecRestoreArray(subCoordinates, &subCoords));
4138: switch (coordinate_type) {
4139: case 0:
4140: PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates));
4141: break;
4142: case 1:
4143: PetscCall(DMSetCellCoordinatesLocal(subdm, subCoordinates));
4144: break;
4145: default:
4146: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "coordinate_type must be {0, 1}");
4147: }
4148: PetscCall(VecDestroy(&subCoordinates));
4149: }
4150: /* Build SF: We need this complexity because subpoints might not be selected on the owning process */
4151: // TODO parallel subcommunicators
4152: PetscMPIInt flag;
4153: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)dm), PetscObjectComm((PetscObject)subdm), &flag));
4154: if (flag == MPI_CONGRUENT || flag == MPI_IDENT) {
4155: PetscSF sfPoint, sfPointSub;
4156: IS subpIS;
4157: const PetscSFNode *remotePoints;
4158: PetscSFNode *sremotePoints = NULL, *newLocalPoints = NULL, *newOwners = NULL;
4159: const PetscInt *localPoints, *subpoints, *rootdegree;
4160: PetscInt *slocalPoints = NULL, *sortedPoints = NULL, *sortedIndices = NULL;
4161: PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl = 0, ll = 0, pStart, pEnd, p;
4162: PetscMPIInt rank, size;
4164: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
4165: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
4166: PetscCall(DMGetPointSF(dm, &sfPoint));
4167: PetscCall(DMGetPointSF(subdm, &sfPointSub));
4168: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
4169: PetscCall(DMPlexGetChart(subdm, NULL, &numSubroots));
4170: PetscCall(DMPlexGetSubpointIS(subdm, &subpIS));
4171: if (subpIS) {
4172: PetscBool sorted = PETSC_TRUE;
4174: PetscCall(ISGetIndices(subpIS, &subpoints));
4175: PetscCall(ISGetLocalSize(subpIS, &numSubpoints));
4176: for (p = 1; p < numSubpoints; ++p) sorted = sorted && (subpoints[p] >= subpoints[p - 1]) ? PETSC_TRUE : PETSC_FALSE;
4177: if (!sorted) {
4178: PetscCall(PetscMalloc2(numSubpoints, &sortedPoints, numSubpoints, &sortedIndices));
4179: for (p = 0; p < numSubpoints; ++p) sortedIndices[p] = p;
4180: PetscCall(PetscArraycpy(sortedPoints, subpoints, numSubpoints));
4181: PetscCall(PetscSortIntWithArray(numSubpoints, sortedPoints, sortedIndices));
4182: }
4183: }
4184: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
4185: if (numRoots >= 0) {
4186: PetscCall(PetscSFComputeDegreeBegin(sfPoint, &rootdegree));
4187: PetscCall(PetscSFComputeDegreeEnd(sfPoint, &rootdegree));
4188: PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners));
4189: for (p = 0; p < pEnd - pStart; ++p) {
4190: newLocalPoints[p].rank = -2;
4191: newLocalPoints[p].index = -2;
4192: }
4193: for (p = pStart; p < pEnd; ++p) {
4194: newOwners[p - pStart].rank = -3;
4195: newOwners[p - pStart].index = -3;
4196: }
4197: if (sanitizeSubmesh) {
4198: /* A subpoint is forced to be owned by a rank that owns */
4199: /* a subcell that contains the subpoint in its closure. */
4200: PetscInt cStart, cEnd, c, clSize, cl;
4201: PetscInt *ownedCells, *closure = NULL;
4203: /* claim ownership */
4204: for (p = 0; p < numSubpoints; ++p) {
4205: const PetscInt point = subpoints[p];
4207: newLocalPoints[point - pStart].rank = rank;
4208: newLocalPoints[point - pStart].index = p;
4209: }
4210: PetscCall(DMGetDimension(subdm, &sdim));
4211: PetscCall(DMPlexGetDepthStratum(dm, sdim, &cStart, &cEnd));
4212: PetscCall(PetscMalloc1(cEnd - cStart, &ownedCells));
4213: for (c = cStart; c < cEnd; ++c) ownedCells[c - cStart] = 0;
4214: for (p = 0; p < numSubpoints; ++p) {
4215: c = subpoints[p];
4216: if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 1;
4217: }
4218: for (l = 0; l < numLeaves; ++l) {
4219: c = localPoints ? localPoints[l] : l;
4220: if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 0;
4221: }
4222: for (c = cStart; c < cEnd; ++c) {
4223: if (ownedCells[c - cStart] == 0) continue;
4224: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
4225: for (cl = 0; cl < clSize * 2; cl += 2) {
4226: p = closure[cl];
4227: if (newLocalPoints[p - pStart].rank < size) newLocalPoints[p - pStart].rank += size;
4228: }
4229: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
4230: }
4231: PetscCall(PetscFree(ownedCells));
4232: for (p = 0; p < numRoots; ++p) {
4233: newOwners[p].rank = newLocalPoints[p].rank;
4234: newOwners[p].index = newLocalPoints[p].index;
4235: }
4236: } else {
4237: /* Set subleaves */
4238: for (l = 0; l < numLeaves; ++l) {
4239: const PetscInt point = localPoints[l];
4240: const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);
4242: if (subpoint < 0) continue;
4243: newLocalPoints[point - pStart].rank = rank;
4244: newLocalPoints[point - pStart].index = subpoint;
4245: ++numSubleaves;
4246: }
4247: /* Must put in owned subpoints */
4248: for (p = 0; p < numSubpoints; ++p) {
4249: /* Hold on to currently owned points */
4250: if (rootdegree[subpoints[p] - pStart]) newOwners[subpoints[p] - pStart].rank = rank + size;
4251: else newOwners[subpoints[p] - pStart].rank = rank;
4252: newOwners[subpoints[p] - pStart].index = p;
4253: }
4254: }
4255: PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC));
4256: PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC));
4257: for (p = pStart; p < pEnd; ++p)
4258: if (newOwners[p - pStart].rank >= size) newOwners[p - pStart].rank -= size;
4259: if (ownershipTransferSF) {
4260: PetscSFNode *iremote1 = NULL, *newOwners1 = NULL;
4261: PetscInt *ilocal1 = NULL;
4262: PetscInt nleaves1, point;
4264: for (p = 0; p < numSubpoints; ++p) {
4265: point = subpoints[p];
4266: newLocalPoints[point - pStart].index = point;
4267: }
4268: PetscCall(PetscMalloc1(numRoots, &newOwners1));
4269: for (p = 0; p < numRoots; ++p) {
4270: newOwners1[p].rank = -1;
4271: newOwners1[p].index = -1;
4272: }
4273: PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC));
4274: PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC));
4275: for (p = 0, nleaves1 = 0; p < numRoots; ++p) {
4276: if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) ++nleaves1;
4277: }
4278: PetscCall(PetscMalloc1(nleaves1, &ilocal1));
4279: PetscCall(PetscMalloc1(nleaves1, &iremote1));
4280: for (p = 0, nleaves1 = 0; p < numRoots; ++p) {
4281: if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) {
4282: ilocal1[nleaves1] = pStart + p;
4283: iremote1[nleaves1].rank = newOwners[p].rank;
4284: iremote1[nleaves1].index = newOwners1[p].index;
4285: ++nleaves1;
4286: }
4287: }
4288: PetscCall(PetscFree(newOwners1));
4289: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), ownershipTransferSF));
4290: PetscCall(PetscSFSetFromOptions(*ownershipTransferSF));
4291: PetscCall(PetscSFSetGraph(*ownershipTransferSF, pEnd - pStart, nleaves1, ilocal1, PETSC_OWN_POINTER, iremote1, PETSC_OWN_POINTER));
4292: }
4293: if (sanitizeSubmesh) {
4294: for (p = pStart; p < pEnd; ++p) {
4295: newLocalPoints[p - pStart].rank = newOwners[p - pStart].rank;
4296: newLocalPoints[p - pStart].index = newOwners[p - pStart].index;
4297: }
4298: }
4299: PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE));
4300: PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE));
4301: if (sanitizeSubmesh) {
4302: for (p = 0; p < numSubpoints; ++p) {
4303: const PetscInt point = subpoints[p];
4305: if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) ++sl;
4306: }
4307: PetscCall(PetscMalloc1(sl, &slocalPoints));
4308: PetscCall(PetscMalloc1(sl, &sremotePoints));
4309: for (p = 0, sl = 0; p < numSubpoints; ++p) {
4310: const PetscInt point = subpoints[p];
4312: if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) {
4313: slocalPoints[sl] = p;
4314: sremotePoints[sl].rank = newLocalPoints[point - pStart].rank;
4315: sremotePoints[sl].index = newLocalPoints[point - pStart].index;
4316: PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank");
4317: PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint");
4318: ++sl;
4319: }
4320: }
4321: } else {
4322: PetscCall(PetscMalloc1(numSubleaves, &slocalPoints));
4323: PetscCall(PetscMalloc1(numSubleaves, &sremotePoints));
4324: for (l = 0; l < numLeaves; ++l) {
4325: const PetscInt point = localPoints[l];
4326: const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices);
4328: if (subpoint < 0) continue;
4329: if (newLocalPoints[point].rank == rank) {
4330: ++ll;
4331: continue;
4332: }
4333: slocalPoints[sl] = subpoint;
4334: sremotePoints[sl].rank = newLocalPoints[point].rank;
4335: sremotePoints[sl].index = newLocalPoints[point].index;
4336: PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point);
4337: PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point);
4338: ++sl;
4339: }
4340: PetscCheck(sl + ll == numSubleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubleaves);
4341: }
4342: PetscCall(PetscFree2(newLocalPoints, newOwners));
4343: PetscCall(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER));
4344: }
4345: if (subpIS) PetscCall(ISRestoreIndices(subpIS, &subpoints));
4346: PetscCall(PetscFree2(sortedPoints, sortedIndices));
4347: }
4348: /* Filter labels */
4349: PetscCall(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm));
4350: /* Cleanup */
4351: for (d = 0; d <= sdim; ++d) {
4352: if (subpointIS[d]) PetscCall(ISRestoreIndices(subpointIS[d], &subpoints[d]));
4353: PetscCall(ISDestroy(&subpointIS[d]));
4354: }
4355: PetscCall(PetscFree4(numSubPoints, firstSubPoint, subpointIS, subpoints));
4356: PetscFunctionReturn(PETSC_SUCCESS);
4357: }
4359: static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm)
4360: {
4361: PetscFunctionBegin;
4362: PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm));
4363: PetscFunctionReturn(PETSC_SUCCESS);
4364: }
4366: /*@
4367: DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label
4369: Input Parameters:
4370: + dm - The original mesh
4371: . vertexLabel - The `DMLabel` marking points contained in the surface
4372: . value - The label value to use
4373: - markedFaces - `PETSC_TRUE` if surface faces are marked in addition to vertices, `PETSC_FALSE` if only vertices are marked
4375: Output Parameter:
4376: . subdm - The surface mesh
4378: Level: developer
4380: Note:
4381: This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`.
4383: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`
4384: @*/
4385: PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm)
4386: {
4387: DMPlexInterpolatedFlag interpolated;
4388: PetscInt dim, cdim;
4390: PetscFunctionBegin;
4392: PetscAssertPointer(subdm, 5);
4393: PetscCall(DMGetDimension(dm, &dim));
4394: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm));
4395: PetscCall(DMSetType(*subdm, DMPLEX));
4396: PetscCall(DMSetDimension(*subdm, dim - 1));
4397: PetscCall(DMGetCoordinateDim(dm, &cdim));
4398: PetscCall(DMSetCoordinateDim(*subdm, cdim));
4399: PetscCall(DMPlexIsInterpolated(dm, &interpolated));
4400: PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes");
4401: if (interpolated) {
4402: PetscCall(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm));
4403: } else {
4404: PetscCall(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm));
4405: }
4406: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm));
4407: PetscFunctionReturn(PETSC_SUCCESS);
4408: }
4410: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm)
4411: {
4412: MPI_Comm comm;
4413: DMLabel subpointMap;
4414: IS subvertexIS;
4415: const PetscInt *subVertices;
4416: PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL;
4417: PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV;
4418: PetscInt c, f;
4420: PetscFunctionBegin;
4421: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
4422: /* Create subpointMap which marks the submesh */
4423: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap));
4424: PetscCall(DMPlexSetSubpointMap(subdm, subpointMap));
4425: PetscCall(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm));
4426: /* Setup chart */
4427: PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices));
4428: PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells));
4429: PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices));
4430: PetscCall(DMPlexSetVTKCellHeight(subdm, 1));
4431: /* Set cone sizes */
4432: firstSubVertex = numSubCells;
4433: firstSubFace = numSubCells + numSubVertices;
4434: newFacePoint = firstSubFace;
4435: PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS));
4436: if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices));
4437: for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1));
4438: for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV));
4439: PetscCall(DMSetUp(subdm));
4440: PetscCall(DMLabelDestroy(&subpointMap));
4441: /* Create face cones */
4442: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
4443: PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface));
4444: for (c = 0; c < numSubCells; ++c) {
4445: const PetscInt cell = subCells[c];
4446: const PetscInt subcell = c;
4447: const PetscInt *cone, *cells;
4448: PetscBool isHybrid = PETSC_FALSE;
4449: PetscInt numCells, subVertex, p, v;
4451: PetscCall(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid));
4452: if (!isHybrid) continue;
4453: PetscCall(DMPlexGetCone(dm, cell, &cone));
4454: for (v = 0; v < nFV; ++v) {
4455: PetscCall(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex));
4456: subface[v] = firstSubVertex + subVertex;
4457: }
4458: PetscCall(DMPlexSetCone(subdm, newFacePoint, subface));
4459: PetscCall(DMPlexSetCone(subdm, subcell, &newFacePoint));
4460: PetscCall(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells));
4461: /* Not true in parallel
4462: PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */
4463: for (p = 0; p < numCells; ++p) {
4464: PetscInt negsubcell;
4465: PetscBool isHybrid = PETSC_FALSE;
4467: PetscCall(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid));
4468: if (isHybrid) continue;
4469: /* I know this is a crap search */
4470: for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) {
4471: if (subCells[negsubcell] == cells[p]) break;
4472: }
4473: PetscCheck(negsubcell != numSubCells, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %" PetscInt_FMT, cell);
4474: PetscCall(DMPlexSetCone(subdm, negsubcell, &newFacePoint));
4475: }
4476: PetscCall(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells));
4477: ++newFacePoint;
4478: }
4479: PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface));
4480: PetscCall(DMPlexSymmetrize(subdm));
4481: PetscCall(DMPlexStratify(subdm));
4482: /* Build coordinates */
4483: {
4484: PetscSection coordSection, subCoordSection;
4485: Vec coordinates, subCoordinates;
4486: PetscScalar *coords, *subCoords;
4487: PetscInt cdim, numComp, coordSize, v;
4488: const char *name;
4490: PetscCall(DMGetCoordinateDim(dm, &cdim));
4491: PetscCall(DMGetCoordinateSection(dm, &coordSection));
4492: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
4493: PetscCall(DMGetCoordinateSection(subdm, &subCoordSection));
4494: PetscCall(PetscSectionSetNumFields(subCoordSection, 1));
4495: PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp));
4496: PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp));
4497: PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices));
4498: for (v = 0; v < numSubVertices; ++v) {
4499: const PetscInt vertex = subVertices[v];
4500: const PetscInt subvertex = firstSubVertex + v;
4501: PetscInt dof;
4503: PetscCall(PetscSectionGetDof(coordSection, vertex, &dof));
4504: PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof));
4505: PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof));
4506: }
4507: PetscCall(PetscSectionSetUp(subCoordSection));
4508: PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize));
4509: PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates));
4510: PetscCall(PetscObjectGetName((PetscObject)coordinates, &name));
4511: PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name));
4512: PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE));
4513: PetscCall(VecSetBlockSize(subCoordinates, cdim));
4514: PetscCall(VecSetType(subCoordinates, VECSTANDARD));
4515: PetscCall(VecGetArray(coordinates, &coords));
4516: PetscCall(VecGetArray(subCoordinates, &subCoords));
4517: for (v = 0; v < numSubVertices; ++v) {
4518: const PetscInt vertex = subVertices[v];
4519: const PetscInt subvertex = firstSubVertex + v;
4520: PetscInt dof, off, sdof, soff, d;
4522: PetscCall(PetscSectionGetDof(coordSection, vertex, &dof));
4523: PetscCall(PetscSectionGetOffset(coordSection, vertex, &off));
4524: PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof));
4525: PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff));
4526: 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);
4527: for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d];
4528: }
4529: PetscCall(VecRestoreArray(coordinates, &coords));
4530: PetscCall(VecRestoreArray(subCoordinates, &subCoords));
4531: PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates));
4532: PetscCall(VecDestroy(&subCoordinates));
4533: }
4534: /* Build SF */
4535: CHKMEMQ;
4536: {
4537: PetscSF sfPoint, sfPointSub;
4538: const PetscSFNode *remotePoints;
4539: PetscSFNode *sremotePoints, *newLocalPoints, *newOwners;
4540: const PetscInt *localPoints;
4541: PetscInt *slocalPoints;
4542: PetscInt numRoots, numLeaves, numSubRoots = numSubCells + numSubFaces + numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd;
4543: PetscMPIInt rank;
4545: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
4546: PetscCall(DMGetPointSF(dm, &sfPoint));
4547: PetscCall(DMGetPointSF(subdm, &sfPointSub));
4548: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
4549: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
4550: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
4551: if (numRoots >= 0) {
4552: /* Only vertices should be shared */
4553: PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners));
4554: for (p = 0; p < pEnd - pStart; ++p) {
4555: newLocalPoints[p].rank = -2;
4556: newLocalPoints[p].index = -2;
4557: }
4558: /* Set subleaves */
4559: for (l = 0; l < numLeaves; ++l) {
4560: const PetscInt point = localPoints[l];
4561: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
4563: PetscCheck(!(point < vStart) || !(point >= vEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %" PetscInt_FMT, point);
4564: if (subPoint < 0) continue;
4565: newLocalPoints[point - pStart].rank = rank;
4566: newLocalPoints[point - pStart].index = subPoint;
4567: ++numSubLeaves;
4568: }
4569: /* Must put in owned subpoints */
4570: for (p = pStart; p < pEnd; ++p) {
4571: const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices);
4573: if (subPoint < 0) {
4574: newOwners[p - pStart].rank = -3;
4575: newOwners[p - pStart].index = -3;
4576: } else {
4577: newOwners[p - pStart].rank = rank;
4578: newOwners[p - pStart].index = subPoint;
4579: }
4580: }
4581: PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC));
4582: PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC));
4583: PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE));
4584: PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE));
4585: PetscCall(PetscMalloc1(numSubLeaves, &slocalPoints));
4586: PetscCall(PetscMalloc1(numSubLeaves, &sremotePoints));
4587: for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) {
4588: const PetscInt point = localPoints[l];
4589: const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices);
4591: if (subPoint < 0) continue;
4592: if (newLocalPoints[point].rank == rank) {
4593: ++ll;
4594: continue;
4595: }
4596: slocalPoints[sl] = subPoint;
4597: sremotePoints[sl].rank = newLocalPoints[point].rank;
4598: sremotePoints[sl].index = newLocalPoints[point].index;
4599: PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point);
4600: PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point);
4601: ++sl;
4602: }
4603: PetscCall(PetscFree2(newLocalPoints, newOwners));
4604: PetscCheck(sl + ll == numSubLeaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubLeaves);
4605: PetscCall(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER));
4606: }
4607: }
4608: CHKMEMQ;
4609: /* Cleanup */
4610: if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices));
4611: PetscCall(ISDestroy(&subvertexIS));
4612: PetscCall(PetscFree(subCells));
4613: PetscFunctionReturn(PETSC_SUCCESS);
4614: }
4616: static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm)
4617: {
4618: DMLabel label = NULL;
4620: PetscFunctionBegin;
4621: if (labelname) PetscCall(DMGetLabel(dm, labelname, &label));
4622: PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm));
4623: PetscFunctionReturn(PETSC_SUCCESS);
4624: }
4626: /*@
4627: 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.
4629: Input Parameters:
4630: + dm - The original mesh
4631: . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells
4632: . label - A label name, or `NULL`
4633: - value - A label value
4635: Output Parameter:
4636: . subdm - The surface mesh
4638: Level: developer
4640: Note:
4641: This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`.
4643: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()`
4644: @*/
4645: PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm)
4646: {
4647: PetscInt dim, cdim, depth;
4649: PetscFunctionBegin;
4651: PetscAssertPointer(subdm, 5);
4652: PetscCall(DMGetDimension(dm, &dim));
4653: PetscCall(DMPlexGetDepth(dm, &depth));
4654: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm));
4655: PetscCall(DMSetType(*subdm, DMPLEX));
4656: PetscCall(DMSetDimension(*subdm, dim - 1));
4657: PetscCall(DMGetCoordinateDim(dm, &cdim));
4658: PetscCall(DMSetCoordinateDim(*subdm, cdim));
4659: if (depth == dim) {
4660: PetscCall(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm));
4661: } else {
4662: PetscCall(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm));
4663: }
4664: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm));
4665: PetscFunctionReturn(PETSC_SUCCESS);
4666: }
4668: /*@
4669: DMPlexReorderCohesiveSupports - Ensure that face supports for cohesive end caps are ordered
4671: Not Collective
4673: Input Parameter:
4674: . dm - The `DM` containing cohesive cells
4676: Level: developer
4678: Note:
4679: For the negative size (first) face, the cohesive cell should be first in the support, and for
4680: the positive side (second) face, the cohesive cell should be second in the support.
4682: .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexCreateCohesiveSubmesh()`
4683: @*/
4684: PetscErrorCode DMPlexReorderCohesiveSupports(DM dm)
4685: {
4686: PetscInt dim, cStart, cEnd;
4688: PetscFunctionBegin;
4689: PetscCall(DMGetDimension(dm, &dim));
4690: PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cStart, &cEnd));
4691: for (PetscInt c = cStart; c < cEnd; ++c) {
4692: const PetscInt *cone;
4693: PetscInt coneSize;
4695: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
4696: PetscCall(DMPlexGetCone(dm, c, &cone));
4697: PetscCheck(coneSize >= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hybrid cell %" PetscInt_FMT " cone size %" PetscInt_FMT " must be >= 2", c, coneSize);
4698: for (PetscInt s = 0; s < 2; ++s) {
4699: const PetscInt *supp;
4700: PetscInt suppSize, newsupp[2];
4702: PetscCall(DMPlexGetSupportSize(dm, cone[s], &suppSize));
4703: PetscCall(DMPlexGetSupport(dm, cone[s], &supp));
4704: if (suppSize == 2) {
4705: /* Reorder hybrid end cap support */
4706: if (supp[s] == c) {
4707: newsupp[0] = supp[1];
4708: newsupp[1] = supp[0];
4709: } else {
4710: newsupp[0] = supp[0];
4711: newsupp[1] = supp[1];
4712: }
4713: 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);
4714: PetscCall(DMPlexSetSupport(dm, cone[s], newsupp));
4715: }
4716: }
4717: }
4718: PetscFunctionReturn(PETSC_SUCCESS);
4719: }
4721: /*@
4722: DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh
4724: Input Parameters:
4725: + dm - The original mesh
4726: . cellLabel - The `DMLabel` marking cells contained in the new mesh
4727: . value - The label value to use
4728: . ignoreLabelHalo - The flag indicating if labeled points that are in the halo are ignored
4729: . 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
4730: - comm - The communicator you want the mesh on (currently supports only a sequential communicator or the same communicator of the original mesh)
4732: Output Parameters:
4733: + ownershipTransferSF - The `PetscSF` representing the ownership transfers between parent local meshes due to submeshing.
4734: - subdm - The new mesh
4736: Level: developer
4738: Note:
4739: This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`.
4741: 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
4742: the remote locations for these leaves are tuples of rank and point that represent the new owners.
4744: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`, `DMPlexCreateSubmesh()`
4745: @*/
4746: PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, MPI_Comm comm, PetscSF *ownershipTransferSF, DM *subdm)
4747: {
4748: PetscInt dim, overlap;
4750: PetscFunctionBegin;
4752: if (ownershipTransferSF) PetscAssertPointer(ownershipTransferSF, 7);
4753: PetscAssertPointer(subdm, 8);
4754: PetscCall(DMGetDimension(dm, &dim));
4755: PetscCall(DMCreate(comm, subdm));
4756: PetscCall(DMSetType(*subdm, DMPLEX));
4757: /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */
4758: PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, ignoreLabelHalo, sanitizeSubmesh, ownershipTransferSF, *subdm));
4759: PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm));
4760: // It is possible to obtain a surface mesh where some faces are in SF
4761: // We should either mark the mesh as having an overlap, or delete these from the SF
4762: PetscCall(DMPlexGetOverlap(dm, &overlap));
4763: if (!overlap) {
4764: PetscSF sf;
4765: const PetscInt *leaves;
4766: PetscInt cStart, cEnd, Nl;
4767: PetscBool hasSubcell = PETSC_FALSE, ghasSubcell;
4769: PetscCall(DMPlexGetHeightStratum(*subdm, 0, &cStart, &cEnd));
4770: PetscCall(DMGetPointSF(*subdm, &sf));
4771: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
4772: for (PetscInt l = 0; l < Nl; ++l) {
4773: const PetscInt point = leaves ? leaves[l] : l;
4775: if (point >= cStart && point < cEnd) {
4776: hasSubcell = PETSC_TRUE;
4777: break;
4778: }
4779: }
4780: PetscCallMPI(MPIU_Allreduce(&hasSubcell, &ghasSubcell, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
4781: if (ghasSubcell) PetscCall(DMPlexSetOverlap(*subdm, NULL, 1));
4782: }
4783: PetscFunctionReturn(PETSC_SUCCESS);
4784: }
4786: /*@
4787: DMPlexGetSubpointMap - Returns a `DMLabel` with point dimension as values
4789: Input Parameter:
4790: . dm - The submesh `DM`
4792: Output Parameter:
4793: . subpointMap - The `DMLabel` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh
4795: Level: developer
4797: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
4798: @*/
4799: PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap)
4800: {
4801: PetscFunctionBegin;
4803: PetscAssertPointer(subpointMap, 2);
4804: *subpointMap = ((DM_Plex *)dm->data)->subpointMap;
4805: PetscFunctionReturn(PETSC_SUCCESS);
4806: }
4808: /*@
4809: DMPlexSetSubpointMap - Sets the `DMLabel` with point dimension as values
4811: Input Parameters:
4812: + dm - The submesh `DM`
4813: - subpointMap - The `DMLabel` of all the points from the original mesh in this submesh
4815: Level: developer
4817: Note:
4818: Should normally not be called by the user, since it is set in `DMPlexCreateSubmesh()`
4820: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()`
4821: @*/
4822: PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap)
4823: {
4824: DM_Plex *mesh = (DM_Plex *)dm->data;
4825: DMLabel tmp;
4827: PetscFunctionBegin;
4829: tmp = mesh->subpointMap;
4830: mesh->subpointMap = subpointMap;
4831: PetscCall(PetscObjectReference((PetscObject)mesh->subpointMap));
4832: PetscCall(DMLabelDestroy(&tmp));
4833: PetscFunctionReturn(PETSC_SUCCESS);
4834: }
4836: static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS)
4837: {
4838: DMLabel spmap;
4839: PetscInt depth, d;
4841: PetscFunctionBegin;
4842: PetscCall(DMPlexGetSubpointMap(dm, &spmap));
4843: PetscCall(DMPlexGetDepth(dm, &depth));
4844: if (spmap && depth >= 0) {
4845: DM_Plex *mesh = (DM_Plex *)dm->data;
4846: PetscInt *points, *depths;
4847: PetscInt pStart, pEnd, p, off;
4849: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
4850: PetscCheck(!pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %" PetscInt_FMT, pStart);
4851: PetscCall(PetscMalloc1(pEnd, &points));
4852: PetscCall(DMGetWorkArray(dm, depth + 1, MPIU_INT, &depths));
4853: depths[0] = depth;
4854: depths[1] = 0;
4855: for (d = 2; d <= depth; ++d) depths[d] = depth + 1 - d;
4856: for (d = 0, off = 0; d <= depth; ++d) {
4857: const PetscInt dep = depths[d];
4858: PetscInt depStart, depEnd, n;
4860: PetscCall(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd));
4861: PetscCall(DMLabelGetStratumSize(spmap, dep, &n));
4862: if ((d < 2 && depth > 1) || d == 1) { /* Only check vertices and cells for now since the map is broken for others */
4863: 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);
4864: } else {
4865: if (!n) {
4866: if (d == 0) {
4867: /* Missing cells */
4868: for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = -1;
4869: } else {
4870: /* Missing faces */
4871: for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = PETSC_INT_MAX;
4872: }
4873: }
4874: }
4875: if (n) {
4876: IS is;
4877: const PetscInt *opoints;
4879: PetscCall(DMLabelGetStratumIS(spmap, dep, &is));
4880: PetscCall(ISGetIndices(is, &opoints));
4881: for (p = 0; p < n; ++p, ++off) points[off] = opoints[p];
4882: PetscCall(ISRestoreIndices(is, &opoints));
4883: PetscCall(ISDestroy(&is));
4884: }
4885: }
4886: PetscCall(DMRestoreWorkArray(dm, depth + 1, MPIU_INT, &depths));
4887: PetscCheck(off == pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " should be %" PetscInt_FMT, off, pEnd);
4888: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS));
4889: PetscCall(PetscObjectStateGet((PetscObject)spmap, &mesh->subpointState));
4890: }
4891: PetscFunctionReturn(PETSC_SUCCESS);
4892: }
4894: /*@
4895: DMPlexGetSubpointIS - Returns an `IS` covering the entire subdm chart with the original points as data
4897: Input Parameter:
4898: . dm - The submesh `DM`
4900: Output Parameter:
4901: . subpointIS - The `IS` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh
4903: Level: developer
4905: Note:
4906: 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.
4908: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()`
4909: @*/
4910: PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS)
4911: {
4912: DM_Plex *mesh = (DM_Plex *)dm->data;
4913: DMLabel spmap;
4914: PetscObjectState state;
4916: PetscFunctionBegin;
4918: PetscAssertPointer(subpointIS, 2);
4919: PetscCall(DMPlexGetSubpointMap(dm, &spmap));
4920: PetscCall(PetscObjectStateGet((PetscObject)spmap, &state));
4921: if (state != mesh->subpointState || !mesh->subpointIS) PetscCall(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS));
4922: *subpointIS = mesh->subpointIS;
4923: PetscFunctionReturn(PETSC_SUCCESS);
4924: }
4926: /*@
4927: DMGetEnclosureRelation - Get the relationship between `dmA` and `dmB`
4929: Input Parameters:
4930: + dmA - The first `DM`
4931: - dmB - The second `DM`
4933: Output Parameter:
4934: . rel - The relation of `dmA` to `dmB`
4936: Level: intermediate
4938: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosurePoint()`
4939: @*/
4940: PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel)
4941: {
4942: DM plexA, plexB, sdm;
4943: DMLabel spmap;
4944: PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB;
4946: PetscFunctionBegin;
4947: PetscAssertPointer(rel, 3);
4948: *rel = DM_ENC_NONE;
4949: if (!dmA || !dmB) PetscFunctionReturn(PETSC_SUCCESS);
4952: if (dmA == dmB) {
4953: *rel = DM_ENC_EQUALITY;
4954: PetscFunctionReturn(PETSC_SUCCESS);
4955: }
4956: PetscCall(DMConvert(dmA, DMPLEX, &plexA));
4957: PetscCall(DMConvert(dmB, DMPLEX, &plexB));
4958: PetscCall(DMPlexGetChart(plexA, &pStartA, &pEndA));
4959: PetscCall(DMPlexGetChart(plexB, &pStartB, &pEndB));
4960: /* Assumption 1: subDMs have smaller charts than the DMs that they originate from
4961: - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */
4962: if (pStartA == pStartB && pEndA == pEndB) {
4963: *rel = DM_ENC_EQUALITY;
4964: goto end;
4965: }
4966: NpA = pEndA - pStartA;
4967: NpB = pEndB - pStartB;
4968: if (NpA == NpB) goto end;
4969: sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */
4970: PetscCall(DMPlexGetSubpointMap(sdm, &spmap));
4971: if (!spmap) goto end;
4972: /* TODO Check the space mapped to by subpointMap is same size as dm */
4973: if (NpA > NpB) {
4974: *rel = DM_ENC_SUPERMESH;
4975: } else {
4976: *rel = DM_ENC_SUBMESH;
4977: }
4978: end:
4979: PetscCall(DMDestroy(&plexA));
4980: PetscCall(DMDestroy(&plexB));
4981: PetscFunctionReturn(PETSC_SUCCESS);
4982: }
4984: /*@
4985: DMGetEnclosurePoint - Get the point `pA` in `dmA` which corresponds to the point `pB` in `dmB`
4987: Input Parameters:
4988: + dmA - The first `DM`
4989: . dmB - The second `DM`
4990: . etype - The type of enclosure relation that `dmA` has to `dmB`
4991: - pB - A point of `dmB`
4993: Output Parameter:
4994: . pA - The corresponding point of `dmA`
4996: Level: intermediate
4998: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosureRelation()`
4999: @*/
5000: PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA)
5001: {
5002: DM sdm;
5003: IS subpointIS;
5004: const PetscInt *subpoints;
5005: PetscInt numSubpoints;
5007: PetscFunctionBegin;
5008: /* TODO Cache the IS, making it look like an index */
5009: switch (etype) {
5010: case DM_ENC_SUPERMESH:
5011: sdm = dmB;
5012: PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS));
5013: PetscCall(ISGetIndices(subpointIS, &subpoints));
5014: *pA = subpoints[pB];
5015: PetscCall(ISRestoreIndices(subpointIS, &subpoints));
5016: break;
5017: case DM_ENC_SUBMESH:
5018: sdm = dmA;
5019: PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS));
5020: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
5021: PetscCall(ISGetIndices(subpointIS, &subpoints));
5022: PetscCall(PetscFindInt(pB, numSubpoints, subpoints, pA));
5023: if (*pA < 0) {
5024: PetscCall(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view"));
5025: PetscCall(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view"));
5026: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB);
5027: }
5028: PetscCall(ISRestoreIndices(subpointIS, &subpoints));
5029: break;
5030: case DM_ENC_EQUALITY:
5031: case DM_ENC_NONE:
5032: *pA = pB;
5033: break;
5034: case DM_ENC_UNKNOWN: {
5035: DMEnclosureType enc;
5037: PetscCall(DMGetEnclosureRelation(dmA, dmB, &enc));
5038: PetscCall(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA));
5039: } break;
5040: default:
5041: SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int)etype);
5042: }
5043: PetscFunctionReturn(PETSC_SUCCESS);
5044: }