1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petscsf.h>
3
4 /*@
5 DMPlexOrientPoint - Act with the given orientation on the cone points of this mesh point, and update its use in the mesh.
6
7 Not Collective
8
9 Input Parameters:
10 + dm - The `DM`
11 . p - The mesh point
12 - o - The orientation
13
14 Level: intermediate
15
16 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`, `DMPlexGetCone()`, `DMPlexGetConeOrientation()`, `DMPlexInterpolate()`, `DMPlexGetChart()`
17 @*/
DMPlexOrientPoint(DM dm,PetscInt p,PetscInt o)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;
24
25 PetscFunctionBegin;
26 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
27 PetscCall(DMPlexGetCellType(dm, p, &ct));
28 arr = DMPolytopeTypeGetArrangement(ct, o);
29 if (!arr) PetscFunctionReturn(PETSC_SUCCESS);
30 PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
31 PetscCall(DMPlexGetCone(dm, p, &cone));
32 PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
33 PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newcone));
34 PetscCall(DMGetWorkArray(dm, coneSize, MPIU_INT, &newornt));
35 for (c = 0; c < coneSize; ++c) {
36 DMPolytopeType ft;
37 PetscInt nO;
38
39 PetscCall(DMPlexGetCellType(dm, cone[c], &ft));
40 nO = DMPolytopeTypeGetNumArrangements(ft) / 2;
41 newcone[c] = cone[arr[c * 2 + 0]];
42 newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], ornt[arr[c * 2 + 0]]);
43 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]);
44 }
45 PetscCall(DMPlexSetCone(dm, p, newcone));
46 PetscCall(DMPlexSetConeOrientation(dm, p, newornt));
47 PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newcone));
48 PetscCall(DMRestoreWorkArray(dm, coneSize, MPIU_INT, &newornt));
49 /* Update orientation of this point in the support points */
50 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize));
51 PetscCall(DMPlexGetSupport(dm, p, &support));
52 for (s = 0; s < supportSize; ++s) {
53 PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize));
54 PetscCall(DMPlexGetCone(dm, support[s], &cone));
55 PetscCall(DMPlexGetConeOrientation(dm, support[s], &ornt));
56 for (c = 0; c < coneSize; ++c) {
57 PetscInt po;
58
59 if (cone[c] != p) continue;
60 /* ornt[c] * 0 = target = po * o so that po = ornt[c] * o^{-1} */
61 po = DMPolytopeTypeComposeOrientationInv(ct, ornt[c], o);
62 PetscCall(DMPlexInsertConeOrientation(dm, support[s], c, po));
63 }
64 }
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
GetPointIndex(PetscInt point,PetscInt pStart,PetscInt pEnd,const PetscInt points[])68 static PetscInt GetPointIndex(PetscInt point, PetscInt pStart, PetscInt pEnd, const PetscInt points[])
69 {
70 if (points) {
71 PetscInt loc;
72
73 PetscCallAbort(PETSC_COMM_SELF, PetscFindInt(point, pEnd - pStart, points, &loc));
74 if (loc >= 0) return loc;
75 } else {
76 if (point >= pStart && point < pEnd) return point - pStart;
77 }
78 return -1;
79 }
80
81 /*
82 - Checks face match
83 - Flips non-matching
84 - Inserts faces of support cells in FIFO
85 */
DMPlexCheckFace_Internal(DM dm,PetscInt * faceFIFO,PetscInt * fTop,PetscInt * fBottom,IS cellIS,IS faceIS,PetscBT seenCells,PetscBT flippedCells,PetscBT seenFaces)86 static PetscErrorCode DMPlexCheckFace_Internal(DM dm, PetscInt *faceFIFO, PetscInt *fTop, PetscInt *fBottom, IS cellIS, IS faceIS, PetscBT seenCells, PetscBT flippedCells, PetscBT seenFaces)
87 {
88 const PetscInt *supp, *coneA, *coneB, *coneOA, *coneOB;
89 PetscInt suppSize, Ns = 0, coneSizeA, coneSizeB, posA = -1, posB = -1;
90 PetscInt face, dim, indC[3], indS[3], seenA, flippedA, seenB, flippedB, mismatch;
91 const PetscInt *cells, *faces;
92 PetscInt cStart, cEnd, fStart, fEnd;
93
94 PetscFunctionBegin;
95 face = faceFIFO[(*fTop)++];
96 PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
97 PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
98 PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &dim));
99 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
100 PetscCall(DMPlexGetSupport(dm, face, &supp));
101 // Filter the support
102 for (PetscInt s = 0; s < suppSize; ++s) {
103 // Filter support
104 indC[Ns] = GetPointIndex(supp[s], cStart, cEnd, cells);
105 indS[Ns] = s;
106 if (indC[Ns] >= 0) ++Ns;
107 }
108 if (Ns < 2) PetscFunctionReturn(PETSC_SUCCESS);
109 PetscCheck(Ns == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, Ns);
110 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]);
111 seenA = PetscBTLookup(seenCells, indC[0]);
112 flippedA = PetscBTLookup(flippedCells, indC[0]) ? 1 : 0;
113 seenB = PetscBTLookup(seenCells, indC[1]);
114 flippedB = PetscBTLookup(flippedCells, indC[1]) ? 1 : 0;
115
116 PetscCall(DMPlexGetConeSize(dm, supp[indS[0]], &coneSizeA));
117 PetscCall(DMPlexGetConeSize(dm, supp[indS[1]], &coneSizeB));
118 PetscCall(DMPlexGetCone(dm, supp[indS[0]], &coneA));
119 PetscCall(DMPlexGetCone(dm, supp[indS[1]], &coneB));
120 PetscCall(DMPlexGetConeOrientation(dm, supp[indS[0]], &coneOA));
121 PetscCall(DMPlexGetConeOrientation(dm, supp[indS[1]], &coneOB));
122 for (PetscInt c = 0; c < coneSizeA; ++c) {
123 const PetscInt indF = GetPointIndex(coneA[c], fStart, fEnd, faces);
124
125 // Filter cone
126 if (indF < 0) continue;
127 if (!PetscBTLookup(seenFaces, indF)) {
128 faceFIFO[(*fBottom)++] = coneA[c];
129 PetscCall(PetscBTSet(seenFaces, indF));
130 }
131 if (coneA[c] == face) posA = c;
132 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);
133 }
134 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]]);
135 for (PetscInt c = 0; c < coneSizeB; ++c) {
136 const PetscInt indF = GetPointIndex(coneB[c], fStart, fEnd, faces);
137
138 // Filter cone
139 if (indF < 0) continue;
140 if (!PetscBTLookup(seenFaces, indF)) {
141 faceFIFO[(*fBottom)++] = coneB[c];
142 PetscCall(PetscBTSet(seenFaces, indF));
143 }
144 if (coneB[c] == face) posB = c;
145 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);
146 }
147 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]]);
148
149 if (dim == 1) {
150 mismatch = posA == posB;
151 } else {
152 mismatch = coneOA[posA] == coneOB[posB];
153 }
154
155 if (mismatch ^ (flippedA ^ flippedB)) {
156 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]]);
157 if (!seenA && !flippedA) {
158 PetscCall(PetscBTSet(flippedCells, indC[0]));
159 } else if (!seenB && !flippedB) {
160 PetscCall(PetscBTSet(flippedCells, indC[1]));
161 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
162 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
163 PetscCall(PetscBTSet(seenCells, indC[0]));
164 PetscCall(PetscBTSet(seenCells, indC[1]));
165 PetscFunctionReturn(PETSC_SUCCESS);
166 }
167
DMPlexCheckFace_Old_Internal(DM dm,PetscInt * faceFIFO,PetscInt * fTop,PetscInt * fBottom,PetscInt cStart,PetscInt fStart,PetscInt fEnd,PetscBT seenCells,PetscBT flippedCells,PetscBT seenFaces)168 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)
169 {
170 const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
171 PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
172 PetscInt face, dim, seenA, flippedA, seenB, flippedB, mismatch, c;
173
174 PetscFunctionBegin;
175 face = faceFIFO[(*fTop)++];
176 PetscCall(DMGetDimension(dm, &dim));
177 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
178 PetscCall(DMPlexGetSupport(dm, face, &support));
179 if (supportSize < 2) PetscFunctionReturn(PETSC_SUCCESS);
180 PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %" PetscInt_FMT, supportSize);
181 seenA = PetscBTLookup(seenCells, support[0] - cStart);
182 flippedA = PetscBTLookup(flippedCells, support[0] - cStart) ? 1 : 0;
183 seenB = PetscBTLookup(seenCells, support[1] - cStart);
184 flippedB = PetscBTLookup(flippedCells, support[1] - cStart) ? 1 : 0;
185
186 PetscCall(DMPlexGetConeSize(dm, support[0], &coneSizeA));
187 PetscCall(DMPlexGetConeSize(dm, support[1], &coneSizeB));
188 PetscCall(DMPlexGetCone(dm, support[0], &coneA));
189 PetscCall(DMPlexGetCone(dm, support[1], &coneB));
190 PetscCall(DMPlexGetConeOrientation(dm, support[0], &coneOA));
191 PetscCall(DMPlexGetConeOrientation(dm, support[1], &coneOB));
192 for (c = 0; c < coneSizeA; ++c) {
193 if (!PetscBTLookup(seenFaces, coneA[c] - fStart)) {
194 faceFIFO[(*fBottom)++] = coneA[c];
195 PetscCall(PetscBTSet(seenFaces, coneA[c] - fStart));
196 }
197 if (coneA[c] == face) posA = c;
198 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);
199 }
200 PetscCheck(posA >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[0]);
201 for (c = 0; c < coneSizeB; ++c) {
202 if (!PetscBTLookup(seenFaces, coneB[c] - fStart)) {
203 faceFIFO[(*fBottom)++] = coneB[c];
204 PetscCall(PetscBTSet(seenFaces, coneB[c] - fStart));
205 }
206 if (coneB[c] == face) posB = c;
207 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);
208 }
209 PetscCheck(posB >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be located in cell %" PetscInt_FMT, face, support[1]);
210
211 if (dim == 1) {
212 mismatch = posA == posB;
213 } else {
214 mismatch = coneOA[posA] == coneOB[posB];
215 }
216
217 if (mismatch ^ (flippedA ^ flippedB)) {
218 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]);
219 if (!seenA && !flippedA) {
220 PetscCall(PetscBTSet(flippedCells, support[0] - cStart));
221 } else if (!seenB && !flippedB) {
222 PetscCall(PetscBTSet(flippedCells, support[1] - cStart));
223 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
224 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
225 PetscCall(PetscBTSet(seenCells, support[0] - cStart));
226 PetscCall(PetscBTSet(seenCells, support[1] - cStart));
227 PetscFunctionReturn(PETSC_SUCCESS);
228 }
229
230 /*
231 DMPlexOrient_Serial - Compute valid orientation for local connected components
232
233 Not collective
234
235 Input Parameters:
236 + dm - The `DM`
237 - cellHeight - The height of k-cells to be oriented
238
239 Output Parameters:
240 + Ncomp - The number of connected component
241 . cellComp - The connected component for each local cell
242 . faceComp - The connected component for each local face
243 - flippedCells - Marked cells should be inverted
244
245 Level: developer
246
247 .seealso: `DMPlexOrient()`
248 */
DMPlexOrient_Serial(DM dm,IS cellIS,IS faceIS,PetscInt * Ncomp,PetscInt cellComp[],PetscInt faceComp[],PetscBT flippedCells)249 static PetscErrorCode DMPlexOrient_Serial(DM dm, IS cellIS, IS faceIS, PetscInt *Ncomp, PetscInt cellComp[], PetscInt faceComp[], PetscBT flippedCells)
250 {
251 PetscBT seenCells, seenFaces;
252 PetscInt *faceFIFO;
253 const PetscInt *cells = NULL, *faces = NULL;
254 PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
255
256 PetscFunctionBegin;
257 if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
258 if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
259 PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
260 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
261 PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
262 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
263 PetscCall(PetscMalloc1(fEnd - fStart, &faceFIFO));
264 *Ncomp = 0;
265 for (PetscInt c = 0; c < cEnd - cStart; ++c) cellComp[c] = -1;
266 do {
267 PetscInt cc, fTop, fBottom;
268
269 // Look for first unmarked cell
270 for (cc = cStart; cc < cEnd; ++cc)
271 if (cellComp[cc - cStart] < 0) break;
272 if (cc >= cEnd) break;
273 // Initialize FIFO with first cell in component
274 {
275 const PetscInt cell = cells ? cells[cc] : cc;
276 const PetscInt *cone;
277 PetscInt coneSize;
278
279 fTop = fBottom = 0;
280 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
281 PetscCall(DMPlexGetCone(dm, cell, &cone));
282 for (PetscInt c = 0; c < coneSize; ++c) {
283 const PetscInt idx = GetPointIndex(cone[c], fStart, fEnd, faces);
284
285 // Cell faces are guaranteed to be in the face set
286 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);
287 faceFIFO[fBottom++] = cone[c];
288 PetscCall(PetscBTSet(seenFaces, idx));
289 }
290 PetscCall(PetscBTSet(seenCells, cc - cStart));
291 }
292 // Consider each face in FIFO
293 while (fTop < fBottom) PetscCall(DMPlexCheckFace_Internal(dm, faceFIFO, &fTop, &fBottom, cellIS, faceIS, seenCells, flippedCells, seenFaces));
294 // Set component for cells and faces
295 for (PetscInt c = 0; c < cEnd - cStart; ++c) {
296 if (PetscBTLookup(seenCells, c)) cellComp[c] = *Ncomp;
297 }
298 for (PetscInt f = 0; f < fEnd - fStart; ++f) {
299 if (PetscBTLookup(seenFaces, f)) faceComp[f] = *Ncomp;
300 }
301 // Wipe seenCells and seenFaces for next component
302 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
303 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
304 ++(*Ncomp);
305 } while (1);
306 PetscCall(PetscBTDestroy(&seenCells));
307 PetscCall(PetscBTDestroy(&seenFaces));
308 PetscCall(PetscFree(faceFIFO));
309 PetscFunctionReturn(PETSC_SUCCESS);
310 }
311
312 /*@
313 DMPlexOrient - Give a consistent orientation to the input mesh
314
315 Input Parameter:
316 . dm - The `DM`
317
318 Note:
319 The orientation data for the `DM` are change in-place.
320
321 This routine will fail for non-orientable surfaces, such as the Moebius strip.
322
323 Level: advanced
324
325 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
326 @*/
DMPlexOrient(DM dm)327 PetscErrorCode DMPlexOrient(DM dm)
328 {
329 #if 0
330 IS cellIS, faceIS;
331
332 PetscFunctionBegin;
333 PetscCall(DMPlexGetAllCells_Internal(dm, &cellIS));
334 PetscCall(DMPlexGetAllFaces_Internal(dm, &faceIS));
335 PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
336 PetscCall(ISDestroy(&cellIS));
337 PetscCall(ISDestroy(&faceIS));
338 PetscFunctionReturn(PETSC_SUCCESS);
339 #else
340 MPI_Comm comm;
341 PetscSF sf;
342 const PetscInt *lpoints;
343 const PetscSFNode *rpoints;
344 PetscSFNode *rorntComp = NULL, *lorntComp = NULL;
345 PetscInt *numNeighbors, **neighbors, *locSupport = NULL;
346 PetscSFNode *nrankComp;
347 PetscBool *match, *flipped;
348 PetscBT seenCells, flippedCells, seenFaces;
349 PetscInt *faceFIFO, fTop, fBottom, *cellComp, *faceComp;
350 PetscInt numLeaves, numRoots, dim, h, cStart, cEnd, c, cell, fStart, fEnd, face, off, totNeighbors = 0;
351 PetscMPIInt rank, size, numComponents, comp = 0;
352 PetscBool flg, flg2;
353 PetscViewer viewer = NULL, selfviewer = NULL;
354
355 PetscFunctionBegin;
356 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
357 PetscCallMPI(MPI_Comm_rank(comm, &rank));
358 PetscCallMPI(MPI_Comm_size(comm, &size));
359 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &flg));
360 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &flg2));
361 PetscCall(DMGetPointSF(dm, &sf));
362 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints));
363 /* Truth Table
364 mismatch flips do action mismatch flipA ^ flipB action
365 F 0 flips no F F F
366 F 1 flip yes F T T
367 F 2 flips no T F T
368 T 0 flips yes T T F
369 T 1 flip no
370 T 2 flips yes
371 */
372 PetscCall(DMGetDimension(dm, &dim));
373 PetscCall(DMPlexGetVTKCellHeight(dm, &h));
374 PetscCall(DMPlexGetHeightStratum(dm, h, &cStart, &cEnd));
375 PetscCall(DMPlexGetHeightStratum(dm, h + 1, &fStart, &fEnd));
376 PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
377 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
378 PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells));
379 PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells));
380 PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
381 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
382 PetscCall(PetscCalloc3(fEnd - fStart, &faceFIFO, cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
383 /*
384 OLD STYLE
385 - Add an integer array over cells and faces (component) for connected component number
386 Foreach component
387 - Mark the initial cell as seen
388 - Process component as usual
389 - Set component for all seenCells
390 - Wipe seenCells and seenFaces (flippedCells can stay)
391 - Generate parallel adjacency for component using SF and seenFaces
392 - Collect numComponents adj data from each proc to 0
393 - Build same serial graph
394 - Use same solver
395 - Use Scatterv to send back flipped flags for each component
396 - Negate flippedCells by component
397
398 NEW STYLE
399 - Create the adj on each process
400 - Bootstrap to complete graph on proc 0
401 */
402 /* Loop over components */
403 for (cell = cStart; cell < cEnd; ++cell) cellComp[cell - cStart] = -1;
404 do {
405 /* Look for first unmarked cell */
406 for (cell = cStart; cell < cEnd; ++cell)
407 if (cellComp[cell - cStart] < 0) break;
408 if (cell >= cEnd) break;
409 /* Initialize FIFO with first cell in component */
410 {
411 const PetscInt *cone;
412 PetscInt coneSize;
413
414 fTop = fBottom = 0;
415 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
416 PetscCall(DMPlexGetCone(dm, cell, &cone));
417 for (c = 0; c < coneSize; ++c) {
418 faceFIFO[fBottom++] = cone[c];
419 PetscCall(PetscBTSet(seenFaces, cone[c] - fStart));
420 }
421 PetscCall(PetscBTSet(seenCells, cell - cStart));
422 }
423 /* Consider each face in FIFO */
424 while (fTop < fBottom) PetscCall(DMPlexCheckFace_Old_Internal(dm, faceFIFO, &fTop, &fBottom, cStart, fStart, fEnd, seenCells, flippedCells, seenFaces));
425 /* Set component for cells and faces */
426 for (cell = 0; cell < cEnd - cStart; ++cell) {
427 if (PetscBTLookup(seenCells, cell)) cellComp[cell] = comp;
428 }
429 for (face = 0; face < fEnd - fStart; ++face) {
430 if (PetscBTLookup(seenFaces, face)) faceComp[face] = comp;
431 }
432 /* Wipe seenCells and seenFaces for next component */
433 PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
434 PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
435 ++comp;
436 } while (1);
437 numComponents = comp;
438 if (flg) {
439 PetscViewer v;
440
441 PetscCall(PetscViewerASCIIGetStdout(comm, &v));
442 PetscCall(PetscViewerASCIIPushSynchronized(v));
443 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
444 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
445 PetscCall(PetscViewerFlush(v));
446 PetscCall(PetscViewerASCIIPopSynchronized(v));
447 }
448 /* Now all subdomains are oriented, but we need a consistent parallel orientation */
449 if (numLeaves >= 0) {
450 PetscInt maxSupportSize, neighbor;
451
452 /* Store orientations of boundary faces*/
453 PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSupportSize));
454 PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSupportSize, &locSupport));
455 for (face = fStart; face < fEnd; ++face) {
456 const PetscInt *cone, *support, *ornt;
457 PetscInt coneSize, supportSize, Ns = 0, s, l;
458
459 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
460 /* Ignore overlapping cells */
461 PetscCall(DMPlexGetSupport(dm, face, &support));
462 for (s = 0; s < supportSize; ++s) {
463 if (lpoints) PetscCall(PetscFindInt(support[s], numLeaves, lpoints, &l));
464 else {
465 if (support[s] >= 0 && support[s] < numLeaves) l = support[s];
466 else l = -1;
467 }
468 if (l >= 0) continue;
469 locSupport[Ns++] = support[s];
470 }
471 if (Ns != 1) continue;
472 neighbor = locSupport[0];
473 PetscCall(DMPlexGetCone(dm, neighbor, &cone));
474 PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize));
475 PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt));
476 for (c = 0; c < coneSize; ++c)
477 if (cone[c] == face) break;
478 if (dim == 1) {
479 /* Use cone position instead, shifted to -1 or 1 */
480 if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = 1 - c * 2;
481 else rorntComp[face].rank = c * 2 - 1;
482 } else {
483 if (PetscBTLookup(flippedCells, neighbor - cStart)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1;
484 else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1;
485 }
486 rorntComp[face].index = faceComp[face - fStart];
487 }
488 /* Communicate boundary edge orientations */
489 PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
490 PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
491 }
492 /* Get process adjacency */
493 PetscCall(PetscMalloc2(numComponents, &numNeighbors, numComponents, &neighbors));
494 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
495 if (flg2) PetscCall(PetscViewerASCIIPushSynchronized(viewer));
496 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
497 for (comp = 0; comp < numComponents; ++comp) {
498 PetscInt l, n;
499
500 numNeighbors[comp] = 0;
501 PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp]));
502 /* I know this is p^2 time in general, but for bounded degree its alright */
503 for (l = 0; l < numLeaves; ++l) {
504 const PetscInt face = lpoints ? lpoints[l] : l;
505
506 /* Find a representative face (edge) separating pairs of procs */
507 if ((face >= fStart) && (face < fEnd) && (faceComp[face - fStart] == comp) && rorntComp[face].rank) {
508 const PetscInt rrank = rpoints[l].rank;
509 const PetscInt rcomp = lorntComp[face].index;
510
511 for (n = 0; n < numNeighbors[comp]; ++n)
512 if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break;
513 if (n >= numNeighbors[comp]) {
514 PetscInt supportSize;
515
516 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize));
517 // We can have internal faces in the SF if we have cells in the SF
518 if (supportSize > 1) continue;
519 if (flg)
520 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,
521 rpoints[l].index, rrank, rcomp, lorntComp[face].rank));
522 neighbors[comp][numNeighbors[comp]++] = l;
523 }
524 }
525 }
526 totNeighbors += numNeighbors[comp];
527 }
528 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
529 if (flg2) PetscCall(PetscViewerASCIIPopSynchronized(viewer));
530 PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match));
531 for (comp = 0, off = 0; comp < numComponents; ++comp) {
532 PetscInt n;
533
534 for (n = 0; n < numNeighbors[comp]; ++n, ++off) {
535 const PetscInt face = lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n];
536 const PetscInt o = rorntComp[face].rank * lorntComp[face].rank;
537
538 if (o < 0) match[off] = PETSC_TRUE;
539 else if (o > 0) match[off] = PETSC_FALSE;
540 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);
541 nrankComp[off].rank = rpoints[neighbors[comp][n]].rank;
542 nrankComp[off].index = lorntComp[lpoints ? lpoints[neighbors[comp][n]] : neighbors[comp][n]].index;
543 }
544 PetscCall(PetscFree(neighbors[comp]));
545 }
546 /* Collect the graph on 0 */
547 if (numLeaves >= 0) {
548 Mat G;
549 PetscBT seenProcs, flippedProcs;
550 PetscInt *procFIFO, pTop, pBottom;
551 PetscInt *N = NULL, *Noff;
552 PetscSFNode *adj = NULL;
553 PetscBool *val = NULL;
554 PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc, p, o, itotNeighbors;
555 PetscMPIInt size = 0;
556
557 PetscCall(PetscCalloc1(numComponents, &flipped));
558 if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
559 PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
560 PetscCallMPI(MPI_Gather(&numComponents, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
561 for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
562 if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
563 PetscCallMPI(MPI_Gatherv(numNeighbors, numComponents, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
564 for (p = 0, o = 0; p < size; ++p) {
565 recvcounts[p] = 0;
566 for (c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
567 displs[p + 1] = displs[p] + recvcounts[p];
568 }
569 if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
570 PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors));
571 PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm));
572 PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPI_C_BOOL, val, recvcounts, displs, MPI_C_BOOL, 0, comm));
573 PetscCall(PetscFree2(numNeighbors, neighbors));
574 if (rank == 0) {
575 for (p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
576 if (flg) {
577 PetscInt n;
578
579 for (p = 0, off = 0; p < size; ++p) {
580 for (c = 0; c < Nc[p]; ++c) {
581 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %d Comp %" PetscInt_FMT ":\n", p, c));
582 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]]));
583 }
584 }
585 }
586 /* Symmetrize the graph */
587 PetscCall(MatCreate(PETSC_COMM_SELF, &G));
588 PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
589 PetscCall(MatSetUp(G));
590 for (p = 0, off = 0; p < size; ++p) {
591 for (c = 0; c < Nc[p]; ++c) {
592 const PetscInt r = Noff[p] + c;
593 PetscInt n;
594
595 for (n = 0; n < N[r]; ++n, ++off) {
596 const PetscInt q = Noff[adj[off].rank] + adj[off].index;
597 const PetscScalar o = val[off] ? 1.0 : 0.0;
598
599 PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
600 PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
601 }
602 }
603 }
604 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
605 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
606
607 PetscCall(PetscBTCreate(Noff[size], &seenProcs));
608 PetscCall(PetscBTMemzero(Noff[size], seenProcs));
609 PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
610 PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
611 PetscCall(PetscMalloc1(Noff[size], &procFIFO));
612 pTop = pBottom = 0;
613 for (p = 0; p < Noff[size]; ++p) {
614 if (PetscBTLookup(seenProcs, p)) continue;
615 /* Initialize FIFO with next proc */
616 procFIFO[pBottom++] = p;
617 PetscCall(PetscBTSet(seenProcs, p));
618 /* Consider each proc in FIFO */
619 while (pTop < pBottom) {
620 const PetscScalar *ornt;
621 const PetscInt *neighbors;
622 PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors, n;
623
624 proc = procFIFO[pTop++];
625 flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
626 PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
627 /* Loop over neighboring procs */
628 for (n = 0; n < numNeighbors; ++n) {
629 nproc = neighbors[n];
630 mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
631 seen = PetscBTLookup(seenProcs, nproc);
632 flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
633
634 if (mismatch ^ (flippedA ^ flippedB)) {
635 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);
636 PetscCheck(!flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
637 PetscCall(PetscBTSet(flippedProcs, nproc));
638 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
639 if (!seen) {
640 procFIFO[pBottom++] = nproc;
641 PetscCall(PetscBTSet(seenProcs, nproc));
642 }
643 }
644 }
645 }
646 PetscCall(PetscFree(procFIFO));
647 PetscCall(MatDestroy(&G));
648 PetscCall(PetscFree2(adj, val));
649 PetscCall(PetscBTDestroy(&seenProcs));
650 }
651 /* Scatter flip flags */
652 {
653 PetscBool *flips = NULL;
654
655 if (rank == 0) {
656 PetscCall(PetscMalloc1(Noff[size], &flips));
657 for (p = 0; p < Noff[size]; ++p) {
658 flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
659 if (flg && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %d:\n", p));
660 }
661 for (p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
662 }
663 PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPI_C_BOOL, flipped, numComponents, MPI_C_BOOL, 0, comm));
664 PetscCall(PetscFree(flips));
665 }
666 if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
667 PetscCall(PetscFree(N));
668 PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
669 PetscCall(PetscFree2(nrankComp, match));
670
671 /* Decide whether to flip cells in each component */
672 for (c = 0; c < cEnd - cStart; ++c) {
673 if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c));
674 }
675 PetscCall(PetscFree(flipped));
676 }
677 if (flg) {
678 PetscViewer v;
679
680 PetscCall(PetscViewerASCIIGetStdout(comm, &v));
681 PetscCall(PetscViewerASCIIPushSynchronized(v));
682 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
683 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
684 PetscCall(PetscViewerFlush(v));
685 PetscCall(PetscViewerASCIIPopSynchronized(v));
686 }
687 /* Reverse flipped cells in the mesh */
688 for (c = cStart; c < cEnd; ++c) {
689 if (PetscBTLookup(flippedCells, c - cStart)) PetscCall(DMPlexOrientPoint(dm, c, -1));
690 }
691 PetscCall(PetscBTDestroy(&seenCells));
692 PetscCall(PetscBTDestroy(&flippedCells));
693 PetscCall(PetscBTDestroy(&seenFaces));
694 PetscCall(PetscFree2(numNeighbors, neighbors));
695 PetscCall(PetscFree3(rorntComp, lorntComp, locSupport));
696 PetscCall(PetscFree3(faceFIFO, cellComp, faceComp));
697 PetscFunctionReturn(PETSC_SUCCESS);
698 #endif
699 }
700
CreateCellAndFaceIS_Private(DM dm,DMLabel label,IS * cellIS,IS * faceIS)701 static PetscErrorCode CreateCellAndFaceIS_Private(DM dm, DMLabel label, IS *cellIS, IS *faceIS)
702 {
703 IS valueIS;
704 const PetscInt *values;
705 PetscInt Nv, depth = 0;
706
707 PetscFunctionBegin;
708 PetscCall(DMLabelGetValueIS(label, &valueIS));
709 PetscCall(ISGetLocalSize(valueIS, &Nv));
710 PetscCall(ISGetIndices(valueIS, &values));
711 for (PetscInt v = 0; v < Nv; ++v) {
712 const PetscInt val = values[v] < 0 || values[v] >= 100 ? 0 : values[v];
713 PetscInt n;
714
715 PetscCall(DMLabelGetStratumSize(label, val, &n));
716 if (!n) continue;
717 depth = PetscMax(val, depth);
718 }
719 PetscCall(ISDestroy(&valueIS));
720 PetscCheck(depth >= 1 || !Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Depth for interface must be at least 1, not %" PetscInt_FMT, depth);
721 PetscCall(DMLabelGetStratumIS(label, depth, cellIS));
722 PetscCall(DMLabelGetStratumIS(label, depth - 1, faceIS));
723 if (!*cellIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cellIS));
724 if (!*faceIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, faceIS));
725 PetscFunctionReturn(PETSC_SUCCESS);
726 }
727
DMPlexOrientLabel(DM dm,DMLabel label)728 PetscErrorCode DMPlexOrientLabel(DM dm, DMLabel label)
729 {
730 IS cellIS, faceIS;
731
732 PetscFunctionBegin;
733 PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
734 PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
735 PetscCall(ISDestroy(&cellIS));
736 PetscCall(ISDestroy(&faceIS));
737 PetscFunctionReturn(PETSC_SUCCESS);
738 }
739
DMPlexOrientCells_Internal(DM dm,IS cellIS,IS faceIS)740 PetscErrorCode DMPlexOrientCells_Internal(DM dm, IS cellIS, IS faceIS)
741 {
742 MPI_Comm comm;
743 PetscSF sf;
744 const PetscInt *lpoints;
745 const PetscSFNode *rpoints;
746 PetscSFNode *rorntComp = NULL, *lorntComp = NULL;
747 PetscInt *numNeighbors, **neighbors, *locSupp = NULL;
748 PetscSFNode *nrankComp;
749 PetscBool *match, *flipped;
750 PetscBT flippedCells;
751 PetscInt *cellComp, *faceComp;
752 const PetscInt *cells = NULL, *faces = NULL;
753 PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
754 PetscInt numLeaves, numRoots, dim, Ncomp, totNeighbors = 0;
755 PetscMPIInt rank, size;
756 PetscBool view, viewSync;
757 PetscViewer viewer = NULL, selfviewer = NULL;
758
759 PetscFunctionBegin;
760 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
761 PetscCallMPI(MPI_Comm_rank(comm, &rank));
762 PetscCallMPI(MPI_Comm_size(comm, &size));
763 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view", &view));
764 PetscCall(PetscOptionsHasName(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-orientation_view_synchronized", &viewSync));
765
766 if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
767 if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
768 PetscCall(DMGetPointSF(dm, &sf));
769 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints));
770 /* Truth Table
771 mismatch flips do action mismatch flipA ^ flipB action
772 F 0 flips no F F F
773 F 1 flip yes F T T
774 F 2 flips no T F T
775 T 0 flips yes T T F
776 T 1 flip no
777 T 2 flips yes
778 */
779 PetscCall(DMGetDimension(dm, &dim));
780 PetscCall(PetscBTCreate(cEnd - cStart, &flippedCells));
781 PetscCall(PetscBTMemzero(cEnd - cStart, flippedCells));
782 PetscCall(PetscCalloc2(cEnd - cStart, &cellComp, fEnd - fStart, &faceComp));
783 /*
784 OLD STYLE
785 - Add an integer array over cells and faces (component) for connected component number
786 Foreach component
787 - Mark the initial cell as seen
788 - Process component as usual
789 - Set component for all seenCells
790 - Wipe seenCells and seenFaces (flippedCells can stay)
791 - Generate parallel adjacency for component using SF and seenFaces
792 - Collect Ncomp adj data from each proc to 0
793 - Build same serial graph
794 - Use same solver
795 - Use Scatterv to send back flipped flags for each component
796 - Negate flippedCells by component
797
798 NEW STYLE
799 - Create the adj on each process
800 - Bootstrap to complete graph on proc 0
801 */
802 PetscCall(DMPlexOrient_Serial(dm, cellIS, faceIS, &Ncomp, cellComp, faceComp, flippedCells));
803 if (view) {
804 PetscViewer v;
805 PetscInt cdepth = -1;
806
807 PetscCall(PetscViewerASCIIGetStdout(comm, &v));
808 PetscCall(PetscViewerASCIIPushSynchronized(v));
809 if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
810 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]New Orientation %" PetscInt_FMT " cells (depth %" PetscInt_FMT ") and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, cdepth, fEnd - fStart));
811 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
812 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
813 PetscCall(PetscViewerFlush(v));
814 PetscCall(PetscViewerASCIIPopSynchronized(v));
815 }
816 /* Now all subdomains are oriented, but we need a consistent parallel orientation */
817 // TODO: This all has to be rewritten to filter cones/supports to the ISes
818 if (numLeaves >= 0) {
819 PetscInt maxSuppSize, neighbor;
820
821 // Store orientations of boundary faces
822 PetscCall(DMPlexGetMaxSizes(dm, NULL, &maxSuppSize));
823 PetscCall(PetscCalloc3(numRoots, &rorntComp, numRoots, &lorntComp, maxSuppSize, &locSupp));
824 for (PetscInt f = fStart; f < fEnd; ++f) {
825 const PetscInt face = faces ? faces[f] : f;
826 const PetscInt *cone, *supp, *ornt;
827 PetscInt coneSize, suppSize, nind, c, Ns = 0;
828
829 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
830 PetscCall(DMPlexGetSupport(dm, face, &supp));
831 for (PetscInt s = 0; s < suppSize; ++s) {
832 PetscInt ind, l;
833
834 // Filter support
835 ind = GetPointIndex(supp[s], cStart, cEnd, cells);
836 if (ind < 0) continue;
837 // Ignore overlapping cells
838 PetscCall(PetscFindInt(supp[s], numLeaves, lpoints, &l));
839 if (l >= 0) continue;
840 locSupp[Ns++] = supp[s];
841 }
842 PetscCheck(Ns < maxSuppSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " exceeds array size %" PetscInt_FMT, Ns, maxSuppSize);
843 if (Ns != 1) continue;
844 neighbor = locSupp[0];
845 nind = GetPointIndex(neighbor, cStart, cEnd, cells);
846 PetscCall(DMPlexGetCone(dm, neighbor, &cone));
847 PetscCall(DMPlexGetConeSize(dm, neighbor, &coneSize));
848 PetscCall(DMPlexGetConeOrientation(dm, neighbor, &ornt));
849 for (c = 0; c < coneSize; ++c)
850 if (cone[c] == face) break;
851 if (dim == 1) {
852 /* Use cone position instead, shifted to -1 or 1 */
853 if (PetscBTLookup(flippedCells, nind)) rorntComp[face].rank = 1 - c * 2;
854 else rorntComp[face].rank = c * 2 - 1;
855 } else {
856 if (PetscBTLookup(flippedCells, nind)) rorntComp[face].rank = ornt[c] < 0 ? -1 : 1;
857 else rorntComp[face].rank = ornt[c] < 0 ? 1 : -1;
858 }
859 rorntComp[face].index = faceComp[GetPointIndex(face, fStart, fEnd, faces)];
860 }
861 // Communicate boundary edge orientations
862 PetscCall(PetscSFBcastBegin(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
863 PetscCall(PetscSFBcastEnd(sf, MPIU_SF_NODE, rorntComp, lorntComp, MPI_REPLACE));
864 }
865 /* Get process adjacency */
866 PetscCall(PetscMalloc2(Ncomp, &numNeighbors, Ncomp, &neighbors));
867 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)dm));
868 if (viewSync) PetscCall(PetscViewerASCIIPushSynchronized(viewer));
869 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
870 for (PetscInt comp = 0; comp < Ncomp; ++comp) {
871 PetscInt n;
872
873 numNeighbors[comp] = 0;
874 PetscCall(PetscMalloc1(PetscMax(numLeaves, 0), &neighbors[comp]));
875 /* I know this is p^2 time in general, but for bounded degree its alright */
876 for (PetscInt l = 0; l < numLeaves; ++l) {
877 const PetscInt face = lpoints[l];
878 PetscInt find;
879
880 /* Find a representative face (edge) separating pairs of procs */
881 find = GetPointIndex(face, fStart, fEnd, faces);
882 if ((find >= 0) && (faceComp[find] == comp) && rorntComp[face].rank) {
883 const PetscInt rrank = rpoints[l].rank;
884 const PetscInt rcomp = lorntComp[face].index;
885
886 for (n = 0; n < numNeighbors[comp]; ++n)
887 if ((rrank == rpoints[neighbors[comp][n]].rank) && (rcomp == lorntComp[lpoints[neighbors[comp][n]]].index)) break;
888 if (n >= numNeighbors[comp]) {
889 const PetscInt *supp;
890 PetscInt suppSize, Ns = 0;
891
892 PetscCall(DMPlexGetSupport(dm, face, &supp));
893 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
894 for (PetscInt s = 0; s < suppSize; ++s) {
895 // Filter support
896 if (GetPointIndex(supp[s], cStart, cEnd, cells) >= 0) ++Ns;
897 }
898 PetscCheck(Ns == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary face %" PetscInt_FMT " should see one cell, not %" PetscInt_FMT, face, Ns);
899 if (view)
900 PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]: component %" PetscInt_FMT ", 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,
901 rpoints[l].index, rrank, rcomp, lorntComp[face].rank));
902 neighbors[comp][numNeighbors[comp]++] = l;
903 }
904 }
905 }
906 totNeighbors += numNeighbors[comp];
907 }
908 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &selfviewer));
909 if (viewSync) PetscCall(PetscViewerASCIIPopSynchronized(viewer));
910 PetscCall(PetscMalloc2(totNeighbors, &nrankComp, totNeighbors, &match));
911 for (PetscInt comp = 0, off = 0; comp < Ncomp; ++comp) {
912 for (PetscInt n = 0; n < numNeighbors[comp]; ++n, ++off) {
913 const PetscInt face = lpoints[neighbors[comp][n]];
914 const PetscInt o = rorntComp[face].rank * lorntComp[face].rank;
915
916 if (o < 0) match[off] = PETSC_TRUE;
917 else if (o > 0) match[off] = PETSC_FALSE;
918 else
919 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid face %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ") neighbor: %" PetscInt_FMT " comp: %" PetscInt_FMT, face, rorntComp[face].rank, lorntComp[face].rank, neighbors[comp][n], comp);
920 nrankComp[off].rank = rpoints[neighbors[comp][n]].rank;
921 nrankComp[off].index = lorntComp[lpoints[neighbors[comp][n]]].index;
922 }
923 PetscCall(PetscFree(neighbors[comp]));
924 }
925 /* Collect the graph on 0 */
926 if (numLeaves >= 0) {
927 Mat G;
928 PetscBT seenProcs, flippedProcs;
929 PetscInt *procFIFO, pTop, pBottom;
930 PetscInt *N = NULL, *Noff;
931 PetscSFNode *adj = NULL;
932 PetscBool *val = NULL;
933 PetscMPIInt *recvcounts = NULL, *displs = NULL, *Nc;
934 PetscMPIInt size = 0, iNcomp, itotNeighbors;
935
936 PetscCall(PetscCalloc1(Ncomp, &flipped));
937 if (rank == 0) PetscCallMPI(MPI_Comm_size(comm, &size));
938 PetscCall(PetscCalloc4(size, &recvcounts, size + 1, &displs, size, &Nc, size + 1, &Noff));
939 PetscCallMPI(MPI_Gather(&Ncomp, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
940 for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
941 if (rank == 0) PetscCall(PetscMalloc1(displs[size], &N));
942 PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
943 PetscCallMPI(MPI_Gatherv(numNeighbors, iNcomp, MPIU_INT, N, Nc, displs, MPIU_INT, 0, comm));
944 for (PetscInt p = 0, o = 0; p < size; ++p) {
945 recvcounts[p] = 0;
946 for (PetscInt c = 0; c < Nc[p]; ++c, ++o) recvcounts[p] += N[o];
947 displs[p + 1] = displs[p] + recvcounts[p];
948 }
949 if (rank == 0) PetscCall(PetscMalloc2(displs[size], &adj, displs[size], &val));
950 PetscCall(PetscMPIIntCast(totNeighbors, &itotNeighbors));
951 PetscCallMPI(MPI_Gatherv(nrankComp, itotNeighbors, MPIU_SF_NODE, adj, recvcounts, displs, MPIU_SF_NODE, 0, comm));
952 PetscCallMPI(MPI_Gatherv(match, itotNeighbors, MPI_C_BOOL, val, recvcounts, displs, MPI_C_BOOL, 0, comm));
953 PetscCall(PetscFree2(numNeighbors, neighbors));
954 if (rank == 0) {
955 for (PetscInt p = 1; p <= size; ++p) Noff[p] = Noff[p - 1] + Nc[p - 1];
956 if (view) {
957 for (PetscInt p = 0, off = 0; p < size; ++p) {
958 for (PetscInt c = 0; c < Nc[p]; ++c) {
959 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Proc %" PetscInt_FMT " Comp %" PetscInt_FMT ":\n", p, c));
960 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]]));
961 }
962 }
963 }
964 /* Symmetrize the graph */
965 PetscCall(MatCreate(PETSC_COMM_SELF, &G));
966 PetscCall(MatSetSizes(G, Noff[size], Noff[size], Noff[size], Noff[size]));
967 PetscCall(MatSetUp(G));
968 for (PetscInt p = 0, off = 0; p < size; ++p) {
969 for (PetscInt c = 0; c < Nc[p]; ++c) {
970 const PetscInt r = Noff[p] + c;
971
972 for (PetscInt n = 0; n < N[r]; ++n, ++off) {
973 const PetscInt q = Noff[adj[off].rank] + adj[off].index;
974 const PetscScalar o = val[off] ? 1.0 : 0.0;
975
976 PetscCall(MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES));
977 PetscCall(MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES));
978 }
979 }
980 }
981 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
982 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
983
984 PetscCall(PetscBTCreate(Noff[size], &seenProcs));
985 PetscCall(PetscBTMemzero(Noff[size], seenProcs));
986 PetscCall(PetscBTCreate(Noff[size], &flippedProcs));
987 PetscCall(PetscBTMemzero(Noff[size], flippedProcs));
988 PetscCall(PetscMalloc1(Noff[size], &procFIFO));
989 pTop = pBottom = 0;
990 for (PetscInt p = 0; p < Noff[size]; ++p) {
991 if (PetscBTLookup(seenProcs, p)) continue;
992 /* Initialize FIFO with next proc */
993 procFIFO[pBottom++] = p;
994 PetscCall(PetscBTSet(seenProcs, p));
995 /* Consider each proc in FIFO */
996 while (pTop < pBottom) {
997 const PetscScalar *ornt;
998 const PetscInt *neighbors;
999 PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors;
1000
1001 proc = procFIFO[pTop++];
1002 flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
1003 PetscCall(MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt));
1004 /* Loop over neighboring procs */
1005 for (PetscInt n = 0; n < numNeighbors; ++n) {
1006 nproc = neighbors[n];
1007 mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
1008 seen = PetscBTLookup(seenProcs, nproc);
1009 flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
1010
1011 if (mismatch ^ (flippedA ^ flippedB)) {
1012 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);
1013 PetscCheck(!flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
1014 PetscCall(PetscBTSet(flippedProcs, nproc));
1015 } else PetscCheck(!mismatch || !flippedA || !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
1016 if (!seen) {
1017 procFIFO[pBottom++] = nproc;
1018 PetscCall(PetscBTSet(seenProcs, nproc));
1019 }
1020 }
1021 }
1022 }
1023 PetscCall(PetscFree(procFIFO));
1024 PetscCall(MatDestroy(&G));
1025 PetscCall(PetscFree2(adj, val));
1026 PetscCall(PetscBTDestroy(&seenProcs));
1027 }
1028 /* Scatter flip flags */
1029 {
1030 PetscBool *flips = NULL;
1031
1032 if (rank == 0) {
1033 PetscCall(PetscMalloc1(Noff[size], &flips));
1034 for (PetscInt p = 0; p < Noff[size]; ++p) {
1035 flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
1036 if (view && flips[p]) PetscCall(PetscPrintf(comm, "Flipping Proc+Comp %" PetscInt_FMT ":\n", p));
1037 }
1038 for (PetscInt p = 0; p < size; ++p) displs[p + 1] = displs[p] + Nc[p];
1039 }
1040 PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
1041 PetscCallMPI(MPI_Scatterv(flips, Nc, displs, MPI_C_BOOL, flipped, iNcomp, MPI_C_BOOL, 0, comm));
1042 PetscCall(PetscFree(flips));
1043 }
1044 if (rank == 0) PetscCall(PetscBTDestroy(&flippedProcs));
1045 PetscCall(PetscFree(N));
1046 PetscCall(PetscFree4(recvcounts, displs, Nc, Noff));
1047 PetscCall(PetscFree2(nrankComp, match));
1048
1049 /* Decide whether to flip cells in each component */
1050 for (PetscInt c = 0; c < cEnd - cStart; ++c) {
1051 if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(flippedCells, c));
1052 }
1053 PetscCall(PetscFree(flipped));
1054 }
1055 if (view) {
1056 PetscViewer v;
1057
1058 PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1059 PetscCall(PetscViewerASCIIPushSynchronized(v));
1060 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
1061 PetscCall(PetscBTView(cEnd - cStart, flippedCells, v));
1062 PetscCall(PetscViewerFlush(v));
1063 PetscCall(PetscViewerASCIIPopSynchronized(v));
1064 }
1065 // Reverse flipped cells in the mesh
1066 PetscViewer v;
1067 const PetscInt *degree = NULL;
1068 PetscInt *points;
1069 PetscInt pStart, pEnd;
1070
1071 if (view) {
1072 PetscCall(PetscViewerASCIIGetStdout(comm, &v));
1073 PetscCall(PetscViewerASCIIPushSynchronized(v));
1074 }
1075 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1076 if (numRoots >= 0) {
1077 PetscCall(PetscSFComputeDegreeBegin(sf, °ree));
1078 PetscCall(PetscSFComputeDegreeEnd(sf, °ree));
1079 }
1080 PetscCall(PetscCalloc1(pEnd - pStart, &points));
1081 for (PetscInt c = cStart; c < cEnd; ++c) {
1082 if (PetscBTLookup(flippedCells, c - cStart)) {
1083 const PetscInt cell = cells ? cells[c] : c;
1084
1085 PetscCall(DMPlexOrientPoint(dm, cell, -1));
1086 if (degree && degree[cell]) points[cell] = 1;
1087 if (view) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT "%s\n", rank, cell, degree && degree[cell] ? " and sending to overlap" : ""));
1088 }
1089 }
1090 // Must propagate flips for cells in the overlap
1091 if (numRoots >= 0) {
1092 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, points, points, MPI_SUM));
1093 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, points, points, MPI_SUM));
1094 }
1095 for (PetscInt c = cStart; c < cEnd; ++c) {
1096 const PetscInt cell = cells ? cells[c] : c;
1097
1098 if (points[cell] && !PetscBTLookup(flippedCells, c - cStart)) {
1099 PetscCall(DMPlexOrientPoint(dm, cell, -1));
1100 if (view) PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]Flipping cell %" PetscInt_FMT " through overlap\n", rank, cell));
1101 }
1102 }
1103 if (view) {
1104 PetscCall(PetscViewerFlush(v));
1105 PetscCall(PetscViewerASCIIPopSynchronized(v));
1106 }
1107 PetscCall(PetscFree(points));
1108 PetscCall(PetscBTDestroy(&flippedCells));
1109 PetscCall(PetscFree2(numNeighbors, neighbors));
1110 PetscCall(PetscFree3(rorntComp, lorntComp, locSupp));
1111 PetscCall(PetscFree2(cellComp, faceComp));
1112 PetscFunctionReturn(PETSC_SUCCESS);
1113 }
1114