Actual source code: plexorient.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
3: #include <petsc/private/hashsetij.h>
5: /*@
6: DMPlexOrientPoint - Act with the given orientation on the cone points of this mesh point, and update its use in the mesh.
8: Not Collective
10: Input Parameters:
11: + dm - The `DM`
12: . p - The mesh point
13: - o - The orientation
15: Level: intermediate
17: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`, `DMPlexGetCone()`, `DMPlexGetConeOrientation()`, `DMPlexInterpolate()`, `DMPlexGetChart()`
18: @*/
19: PetscErrorCode DMPlexOrientPoint(DM dm, PetscInt p, PetscInt o)
20: {
21: DMPolytopeType ct;
22: const PetscInt *arr, *cone, *ornt, *support;
23: PetscInt *newcone, *newornt;
24: PetscInt coneSize, c, supportSize, s;
26: PetscFunctionBegin;
28: PetscCall(DMPlexGetCellType(dm, p, &ct));
29: arr = DMPolytopeTypeGetArrangement(ct, o);
30: if (!arr) PetscFunctionReturn(PETSC_SUCCESS);
31: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
32: PetscCall(DMPlexGetCone(dm, p, &cone));
33: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
34: PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newcone));
35: PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newornt));
36: for (c = 0; c < coneSize; ++c) {
37: DMPolytopeType ft;
38: PetscInt nO;
40: PetscCall(DMPlexGetCellType(dm, cone[c], &ft));
41: nO = DMPolytopeTypeGetNumArrangements(ft) / 2;
42: newcone[c] = cone[arr[c * 2 + 0]];
43: newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], ornt[arr[c * 2 + 0]]);
44: PetscCheck(!newornt[c] || !(newornt[c] >= nO || newornt[c] < -nO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid orientation %" PetscInt_FMT " not in [%" PetscInt_FMT ",%" PetscInt_FMT ") for %s %" PetscInt_FMT, newornt[c], -nO, nO, DMPolytopeTypes[ft], cone[c]);
45: }
46: PetscCall(DMPlexSetCone(dm, p, newcone));
47: PetscCall(DMPlexSetConeOrientation(dm, p, newornt));
48: PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newcone));
49: PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newornt));
50: /* Update orientation of this point in the support points */
51: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
52: PetscCall(DMPlexGetSupport(dm, p, &support));
53: for (s = 0; s < supportSize; ++s) {
54: PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize));
55: PetscCall(DMPlexGetCone(dm, support[s], &cone));
56: PetscCall(DMPlexGetConeOrientation(dm, support[s], &ornt));
57: for (c = 0; c < coneSize; ++c) {
58: PetscInt po;
60: if (cone[c] != p) continue;
61: /* ornt[c] * 0 = target = po * o so that po = ornt[c] * o^{-1} */
62: po = DMPolytopeTypeComposeOrientationInv(ct, ornt[c], o);
63: PetscCall(DMPlexInsertConeOrientation(dm, support[s], c, po));
64: }
65: }
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscInt GetPointIndex(PetscInt point, PetscInt pStart, PetscInt pEnd, const PetscInt points[])
70: {
71: if (points) {
72: PetscInt loc;
74: PetscCallAbort(PETSC_COMM_SELF, PetscFindInt(point, pEnd - pStart, points, &loc));
75: if (loc >= 0) return loc;
76: } else {
77: if (point >= pStart && point < pEnd) return point - pStart;
78: }
79: return -1;
80: }
82: /*
83: - Checks face match
84: - Flips non-matching
85: - Inserts faces of support cells in FIFO
86: */
87: static PetscErrorCode DMPlexCheckFace_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, IS cellIS, IS faceIS, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces)
88: {
89: const PetscInt *supp, *coneA, *coneB, *coneOA, *coneOB;
90: PetscInt suppSize, Ns = 0, coneSizeA, coneSizeB, posA = -1, posB = -1;
91: PetscInt face, dim, indC[3], indS[3], seenA, flippedA, seenB, flippedB, mismatch;
92: const PetscInt *cells, *faces;
93: PetscInt cStart, cEnd, fStart, fEnd;
95: PetscFunctionBegin;
96: face = faceFIFO[(*fTop)++];
97: PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
98: PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
99: PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &dim));
100: PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
101: PetscCall(DMPlexGetSupport(dm, face, &supp));
102: // Filter the support
103: for (PetscInt s = 0; s < suppSize; ++s) {
104: // Filter support
105: indC[Ns] = GetPointIndex(supp[s], cStart, cEnd, cells);
106: indS[Ns] = s;
107: if (indC[Ns] >= 0) ++Ns;
108: }
109: if (Ns < 2) PetscFunctionReturn(PETSC_SUCCESS);
110: PetscCheck(Ns == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, Ns);
111: PetscCheck(indC[0] >= 0 && indC[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Support cells %" PetscInt_FMT " (%" PetscInt_FMT ") and %" PetscInt_FMT " (%" PetscInt_FMT ") are not both valid", supp[0], indC[0], supp[1], indC[1]);
112: seenA = PetscBTLookup(seenCells, indC[0]);
113: flippedA = PetscBTLookup(flippedCells, indC[0]) ? 1 : 0;
114: seenB = PetscBTLookup(seenCells, indC[1]);
115: flippedB = PetscBTLookup(flippedCells, indC[1]) ? 1 : 0;
117: PetscCall(DMPlexGetConeSize(dm, supp[indS[0]], &coneSizeA));
118: PetscCall(DMPlexGetConeSize(dm, supp[indS[1]], &coneSizeB));
119: PetscCall(DMPlexGetCone(dm, supp[indS[0]], &coneA));
120: PetscCall(DMPlexGetCone(dm, supp[indS[1]], &coneB));
121: PetscCall(DMPlexGetConeOrientation(dm, supp[indS[0]], &coneOA));
122: PetscCall(DMPlexGetConeOrientation(dm, supp[indS[1]], &coneOB));
123: for (PetscInt c = 0; c < coneSizeA; ++c) {
124: const PetscInt indF = GetPointIndex(coneA[c], fStart, fEnd, faces);
126: // Filter cone
127: if (indF < 0) continue;
128: if (!PetscBTLookup(seenFaces, indF)) {
129: faceFIFO[(*fBottom)++] = coneA[c];
130: PetscCall(PetscBTSet(seenFaces, indF));
131: }
132: if (coneA[c] == face) posA = c;
133: PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart);
134: }
135: PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, supp[indS[0]]);
136: for (PetscInt c = 0; c < coneSizeB; ++c) {
137: const PetscInt indF = GetPointIndex(coneB[c], fStart, fEnd, faces);
139: // Filter cone
140: if (indF < 0) continue;
141: if (!PetscBTLookup(seenFaces, indF)) {
142: faceFIFO[(*fBottom)++] = coneB[c];
143: PetscCall(PetscBTSet(seenFaces, indF));
144: }
145: if (coneB[c] == face) posB = c;
146: PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart);
147: }
148: PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, supp[indS[1]]);
150: if (dim == 1) {
151: mismatch = posA == posB;
152: } else {
153: mismatch = coneOA[posA] == coneOB[posB];
154: }
156: if (mismatch ^ (flippedA ^ flippedB)) {
157: PetscCheck(!seenA || !seenB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", supp[indS[0]], supp[indS[1]]);
158: if (!seenA && !flippedA) PetscCall(PetscBTSet(flippedCells, indC[0]));
159: else {
160: PetscCheck(!seenB && !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
161: PetscCall(PetscBTSet(flippedCells, indC[1]));
162: }
163: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
164: PetscCall(PetscBTSet(seenCells, indC[0]));
165: PetscCall(PetscBTSet(seenCells, indC[1]));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode DMPlexCheckFace_Old_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, PetscInt cStart, PetscInt fStart, PetscInt fEnd, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces)
170: {
171: const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
172: PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
173: PetscInt face, dim, seenA, flippedA, seenB, flippedB, mismatch, c;
175: PetscFunctionBegin;
176: face = faceFIFO[(*fTop)++];
177: PetscCall(DMGetDimension(dm, &dim));
178: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
179: PetscCall(DMPlexGetSupport(dm, face, &support));
180: if (supportSize < 2) PetscFunctionReturn(PETSC_SUCCESS);
181: PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, supportSize);
182: seenA = PetscBTLookup(seenCells, support[0] - cStart);
183: flippedA = PetscBTLookup(flippedCells, support[0] - cStart) ? 1 : 0;
184: seenB = PetscBTLookup(seenCells, support[1] - cStart);
185: flippedB = PetscBTLookup(flippedCells, support[1] - cStart) ? 1 : 0;
187: PetscCall(DMPlexGetConeSize(dm, support[0], &coneSizeA));
188: PetscCall(DMPlexGetConeSize(dm, support[1], &coneSizeB));
189: PetscCall(DMPlexGetCone(dm, support[0], &coneA));
190: PetscCall(DMPlexGetCone(dm, support[1], &coneB));
191: PetscCall(DMPlexGetConeOrientation(dm, support[0], &coneOA));
192: PetscCall(DMPlexGetConeOrientation(dm, support[1], &coneOB));
193: for (c = 0; c < coneSizeA; ++c) {
194: if (!PetscBTLookup(seenFaces, coneA[c] - fStart)) {
195: faceFIFO[(*fBottom)++] = coneA[c];
196: PetscCall(PetscBTSet(seenFaces, coneA[c] - fStart));
197: }
198: if (coneA[c] == face) posA = c;
199: PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart);
200: }
201: PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[0]);
202: for (c = 0; c < coneSizeB; ++c) {
203: if (!PetscBTLookup(seenFaces, coneB[c] - fStart)) {
204: faceFIFO[(*fBottom)++] = coneB[c];
205: PetscCall(PetscBTSet(seenFaces, coneB[c] - fStart));
206: }
207: if (coneB[c] == face) posB = c;
208: PetscCheck(*fBottom <= fEnd - fStart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %" PetscInt_FMT " was pushed exceeding capacity %" PetscInt_FMT " > %" PetscInt_FMT, coneA[c], *fBottom, fEnd - fStart);
209: }
210: PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[1]);
212: if (dim == 1) {
213: mismatch = posA == posB;
214: } else {
215: mismatch = coneOA[posA] == coneOB[posB];
216: }
218: if (mismatch ^ (flippedA ^ flippedB)) {
219: PetscCheck(!seenA || !seenB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", support[0], support[1]);
220: if (!seenA && !flippedA) {
221: PetscCall(PetscBTSet(flippedCells, support[0] - cStart));
222: } else if (!seenB && !flippedB) {
223: PetscCall(PetscBTSet(flippedCells, support[1] - cStart));
224: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
225: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
226: PetscCall(PetscBTSet(seenCells, support[0] - cStart));
227: PetscCall(PetscBTSet(seenCells, support[1] - cStart));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*
232: DMPlexOrient_Serial - Compute valid orientation for local connected components
234: Not collective
236: Input Parameters:
237: + dm - The `DM`
238: - cellHeight - The height of k-cells to be oriented
240: Output Parameters:
241: + Ncomp - The number of connected component
242: . cellComp - The connected component for each local cell
243: . faceComp - The connected component for each local face
244: - flippedCells - Marked cells should be inverted
246: Level: developer
248: .seealso: `DMPlexOrient()`
249: */
250: static PetscErrorCode DMPlexOrient_Serial(DM dm, IS cellIS, IS faceIS, PetscInt *Ncomp, PetscInt cellComp[], PetscInt faceComp[], PetscBT flippedCells)
251: {
252: PetscBT seenCells, seenFaces;
253: PetscInt *faceFIFO;
254: const PetscInt *cells = NULL, *faces = NULL;
255: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
257: PetscFunctionBegin;
258: /* Truth Table
259: mismatch flips do action mismatch flipA ^ flipB action
260: F 0 flips no F F F
261: F 1 flip yes F T T
262: F 2 flips no T F T
263: T 0 flips yes T T F
264: T 1 flip no
265: T 2 flips yes
266: */
267: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
268: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
269: PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
270: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
271: PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
272: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
273: PetscCall(PetscMalloc1(fEnd - fStart, &faceFIFO));
274: *Ncomp = 0;
275: for (PetscInt c = 0; c < cEnd - cStart; ++c) cellComp[c] = -1;
276: do {
277: PetscInt cc, fTop, fBottom;
279: // Look for first unmarked cell
280: for (cc = cStart; cc < cEnd; ++cc)
281: if (cellComp[cc - cStart] < 0) break;
282: if (cc >= cEnd) break;
283: // Initialize FIFO with first cell in component
284: {
285: const PetscInt cell = cells ? cells[cc] : cc;
286: const PetscInt *cone;
287: PetscInt coneSize;
289: fTop = fBottom = 0;
290: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
291: PetscCall(DMPlexGetCone(dm, cell, &cone));
292: for (PetscInt c = 0; c < coneSize; ++c) {
293: const PetscInt idx = GetPointIndex(cone[c], fStart, fEnd, faces);
295: // Cell faces are guaranteed to be in the face set
296: PetscCheck(idx >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " of cell %" PetscInt_FMT " is not present in the label", cone[c], cell);
297: faceFIFO[fBottom++] = cone[c];
298: PetscCall(PetscBTSet(seenFaces, idx));
299: }
300: PetscCall(PetscBTSet(seenCells, cc - cStart));
301: }
302: // Consider each face in FIFO
303: while (fTop < fBottom) PetscCall(DMPlexCheckFace_Internal(dm, faceFIFO, &fTop, &fBottom, cellIS, faceIS, seenCells, flippedCells, seenFaces));
304: // Set component for cells and faces
305: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
306: if (PetscBTLookup(seenCells, c)) cellComp[c] = *Ncomp;
307: }
308: for (PetscInt f = 0; f < fEnd - fStart; ++f) {
309: if (PetscBTLookup(seenFaces, f)) faceComp[f] = *Ncomp;
310: }
311: // Wipe seenCells and seenFaces for next component
312: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
313: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
314: ++(*Ncomp);
315: } while (1);
316: PetscCall(PetscBTDestroy(&seenCells));
317: PetscCall(PetscBTDestroy(&seenFaces));
318: PetscCall(PetscFree(faceFIFO));
319: for (PetscInt c = 0; c < cEnd - cStart; ++c)
320: PetscCheck(0 <= cellComp[c] && cellComp[c] < *Ncomp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT " for cell %" PetscInt_FMT " (%" PetscInt_FMT ")", cellComp[c], cells ? cells[c] : c, c);
321: for (PetscInt f = 0; f < fEnd - fStart; ++f)
322: PetscCheck(0 <= faceComp[f] && faceComp[f] < *Ncomp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT " for face %" PetscInt_FMT " (%" PetscInt_FMT ")", faceComp[f], faces ? faces[f] : f, f);
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /*@
327: DMPlexOrient - Give a consistent orientation to the input mesh
329: Input Parameter:
330: . dm - The `DM`
332: Notes:
333: The orientation data for the `DM` are changed in-place.
335: This routine will fail for non-orientable surfaces, such as the Moebius strip.
337: Level: advanced
339: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexOrientLabel()`
340: @*/
341: PetscErrorCode DMPlexOrient(DM dm)
342: {
343: #if 0
344: IS cellIS, faceIS;
346: PetscFunctionBegin;
347: PetscCall(DMPlexGetAllCells_Internal(dm, &cellIS));
348: PetscCall(DMPlexGetAllFaces_Internal(dm, &faceIS));
349: PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
350: PetscCall(ISDestroy(&cellIS));
351: PetscCall(ISDestroy(&faceIS));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: #else
354: MPI_Comm comm;
355: PetscSF sf;
356: const PetscInt *lpoints;
357: const PetscSFNode *rpoints;
358: PetscSFNode *rorntComp = NULL, *lorntComp = NULL;
359: PetscInt *numNeighbors, **neighbors, *locSupport = NULL;
360: PetscSFNode *nrankComp;
361: PetscBool *match, *flipped;
362: PetscBT seenCells, flippedCells, seenFaces;
363: PetscInt *faceFIFO, fTop, fBottom, *cellComp, *faceComp;
364: PetscInt numLeaves, numRoots, dim, h, cStart, cEnd, c, cell, fStart, fEnd, face, off, totNeighbors = 0;
365: PetscMPIInt rank, size, numComponents, comp = 0;
366: PetscBool flg, flg2;
367: PetscViewer viewer = NULL, selfviewer = NULL;
369: PetscFunctionBegin;
370: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
371: PetscCallMPI(MPI_Comm_rank(comm, &rank));
372: PetscCallMPI(MPI_Comm_size(comm, &size));
373: PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &flg));
374: PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &flg2));
375: PetscCall(DMGetPointSF(dm, &sf));
376: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints));
377: /* Truth Table
378: mismatch flips do action mismatch flipA ^ flipB action
379: F 0 flips no F F F
380: F 1 flip yes F T T
381: F 2 flips no T F T
382: T 0 flips yes T T F
383: T 1 flip no
384: T 2 flips yes
385: */
386: PetscCall(DMGetDimension(dm, &dim));
387: PetscCall(DMPlexGetVTKCellHeight(dm, &h));
388: PetscCall(DMPlexGetHeightStratum(dm, h, &cStart, &cEnd));
389: PetscCall(DMPlexGetHeightStratum(dm, h + 1, &fStart, &fEnd));
390: PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
391: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
392: PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells));
393: PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells));
394: PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
395: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
396: PetscCall(PetscCalloc3(fEnd - fStart, &faceFIFO, cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
397: /*
398: OLD STYLE
399: - Add an integer array over cells and faces (component) for connected component number
400: Foreach component
401: - Mark the initial cell as seen
402: - Process component as usual
403: - Set component for all seenCells
404: - Wipe seenCells and seenFaces (flippedCells can stay)
405: - Generate parallel adjacency for component using SF and seenFaces
406: - Collect numComponents adj data from each proc to 0
407: - Build same serial graph
408: - Use same solver
409: - Use Scatterv to send back flipped flags for each component
410: - Negate flippedCells by component
412: NEW STYLE
413: - Create the adj on each process
414: - Bootstrap to complete graph on proc 0
415: */
416: /* Loop over components */
417: for (cell = cStart; cell < cEnd; ++cell) cellComp[cell - cStart] = -1;
418: do {
419: /* Look for first unmarked cell */
420: for (cell = cStart; cell < cEnd; ++cell)
421: if (cellComp[cell - cStart] < 0) break;
422: if (cell >= cEnd) break;
423: /* Initialize FIFO with first cell in component */
424: {
425: const PetscInt *cone;
426: PetscInt coneSize;
428: fTop = fBottom = 0;
429: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
430: PetscCall(DMPlexGetCone(dm, cell, &cone));
431: for (c = 0; c < coneSize; ++c) {
432: faceFIFO[fBottom++] = cone[c];
433: PetscCall(PetscBTSet(seenFaces, cone[c] - fStart));
434: }
435: PetscCall(PetscBTSet(seenCells, cell - cStart));
436: }
437: /* Consider each face in FIFO */
438: while (fTop < fBottom) PetscCall(DMPlexCheckFace_Old_Internal(dm, faceFIFO, &fTop, &fBottom, cStart, fStart, fEnd, seenCells, flippedCells, seenFaces));
439: /* Set component for cells and faces */
440: for (cell = 0; cell < cEnd - cStart; ++cell) {
441: if (PetscBTLookup(seenCells, cell)) cellComp[cell] = comp;
442: }
443: for (face = 0; face < fEnd - fStart; ++face) {
444: if (PetscBTLookup(seenFaces, face)) faceComp[face] = comp;
445: }
446: /* Wipe seenCells and seenFaces for next component */
447: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
448: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
449: ++comp;
450: } while (1);
451: numComponents = comp;
452: if (flg) {
453: PetscViewer v;
455: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
456: PetscCall(PetscViewerASCIIPushSynchronized(v));
457: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
458: PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
459: PetscCall(PetscViewerFlush(v));
460: PetscCall(PetscViewerASCIIPopSynchronized(v));
461: }
462: /* Now all subdomains are oriented, but we need a consistent parallel orientation */
463: if (numLeaves >= 0) {
464: PetscInt maxSupportSize, neighbor;
466: /* Store orientations of boundary faces*/
467: PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSupportSize));
468: PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSupportSize, &locSupport));
469: for (face = fStart; face < fEnd; ++face) {
470: const PetscInt *cone, *support, *ornt;
471: PetscInt coneSize, supportSize, Ns = 0, s, l;
473: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
474: /* Ignore overlapping cells */
475: PetscCall(DMPlexGetSupport(dm, face, &support));
476: for (s = 0; s < supportSize; ++s) {
477: if (lpoints) PetscCall(PetscFindInt(support[s], numLeaves, lpoints, &l));
478: else {
479: if (support[s] >= 0 && support[s] < numLeaves) l = support[s];
480: else l = -1;
481: }
482: if (l >= 0) continue;
483: locSupport[Ns++] = support[s];
484: }
485: if (Ns != 1) continue;
486: neighbor = locSupport[0];
487: PetscCall(DMPlexGetCone(dm, neighbor, &cone));
488: PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize));
489: PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt));
490: for (c = 0; c < coneSize; ++c)
491: if (cone[c] == face) break;
492: if (dim == 1) {
493: /* Use cone position instead, shifted to -1 or 1 */
494: if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = 1 - c * 2;
495: else rorntComp[face].rank = c * 2 - 1;
496: } else {
497: if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1;
498: else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1;
499: }
500: rorntComp[face].index = faceComp[face - fStart];
501: }
502: /* Communicate boundary edge orientations */
503: PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
504: PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
505: }
506: /* Get process adjacency */
507: PetscCall(PetscMalloc2(numComponents, &numNeighbors, numComponents, &neighbors));
508: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
509: if (flg2) PetscCall(PetscViewerASCIIPushSynchronized(viewer));
510: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
511: for (comp = 0; comp < numComponents; ++comp) {
512: PetscInt n;
514: numNeighbors[comp] = 0;
515: PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp]));
516: /* I know this is p^2 time in general, but for bounded degree its alright */
517: for (PetscInt l = 0; l < numLeaves; ++l) {
518: const PetscInt face = lpoints ? lpoints[l] : l;
520: /* Find a representative face (edge) separating pairs of procs */
521: if ((face >= fStart) && (face < fEnd) && (faceComp[face - fStart] == comp) && rorntComp[face].rank) {
522: const PetscInt rrank = rpoints[l].rank;
523: const PetscInt rcomp = lorntComp[face].index;
525: for (n = 0; n < numNeighbors[comp]; ++n)
526: if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break;
527: if (n >= numNeighbors[comp]) {
528: PetscInt supportSize;
530: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
531: // We can have internal faces in the SF if we have cells in the SF
532: if (supportSize > 1) continue;
533: if (flg)
534: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %d, Found representative leaf %" PetscInt_FMT " (face %" PetscInt_FMT ") connecting to face %" PetscInt_FMT " on (%" PetscInt_FMT ", %" PetscInt_FMT ") with orientation %" PetscInt_FMT "\n", rank, comp, l, face,
535: rpoints[l].index, rrank, rcomp, lorntComp[face].rank));
536: neighbors[comp][numNeighbors[comp]++] = l;
537: }
538: }
539: }
540: totNeighbors += numNeighbors[comp];
541: }
542: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
543: if (flg2) PetscCall(PetscViewerASCIIPopSynchronized(viewer));
544: PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match));
545: for (comp = 0, off = 0; comp < numComponents; ++comp) {
546: PetscInt n;
548: for (n = 0; n < numNeighbors[comp]; ++n, ++off) {
549: const PetscInt face = lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n];
550: const PetscInt o = rorntComp[face].rank * lorntComp[face].rank;
552: if (o < 0) match[off] = PETSC_TRUE;
553: else if (o > 0) match[off] = PETSC_FALSE;
554: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid face %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ") neighbor: %" PetscInt_FMT " comp: %d", face, rorntComp[face].rank, lorntComp[face].rank, neighbors[comp][n], comp);
555: nrankComp[off].rank = rpoints[neighbors[comp][n]].rank;
556: nrankComp[off].index = lorntComp[lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n]].index;
557: }
558: PetscCall(PetscFree(neighbors[comp]));
559: }
560: /* Collect the graph on 0 */
561: if (numLeaves >= 0) {
562: Mat G;
563: PetscBT seenProcs, flippedProcs;
564: PetscInt *procFIFO, pTop, pBottom;
565: PetscInt *N = NULL, *Noff;
566: PetscSFNode *adj = NULL;
567: PetscBool *val = NULL;
568: PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc, p, o, itotNeighbors;
569: PetscMPIInt size = 0;
571: PetscCall(PetscCalloc1(numComponents, &flipped));
572: if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
573: PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
574: PetscCallMPI(MPI_Gather(&numComponents, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
575: for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
576: if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
577: PetscCallMPI(MPI_Gatherv(numNeighbors, numComponents, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
578: for (p = 0, o = 0; p < size; ++p) {
579: recvcounts[p] = 0;
580: for (c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
581: displs[p + 1] = displs[p] + recvcounts[p];
582: }
583: if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
584: PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors));
585: PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm));
586: PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPI_C_BOOL, val, recvcounts, displs, MPI_C_BOOL, 0, comm));
587: PetscCall(PetscFree2(numNeighbors, neighbors));
588: if (rank == 0) {
589: for (p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
590: if (flg) {
591: PetscInt n;
593: for (p = 0, off = 0; p < size; ++p) {
594: for (c = 0; c < Nc[p]; ++c) {
595: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %d Comp %" PetscInt_FMT ":\n", p, c));
596: for (n = 0; n < N[Noff[p] + c]; ++n, ++off) PetscCall(PetscPrintf(PETSC_COMM_SELF, " edge (%" PetscInt_FMT ", %" PetscInt_FMT ") (%s):\n", adj[off].rank, adj[off].index, PetscBools[val[off]]));
597: }
598: }
599: }
600: /* Symmetrize the graph */
601: PetscCall(MatCreate(PETSC_COMM_SELF, &G));
602: PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
603: PetscCall(MatSetUp(G));
604: for (p = 0, off = 0; p < size; ++p) {
605: for (c = 0; c < Nc[p]; ++c) {
606: const PetscInt r = Noff[p] + c;
608: for (PetscInt n = 0; n < N[r]; ++n, ++off) {
609: const PetscInt q = Noff[adj[off].rank] + adj[off].index;
610: const PetscScalar o = val[off] ? 1.0 : 0.0;
612: PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
613: PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
614: }
615: }
616: }
617: PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
618: PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
620: PetscCall(PetscBTCreate(Noff[size], &seenProcs));
621: PetscCall(PetscBTMemzero(Noff[size], seenProcs));
622: PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
623: PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
624: PetscCall(PetscMalloc1(Noff[size], &procFIFO));
625: pTop = pBottom = 0;
626: for (p = 0; p < Noff[size]; ++p) {
627: if (PetscBTLookup(seenProcs, p)) continue;
628: /* Initialize FIFO with next proc */
629: procFIFO[pBottom++] = p;
630: PetscCall(PetscBTSet(seenProcs, p));
631: /* Consider each proc in FIFO */
632: while (pTop < pBottom) {
633: const PetscScalar *ornt;
634: const PetscInt *neighbors;
635: PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors, n;
637: proc = procFIFO[pTop++];
638: flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
639: PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
640: /* Loop over neighboring procs */
641: for (n = 0; n < numNeighbors; ++n) {
642: nproc = neighbors[n];
643: mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
644: seen = PetscBTLookup(seenProcs, nproc);
645: flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
647: if (mismatch ^ (flippedA ^ flippedB)) {
648: PetscCheck(!seen, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", proc, nproc);
649: PetscCheck(!flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
650: PetscCall(PetscBTSet(flippedProcs, nproc));
651: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
652: if (!seen) {
653: procFIFO[pBottom++] = nproc;
654: PetscCall(PetscBTSet(seenProcs, nproc));
655: }
656: }
657: }
658: }
659: PetscCall(PetscFree(procFIFO));
660: PetscCall(MatDestroy(&G));
661: PetscCall(PetscFree2(adj, val));
662: PetscCall(PetscBTDestroy(&seenProcs));
663: }
664: /* Scatter flip flags */
665: {
666: PetscBool *flips = NULL;
668: if (rank == 0) {
669: PetscCall(PetscMalloc1(Noff[size], &flips));
670: for (p = 0; p < Noff[size]; ++p) {
671: flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
672: if (flg && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %d:\n", p));
673: }
674: for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
675: }
676: PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPI_C_BOOL, flipped, numComponents, MPI_C_BOOL, 0, comm));
677: PetscCall(PetscFree(flips));
678: }
679: if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
680: PetscCall(PetscFree(N));
681: PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
682: PetscCall(PetscFree2(nrankComp, match));
684: /* Decide whether to flip cells in each component */
685: for (c = 0; c < cEnd - cStart; ++c) {
686: if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c));
687: }
688: PetscCall(PetscFree(flipped));
689: }
690: if (flg) {
691: PetscViewer v;
693: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
694: PetscCall(PetscViewerASCIIPushSynchronized(v));
695: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
696: PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
697: PetscCall(PetscViewerFlush(v));
698: PetscCall(PetscViewerASCIIPopSynchronized(v));
699: }
700: /* Reverse flipped cells in the mesh */
701: for (c = cStart; c < cEnd; ++c) {
702: if (PetscBTLookup(flippedCells, c - cStart)) PetscCall(DMPlexOrientPoint(dm, c, -1));
703: }
704: PetscCall(PetscBTDestroy(&seenCells));
705: PetscCall(PetscBTDestroy(&flippedCells));
706: PetscCall(PetscBTDestroy(&seenFaces));
707: PetscCall(PetscFree2(numNeighbors, neighbors));
708: PetscCall(PetscFree3(rorntComp, lorntComp, locSupport));
709: PetscCall(PetscFree3(faceFIFO, cellComp, faceComp));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: #endif
712: }
714: static PetscErrorCode CreateCellAndFaceIS_Private(DM dm, DMLabel label, IS *cellIS, IS *faceIS)
715: {
716: IS valueIS;
717: const PetscInt *values;
718: PetscInt Nv, depth = 0;
720: PetscFunctionBegin;
721: PetscCall(DMLabelGetValueIS(label, &valueIS));
722: PetscCall(ISGetLocalSize(valueIS, &Nv));
723: PetscCall(ISGetIndices(valueIS, &values));
724: for (PetscInt v = 0; v < Nv; ++v) {
725: const PetscInt val = values[v] < 0 || values[v] >= 100 ? 0 : values[v];
726: PetscInt n;
728: PetscCall(DMLabelGetStratumSize(label, val, &n));
729: if (!n) continue;
730: depth = PetscMax(val, depth);
731: }
732: PetscCall(ISDestroy(&valueIS));
733: PetscCheck(depth >= 1 || !Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Depth for interface must be at least 1, not %" PetscInt_FMT, depth);
734: PetscCall(DMLabelGetStratumIS(label, depth, cellIS));
735: PetscCall(DMLabelGetStratumIS(label, depth - 1, faceIS));
736: if (!*cellIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cellIS));
737: if (!*faceIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, faceIS));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: /*@
742: DMPlexOrientLabel - Give a consistent orientation to the hypersurface marked by the `DMLabel` in the input mesh
744: Collective on dm
746: Input Parameters:
747: + dm - The `DM`
748: - label - The `DMLabel`
750: Notes:
751: The orientation data for the `DM` are changed in-place.
753: This routine will fail for non-orientable surfaces, such as the Moebius strip.
755: Level: advanced
757: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexOrient()`
758: @*/
759: PetscErrorCode DMPlexOrientLabel(DM dm, DMLabel label)
760: {
761: IS cellIS, faceIS;
763: PetscFunctionBegin;
764: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
765: PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
766: PetscCall(ISDestroy(&cellIS));
767: PetscCall(ISDestroy(&faceIS));
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: typedef struct {
772: PetscInt comp; // connected component number
773: PetscInt cell; // cell determining face orientation
774: PetscInt ornt; // face orientation
775: } SharedFace;
777: static PetscErrorCode DMPlexOrientCreateSharedFaces_Internal(DM dm, IS cellIS, IS faceIS, PetscInt Ncomp, const PetscInt faceComp[], PetscBT cellFlip, SharedFace **localFace, SharedFace **remoteFace)
778: {
779: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
780: const PetscInt *cells = NULL, *faces = NULL;
781: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
782: PetscSF sf;
783: const PetscInt *lpoints, *rootdegree;
784: const PetscSFNode *rpoints;
785: PetscInt Nr, Nl;
786: PetscInt depth, fdepth;
787: PetscBool faceIsVertex = PETSC_FALSE;
788: MPI_Datatype MPIU_3INT;
789: PetscViewer viewer = NULL, selfviewer = NULL;
790: PetscMPIInt rank;
792: PetscFunctionBegin;
793: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
794: PetscCall(DMPlexGetDepth(dm, &depth));
795: PetscCall(DMGetPointSF(dm, &sf));
796: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &lpoints, &rpoints));
797: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
798: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
799: PetscCall(PetscMalloc2(Nr, localFace, Nr, remoteFace));
800: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
801: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
802: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &fdepth));
803: if (!fdepth) faceIsVertex = PETSC_TRUE;
804: for (PetscInt r = 0; r < Nr; ++r) {
805: (*localFace)[r].comp = -1;
806: (*remoteFace)[r].comp = -1;
807: }
808: // Get information for shared faces
809: if (debug) {
810: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
811: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
812: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
813: }
814: for (PetscInt f = fStart; f < fEnd; ++f) {
815: const PetscInt face = faces ? faces[f] : f;
816: const PetscInt *supp, *cone, *ornt;
817: PetscInt lf, sS, cS, cind = -1, c;
819: (*localFace)[face].comp = -1;
820: (*localFace)[face].cell = -1;
821: (*localFace)[face].ornt = 0;
822: PetscCall(PetscFindInt(face, Nl, lpoints, &lf));
823: if (!rootdegree[face] && lf < 0) continue;
824: PetscCall(DMPlexGetSupportSize(dm, face, &sS));
825: PetscCall(DMPlexGetSupport(dm, face, &supp));
826: for (PetscInt s = 0; s < sS; ++s) {
827: PetscInt pdepth, l;
829: // Filter support
830: cind = GetPointIndex(supp[s], cStart, cEnd, cells);
831: if (cind < 0) continue;
832: // Ignore overlapping cells, but not for embedded manifolds
833: PetscCall(DMPlexGetPointDepth(dm, supp[s], &pdepth));
834: PetscCall(PetscFindInt(supp[s], Nl, lpoints, &l));
835: if (pdepth == depth && l >= 0) continue;
836: (*localFace)[face].cell = supp[s];
837: break;
838: }
839: (*localFace)[face].comp = faceComp[f - fStart];
840: // Cannot determine orientation without a cell
841: if (cind < 0) continue;
842: PetscCheck(0 <= faceComp[f - fStart] && faceComp[f - fStart] < Ncomp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT " for face %" PetscInt_FMT " (%" PetscInt_FMT ")", faceComp[f - fStart], face, f);
843: PetscCall(DMPlexGetOrientedCone(dm, (*localFace)[face].cell, &cone, &ornt));
844: PetscCall(DMPlexGetConeSize(dm, (*localFace)[face].cell, &cS));
845: for (c = 0; c < cS; ++c)
846: if (cone[c] == face) break;
847: PetscCheck(c < cS, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " not found in cone of cell %" PetscInt_FMT, face, (*localFace)[face].cell);
848: if (faceIsVertex) {
849: // Use cone position, shifted to -1 or 1
850: if (PetscBTLookup(cellFlip, cind)) (*localFace)[face].ornt = 1 - c * 2;
851: else (*localFace)[face].ornt = c * 2 - 1;
852: } else {
853: // Use orientation sense
854: if (PetscBTLookup(cellFlip, cind)) (*localFace)[face].ornt = ornt[c] < 0 ? -1 : 1;
855: else (*localFace)[face].ornt = ornt[c] < 0 ? 1 : -1;
856: }
857: PetscCall(DMPlexRestoreOrientedCone(dm, (*localFace)[face].cell, &cone, &ornt));
858: if (debug)
859: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: Local shared face %" PetscInt_FMT " component %" PetscInt_FMT " cell %" PetscInt_FMT " orientation %" PetscInt_FMT "\n", rank, face, (*localFace)[face].comp, (*localFace)[face].cell,
860: (*localFace)[face].ornt));
861: }
862: if (debug) {
863: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
864: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
865: }
866: // Get information from owners
867: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &MPIU_3INT));
868: PetscCallMPI(MPI_Type_commit(&MPIU_3INT));
869: PetscCall(PetscSFBcastBegin(sf, MPIU_3INT, *localFace, *remoteFace, MPI_REPLACE));
870: PetscCall(PetscSFBcastEnd(sf, MPIU_3INT, *localFace, *remoteFace, MPI_REPLACE));
871: PetscCallMPI(MPI_Type_free(&MPIU_3INT));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: typedef struct {
876: PetscInt lface; // local shared face
877: PetscInt lcell; // local cell determining face orientation
878: PetscInt lornt; // local face orientation
879: PetscInt rrank; // remote rank
880: PetscInt rcomp; // remote connected component
881: PetscInt rface; // remote shared face
882: PetscInt rcell; // remote cell determining face orientation
883: PetscInt rornt; // remote face orientation
884: } Neighbor;
886: static PetscErrorCode DMPlexOrientCreateNeighbors_Internal(DM dm, IS cellIS, IS faceIS, PetscInt Ncomp, const PetscInt faceComp[], const SharedFace localFace[], const SharedFace remoteFace[], PetscInt **Nneigh, Neighbor ***neighbors)
887: {
888: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
889: const PetscInt *cells = NULL, *faces = NULL;
890: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
891: PetscSF sf;
892: const PetscInt *lpoints;
893: const PetscSFNode *rpoints;
894: PetscInt Nl;
895: PetscHSetIJ *ht;
896: PetscInt *counts;
897: PetscViewer viewer = NULL, selfviewer = NULL;
898: PetscMPIInt rank;
900: PetscFunctionBegin;
901: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
902: if (debug) {
903: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
904: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
905: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
906: }
907: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
908: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
909: PetscCall(DMGetPointSF(dm, &sf));
910: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
911: PetscCall(PetscCalloc2(Ncomp, &ht, Ncomp, &counts));
912: for (PetscInt c = 0; c < Ncomp; ++c) PetscCall(PetscHSetIJCreate(&ht[c]));
913: PetscCall(PetscCalloc2(Ncomp, Nneigh, Ncomp, neighbors));
914: // Count the number of unique connections between process/component
915: for (PetscInt f = fStart; f < fEnd; ++f) {
916: const PetscInt face = faces ? faces[f] : f;
917: PetscHashIJKey key;
918: PetscInt l;
920: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
921: if (l < 0) continue;
922: // Remote face is not part of the label
923: if (remoteFace[face].comp < 0) continue;
924: // Remote face is isolated from any surface cell, so cannot get orientation
925: if (!remoteFace[face].ornt) continue;
926: key.i = rpoints[l].rank;
927: key.j = remoteFace[face].comp;
928: PetscCheck(0 <= remoteFace[face].comp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT " on remote face %" PetscInt_FMT " from rank %" PetscInt_FMT " for local face %" PetscInt_FMT, remoteFace[face].comp, rpoints[l].index,
929: rpoints[l].rank, face);
930: PetscCall(PetscHSetIJAdd(ht[faceComp[f - fStart]], key));
931: }
932: // Get the sizes from each hash
933: for (PetscInt c = 0; c < Ncomp; ++c) {
934: PetscCall(PetscHSetIJGetSize(ht[c], &(*Nneigh)[c]));
935: PetscCall(PetscMalloc1((*Nneigh)[c], &(*neighbors)[c]));
936: PetscCall(PetscHSetIJClear(ht[c]));
937: if (debug) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %" PetscInt_FMT ", Found %" PetscInt_FMT " connections\n", rank, c, (*Nneigh)[c]));
938: }
939: // Gather the neighbor information
940: for (PetscInt f = fStart; f < fEnd; ++f) {
941: const PetscInt face = faces ? faces[f] : f;
942: const PetscInt comp = faceComp[f - fStart];
943: const PetscInt ind = counts[comp];
944: PetscHashIJKey key;
945: PetscBool missing;
946: PetscInt l;
948: PetscCheck(0 <= comp && comp < Ncomp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT, comp);
949: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
950: if (l < 0) continue;
951: if (remoteFace[face].comp < 0) continue;
952: if (!remoteFace[face].ornt) continue;
953: key.i = rpoints[l].rank;
954: key.j = remoteFace[face].comp;
955: PetscCall(PetscHSetIJQueryAdd(ht[comp], key, &missing));
956: if (!missing) continue;
957: PetscCheck(localFace[face].ornt != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid local face %" PetscInt_FMT " orientation: %" PetscInt_FMT, face, localFace[face].ornt);
958: (*neighbors)[comp][ind].lface = face;
959: (*neighbors)[comp][ind].lcell = localFace[face].cell;
960: (*neighbors)[comp][ind].lornt = localFace[face].ornt;
961: (*neighbors)[comp][ind].rrank = rpoints[l].rank;
962: (*neighbors)[comp][ind].rcomp = remoteFace[face].comp;
963: (*neighbors)[comp][ind].rface = rpoints[l].index;
964: (*neighbors)[comp][ind].rcell = remoteFace[face].cell;
965: (*neighbors)[comp][ind].rornt = remoteFace[face].ornt;
966: if (debug)
967: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %" PetscInt_FMT ", Found representative %" PetscInt_FMT " leaf %" PetscInt_FMT " (face %" PetscInt_FMT ") connecting to face %" PetscInt_FMT " on (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with orientation %" PetscInt_FMT "\n", rank,
968: remoteFace[face].comp, ind, l, face, (*neighbors)[comp][ind].rface, (*neighbors)[comp][ind].rrank, (*neighbors)[comp][ind].rcomp, (*neighbors)[comp][ind].rcell, (*neighbors)[comp][ind].rornt));
969: ++counts[comp];
970: }
971: // Cleanup
972: if (debug) {
973: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
974: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
975: }
976: for (PetscInt c = 0; c < Ncomp; ++c) {
977: PetscCheck(counts[c] == (*Nneigh)[c], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Neigh count for component %" PetscInt_FMT ": %" PetscInt_FMT " != %" PetscInt_FMT " allocated size", c, counts[c], (*Nneigh)[c]);
978: PetscCall(PetscHSetIJDestroy(&ht[c]));
979: }
980: PetscCall(PetscFree2(ht, counts));
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: static PetscErrorCode DMPlexOrientCreateProcessGraph_Internal(DM dm, IS faceIS, PetscInt Ncomp, const PetscInt Nneigh[], Neighbor *neighbors[], PetscSFNode **neighborAdj, PetscBool **neighborVal)
985: {
986: const PetscInt *faces = NULL;
987: PetscInt fStart = 0, fEnd = 0;
988: PetscSFNode *nrankComp; // The (rank, comp) of each neighbor
989: PetscBool *match; // Whether neighbors currently match
990: PetscSF sf;
991: const PetscInt *lpoints;
992: const PetscSFNode *rpoints;
993: PetscInt totNeighbors = 0, Nl, fdepth;
994: PetscBool faceIsVertex = PETSC_FALSE;
996: PetscFunctionBegin;
997: PetscCall(DMGetPointSF(dm, &sf));
998: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
999: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
1000: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &fdepth));
1001: if (!fdepth) faceIsVertex = PETSC_TRUE;
1002: for (PetscInt c = 0; c < Ncomp; ++c) totNeighbors += Nneigh[c];
1003: PetscCall(PetscMalloc2(totNeighbors, neighborAdj, totNeighbors, neighborVal));
1004: nrankComp = *neighborAdj;
1005: match = *neighborVal;
1006: for (PetscInt c = 0, off = 0; c < Ncomp; ++c) {
1007: for (PetscInt n = 0; n < Nneigh[c]; ++n, ++off) {
1008: PetscInt l;
1010: if (faceIsVertex) {
1011: match[off] = neighbors[c][n].lornt != neighbors[c][n].rornt ? PETSC_TRUE : PETSC_FALSE;
1012: } else {
1013: const PetscInt o = neighbors[c][n].lornt * neighbors[c][n].rornt;
1015: if (o < 0) match[off] = PETSC_TRUE;
1016: else match[off] = PETSC_FALSE;
1017: }
1018: // Flip sense if we are matching from an unowned cell
1019: PetscCall(PetscFindInt(neighbors[c][n].lcell, Nl, lpoints, &l));
1020: if (l >= 0) {
1021: if (rpoints[l].rank == neighbors[c][n].rrank && rpoints[l].index == neighbors[c][n].rcell) match[off] = match[off] ? PETSC_FALSE : PETSC_TRUE;
1022: }
1023: nrankComp[off].rank = neighbors[c][n].rrank;
1024: nrankComp[off].index = neighbors[c][n].rcomp;
1025: }
1026: }
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: static PetscErrorCode DMPlexOrientSolveProcessGraph_Internal(DM dm, IS cellIS, PetscInt Ncomp, const PetscInt cellComp[], const PetscInt Nneigh[], const PetscSFNode nrankComp[], const PetscBool match[], PetscBT cellFlip)
1031: {
1032: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1033: const PetscInt *cells = NULL;
1034: PetscInt cStart = 0, cEnd = 0;
1035: PetscSF sf;
1036: const PetscInt *lpoints;
1037: const PetscSFNode *rpoints;
1038: PetscInt Nl, totNeighbors = 0;
1039: PetscBool *flipped;
1040: MPI_Comm comm;
1041: PetscMPIInt rank;
1043: PetscFunctionBegin;
1044: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1045: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1046: PetscCall(DMGetPointSF(dm, &sf));
1047: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
1048: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1049: /* Collect the graph on 0 */
1050: for (PetscInt c = 0; c < Ncomp; ++c) totNeighbors += Nneigh[c];
1051: if (Nl >= 0) {
1052: Mat G;
1053: PetscBT seenProcs, flippedProcs;
1054: PetscInt *procFIFO, pTop, pBottom;
1055: PetscInt *N = NULL, *Noff;
1056: PetscSFNode *adj = NULL;
1057: PetscBool *val = NULL;
1058: PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc;
1059: PetscMPIInt size = 0, iNcomp, itotNeighbors;
1061: PetscCall(PetscCalloc1(Ncomp, &flipped));
1062: if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
1063: PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
1064: PetscCallMPI(MPI_Gather(&Ncomp, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
1065: for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
1066: if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
1067: PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
1068: PetscCallMPI(MPI_Gatherv(Nneigh, iNcomp, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
1069: for (PetscInt p = 0, o = 0; p < size; ++p) {
1070: recvcounts[p] = 0;
1071: for (PetscInt c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
1072: displs[p + 1] = displs[p] + recvcounts[p];
1073: }
1074: if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
1075: PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors));
1076: PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm));
1077: PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPI_C_BOOL, val, recvcounts, displs, MPI_C_BOOL, 0, comm));
1078: if (rank == 0) {
1079: for (PetscInt p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
1080: if (debug) {
1081: for (PetscInt p = 0, off = 0; p < size; ++p) {
1082: for (PetscInt c = 0; c < Nc[p]; ++c) {
1083: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %" PetscInt_FMT " Comp %" PetscInt_FMT ":\n", p, c));
1084: for (PetscInt n = 0; n < N[Noff[p] + c]; ++n, ++off) PetscCall(PetscPrintf(PETSC_COMM_SELF, " edge (%" PetscInt_FMT ", %" PetscInt_FMT ") (%s):\n", adj[off].rank, adj[off].index, PetscBools[val[off]]));
1085: }
1086: }
1087: }
1088: /* Symmetrize the graph */
1089: PetscCall(MatCreate(PETSC_COMM_SELF, &G));
1090: PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
1091: PetscCall(MatSetUp(G));
1092: for (PetscInt p = 0, off = 0; p < size; ++p) {
1093: for (PetscInt c = 0; c < Nc[p]; ++c) {
1094: const PetscInt r = Noff[p] + c;
1096: for (PetscInt n = 0; n < N[r]; ++n, ++off) {
1097: const PetscInt q = Noff[adj[off].rank] + adj[off].index;
1098: const PetscScalar o = val[off] ? 1.0 : 0.0;
1100: // Do not set values for processes that have no face to orient
1101: if (!Nc[adj[off].rank]) continue;
1102: PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
1103: PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
1104: }
1105: }
1106: }
1107: PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
1108: PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
1110: PetscCall(PetscBTCreate(Noff[size], &seenProcs));
1111: PetscCall(PetscBTMemzero(Noff[size], seenProcs));
1112: PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
1113: PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
1114: PetscCall(PetscMalloc1(Noff[size], &procFIFO));
1115: pTop = pBottom = 0;
1116: for (PetscInt p = 0; p < Noff[size]; ++p) {
1117: if (PetscBTLookup(seenProcs, p)) continue;
1118: /* Initialize FIFO with next proc */
1119: procFIFO[pBottom++] = p;
1120: PetscCall(PetscBTSet(seenProcs, p));
1121: /* Consider each proc in FIFO */
1122: while (pTop < pBottom) {
1123: const PetscScalar *ornt;
1124: const PetscInt *neighbors;
1125: PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors;
1127: proc = procFIFO[pTop++];
1128: flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
1129: PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
1130: /* Loop over neighboring procs */
1131: for (PetscInt n = 0; n < numNeighbors; ++n) {
1132: nproc = neighbors[n];
1133: mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
1134: seen = PetscBTLookup(seenProcs, nproc);
1135: flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
1137: if (mismatch ^ (flippedA ^ flippedB)) {
1138: PetscCheck(!seen, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %" PetscInt_FMT " and %" PetscInt_FMT " do not match: Fault mesh is non-orientable", proc, nproc);
1139: PetscCheck(!flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
1140: PetscCall(PetscBTSet(flippedProcs, nproc));
1141: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
1142: if (!seen) {
1143: procFIFO[pBottom++] = nproc;
1144: PetscCall(PetscBTSet(seenProcs, nproc));
1145: }
1146: }
1147: }
1148: }
1149: PetscCall(PetscFree(procFIFO));
1150: PetscCall(MatDestroy(&G));
1151: PetscCall(PetscFree2(adj, val));
1152: PetscCall(PetscBTDestroy(&seenProcs));
1153: }
1154: /* Scatter flip flags */
1155: {
1156: PetscBool *flips = NULL;
1158: if (rank == 0) {
1159: PetscCall(PetscMalloc1(Noff[size], &flips));
1160: for (PetscInt p = 0; p < Noff[size]; ++p) {
1161: flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
1162: if (debug && flips[p]) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flipping Proc+Comp %" PetscInt_FMT ":\n", p));
1163: }
1164: for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
1165: }
1166: PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
1167: PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPI_C_BOOL, flipped, iNcomp, MPI_C_BOOL, 0, comm));
1168: PetscCall(PetscFree(flips));
1169: }
1170: if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
1171: PetscCall(PetscFree(N));
1172: PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
1174: /* Decide whether to flip cells in each component */
1175: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
1176: if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(cellFlip, c));
1177: }
1178: PetscCall(PetscFree(flipped));
1179: }
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1183: // Reverse flipped cells in the mesh
1184: static PetscErrorCode DMPlexOrientPlex_Internal(DM dm, IS cellIS, PetscBT cellFlip)
1185: {
1186: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1187: const PetscInt *cells = NULL;
1188: PetscInt cStart = 0, cEnd = 0;
1189: PetscSF sf;
1190: const PetscInt *rootdegree = NULL;
1191: PetscInt *points = NULL;
1192: PetscInt pStart, pEnd, Nr, depth, cdepth = -1;
1193: PetscViewer v;
1194: PetscMPIInt rank;
1196: PetscFunctionBegin;
1197: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1198: if (debug) {
1199: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &v));
1200: PetscCall(PetscViewerASCIIPushSynchronized(v));
1201: }
1202: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1203: PetscCall(DMPlexGetDepth(dm, &depth));
1204: if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
1205: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &cdepth, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
1206: PetscCall(DMGetPointSF(dm, &sf));
1207: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
1208: if (Nr >= 0) {
1209: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1210: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1211: }
1212: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1213: if (cdepth == depth) PetscCall(PetscCalloc1(pEnd - pStart, &points));
1214: for (PetscInt c = cStart; c < cEnd; ++c) {
1215: if (PetscBTLookup(cellFlip, c - cStart)) {
1216: const PetscInt cell = cells ? cells[c] : c;
1218: PetscCall(DMPlexOrientPoint(dm, cell, -1));
1219: if (points && rootdegree && rootdegree[cell]) points[cell] = 1;
1220: if (debug) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT "%s\n", rank, cell, points && rootdegree && rootdegree[cell] ? " and sending to overlap" : ""));
1221: }
1222: }
1223: // Propagate flips for volumetric cells in the overlap
1224: if (cdepth == depth) {
1225: if (Nr >= 0) {
1226: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, points, points, MPI_SUM));
1227: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, points, points, MPI_SUM));
1228: }
1229: for (PetscInt c = cStart; c < cEnd; ++c) {
1230: const PetscInt cell = cells ? cells[c] : c;
1232: if (points[cell] && !PetscBTLookup(cellFlip, c - cStart)) {
1233: PetscCall(DMPlexOrientPoint(dm, cell, -1));
1234: if (debug) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT " through overlap\n", rank, cell));
1235: }
1236: }
1237: PetscCall(PetscFree(points));
1238: }
1239: if (debug) {
1240: PetscCall(PetscViewerFlush(v));
1241: PetscCall(PetscViewerASCIIPopSynchronized(v));
1242: }
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: /*
1247: The orientation operation needs to perform two functions:
1249: 1) Orient a bulk mesh. Here we assume that faces are in the correct order, but may
1250: be reversed.
1252: 2) Orient a manifold embedded in a given mesh. This introduces many complications, since
1253: lower-dimensional pieces of the manifold might be shared with other processes.
1255: We will divide the operation into phases
1257: 1) Each process orients the local mesh or manifold, and also counts the connected components.
1259: 2) Send information about all neighboring faces to leaves
1261: 3) Determine which processes are connected to each component, and pick a representative shared face connecting them
1263: For volumetric meshes (ignoring overlap cells), shared faces connect a single cell on each process.
1265: For manifolds, interior faces can be shared, so multiple cells can be connected. We must determine whether the cells
1266: we are using to compute orientation are identified.
1268: 4) Determine whether the orientations match across each process boundary, to make the local piece of the process graph
1270: For manifolds, if the assoicated cells are identified, we flip the sense of the comparison.
1272: 5) Gather the process graph to process 0, compute a satisfying orientation, communicate the orientation
1274: 6) Flip the relevant cells and for volumetric meshes, propagate reordering to overlap cells
1275: */
1276: PetscErrorCode DMPlexOrientCells_Internal(DM dm, IS cellIS, IS faceIS)
1277: {
1278: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1279: const PetscInt *cells = NULL, *faces = NULL;
1280: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
1281: PetscBT cellFlip; // The bit is true if a cell should have its orientation reversed
1282: PetscInt *cellComp; // The connected component number of each cell
1283: PetscInt *faceComp; // The connected component number of each face
1284: PetscInt Ncomp = 0; // The number of local connected components
1285: SharedFace *localFace; // Holds local information for owned and ghost shared faces
1286: SharedFace *remoteFace; // Holds remote information from owners of shared faces
1287: PetscInt *Nneigh; // The number of neighboring processes for each component
1288: Neighbor **neighbors; // The neighbor for each component
1289: PetscSFNode *neighborAdj; // The (rank, comp) of each neighbor
1290: PetscBool *neighborVal; // Whether neighbors currently match
1291: MPI_Comm comm;
1292: PetscMPIInt rank, size;
1294: PetscFunctionBegin;
1295: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1296: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1297: PetscCallMPI(MPI_Comm_size(comm, &size));
1298: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1299: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
1300: PetscCall(PetscBTCreate(cEnd - cStart, &cellFlip));
1301: PetscCall(PetscBTMemzero(cEnd - cStart, cellFlip));
1302: PetscCall(PetscCalloc2(cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
1303: // Phase 1: Serial Orientation
1304: PetscCall(DMPlexOrient_Serial(dm, cellIS, faceIS, &Ncomp, cellComp, faceComp, cellFlip));
1305: if (debug) {
1306: PetscViewer v;
1307: PetscInt cdepth = -1;
1309: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1310: PetscCall(PetscViewerASCIIPushSynchronized(v));
1311: if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
1312: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]New Orientation %" PetscInt_FMT " cells (depth %" PetscInt_FMT ") and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, cdepth, fEnd - fStart));
1313: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
1314: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
1315: PetscCall(PetscViewerFlush(v));
1316: PetscCall(PetscViewerASCIIPopSynchronized(v));
1317: }
1318: if (size == 1) goto end;
1319: // Phase 2
1320: PetscCall(DMPlexOrientCreateSharedFaces_Internal(dm, cellIS, faceIS, Ncomp, faceComp, cellFlip, &localFace, &remoteFace));
1321: // Phase 3
1322: PetscCall(DMPlexOrientCreateNeighbors_Internal(dm, cellIS, faceIS, Ncomp, faceComp, localFace, remoteFace, &Nneigh, &neighbors));
1323: PetscCall(PetscFree2(localFace, remoteFace));
1324: // Phase 4
1325: PetscCall(DMPlexOrientCreateProcessGraph_Internal(dm, faceIS, Ncomp, Nneigh, neighbors, &neighborAdj, &neighborVal));
1326: // Phase 5
1327: PetscCall(DMPlexOrientSolveProcessGraph_Internal(dm, cellIS, Ncomp, cellComp, Nneigh, neighborAdj, neighborVal, cellFlip));
1328: if (debug) {
1329: PetscViewer v;
1331: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1332: PetscCall(PetscViewerASCIIPushSynchronized(v));
1333: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
1334: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
1335: PetscCall(PetscViewerFlush(v));
1336: PetscCall(PetscViewerASCIIPopSynchronized(v));
1337: }
1338: for (PetscInt c = 0; c < Ncomp; ++c) PetscCall(PetscFree(neighbors[c]));
1339: PetscCall(PetscFree2(Nneigh, neighbors));
1340: PetscCall(PetscFree2(neighborAdj, neighborVal));
1341: end:
1342: // Phase 6
1343: PetscCall(DMPlexOrientPlex_Internal(dm, cellIS, cellFlip));
1344: PetscCall(PetscBTDestroy(&cellFlip));
1345: PetscCall(PetscFree2(cellComp, faceComp));
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: static PetscErrorCode DMPlexCheckOrientation_Internal(DM dm, IS cellIS, IS faceIS)
1350: {
1351: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1352: PetscViewer viewer = NULL, selfviewer = NULL;
1353: PetscSF sf;
1354: const PetscInt *lpoints, *rootdegree = NULL;
1355: const PetscSFNode *rpoints;
1356: const PetscInt *cells = NULL, *faces = NULL;
1357: PetscSFNode *oornt, *rornt, *lornt;
1358: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
1359: PetscInt pdepth, Nr, Nl;
1360: PetscBool faceIsVertex = PETSC_FALSE, valid = PETSC_TRUE;
1361: MPI_Comm comm;
1362: PetscMPIInt size, rank;
1364: PetscFunctionBegin;
1365: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1366: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1367: PetscCallMPI(MPI_Comm_size(comm, &size));
1368: if (debug) {
1369: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
1370: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1371: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
1372: }
1374: PetscCall(DMGetPointSF(dm, &sf));
1375: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &lpoints, &rpoints));
1376: if (Nr >= 0) {
1377: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1378: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1379: } else {
1380: sf = NULL;
1381: Nr = 0;
1382: Nl = 0;
1383: }
1384: PetscCall(PetscCalloc3(Nr, &oornt, Nr, &rornt, Nl, &lornt));
1385: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1386: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
1387: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &pdepth));
1388: if (!pdepth) faceIsVertex = PETSC_TRUE;
1389: if (debug) {
1390: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Checking orientation of %" PetscInt_FMT " cells and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, fEnd - fStart));
1391: PetscCall(PetscViewerASCIIPushTab(selfviewer));
1392: }
1393: for (PetscInt f = fStart; f < fEnd; ++f) {
1394: const PetscInt face = faces ? faces[f] : f;
1395: const PetscBool owner = rootdegree && rootdegree[face] ? PETSC_TRUE : PETSC_FALSE;
1396: const PetscInt *supp;
1397: PetscInt neighbors[2], o[2] = {0, 0};
1398: PetscInt lsS = 0, sS;
1400: PetscCall(DMPlexGetSupport(dm, face, &supp));
1401: PetscCall(DMPlexGetSupportSize(dm, face, &sS));
1402: // Filter support for local cells
1403: for (PetscInt s = 0; s < sS; ++s) {
1404: PetscInt ind;
1406: ind = GetPointIndex(supp[s], cStart, cEnd, cells);
1407: if (ind >= 0) {
1408: neighbors[PetscMin(lsS, 1)] = supp[s];
1409: ++lsS;
1410: }
1411: }
1412: PetscCheck(lsS < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " > 2", face, lsS);
1413: // Extract orientations
1414: for (PetscInt s = 0; s < lsS; ++s) {
1415: const PetscInt *cone, *ornt;
1416: PetscInt cS;
1418: PetscCall(DMPlexGetConeSize(dm, neighbors[s], &cS));
1419: PetscCall(DMPlexGetOrientedCone(dm, neighbors[s], &cone, &ornt));
1420: for (PetscInt c = 0; c < cS; ++c) {
1421: if (cone[c] == face) {
1422: if (faceIsVertex) o[s] = c * 2 - 1;
1423: else o[s] = ornt[c] < 0 ? -1 : 1;
1424: break;
1425: }
1426: }
1427: PetscCall(DMPlexRestoreOrientedCone(dm, neighbors[s], &cone, &ornt));
1428: }
1429: if (lsS == 2) {
1430: // Check internal face
1431: if (o[0] * o[1] >= 0) {
1432: valid = PETSC_FALSE;
1433: if (debug)
1434: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Internal Face %" PetscInt_FMT " is mismatched: cell %" PetscInt_FMT " (%" PetscInt_FMT ") ~ cell %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, face, neighbors[0], o[0], neighbors[1], o[1]));
1435: } else if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Internal Face %" PetscInt_FMT " valid\n", rank, face));
1436: } else {
1437: // Check shared and boundary faces
1438: PetscInt l;
1440: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
1441: // Boundary face
1442: if (l < 0 && !owner) {
1443: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Boundary Face %" PetscInt_FMT " valid\n", rank, face));
1444: continue;
1445: }
1446: if (l >= 0) {
1447: lornt[l].index = neighbors[0];
1448: lornt[l].rank = o[0];
1449: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Ghost Face %" PetscInt_FMT " (%" PetscInt_FMT ") was stored\n", rank, face, o[0]));
1450: }
1451: }
1452: if (owner) {
1453: // Store shared face orientation and cell from owner
1454: oornt[face].index = neighbors[0];
1455: oornt[face].rank = o[0];
1456: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Owned Face %" PetscInt_FMT " (%" PetscInt_FMT ") was stored\n", rank, face, o[0]));
1457: }
1458: }
1459: // Communicate shared face orientations from owner
1460: if (sf) {
1461: PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, oornt, rornt, MPI_REPLACE));
1462: PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, oornt, rornt, MPI_REPLACE));
1463: }
1464: // Check unowned shared faces
1465: for (PetscInt l = 0; l < Nl; ++l) {
1466: const PetscInt face = lpoints ? lpoints[l] : l;
1467: PetscBool flip = PETSC_FALSE;
1468: PetscInt ind, cl, o[2] = {0, 0};
1470: // Filter for local faces
1471: ind = GetPointIndex(face, fStart, fEnd, faces);
1472: if (ind < 0) continue;
1473: // Check for shared cell
1474: PetscCall(PetscFindInt(lornt[l].index, Nl, lpoints, &cl));
1475: if (cl >= 0) {
1476: if (rpoints[cl].index == rornt[face].index) {
1477: flip = PETSC_TRUE;
1478: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "Shared cell %" PetscInt_FMT " maps to %" PetscInt_FMT " (%" PetscInt_FMT ") so we reverse the orientation check\n", lornt[l].index, rpoints[cl].index, rpoints[cl].rank));
1479: } else {
1480: if (debug > 1)
1481: PetscCall(PetscViewerASCIIPrintf(selfviewer, "Shared cell %" PetscInt_FMT " maps to %" PetscInt_FMT " (%" PetscInt_FMT ") instead of %" PetscInt_FMT " so we do not reverse the orientation check\n", lornt[l].index, rpoints[cl].index,
1482: rpoints[cl].rank, rornt[face].index));
1483: }
1484: }
1485: o[0] = lornt[l].rank;
1486: o[1] = rornt[face].rank;
1487: // This was not a shared face
1488: if (!o[0]) continue;
1489: if ((o[0] * o[1] >= 0 && !flip) || (o[0] * o[1] < 0 && flip)) {
1490: valid = PETSC_FALSE;
1491: if (debug)
1492: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Ghost Face %" PetscInt_FMT " (%" PetscInt_FMT ") does not match face %" PetscInt_FMT " rank %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, face, o[0], rpoints[l].index, rpoints[l].rank, o[1]));
1493: } else if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Ghost Face %" PetscInt_FMT " matches %" PetscInt_FMT " (%" PetscInt_FMT ") valid\n", rank, face, rpoints[l].index, rpoints[l].rank));
1494: }
1495: // Cleanup
1496: PetscCall(PetscFree3(oornt, rornt, lornt));
1497: if (debug) {
1498: PetscCall(PetscViewerASCIIPopTab(selfviewer));
1499: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
1500: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1501: }
1502: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &valid, 1, MPI_C_BOOL, MPI_LAND, comm));
1503: PetscCheck(valid, comm, PETSC_ERR_ARG_WRONGSTATE, "Mesh was not properly oriented");
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1507: /*@
1508: DMPlexCheckOrientationLabel - Check that the surface defined by the given `DMLabel` is oriented
1510: Collective
1512: Input Parameters:
1513: + dm - The `DM`
1514: - label - The `DMLabel` defining the embedded surface
1516: Level: advanced
1518: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`
1519: @*/
1520: PetscErrorCode DMPlexCheckOrientationLabel(DM dm, DMLabel label)
1521: {
1522: IS cellIS, faceIS;
1524: PetscFunctionBegin;
1525: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
1526: PetscCall(DMPlexCheckOrientation_Internal(dm, cellIS, faceIS));
1527: PetscCall(ISDestroy(&cellIS));
1528: PetscCall(ISDestroy(&faceIS));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }