Actual source code: plexorient.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: /*@
5: DMPlexOrientPoint - Act with the given orientation on the cone points of this mesh point, and update its use in the mesh.
7: Not Collective
9: Input Parameters:
10: + dm - The `DM`
11: . p - The mesh point
12: - o - The orientation
14: Level: intermediate
16: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`, `DMPlexGetCone()`, `DMPlexGetConeOrientation()`, `DMPlexInterpolate()`, `DMPlexGetChart()`
17: @*/
18: PetscErrorCode DMPlexOrientPoint(DM dm, PetscInt p, PetscInt o)
19: {
20: DMPolytopeType ct;
21: const PetscInt *arr, *cone, *ornt, *support;
22: PetscInt *newcone, *newornt;
23: PetscInt coneSize, c, supportSize, s;
25: PetscFunctionBegin;
27: PetscCall(DMPlexGetCellType(dm, p, &ct));
28: arr = DMPolytopeTypeGetArrangment(ct, o);
29: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
30: PetscCall(DMPlexGetCone(dm, p, &cone));
31: PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
32: PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newcone));
33: PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newornt));
34: for (c = 0; c < coneSize; ++c) {
35: DMPolytopeType ft;
36: PetscInt nO;
38: PetscCall(DMPlexGetCellType(dm, cone[c], &ft));
39: nO = DMPolytopeTypeGetNumArrangments(ft) / 2;
40: newcone[c] = cone[arr[c * 2 + 0]];
41: newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], ornt[arr[c * 2 + 0]]);
42: 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]);
43: }
44: PetscCall(DMPlexSetCone(dm, p, newcone));
45: PetscCall(DMPlexSetConeOrientation(dm, p, newornt));
46: PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newcone));
47: PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newornt));
48: /* Update orientation of this point in the support points */
49: PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
50: PetscCall(DMPlexGetSupport(dm, p, &support));
51: for (s = 0; s < supportSize; ++s) {
52: PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize));
53: PetscCall(DMPlexGetCone(dm, support[s], &cone));
54: PetscCall(DMPlexGetConeOrientation(dm, support[s], &ornt));
55: for (c = 0; c < coneSize; ++c) {
56: PetscInt po;
58: if (cone[c] != p) continue;
59: /* ornt[c] * 0 = target = po * o so that po = ornt[c] * o^{-1} */
60: po = DMPolytopeTypeComposeOrientationInv(ct, ornt[c], o);
61: PetscCall(DMPlexInsertConeOrientation(dm, support[s], c, po));
62: }
63: }
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*
68: - Checks face match
69: - Flips non-matching
70: - Inserts faces of support cells in FIFO
71: */
72: static PetscErrorCode DMPlexCheckFace_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, PetscInt cStart, PetscInt fStart, PetscInt fEnd, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces)
73: {
74: const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
75: PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
76: PetscInt face, dim, seenA, flippedA, seenB, flippedB, mismatch, c;
78: PetscFunctionBegin;
79: face = faceFIFO[(*fTop)++];
80: PetscCall(DMGetDimension(dm, &dim));
81: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
82: PetscCall(DMPlexGetSupport(dm, face, &support));
83: if (supportSize < 2) PetscFunctionReturn(PETSC_SUCCESS);
84: PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, supportSize);
85: seenA = PetscBTLookup(seenCells, support[0] - cStart);
86: flippedA = PetscBTLookup(flippedCells, support[0] - cStart) ? 1 : 0;
87: seenB = PetscBTLookup(seenCells, support[1] - cStart);
88: flippedB = PetscBTLookup(flippedCells, support[1] - cStart) ? 1 : 0;
90: PetscCall(DMPlexGetConeSize(dm, support[0], &coneSizeA));
91: PetscCall(DMPlexGetConeSize(dm, support[1], &coneSizeB));
92: PetscCall(DMPlexGetCone(dm, support[0], &coneA));
93: PetscCall(DMPlexGetCone(dm, support[1], &coneB));
94: PetscCall(DMPlexGetConeOrientation(dm, support[0], &coneOA));
95: PetscCall(DMPlexGetConeOrientation(dm, support[1], &coneOB));
96: for (c = 0; c < coneSizeA; ++c) {
97: if (!PetscBTLookup(seenFaces, coneA[c] - fStart)) {
98: faceFIFO[(*fBottom)++] = coneA[c];
99: PetscCall(PetscBTSet(seenFaces, coneA[c] - fStart));
100: }
101: if (coneA[c] == face) posA = c;
102: 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);
103: }
104: PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[0]);
105: for (c = 0; c < coneSizeB; ++c) {
106: if (!PetscBTLookup(seenFaces, coneB[c] - fStart)) {
107: faceFIFO[(*fBottom)++] = coneB[c];
108: PetscCall(PetscBTSet(seenFaces, coneB[c] - fStart));
109: }
110: if (coneB[c] == face) posB = c;
111: 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);
112: }
113: PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[1]);
115: if (dim == 1) {
116: mismatch = posA == posB;
117: } else {
118: mismatch = coneOA[posA] == coneOB[posB];
119: }
121: if (mismatch ^ (flippedA ^ flippedB)) {
122: 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]);
123: if (!seenA && !flippedA) {
124: PetscCall(PetscBTSet(flippedCells, support[0] - cStart));
125: } else if (!seenB && !flippedB) {
126: PetscCall(PetscBTSet(flippedCells, support[1] - cStart));
127: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
128: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
129: PetscCall(PetscBTSet(seenCells, support[0] - cStart));
130: PetscCall(PetscBTSet(seenCells, support[1] - cStart));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: DMPlexOrient - Give a consistent orientation to the input mesh
137: Input Parameter:
138: . dm - The `DM`
140: Note:
141: The orientation data for the `DM` are change in-place.
143: This routine will fail for non-orientable surfaces, such as the Moebius strip.
145: Level: advanced
147: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
148: @*/
149: PetscErrorCode DMPlexOrient(DM dm)
150: {
151: MPI_Comm comm;
152: PetscSF sf;
153: const PetscInt *lpoints;
154: const PetscSFNode *rpoints;
155: PetscSFNode *rorntComp = NULL, *lorntComp = NULL;
156: PetscInt *numNeighbors, **neighbors, *locSupport = NULL;
157: PetscSFNode *nrankComp;
158: PetscBool *match, *flipped;
159: PetscBT seenCells, flippedCells, seenFaces;
160: PetscInt *faceFIFO, fTop, fBottom, *cellComp, *faceComp;
161: PetscInt numLeaves, numRoots, dim, h, cStart, cEnd, c, cell, fStart, fEnd, face, off, totNeighbors = 0;
162: PetscMPIInt rank, size, numComponents, comp = 0;
163: PetscBool flg, flg2;
164: PetscViewer viewer = NULL, selfviewer = NULL;
166: PetscFunctionBegin;
167: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
168: PetscCallMPI(MPI_Comm_rank(comm, &rank));
169: PetscCallMPI(MPI_Comm_size(comm, &size));
170: PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &flg));
171: PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &flg2));
172: PetscCall(DMGetPointSF(dm, &sf));
173: PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints));
174: /* Truth Table
175: mismatch flips do action mismatch flipA ^ flipB action
176: F 0 flips no F F F
177: F 1 flip yes F T T
178: F 2 flips no T F T
179: T 0 flips yes T T F
180: T 1 flip no
181: T 2 flips yes
182: */
183: PetscCall(DMGetDimension(dm, &dim));
184: PetscCall(DMPlexGetVTKCellHeight(dm, &h));
185: PetscCall(DMPlexGetHeightStratum(dm, h, &cStart, &cEnd));
186: PetscCall(DMPlexGetHeightStratum(dm, h + 1, &fStart, &fEnd));
187: PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
188: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
189: PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells));
190: PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells));
191: PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
192: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
193: PetscCall(PetscCalloc3(fEnd - fStart, &faceFIFO, cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
194: /*
195: OLD STYLE
196: - Add an integer array over cells and faces (component) for connected component number
197: Foreach component
198: - Mark the initial cell as seen
199: - Process component as usual
200: - Set component for all seenCells
201: - Wipe seenCells and seenFaces (flippedCells can stay)
202: - Generate parallel adjacency for component using SF and seenFaces
203: - Collect numComponents adj data from each proc to 0
204: - Build same serial graph
205: - Use same solver
206: - Use Scatterv to to send back flipped flags for each component
207: - Negate flippedCells by component
209: NEW STYLE
210: - Create the adj on each process
211: - Bootstrap to complete graph on proc 0
212: */
213: /* Loop over components */
214: for (cell = cStart; cell < cEnd; ++cell) cellComp[cell - cStart] = -1;
215: do {
216: /* Look for first unmarked cell */
217: for (cell = cStart; cell < cEnd; ++cell)
218: if (cellComp[cell - cStart] < 0) break;
219: if (cell >= cEnd) break;
220: /* Initialize FIFO with first cell in component */
221: {
222: const PetscInt *cone;
223: PetscInt coneSize;
225: fTop = fBottom = 0;
226: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
227: PetscCall(DMPlexGetCone(dm, cell, &cone));
228: for (c = 0; c < coneSize; ++c) {
229: faceFIFO[fBottom++] = cone[c];
230: PetscCall(PetscBTSet(seenFaces, cone[c] - fStart));
231: }
232: PetscCall(PetscBTSet(seenCells, cell - cStart));
233: }
234: /* Consider each face in FIFO */
235: while (fTop < fBottom) PetscCall(DMPlexCheckFace_Internal(dm, faceFIFO, &fTop, &fBottom, cStart, fStart, fEnd, seenCells, flippedCells, seenFaces));
236: /* Set component for cells and faces */
237: for (cell = 0; cell < cEnd - cStart; ++cell) {
238: if (PetscBTLookup(seenCells, cell)) cellComp[cell] = comp;
239: }
240: for (face = 0; face < fEnd - fStart; ++face) {
241: if (PetscBTLookup(seenFaces, face)) faceComp[face] = comp;
242: }
243: /* Wipe seenCells and seenFaces for next component */
244: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
245: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
246: ++comp;
247: } while (1);
248: numComponents = comp;
249: if (flg) {
250: PetscViewer v;
252: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
253: PetscCall(PetscViewerASCIIPushSynchronized(v));
254: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
255: PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
256: PetscCall(PetscViewerFlush(v));
257: PetscCall(PetscViewerASCIIPopSynchronized(v));
258: }
259: /* Now all subdomains are oriented, but we need a consistent parallel orientation */
260: if (numLeaves >= 0) {
261: PetscInt maxSupportSize, neighbor;
263: /* Store orientations of boundary faces*/
264: PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSupportSize));
265: PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSupportSize, &locSupport));
266: for (face = fStart; face < fEnd; ++face) {
267: const PetscInt *cone, *support, *ornt;
268: PetscInt coneSize, supportSize, Ns = 0, s, l;
270: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
271: /* Ignore overlapping cells */
272: PetscCall(DMPlexGetSupport(dm, face, &support));
273: for (s = 0; s < supportSize; ++s) {
274: PetscCall(PetscFindInt(support[s], numLeaves, lpoints, &l));
275: if (l >= 0) continue;
276: locSupport[Ns++] = support[s];
277: }
278: if (Ns != 1) continue;
279: neighbor = locSupport[0];
280: PetscCall(DMPlexGetCone(dm, neighbor, &cone));
281: PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize));
282: PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt));
283: for (c = 0; c < coneSize; ++c)
284: if (cone[c] == face) break;
285: if (dim == 1) {
286: /* Use cone position instead, shifted to -1 or 1 */
287: if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = 1 - c * 2;
288: else rorntComp[face].rank = c * 2 - 1;
289: } else {
290: if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1;
291: else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1;
292: }
293: rorntComp[face].index = faceComp[face - fStart];
294: }
295: /* Communicate boundary edge orientations */
296: PetscCall(PetscSFBcastBegin(sf, MPIU_2INT, rorntComp, lorntComp, MPI_REPLACE));
297: PetscCall(PetscSFBcastEnd(sf, MPIU_2INT, rorntComp, lorntComp, MPI_REPLACE));
298: }
299: /* Get process adjacency */
300: PetscCall(PetscMalloc2(numComponents, &numNeighbors, numComponents, &neighbors));
301: viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
302: if (flg2) PetscCall(PetscViewerASCIIPushSynchronized(viewer));
303: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
304: for (comp = 0; comp < numComponents; ++comp) {
305: PetscInt l, n;
307: numNeighbors[comp] = 0;
308: PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp]));
309: /* I know this is p^2 time in general, but for bounded degree its alright */
310: for (l = 0; l < numLeaves; ++l) {
311: const PetscInt face = lpoints[l];
313: /* Find a representative face (edge) separating pairs of procs */
314: if ((face >= fStart) && (face < fEnd) && (faceComp[face - fStart] == comp) && rorntComp[face].rank) {
315: const PetscInt rrank = rpoints[l].rank;
316: const PetscInt rcomp = lorntComp[face].index;
318: for (n = 0; n < numNeighbors[comp]; ++n)
319: if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break;
320: if (n >= numNeighbors[comp]) {
321: PetscInt supportSize;
323: PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
324: PetscCheck(supportSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary faces should see one cell, not %" PetscInt_FMT, supportSize);
325: if (flg)
326: 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,
327: rpoints[l].index, rrank, rcomp, lorntComp[face].rank));
328: neighbors[comp][numNeighbors[comp]++] = l;
329: }
330: }
331: }
332: totNeighbors += numNeighbors[comp];
333: }
334: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
335: PetscCall(PetscViewerFlush(viewer));
336: if (flg2) PetscCall(PetscViewerASCIIPopSynchronized(viewer));
337: PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match));
338: for (comp = 0, off = 0; comp < numComponents; ++comp) {
339: PetscInt n;
341: for (n = 0; n < numNeighbors[comp]; ++n, ++off) {
342: const PetscInt face = lpoints[neighbors[comp][n]];
343: const PetscInt o = rorntComp[face].rank * lorntComp[face].rank;
345: if (o < 0) match[off] = PETSC_TRUE;
346: else if (o > 0) match[off] = PETSC_FALSE;
347: 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);
348: nrankComp[off].rank = rpoints[neighbors[comp][n]].rank;
349: nrankComp[off].index = lorntComp[lpoints[neighbors[comp][n]]].index;
350: }
351: PetscCall(PetscFree(neighbors[comp]));
352: }
353: /* Collect the graph on 0 */
354: if (numLeaves >= 0) {
355: Mat G;
356: PetscBT seenProcs, flippedProcs;
357: PetscInt *procFIFO, pTop, pBottom;
358: PetscInt *N = NULL, *Noff;
359: PetscSFNode *adj = NULL;
360: PetscBool *val = NULL;
361: PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc, p, o;
362: PetscMPIInt size = 0;
364: PetscCall(PetscCalloc1(numComponents, &flipped));
365: if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
366: PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
367: PetscCallMPI(MPI_Gather(&numComponents, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
368: for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
369: if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
370: PetscCallMPI(MPI_Gatherv(numNeighbors, numComponents, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
371: for (p = 0, o = 0; p < size; ++p) {
372: recvcounts[p] = 0;
373: for (c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
374: displs[p + 1] = displs[p] + recvcounts[p];
375: }
376: if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
377: PetscCallMPI(MPI_Gatherv(nrankComp, totNeighbors, MPIU_2INT, adj, recvcounts, displs, MPIU_2INT, 0, comm));
378: PetscCallMPI(MPI_Gatherv(match, totNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm));
379: PetscCall(PetscFree2(numNeighbors, neighbors));
380: if (rank == 0) {
381: for (p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
382: if (flg) {
383: PetscInt n;
385: for (p = 0, off = 0; p < size; ++p) {
386: for (c = 0; c < Nc[p]; ++c) {
387: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %d Comp %" PetscInt_FMT ":\n", p, c));
388: 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]]));
389: }
390: }
391: }
392: /* Symmetrize the graph */
393: PetscCall(MatCreate(PETSC_COMM_SELF, &G));
394: PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
395: PetscCall(MatSetUp(G));
396: for (p = 0, off = 0; p < size; ++p) {
397: for (c = 0; c < Nc[p]; ++c) {
398: const PetscInt r = Noff[p] + c;
399: PetscInt n;
401: for (n = 0; n < N[r]; ++n, ++off) {
402: const PetscInt q = Noff[adj[off].rank] + adj[off].index;
403: const PetscScalar o = val[off] ? 1.0 : 0.0;
405: PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
406: PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
407: }
408: }
409: }
410: PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
411: PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
413: PetscCall(PetscBTCreate(Noff[size], &seenProcs));
414: PetscCall(PetscBTMemzero(Noff[size], seenProcs));
415: PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
416: PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
417: PetscCall(PetscMalloc1(Noff[size], &procFIFO));
418: pTop = pBottom = 0;
419: for (p = 0; p < Noff[size]; ++p) {
420: if (PetscBTLookup(seenProcs, p)) continue;
421: /* Initialize FIFO with next proc */
422: procFIFO[pBottom++] = p;
423: PetscCall(PetscBTSet(seenProcs, p));
424: /* Consider each proc in FIFO */
425: while (pTop < pBottom) {
426: const PetscScalar *ornt;
427: const PetscInt *neighbors;
428: PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors, n;
430: proc = procFIFO[pTop++];
431: flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
432: PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
433: /* Loop over neighboring procs */
434: for (n = 0; n < numNeighbors; ++n) {
435: nproc = neighbors[n];
436: mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
437: seen = PetscBTLookup(seenProcs, nproc);
438: flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
440: if (mismatch ^ (flippedA ^ flippedB)) {
441: 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);
442: if (!flippedB) {
443: PetscCall(PetscBTSet(flippedProcs, nproc));
444: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
445: } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
446: if (!seen) {
447: procFIFO[pBottom++] = nproc;
448: PetscCall(PetscBTSet(seenProcs, nproc));
449: }
450: }
451: }
452: }
453: PetscCall(PetscFree(procFIFO));
454: PetscCall(MatDestroy(&G));
455: PetscCall(PetscFree2(adj, val));
456: PetscCall(PetscBTDestroy(&seenProcs));
457: }
458: /* Scatter flip flags */
459: {
460: PetscBool *flips = NULL;
462: if (rank == 0) {
463: PetscCall(PetscMalloc1(Noff[size], &flips));
464: for (p = 0; p < Noff[size]; ++p) {
465: flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
466: if (flg && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %d:\n", p));
467: }
468: for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
469: }
470: PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPIU_BOOL, flipped, numComponents, MPIU_BOOL, 0, comm));
471: PetscCall(PetscFree(flips));
472: }
473: if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
474: PetscCall(PetscFree(N));
475: PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
476: PetscCall(PetscFree2(nrankComp, match));
478: /* Decide whether to flip cells in each component */
479: for (c = 0; c < cEnd - cStart; ++c) {
480: if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c));
481: }
482: PetscCall(PetscFree(flipped));
483: }
484: if (flg) {
485: PetscViewer v;
487: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
488: PetscCall(PetscViewerASCIIPushSynchronized(v));
489: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
490: PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
491: PetscCall(PetscViewerFlush(v));
492: PetscCall(PetscViewerASCIIPopSynchronized(v));
493: }
494: /* Reverse flipped cells in the mesh */
495: for (c = cStart; c < cEnd; ++c) {
496: if (PetscBTLookup(flippedCells, c - cStart)) PetscCall(DMPlexOrientPoint(dm, c, -1));
497: }
498: PetscCall(PetscBTDestroy(&seenCells));
499: PetscCall(PetscBTDestroy(&flippedCells));
500: PetscCall(PetscBTDestroy(&seenFaces));
501: PetscCall(PetscFree2(numNeighbors, neighbors));
502: PetscCall(PetscFree3(rorntComp, lorntComp, locSupport));
503: PetscCall(PetscFree3(faceFIFO, cellComp, faceComp));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }