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 = 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;
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;
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: }
68: static PetscInt GetPointIndex(PetscInt point, PetscInt pStart, PetscInt pEnd, const PetscInt points[])
69: {
70: if (points) {
71: PetscInt loc;
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: }
81: /*
82: - Checks face match
83: - Flips non-matching
84: - Inserts faces of support cells in FIFO
85: */
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;
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;
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);
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);
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]]);
149: if (dim == 1) {
150: mismatch = posA == posB;
151: } else {
152: mismatch = coneOA[posA] == coneOB[posB];
153: }
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) PetscCall(PetscBTSet(flippedCells, indC[0]));
158: else {
159: PetscCheck(!seenB && !flippedB, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
160: PetscCall(PetscBTSet(flippedCells, indC[1]));
161: }
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: }
168: /*
169: DMPlexOrient_Serial - Compute valid orientation for local connected components
171: Not collective
173: Input Parameters:
174: + dm - The `DM`
175: . cellIS - The cells to orient
176: - faceIS - The faces between the cells
178: Output Parameters:
179: + Ncomp - The number of connected component
180: . cellComp - The connected component for each local cell
181: - flippedCells - Marked cells should be inverted
183: Level: developer
185: .seealso: `DMPlexOrient()`
186: */
187: static PetscErrorCode DMPlexOrient_Serial(DM dm, IS cellIS, IS faceIS, PetscInt *Ncomp, PetscInt cellComp[], PetscBT flippedCells)
188: {
189: PetscBT seenCells, seenFaces;
190: PetscInt *faceFIFO;
191: const PetscInt *cells = NULL, *faces = NULL;
192: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
194: PetscFunctionBegin;
195: /* Truth Table
196: mismatch flips do action mismatch flipA ^ flipB action
197: F 0 flips no F F F
198: F 1 flip yes F T T
199: F 2 flips no T F T
200: T 0 flips yes T T F
201: T 1 flip no
202: T 2 flips yes
203: */
204: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
205: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
206: PetscCall(PetscBTCreate(cEnd - cStart, &seenCells));
207: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
208: PetscCall(PetscBTCreate(fEnd - fStart, &seenFaces));
209: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
210: PetscCall(PetscMalloc1(fEnd - fStart, &faceFIFO));
211: *Ncomp = 0;
212: for (PetscInt c = 0; c < cEnd - cStart; ++c) cellComp[c] = -1;
213: do {
214: PetscInt cc, fTop, fBottom;
216: // Look for first unmarked cell
217: for (cc = cStart; cc < cEnd; ++cc)
218: if (cellComp[cc - cStart] < 0) break;
219: if (cc >= cEnd) break;
220: // Initialize FIFO with first cell in component
221: {
222: const PetscInt cell = cells ? cells[cc] : cc;
223: const PetscInt *cone;
224: PetscInt coneSize;
226: fTop = fBottom = 0;
227: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
228: PetscCall(DMPlexGetCone(dm, cell, &cone));
229: for (PetscInt c = 0; c < coneSize; ++c) {
230: const PetscInt idx = GetPointIndex(cone[c], fStart, fEnd, faces);
232: // Cell faces are guaranteed to be in the face set
233: 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);
234: faceFIFO[fBottom++] = cone[c];
235: PetscCall(PetscBTSet(seenFaces, idx));
236: }
237: PetscCall(PetscBTSet(seenCells, cc - cStart));
238: }
239: // Consider each face in FIFO
240: while (fTop < fBottom) PetscCall(DMPlexCheckFace_Internal(dm, faceFIFO, &fTop, &fBottom, cellIS, faceIS, seenCells, flippedCells, seenFaces));
241: // Set component for cells
242: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
243: if (PetscBTLookup(seenCells, c)) cellComp[c] = *Ncomp;
244: }
245: // Wipe seenCells and seenFaces for next component
246: PetscCall(PetscBTMemzero(fEnd - fStart, seenFaces));
247: PetscCall(PetscBTMemzero(cEnd - cStart, seenCells));
248: ++(*Ncomp);
249: } while (1);
250: PetscCall(PetscBTDestroy(&seenCells));
251: PetscCall(PetscBTDestroy(&seenFaces));
252: PetscCall(PetscFree(faceFIFO));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: // A local cell next to a face, with the orientation that it induces on the face
257: typedef struct {
258: PetscInt rank; // process holding the cell
259: PetscInt comp; // connected component of the cell on that process
260: PetscInt ornt; // orientation of the face in the cell, 1 or -1, or 0 if there is no cell
261: PetscSFNode owner; // owner of the cell, which identifies the copies of the cell on other processes
262: } FaceSample;
264: // Copies of the same cell induce the same orientation on a face, and two different cells induce opposite orientations
265: static PetscBool FaceSamplesMatch(const FaceSample *a, const FaceSample *b)
266: {
267: const PetscBool same = a->owner.rank == b->owner.rank && a->owner.index == b->owner.index ? PETSC_TRUE : PETSC_FALSE;
269: return (a->ornt == b->ornt) == same ? PETSC_TRUE : PETSC_FALSE;
270: }
272: /*
273: DMPlexGetFaceSamples_Private - Collect on the owner of each face the samples from all processes that hold the face
275: Collective
277: Each process samples the local cells of cellIS next to each face of faceIS, with the orientation of cellFlip. For each
278: face faces[f] that this process owns, the samples are samples[off[f]] to samples[off[f + 1]]. With overlap, a cell can
279: have copies on several processes, so a face can have more than two samples.
280: */
281: static PetscErrorCode DMPlexGetFaceSamples_Private(DM dm, IS cellIS, IS faceIS, const PetscInt cellComp[], PetscBT cellFlip, PetscInt **off, FaceSample **samples)
282: {
283: PetscSF sf;
284: const PetscInt *lpoints, *rootdegree = NULL, *cells = NULL, *faces = NULL;
285: const PetscSFNode *rpoints;
286: FaceSample *local, *remote = NULL;
287: PetscInt *roff;
288: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0, pEnd, Nr, Nl;
289: PetscMPIInt rank;
291: PetscFunctionBegin;
292: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
293: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
294: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
295: PetscCall(DMPlexGetChart(dm, NULL, &pEnd));
296: PetscCall(DMGetPointSF(dm, &sf));
297: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &lpoints, &rpoints));
298: if (Nr < 0) {
299: sf = NULL;
300: Nl = 0;
301: }
302: PetscCall(PetscCalloc2(2 * pEnd, &local, pEnd + 1, &roff));
303: for (PetscInt f = fStart; f < fEnd; ++f) {
304: const PetscInt face = faces ? faces[f] : f;
305: const PetscInt *supp;
306: PetscInt sS, depth, n = 0;
308: PetscCall(DMPlexGetPointDepth(dm, face, &depth));
309: PetscCall(DMPlexGetSupportSize(dm, face, &sS));
310: PetscCall(DMPlexGetSupport(dm, face, &supp));
311: for (PetscInt s = 0; s < sS; ++s) {
312: const PetscInt cind = GetPointIndex(supp[s], cStart, cEnd, cells);
313: const PetscInt l = GetPointIndex(supp[s], 0, Nl, lpoints);
314: FaceSample *fs = &local[2 * face + n];
315: const PetscInt *cone, *ornt;
316: PetscInt cS, c;
318: if (cind < 0) continue;
319: PetscCheck(n < 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " separates more than two cells", face);
320: PetscCall(DMPlexGetConeSize(dm, supp[s], &cS));
321: PetscCall(DMPlexGetOrientedCone(dm, supp[s], &cone, &ornt));
322: for (c = 0; c < cS; ++c)
323: if (cone[c] == face) break;
324: // A vertex is oriented by its position in the cone of an edge
325: fs->ornt = depth ? (ornt[c] < 0 ? -1 : 1) : 2 * c - 1;
326: PetscCall(DMPlexRestoreOrientedCone(dm, supp[s], &cone, &ornt));
327: if (cellFlip && PetscBTLookup(cellFlip, cind)) fs->ornt = -fs->ornt;
328: fs->rank = rank;
329: fs->comp = cellComp ? cellComp[cind] : 0;
330: if (l >= 0) fs->owner = rpoints[l];
331: else {
332: fs->owner.rank = rank;
333: fs->owner.index = supp[s];
334: }
335: ++n;
336: }
337: }
338: if (sf) {
339: MPI_Datatype MPIU_10INT;
341: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
342: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
343: for (PetscInt p = 0; p < pEnd; ++p) roff[p + 1] = roff[p] + (p < Nr ? rootdegree[p] : 0);
344: PetscCall(PetscMalloc1(2 * roff[pEnd], &remote));
345: PetscCallMPI(MPI_Type_contiguous(10, MPIU_INT, &MPIU_10INT));
346: PetscCallMPI(MPI_Type_commit(&MPIU_10INT));
347: PetscCall(PetscSFGatherBegin(sf, MPIU_10INT, local, remote));
348: PetscCall(PetscSFGatherEnd(sf, MPIU_10INT, local, remote));
349: PetscCallMPI(MPI_Type_free(&MPIU_10INT));
350: }
351: // An owned face has its local samples, and two from each process holding a copy of it
352: PetscCall(PetscMalloc1(fEnd - fStart + 1, off));
353: (*off)[0] = 0;
354: for (PetscInt f = fStart; f < fEnd; ++f) {
355: const PetscInt face = faces ? faces[f] : f;
357: (*off)[f - fStart + 1] = (*off)[f - fStart] + (GetPointIndex(face, 0, Nl, lpoints) < 0 ? 2 * (1 + roff[face + 1] - roff[face]) : 0);
358: }
359: PetscCall(PetscMalloc1((*off)[fEnd - fStart], samples));
360: for (PetscInt f = fStart; f < fEnd; ++f) {
361: const PetscInt face = faces ? faces[f] : f;
362: FaceSample *fs;
364: if ((*off)[f - fStart + 1] == (*off)[f - fStart]) continue;
365: fs = PetscSafePointerPlusOffset(*samples, (*off)[f - fStart]);
366: PetscCall(PetscArraycpy(fs, PetscSafePointerPlusOffset(local, 2 * face), 2));
367: if (remote) PetscCall(PetscArraycpy(&fs[2], PetscSafePointerPlusOffset(remote, 2 * roff[face]), 2 * (roff[face + 1] - roff[face])));
368: }
369: PetscCall(PetscFree(remote));
370: PetscCall(PetscFree2(local, roff));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*
375: DMPlexOrientSolve_Private - Choose which components to flip, so that all edges of the process graph match
377: Collective
379: Each edge links component (rank, comp) with component (rank, comp) and says whether their orientations match now.
380: Process 0 gathers the graph and flips the components breadth first, and flipped[c] is true if the local component c
381: must be flipped.
382: */
383: static PetscErrorCode DMPlexOrientSolve_Private(DM dm, PetscInt Ncomp, PetscInt Ne, const PetscInt edges[], PetscBool flipped[])
384: {
385: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
386: MPI_Comm comm;
387: PetscMPIInt rank, size, iNcomp, iNe;
388: PetscMPIInt *Nc = NULL, *coff = NULL, *ecounts = NULL, *eoff = NULL;
389: PetscInt *allEdges = NULL, *adj = NULL, *aoff = NULL, *queue = NULL;
390: PetscBool *flips = NULL, *seen = NULL, orientable = PETSC_TRUE;
392: PetscFunctionBegin;
393: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
394: PetscCallMPI(MPI_Comm_rank(comm, &rank));
395: PetscCallMPI(MPI_Comm_size(comm, &size));
396: PetscCall(PetscMPIIntCast(Ncomp, &iNcomp));
397: PetscCall(PetscMPIIntCast(5 * Ne, &iNe));
398: if (rank == 0) PetscCall(PetscCalloc4(size, &Nc, size + 1, &coff, size, &ecounts, size + 1, &eoff));
399: PetscCallMPI(MPI_Gather(&iNcomp, 1, MPI_INT, Nc, 1, MPI_INT, 0, comm));
400: PetscCallMPI(MPI_Gather(&iNe, 1, MPI_INT, ecounts, 1, MPI_INT, 0, comm));
401: if (rank == 0) {
402: for (PetscMPIInt p = 0; p < size; ++p) {
403: coff[p + 1] = coff[p] + Nc[p];
404: eoff[p + 1] = eoff[p] + ecounts[p];
405: }
406: PetscCall(PetscMalloc1(eoff[size], &allEdges));
407: }
408: PetscCallMPI(MPI_Gatherv(edges, iNe, MPIU_INT, allEdges, ecounts, eoff, MPIU_INT, 0, comm));
409: if (rank == 0) {
410: const PetscInt N = coff[size], E = eoff[size] / 5;
411: PetscInt qTop = 0, qBottom = 0;
413: // Store both directions of each edge, as the neighbor and whether it matches
414: PetscCall(PetscCalloc5(N + 1, &aoff, 4 * E, &adj, N, &queue, N, &flips, N, &seen));
415: for (PetscInt e = 0; e < E; ++e) {
416: ++aoff[coff[allEdges[5 * e + 0]] + allEdges[5 * e + 1] + 1];
417: ++aoff[coff[allEdges[5 * e + 2]] + allEdges[5 * e + 3] + 1];
418: }
419: for (PetscInt n = 0; n < N; ++n) aoff[n + 1] += aoff[n];
420: for (PetscInt e = 0; e < E; ++e) {
421: const PetscInt a = coff[allEdges[5 * e + 0]] + allEdges[5 * e + 1], b = coff[allEdges[5 * e + 2]] + allEdges[5 * e + 3];
423: if (debug)
424: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Edge (%" PetscInt_FMT ", %" PetscInt_FMT ") ~ (%" PetscInt_FMT ", %" PetscInt_FMT ") (%s)\n", allEdges[5 * e + 0], allEdges[5 * e + 1], allEdges[5 * e + 2], allEdges[5 * e + 3], PetscBools[allEdges[5 * e + 4]]));
425: adj[2 * aoff[a]] = b;
426: adj[2 * aoff[a] + 1] = allEdges[5 * e + 4];
427: ++aoff[a];
428: adj[2 * aoff[b]] = a;
429: adj[2 * aoff[b] + 1] = allEdges[5 * e + 4];
430: ++aoff[b];
431: }
432: for (PetscInt n = N; n > 0; --n) aoff[n] = aoff[n - 1];
433: aoff[0] = 0;
434: // A component is flipped relative to its neighbor if their orientations do not match
435: for (PetscInt n = 0; n < N; ++n) {
436: if (seen[n]) continue;
437: seen[n] = PETSC_TRUE;
438: queue[qBottom++] = n;
439: while (qTop < qBottom) {
440: const PetscInt a = queue[qTop++];
442: for (PetscInt i = aoff[a]; i < aoff[a + 1]; ++i) {
443: const PetscInt b = adj[2 * i];
444: const PetscBool flip = adj[2 * i + 1] ? flips[a] : (PetscBool)!flips[a];
446: if (!seen[b]) {
447: seen[b] = PETSC_TRUE;
448: flips[b] = flip;
449: queue[qBottom++] = b;
450: } else if (flips[b] != flip) orientable = PETSC_FALSE;
451: }
452: }
453: }
454: }
455: PetscCallMPI(MPI_Bcast(&orientable, 1, MPI_C_BOOL, 0, comm));
456: PetscCheck(orientable, comm, PETSC_ERR_ARG_WRONG, "Mesh is non-orientable");
457: PetscCallMPI(MPI_Scatterv(flips, Nc, coff, MPI_C_BOOL, flipped, iNcomp, MPI_C_BOOL, 0, comm));
458: if (rank == 0) {
459: PetscCall(PetscFree5(aoff, adj, queue, flips, seen));
460: PetscCall(PetscFree(allEdges));
461: PetscCall(PetscFree4(Nc, coff, ecounts, eoff));
462: }
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /*
467: The orientation operation needs to orient both a bulk mesh and a manifold embedded in a mesh, whose faces can be
468: shared by several processes, and whose cells can have copies on several processes when the mesh has overlap.
470: 1) Each process orients its local cells, and counts the connected components.
472: 2) The owner of each face gathers the orientation that every copy of every cell next to the face induces on it.
474: 3) The owner links the components of these cells into a process graph, whose edges say whether two components match.
476: 4) Process 0 gathers the graph, finds which components to flip, and returns the result to each process.
478: 5) Each process flips the cells of the flipped components.
479: */
480: PetscErrorCode DMPlexOrientCells_Internal(DM dm, IS cellIS, IS faceIS)
481: {
482: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
483: const PetscInt *cells = NULL, *faces = NULL;
484: PetscInt cStart = 0, cEnd = 0, fStart = 0, fEnd = 0;
485: PetscBT cellFlip; // The bit is true if a cell should have its orientation reversed
486: PetscInt *cellComp; // The connected component number of each cell
487: PetscInt Ncomp = 0; // The number of local connected components
488: PetscInt *off, *edges, Ne = 0;
489: FaceSample *samples;
490: PetscBool *flipped;
491: MPI_Comm comm;
492: PetscMPIInt rank;
494: PetscFunctionBegin;
495: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
496: PetscCallMPI(MPI_Comm_rank(comm, &rank));
497: if (cellIS) PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
498: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
499: PetscCall(PetscBTCreate(cEnd - cStart, &cellFlip));
500: PetscCall(PetscBTMemzero(cEnd - cStart, cellFlip));
501: PetscCall(PetscMalloc1(cEnd - cStart, &cellComp));
502: // Phase 1: Serial Orientation
503: PetscCall(DMPlexOrient_Serial(dm, cellIS, faceIS, &Ncomp, cellComp, cellFlip));
504: if (debug) {
505: PetscViewer v;
506: PetscInt cdepth = -1;
508: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
509: PetscCall(PetscViewerASCIIPushSynchronized(v));
510: if (cEnd > cStart) PetscCall(DMPlexGetPointDepth(dm, cells ? cells[cStart] : cStart, &cdepth));
511: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]New Orientation %" PetscInt_FMT " cells (depth %" PetscInt_FMT ") and %" PetscInt_FMT " faces\n", rank, cEnd - cStart, cdepth, fEnd - fStart));
512: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for serial flipped cells:\n", rank));
513: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
514: PetscCall(PetscViewerFlush(v));
515: PetscCall(PetscViewerASCIIPopSynchronized(v));
516: }
517: // Phase 2
518: PetscCall(DMPlexGetFaceSamples_Private(dm, cellIS, faceIS, cellComp, cellFlip, &off, &samples));
519: // Phase 3: Link the first sample of each face to the samples from other components
520: for (PetscInt pass = 0; pass < 2; ++pass) {
521: PetscInt e = 0;
523: for (PetscInt f = 0; f < fEnd - fStart; ++f) {
524: const FaceSample *first = NULL;
526: for (PetscInt i = off[f]; i < off[f + 1]; ++i) {
527: const FaceSample *s = &samples[i];
529: if (!s->ornt) continue;
530: if (!first) first = s;
531: else if (s->rank != first->rank || s->comp != first->comp) {
532: if (pass) {
533: edges[5 * e + 0] = first->rank;
534: edges[5 * e + 1] = first->comp;
535: edges[5 * e + 2] = s->rank;
536: edges[5 * e + 3] = s->comp;
537: edges[5 * e + 4] = FaceSamplesMatch(first, s);
538: }
539: ++e;
540: }
541: }
542: }
543: if (!pass) {
544: Ne = e;
545: PetscCall(PetscMalloc1(5 * Ne, &edges));
546: }
547: }
548: // Phase 4
549: PetscCall(PetscMalloc1(Ncomp, &flipped));
550: PetscCall(DMPlexOrientSolve_Private(dm, Ncomp, Ne, edges, flipped));
551: // Phase 5
552: for (PetscInt c = 0; c < cEnd - cStart; ++c) {
553: if (flipped[cellComp[c]]) PetscCall(PetscBTNegate(cellFlip, c));
554: if (PetscBTLookup(cellFlip, c)) PetscCall(DMPlexOrientPoint(dm, cells ? cells[cStart + c] : cStart + c, -1));
555: }
556: if (debug) {
557: PetscViewer v;
559: PetscCall(PetscViewerASCIIGetStdout(comm, &v));
560: PetscCall(PetscViewerASCIIPushSynchronized(v));
561: PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d]BT for parallel flipped cells:\n", rank));
562: PetscCall(PetscBTView(cEnd - cStart, cellFlip, v));
563: PetscCall(PetscViewerFlush(v));
564: PetscCall(PetscViewerASCIIPopSynchronized(v));
565: }
566: PetscCall(PetscFree(flipped));
567: PetscCall(PetscFree(edges));
568: PetscCall(PetscFree(off));
569: PetscCall(PetscFree(samples));
570: PetscCall(PetscBTDestroy(&cellFlip));
571: PetscCall(PetscFree(cellComp));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*@
576: DMPlexOrient - Give a consistent orientation to the input mesh
578: Input Parameter:
579: . dm - The `DM`
581: Notes:
582: The orientation data for the `DM` are changed in-place.
584: This routine will fail for non-orientable surfaces, such as the Moebius strip.
586: Level: advanced
588: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexOrientLabel()`
589: @*/
590: PetscErrorCode DMPlexOrient(DM dm)
591: {
592: IS cellIS, faceIS;
593: PetscInt h, cStart, cEnd, fStart, fEnd;
595: PetscFunctionBegin;
596: PetscCall(DMPlexGetVTKCellHeight(dm, &h));
597: PetscCall(DMPlexGetHeightStratum(dm, h, &cStart, &cEnd));
598: PetscCall(DMPlexGetHeightStratum(dm, h + 1, &fStart, &fEnd));
599: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cStart, cStart, 1, &cellIS));
600: PetscCall(ISCreateStride(PETSC_COMM_SELF, fEnd - fStart, fStart, 1, &faceIS));
601: PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
602: PetscCall(ISDestroy(&cellIS));
603: PetscCall(ISDestroy(&faceIS));
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: // The depth of the surface is the largest depth in the label on any process, since a process can hold only part of its closure
608: static PetscErrorCode CreateCellAndFaceIS_Private(DM dm, DMLabel label, IS *cellIS, IS *faceIS)
609: {
610: IS valueIS;
611: const PetscInt *values;
612: PetscInt Nv, depth;
614: PetscFunctionBegin;
615: PetscCall(DMLabelGetValueIS(label, &valueIS));
616: PetscCall(ISGetLocalSize(valueIS, &Nv));
617: PetscCall(ISGetIndices(valueIS, &values));
618: depth = Nv ? 0 : -1;
619: for (PetscInt v = 0; v < Nv; ++v) {
620: const PetscInt val = values[v] < 0 || values[v] >= 100 ? 0 : values[v];
621: PetscInt n;
623: PetscCall(DMLabelGetStratumSize(label, val, &n));
624: if (!n) continue;
625: depth = PetscMax(val, depth);
626: }
627: PetscCall(ISRestoreIndices(valueIS, &values));
628: PetscCall(ISDestroy(&valueIS));
629: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &depth, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
630: PetscCheck(depth, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Depth for interface must be at least 1, not %" PetscInt_FMT, depth);
631: *cellIS = *faceIS = NULL;
632: if (depth > 0) {
633: PetscCall(DMLabelGetStratumIS(label, depth, cellIS));
634: PetscCall(DMLabelGetStratumIS(label, depth - 1, faceIS));
635: }
636: if (!*cellIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cellIS));
637: if (!*faceIS) PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, faceIS));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@
642: DMPlexOrientLabel - Give a consistent orientation to the hypersurface marked by the `DMLabel` in the input mesh
644: Collective on dm
646: Input Parameters:
647: + dm - The `DM`
648: - label - The `DMLabel`
650: Notes:
651: The orientation data for the `DM` are changed in-place.
653: This routine will fail for non-orientable surfaces, such as the Moebius strip.
655: Level: advanced
657: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexOrient()`
658: @*/
659: PetscErrorCode DMPlexOrientLabel(DM dm, DMLabel label)
660: {
661: IS cellIS, faceIS;
663: PetscFunctionBegin;
664: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
665: PetscCall(DMPlexOrientCells_Internal(dm, cellIS, faceIS));
666: PetscCall(ISDestroy(&cellIS));
667: PetscCall(ISDestroy(&faceIS));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: // Every pair of samples of a face must match
672: static PetscErrorCode DMPlexCheckOrientation_Internal(DM dm, IS cellIS, IS faceIS)
673: {
674: const PetscInt debug = ((DM_Plex *)dm->data)->printOrient;
675: const PetscInt *faces = NULL;
676: PetscInt fStart = 0, fEnd = 0, *off;
677: FaceSample *samples;
678: PetscBool valid = PETSC_TRUE;
680: PetscFunctionBegin;
681: if (faceIS) PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
682: PetscCall(DMPlexGetFaceSamples_Private(dm, cellIS, faceIS, NULL, NULL, &off, &samples));
683: for (PetscInt f = 0; f < fEnd - fStart; ++f) {
684: for (PetscInt i = off[f]; i < off[f + 1]; ++i) {
685: for (PetscInt j = i + 1; j < off[f + 1]; ++j) {
686: if (!samples[i].ornt || !samples[j].ornt || FaceSamplesMatch(&samples[i], &samples[j])) continue;
687: valid = PETSC_FALSE;
688: if (debug)
689: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Face %" PetscInt_FMT " is mismatched: cell (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ") ~ cell (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ")\n", faces ? faces[fStart + f] : fStart + f,
690: samples[i].owner.rank, samples[i].owner.index, samples[i].ornt, samples[j].owner.rank, samples[j].owner.index, samples[j].ornt));
691: }
692: }
693: }
694: PetscCall(PetscFree(off));
695: PetscCall(PetscFree(samples));
696: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &valid, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)dm)));
697: PetscCheck(valid, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Mesh was not properly oriented");
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: /*@
702: DMPlexCheckOrientationLabel - Check that the surface defined by the given `DMLabel` is oriented
704: Collective
706: Input Parameters:
707: + dm - The `DM`
708: - label - The `DMLabel` defining the embedded surface
710: Level: advanced
712: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexOrient()`
713: @*/
714: PetscErrorCode DMPlexCheckOrientationLabel(DM dm, DMLabel label)
715: {
716: IS cellIS, faceIS;
718: PetscFunctionBegin;
719: PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS));
720: PetscCall(DMPlexCheckOrientation_Internal(dm, cellIS, faceIS));
721: PetscCall(ISDestroy(&cellIS));
722: PetscCall(ISDestroy(&faceIS));
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }