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: Note:
333: The orientation data for the `DM` are change 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()`
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: PetscErrorCode DMPlexOrientLabel(DM dm, DMLabel label)
742: {
743: IS cellIS, faceIS;
745: PetscFunctionBegin;
746: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
747: PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
748: PetscCall(ISDestroy(&cellIS));
749: PetscCall(ISDestroy(&faceIS));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: typedef struct {
754: PetscInt comp; // connected component number
755: PetscInt cell; // cell determining face orientation
756: PetscInt ornt; // face orientation
757: } SharedFace;
759: static PetscErrorCode DMPlexOrientCreateSharedFaces_Internal(DM dm, IS cellIS, IS faceIS, PetscInt Ncomp, const PetscInt faceComp[], PetscBT cellFlip, SharedFace **localFace, SharedFace **remoteFace)
760: {
761: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
762: const PetscInt *cells = NULL, *faces = NULL;
763: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
764: PetscSF sf;
765: const PetscInt *lpoints, *rootdegree;
766: const PetscSFNode *rpoints;
767: PetscInt Nr, Nl;
768: PetscInt depth, fdepth;
769: PetscBool faceIsVertex = PETSC_FALSE;
770: MPI_Datatype MPIU_3INT;
771: PetscViewer viewer = NULL, selfviewer = NULL;
772: PetscMPIInt rank;
774: PetscFunctionBegin;
775: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
776: PetscCall(DMPlexGetDepth(dm, &depth));
777: PetscCall(DMGetPointSF(dm, &sf));
778: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &lpoints, &rpoints));
779: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
780: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
781: PetscCall(PetscMalloc2(Nr, localFace, Nr, remoteFace));
782: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
783: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
784: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &fdepth));
785: if (!fdepth) faceIsVertex = PETSC_TRUE;
786: for (PetscInt r = 0; r < Nr; ++r) {
787: (*localFace)[r].comp = -1;
788: (*remoteFace)[r].comp = -1;
789: }
790: // Get information for shared faces
791: if (debug) {
792: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
793: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
794: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
795: }
796: for (PetscInt f = fStart; f < fEnd; ++f) {
797: const PetscInt face = faces ? faces[f] : f;
798: const PetscInt *supp, *cone, *ornt;
799: PetscInt lf, sS, cS, cind = -1, c;
801: (*localFace)[face].comp = -1;
802: (*localFace)[face].cell = -1;
803: (*localFace)[face].ornt = 0;
804: PetscCall(PetscFindInt(face, Nl, lpoints, &lf));
805: if (!rootdegree[face] && lf < 0) continue;
806: PetscCall(DMPlexGetSupportSize(dm, face, &sS));
807: PetscCall(DMPlexGetSupport(dm, face, &supp));
808: for (PetscInt s = 0; s < sS; ++s) {
809: PetscInt pdepth, l;
811: // Filter support
812: cind = GetPointIndex(supp[s], cStart, cEnd, cells);
813: if (cind < 0) continue;
814: // Ignore overlapping cells, but not for embedded manifolds
815: PetscCall(DMPlexGetPointDepth(dm, supp[s], &pdepth));
816: PetscCall(PetscFindInt(supp[s], Nl, lpoints, &l));
817: if (pdepth == depth && l >= 0) continue;
818: (*localFace)[face].cell = supp[s];
819: break;
820: }
821: (*localFace)[face].comp = faceComp[f - fStart];
822: // Cannot determine orientation without a cell
823: if (cind < 0) continue;
824: 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);
825: PetscCall(DMPlexGetOrientedCone(dm, (*localFace)[face].cell, &cone, &ornt));
826: PetscCall(DMPlexGetConeSize(dm, (*localFace)[face].cell, &cS));
827: for (c = 0; c < cS; ++c)
828: if (cone[c] == face) break;
829: 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);
830: if (faceIsVertex) {
831: // Use cone position, shifted to -1 or 1
832: if (PetscBTLookup(cellFlip, cind)) (*localFace)[face].ornt = 1 - c * 2;
833: else (*localFace)[face].ornt = c * 2 - 1;
834: } else {
835: // Use orientation sense
836: if (PetscBTLookup(cellFlip, cind)) (*localFace)[face].ornt = ornt[c] < 0 ? -1 : 1;
837: else (*localFace)[face].ornt = ornt[c] < 0 ? 1 : -1;
838: }
839: PetscCall(DMPlexRestoreOrientedCone(dm, (*localFace)[face].cell, &cone, &ornt));
840: if (debug)
841: 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,
842: (*localFace)[face].ornt));
843: }
844: if (debug) {
845: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
846: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
847: }
848: // Get information from owners
849: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &MPIU_3INT));
850: PetscCallMPI(MPI_Type_commit(&MPIU_3INT));
851: PetscCall(PetscSFBcastBegin(sf, MPIU_3INT, *localFace, *remoteFace, MPI_REPLACE));
852: PetscCall(PetscSFBcastEnd(sf, MPIU_3INT, *localFace, *remoteFace, MPI_REPLACE));
853: PetscCallMPI(MPI_Type_free(&MPIU_3INT));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: typedef struct {
858: PetscInt lface; // local shared face
859: PetscInt lcell; // local cell determining face orientation
860: PetscInt lornt; // local face orientation
861: PetscInt rrank; // remote rank
862: PetscInt rcomp; // remote connected component
863: PetscInt rface; // remote shared face
864: PetscInt rcell; // remote cell determining face orientation
865: PetscInt rornt; // remote face orientation
866: } Neighbor;
868: 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)
869: {
870: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
871: const PetscInt *cells = NULL, *faces = NULL;
872: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
873: PetscSF sf;
874: const PetscInt *lpoints;
875: const PetscSFNode *rpoints;
876: PetscInt Nl;
877: PetscHSetIJ *ht;
878: PetscInt *counts;
879: PetscViewer viewer = NULL, selfviewer = NULL;
880: PetscMPIInt rank;
882: PetscFunctionBegin;
883: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
884: if (debug) {
885: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
886: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
887: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
888: }
889: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
890: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
891: PetscCall(DMGetPointSF(dm, &sf));
892: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
893: PetscCall(PetscCalloc2(Ncomp, &ht, Ncomp, &counts));
894: for (PetscInt c = 0; c < Ncomp; ++c) PetscCall(PetscHSetIJCreate(&ht[c]));
895: PetscCall(PetscCalloc2(Ncomp, Nneigh, Ncomp, neighbors));
896: // Count the number of unique connections between process/component
897: for (PetscInt f = fStart; f < fEnd; ++f) {
898: const PetscInt face = faces ? faces[f] : f;
899: PetscHashIJKey key;
900: PetscInt l;
902: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
903: if (l < 0) continue;
904: // Remote face is not part of the label
905: if (remoteFace[face].comp < 0) continue;
906: // Remote face is isolated from any surface cell, so cannot get orientation
907: if (!remoteFace[face].ornt) continue;
908: key.i = rpoints[l].rank;
909: key.j = remoteFace[face].comp;
910: 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,
911: rpoints[l].rank, face);
912: PetscCall(PetscHSetIJAdd(ht[faceComp[f - fStart]], key));
913: }
914: // Get the sizes from each hash
915: for (PetscInt c = 0; c < Ncomp; ++c) {
916: PetscCall(PetscHSetIJGetSize(ht[c], &(*Nneigh)[c]));
917: PetscCall(PetscMalloc1((*Nneigh)[c], &(*neighbors)[c]));
918: PetscCall(PetscHSetIJClear(ht[c]));
919: if (debug) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %" PetscInt_FMT ", Found %" PetscInt_FMT " connections\n", rank, c, (*Nneigh)[c]));
920: }
921: // Gather the neighbor information
922: for (PetscInt f = fStart; f < fEnd; ++f) {
923: const PetscInt face = faces ? faces[f] : f;
924: const PetscInt comp = faceComp[f - fStart];
925: const PetscInt ind = counts[comp];
926: PetscHashIJKey key;
927: PetscBool missing;
928: PetscInt l;
930: PetscCheck(0 <= comp && comp < Ncomp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid component %" PetscInt_FMT, comp);
931: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
932: if (l < 0) continue;
933: if (remoteFace[face].comp < 0) continue;
934: if (!remoteFace[face].ornt) continue;
935: key.i = rpoints[l].rank;
936: key.j = remoteFace[face].comp;
937: PetscCall(PetscHSetIJQueryAdd(ht[comp], key, &missing));
938: if (!missing) continue;
939: PetscCheck(localFace[face].ornt != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid local face %" PetscInt_FMT " orientation: %" PetscInt_FMT, face, localFace[face].ornt);
940: (*neighbors)[comp][ind].lface = face;
941: (*neighbors)[comp][ind].lcell = localFace[face].cell;
942: (*neighbors)[comp][ind].lornt = localFace[face].ornt;
943: (*neighbors)[comp][ind].rrank = rpoints[l].rank;
944: (*neighbors)[comp][ind].rcomp = remoteFace[face].comp;
945: (*neighbors)[comp][ind].rface = rpoints[l].index;
946: (*neighbors)[comp][ind].rcell = remoteFace[face].cell;
947: (*neighbors)[comp][ind].rornt = remoteFace[face].ornt;
948: if (debug)
949: 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,
950: 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));
951: ++counts[comp];
952: }
953: // Cleanup
954: if (debug) {
955: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
956: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
957: }
958: for (PetscInt c = 0; c < Ncomp; ++c) {
959: 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]);
960: PetscCall(PetscHSetIJDestroy(&ht[c]));
961: }
962: PetscCall(PetscFree2(ht, counts));
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: static PetscErrorCode DMPlexOrientCreateProcessGraph_Internal(DM dm, IS faceIS, PetscInt Ncomp, const PetscInt Nneigh[], Neighbor *neighbors[], PetscSFNode **neighborAdj, PetscBool **neighborVal)
967: {
968: const PetscInt *faces = NULL;
969: PetscInt fStart = 0, fEnd = 0;
970: PetscSFNode *nrankComp; // The (rank, comp) of each neighbor
971: PetscBool *match; // Whether neighbors currently match
972: PetscSF sf;
973: const PetscInt *lpoints;
974: const PetscSFNode *rpoints;
975: PetscInt totNeighbors = 0, Nl, fdepth;
976: PetscBool faceIsVertex = PETSC_FALSE;
978: PetscFunctionBegin;
979: PetscCall(DMGetPointSF(dm, &sf));
980: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
981: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
982: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &fdepth));
983: if (!fdepth) faceIsVertex = PETSC_TRUE;
984: for (PetscInt c = 0; c < Ncomp; ++c) totNeighbors += Nneigh[c];
985: PetscCall(PetscMalloc2(totNeighbors, neighborAdj, totNeighbors, neighborVal));
986: nrankComp = *neighborAdj;
987: match = *neighborVal;
988: for (PetscInt c = 0, off = 0; c < Ncomp; ++c) {
989: for (PetscInt n = 0; n < Nneigh[c]; ++n, ++off) {
990: PetscInt l;
992: if (faceIsVertex) {
993: match[off] = neighbors[c][n].lornt != neighbors[c][n].rornt ? PETSC_TRUE : PETSC_FALSE;
994: } else {
995: const PetscInt o = neighbors[c][n].lornt * neighbors[c][n].rornt;
997: if (o < 0) match[off] = PETSC_TRUE;
998: else match[off] = PETSC_FALSE;
999: }
1000: // Flip sense if we are matching from an unowned cell
1001: PetscCall(PetscFindInt(neighbors[c][n].lcell, Nl, lpoints, &l));
1002: if (l >= 0) {
1003: if (rpoints[l].rank == neighbors[c][n].rrank && rpoints[l].index == neighbors[c][n].rcell) match[off] = match[off] ? PETSC_FALSE : PETSC_TRUE;
1004: }
1005: nrankComp[off].rank = neighbors[c][n].rrank;
1006: nrankComp[off].index = neighbors[c][n].rcomp;
1007: }
1008: }
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: static PetscErrorCode DMPlexOrientSolveProcessGraph_Internal(DM dm, IS cellIS, PetscInt Ncomp, const PetscInt cellComp[], const PetscInt Nneigh[], const PetscSFNode nrankComp[], const PetscBool match[], PetscBT cellFlip)
1013: {
1014: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1015: const PetscInt *cells = NULL;
1016: PetscInt cStart = 0, cEnd = 0;
1017: PetscSF sf;
1018: const PetscInt *lpoints;
1019: const PetscSFNode *rpoints;
1020: PetscInt Nl, totNeighbors = 0;
1021: PetscBool *flipped;
1022: MPI_Comm comm;
1023: PetscMPIInt rank;
1025: PetscFunctionBegin;
1026: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1027: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1028: PetscCall(DMGetPointSF(dm, &sf));
1029: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &lpoints, &rpoints));
1030: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1031: /* Collect the graph on 0 */
1032: for (PetscInt c = 0; c < Ncomp; ++c) totNeighbors += Nneigh[c];
1033: if (Nl >= 0) {
1034: Mat G;
1035: PetscBT seenProcs, flippedProcs;
1036: PetscInt *procFIFO, pTop, pBottom;
1037: PetscInt *N = NULL, *Noff;
1038: PetscSFNode *adj = NULL;
1039: PetscBool *val = NULL;
1040: PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc;
1041: PetscMPIInt size = 0, iNcomp, itotNeighbors;
1043: PetscCall(PetscCalloc1(Ncomp, &flipped));
1044: if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
1045: PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
1046: PetscCallMPI(MPI_Gather(&Ncomp, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
1047: for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
1048: if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
1049: PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
1050: PetscCallMPI(MPI_Gatherv(Nneigh, iNcomp, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
1051: for (PetscInt p = 0, o = 0; p < size; ++p) {
1052: recvcounts[p] = 0;
1053: for (PetscInt c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
1054: displs[p + 1] = displs[p] + recvcounts[p];
1055: }
1056: if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
1057: PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors));
1058: PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm));
1059: PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPI_C_BOOL, val, recvcounts, displs, MPI_C_BOOL, 0, comm));
1060: if (rank == 0) {
1061: for (PetscInt p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
1062: if (debug) {
1063: for (PetscInt p = 0, off = 0; p < size; ++p) {
1064: for (PetscInt c = 0; c < Nc[p]; ++c) {
1065: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %" PetscInt_FMT " Comp %" PetscInt_FMT ":\n", p, c));
1066: 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]]));
1067: }
1068: }
1069: }
1070: /* Symmetrize the graph */
1071: PetscCall(MatCreate(PETSC_COMM_SELF, &G));
1072: PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
1073: PetscCall(MatSetUp(G));
1074: for (PetscInt p = 0, off = 0; p < size; ++p) {
1075: for (PetscInt c = 0; c < Nc[p]; ++c) {
1076: const PetscInt r = Noff[p] + c;
1078: for (PetscInt n = 0; n < N[r]; ++n, ++off) {
1079: const PetscInt q = Noff[adj[off].rank] + adj[off].index;
1080: const PetscScalar o = val[off] ? 1.0 : 0.0;
1082: // Do not set values for processes that have no face to orient
1083: if (!Nc[adj[off].rank]) continue;
1084: PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
1085: PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
1086: }
1087: }
1088: }
1089: PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
1090: PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
1092: PetscCall(PetscBTCreate(Noff[size], &seenProcs));
1093: PetscCall(PetscBTMemzero(Noff[size], seenProcs));
1094: PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
1095: PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
1096: PetscCall(PetscMalloc1(Noff[size], &procFIFO));
1097: pTop = pBottom = 0;
1098: for (PetscInt p = 0; p < Noff[size]; ++p) {
1099: if (PetscBTLookup(seenProcs, p)) continue;
1100: /* Initialize FIFO with next proc */
1101: procFIFO[pBottom++] = p;
1102: PetscCall(PetscBTSet(seenProcs, p));
1103: /* Consider each proc in FIFO */
1104: while (pTop < pBottom) {
1105: const PetscScalar *ornt;
1106: const PetscInt *neighbors;
1107: PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors;
1109: proc = procFIFO[pTop++];
1110: flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
1111: PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
1112: /* Loop over neighboring procs */
1113: for (PetscInt n = 0; n < numNeighbors; ++n) {
1114: nproc = neighbors[n];
1115: mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
1116: seen = PetscBTLookup(seenProcs, nproc);
1117: flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
1119: if (mismatch ^ (flippedA ^ flippedB)) {
1120: 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);
1121: PetscCheck(!flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
1122: PetscCall(PetscBTSet(flippedProcs, nproc));
1123: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
1124: if (!seen) {
1125: procFIFO[pBottom++] = nproc;
1126: PetscCall(PetscBTSet(seenProcs, nproc));
1127: }
1128: }
1129: }
1130: }
1131: PetscCall(PetscFree(procFIFO));
1132: PetscCall(MatDestroy(&G));
1133: PetscCall(PetscFree2(adj, val));
1134: PetscCall(PetscBTDestroy(&seenProcs));
1135: }
1136: /* Scatter flip flags */
1137: {
1138: PetscBool *flips = NULL;
1140: if (rank == 0) {
1141: PetscCall(PetscMalloc1(Noff[size], &flips));
1142: for (PetscInt p = 0; p < Noff[size]; ++p) {
1143: flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
1144: if (debug && flips[p]) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flipping Proc+Comp %" PetscInt_FMT ":\n", p));
1145: }
1146: for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
1147: }
1148: PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
1149: PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPI_C_BOOL, flipped, iNcomp, MPI_C_BOOL, 0, comm));
1150: PetscCall(PetscFree(flips));
1151: }
1152: if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
1153: PetscCall(PetscFree(N));
1154: PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
1156: /* Decide whether to flip cells in each component */
1157: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
1158: if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(cellFlip, c));
1159: }
1160: PetscCall(PetscFree(flipped));
1161: }
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: // Reverse flipped cells in the mesh
1166: static PetscErrorCode DMPlexOrientPlex_Internal(DM dm, IS cellIS, PetscBT cellFlip)
1167: {
1168: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1169: const PetscInt *cells = NULL;
1170: PetscInt cStart = 0, cEnd = 0;
1171: PetscSF sf;
1172: const PetscInt *rootdegree = NULL;
1173: PetscInt *points = NULL;
1174: PetscInt pStart, pEnd, Nr, depth, cdepth = -1;
1175: PetscViewer v;
1176: PetscMPIInt rank;
1178: PetscFunctionBegin;
1179: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1180: if (debug) {
1181: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm), &v));
1182: PetscCall(PetscViewerASCIIPushSynchronized(v));
1183: }
1184: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1185: PetscCall(DMPlexGetDepth(dm, &depth));
1186: if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
1187: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &cdepth, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
1188: PetscCall(DMGetPointSF(dm, &sf));
1189: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
1190: if (Nr >= 0) {
1191: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1192: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1193: }
1194: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1195: if (cdepth == depth) PetscCall(PetscCalloc1(pEnd - pStart, &points));
1196: for (PetscInt c = cStart; c < cEnd; ++c) {
1197: if (PetscBTLookup(cellFlip, c - cStart)) {
1198: const PetscInt cell = cells ? cells[c] : c;
1200: PetscCall(DMPlexOrientPoint(dm, cell, -1));
1201: if (points && rootdegree && rootdegree[cell]) points[cell] = 1;
1202: if (debug) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT "%s\n", rank, cell, points && rootdegree && rootdegree[cell] ? " and sending to overlap" : ""));
1203: }
1204: }
1205: // Propagate flips for volumetric cells in the overlap
1206: if (cdepth == depth) {
1207: if (Nr >= 0) {
1208: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, points, points, MPI_SUM));
1209: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, points, points, MPI_SUM));
1210: }
1211: for (PetscInt c = cStart; c < cEnd; ++c) {
1212: const PetscInt cell = cells ? cells[c] : c;
1214: if (points[cell] && !PetscBTLookup(cellFlip, c - cStart)) {
1215: PetscCall(DMPlexOrientPoint(dm, cell, -1));
1216: if (debug) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT " through overlap\n", rank, cell));
1217: }
1218: }
1219: PetscCall(PetscFree(points));
1220: }
1221: if (debug) {
1222: PetscCall(PetscViewerFlush(v));
1223: PetscCall(PetscViewerASCIIPopSynchronized(v));
1224: }
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: /*
1229: The orientation operation needs to perform two functions:
1231: 1) Orient a bulk mesh. Here we assume that faces are in the correct order, but may
1232: be reversed.
1234: 2) Orient a manifold embedded in a given mesh. This introduces many complications, since
1235: lower-dimensional pieces of the manifold might be shared with other processes.
1237: We will divide the operation into phases
1239: 1) Each process orients the local mesh or manifold, and also counts the connected components.
1241: 2) Send information about all neighboring faces to leaves
1243: 3) Determine which processes are connected to each component, and pick a representative shared face connecting them
1245: For volumetric meshes (ignoring overlap cells), shared faces connect a single cell on each process.
1247: For manifolds, interior faces can be shared, so multiple cells can be connected. We must determine whether the cells
1248: we are using to compute orientation are identified.
1250: 4) Determine whether the orientations match across each process boundary, to make the local piece of the process graph
1252: For manifolds, if the assoicated cells are identified, we flip the sense of the comparison.
1254: 5) Gather the process graph to process 0, compute a satisfying orientation, communicate the orientation
1256: 6) Flip the relevant cells and for volumetric meshes, propagate reordering to overlap cells
1257: */
1258: PetscErrorCode DMPlexOrientCells_Internal(DM dm, IS cellIS, IS faceIS)
1259: {
1260: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1261: const PetscInt *cells = NULL, *faces = NULL;
1262: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
1263: PetscBT cellFlip; // The bit is true if a cell should have its orientation reversed
1264: PetscInt *cellComp; // The connected component number of each cell
1265: PetscInt *faceComp; // The connected component number of each face
1266: PetscInt Ncomp = 0; // The number of local connected components
1267: SharedFace *localFace; // Holds local information for owned and ghost shared faces
1268: SharedFace *remoteFace; // Holds remote information from owners of shared faces
1269: PetscInt *Nneigh; // The number of neighboring processes for each component
1270: Neighbor **neighbors; // The neighbor for each component
1271: PetscSFNode *neighborAdj; // The (rank, comp) of each neighbor
1272: PetscBool *neighborVal; // Whether neighbors currently match
1273: MPI_Comm comm;
1274: PetscMPIInt rank, size;
1276: PetscFunctionBegin;
1277: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1278: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1279: PetscCallMPI(MPI_Comm_size(comm, &size));
1280: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1281: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
1282: PetscCall(PetscBTCreate(cEnd - cStart, &cellFlip));
1283: PetscCall(PetscBTMemzero(cEnd - cStart, cellFlip));
1284: PetscCall(PetscCalloc2(cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
1285: // Phase 1: Serial Orientation
1286: PetscCall(DMPlexOrient_Serial(dm, cellIS, faceIS, &Ncomp, cellComp, faceComp, cellFlip));
1287: if (debug) {
1288: PetscViewer v;
1289: PetscInt cdepth = -1;
1291: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1292: PetscCall(PetscViewerASCIIPushSynchronized(v));
1293: if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
1294: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]New Orientation %" PetscInt_FMT " cells (depth %" PetscInt_FMT ") and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, cdepth, fEnd - fStart));
1295: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
1296: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
1297: PetscCall(PetscViewerFlush(v));
1298: PetscCall(PetscViewerASCIIPopSynchronized(v));
1299: }
1300: if (size == 1) goto end;
1301: // Phase 2
1302: PetscCall(DMPlexOrientCreateSharedFaces_Internal(dm, cellIS, faceIS, Ncomp, faceComp, cellFlip, &localFace, &remoteFace));
1303: // Phase 3
1304: PetscCall(DMPlexOrientCreateNeighbors_Internal(dm, cellIS, faceIS, Ncomp, faceComp, localFace, remoteFace, &Nneigh, &neighbors));
1305: PetscCall(PetscFree2(localFace, remoteFace));
1306: // Phase 4
1307: PetscCall(DMPlexOrientCreateProcessGraph_Internal(dm, faceIS, Ncomp, Nneigh, neighbors, &neighborAdj, &neighborVal));
1308: // Phase 5
1309: PetscCall(DMPlexOrientSolveProcessGraph_Internal(dm, cellIS, Ncomp, cellComp, Nneigh, neighborAdj, neighborVal, cellFlip));
1310: if (debug) {
1311: PetscViewer v;
1313: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1314: PetscCall(PetscViewerASCIIPushSynchronized(v));
1315: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
1316: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
1317: PetscCall(PetscViewerFlush(v));
1318: PetscCall(PetscViewerASCIIPopSynchronized(v));
1319: }
1320: for (PetscInt c = 0; c < Ncomp; ++c) PetscCall(PetscFree(neighbors[c]));
1321: PetscCall(PetscFree2(Nneigh, neighbors));
1322: PetscCall(PetscFree2(neighborAdj, neighborVal));
1323: end:
1324: // Phase 6
1325: PetscCall(DMPlexOrientPlex_Internal(dm, cellIS, cellFlip));
1326: PetscCall(PetscBTDestroy(&cellFlip));
1327: PetscCall(PetscFree2(cellComp, faceComp));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: static PetscErrorCode DMPlexCheckOrientation_Internal(DM dm, IS cellIS, IS faceIS)
1332: {
1333: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
1334: PetscViewer viewer = NULL, selfviewer = NULL;
1335: PetscSF sf;
1336: const PetscInt *lpoints, *rootdegree = NULL;
1337: const PetscSFNode *rpoints;
1338: const PetscInt *cells = NULL, *faces = NULL;
1339: PetscSFNode *oornt, *rornt, *lornt;
1340: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
1341: PetscInt pdepth, Nr, Nl;
1342: PetscBool faceIsVertex = PETSC_FALSE, valid = PETSC_TRUE;
1343: MPI_Comm comm;
1344: PetscMPIInt size, rank;
1346: PetscFunctionBegin;
1347: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1348: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1349: PetscCallMPI(MPI_Comm_size(comm, &size));
1350: if (debug) {
1351: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
1352: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1353: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
1354: }
1356: PetscCall(DMGetPointSF(dm, &sf));
1357: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &lpoints, &rpoints));
1358: if (Nr >= 0) {
1359: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
1360: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
1361: } else {
1362: sf = NULL;
1363: Nr = 0;
1364: Nl = 0;
1365: }
1366: PetscCall(PetscCalloc3(Nr, &oornt, Nr, &rornt, Nl, &lornt));
1367: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
1368: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
1369: PetscCall(DMPlexGetPointDepth(dm, faces ? faces[fStart] : fStart, &pdepth));
1370: if (!pdepth) faceIsVertex = PETSC_TRUE;
1371: if (debug) {
1372: PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Checking orientation of %" PetscInt_FMT " cells and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, fEnd - fStart));
1373: PetscCall(PetscViewerASCIIPushTab(selfviewer));
1374: }
1375: for (PetscInt f = fStart; f < fEnd; ++f) {
1376: const PetscInt face = faces ? faces[f] : f;
1377: const PetscBool owner = rootdegree && rootdegree[face] ? PETSC_TRUE : PETSC_FALSE;
1378: const PetscInt *supp;
1379: PetscInt neighbors[2], o[2] = {0, 0};
1380: PetscInt lsS = 0, sS;
1382: PetscCall(DMPlexGetSupport(dm, face, &supp));
1383: PetscCall(DMPlexGetSupportSize(dm, face, &sS));
1384: // Filter support for local cells
1385: for (PetscInt s = 0; s < sS; ++s) {
1386: PetscInt ind;
1388: ind = GetPointIndex(supp[s], cStart, cEnd, cells);
1389: if (ind >= 0) {
1390: neighbors[PetscMin(lsS, 1)] = supp[s];
1391: ++lsS;
1392: }
1393: }
1394: PetscCheck(lsS < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " > 2", face, lsS);
1395: // Extract orientations
1396: for (PetscInt s = 0; s < lsS; ++s) {
1397: const PetscInt *cone, *ornt;
1398: PetscInt cS;
1400: PetscCall(DMPlexGetConeSize(dm, neighbors[s], &cS));
1401: PetscCall(DMPlexGetOrientedCone(dm, neighbors[s], &cone, &ornt));
1402: for (PetscInt c = 0; c < cS; ++c) {
1403: if (cone[c] == face) {
1404: if (faceIsVertex) o[s] = c * 2 - 1;
1405: else o[s] = ornt[c] < 0 ? -1 : 1;
1406: break;
1407: }
1408: }
1409: PetscCall(DMPlexRestoreOrientedCone(dm, neighbors[s], &cone, &ornt));
1410: }
1411: if (lsS == 2) {
1412: // Check internal face
1413: if (o[0] * o[1] >= 0) {
1414: valid = PETSC_FALSE;
1415: if (debug)
1416: 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]));
1417: } else if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Internal Face %" PetscInt_FMT " valid\n", rank, face));
1418: } else {
1419: // Check shared and boundary faces
1420: PetscInt l;
1422: PetscCall(PetscFindInt(face, Nl, lpoints, &l));
1423: // Boundary face
1424: if (l < 0 && !owner) {
1425: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Boundary Face %" PetscInt_FMT " valid\n", rank, face));
1426: continue;
1427: }
1428: if (l >= 0) {
1429: lornt[l].index = neighbors[0];
1430: lornt[l].rank = o[0];
1431: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Ghost Face %" PetscInt_FMT " (%" PetscInt_FMT ") was stored\n", rank, face, o[0]));
1432: }
1433: }
1434: if (owner) {
1435: // Store shared face orientation and cell from owner
1436: oornt[face].index = neighbors[0];
1437: oornt[face].rank = o[0];
1438: if (debug > 1) PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]Owned Face %" PetscInt_FMT " (%" PetscInt_FMT ") was stored\n", rank, face, o[0]));
1439: }
1440: }
1441: // Communicate shared face orientations from owner
1442: if (sf) {
1443: PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, oornt, rornt, MPI_REPLACE));
1444: PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, oornt, rornt, MPI_REPLACE));
1445: }
1446: // Check unowned shared faces
1447: for (PetscInt l = 0; l < Nl; ++l) {
1448: const PetscInt face = lpoints ? lpoints[l] : l;
1449: PetscBool flip = PETSC_FALSE;
1450: PetscInt ind, cl, o[2] = {0, 0};
1452: // Filter for local faces
1453: ind = GetPointIndex(face, fStart, fEnd, faces);
1454: if (ind < 0) continue;
1455: // Check for shared cell
1456: PetscCall(PetscFindInt(lornt[l].index, Nl, lpoints, &cl));
1457: if (cl >= 0) {
1458: if (rpoints[cl].index == rornt[face].index) {
1459: flip = PETSC_TRUE;
1460: 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));
1461: } else {
1462: if (debug > 1)
1463: 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,
1464: rpoints[cl].rank, rornt[face].index));
1465: }
1466: }
1467: o[0] = lornt[l].rank;
1468: o[1] = rornt[face].rank;
1469: // This was not a shared face
1470: if (!o[0]) continue;
1471: if ((o[0] * o[1] >= 0 && !flip) || (o[0] * o[1] < 0 && flip)) {
1472: valid = PETSC_FALSE;
1473: if (debug)
1474: 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]));
1475: } 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));
1476: }
1477: // Cleanup
1478: PetscCall(PetscFree3(oornt, rornt, lornt));
1479: if (debug) {
1480: PetscCall(PetscViewerASCIIPopTab(selfviewer));
1481: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
1482: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1483: }
1484: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &valid, 1, MPI_C_BOOL, MPI_LAND, comm));
1485: PetscCheck(valid, comm, PETSC_ERR_ARG_WRONGSTATE, "Mesh was not properly oriented");
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: /*@
1490: DMPlexCheckOrientationLabel - Check that the surface defined by the given `DMLabel` is oriented
1492: Collective
1494: Input Parameters:
1495: + dm - The `DM`
1496: - label - The `DMLabel` defining the embedded surface
1498: Level: advanced
1500: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`
1501: @*/
1502: PetscErrorCode DMPlexCheckOrientationLabel(DM dm, DMLabel label)
1503: {
1504: IS cellIS, faceIS;
1506: PetscFunctionBegin;
1507: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
1508: PetscCall(DMPlexCheckOrientation_Internal(dm, cellIS, faceIS));
1509: PetscCall(ISDestroy(&cellIS));
1510: PetscCall(ISDestroy(&faceIS));
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }