From f3a1137d2577956f5401af311a335e2297e8301a Mon Sep 17 00:00:00 2001 From: Mark Olesen <Mark.Olesen@esi-group.com> Date: Fri, 1 Oct 2021 18:26:49 +0200 Subject: [PATCH] BUG: parallel construct finiteArea fails with arbitrary connections (#2152) - the case of 'fan-like' processor was previously assumed to be rare (see merge-request !487 and issue #2084). However, Vaggelis Papoutsis noticed that even fairly normal geometries can trigger problems. - replaced the old patch/patch matching style with a more general edge-based synchronisation and matching that appears to handle the corner cases inherently. The internal communication overhead is essentially unchanged, and the logic is simpler. --- src/finiteArea/faMesh/faMesh.H | 20 +- src/finiteArea/faMesh/faMeshPatches.C | 1285 ++++++++++------- .../faMesh/faPatches/faPatch/faPatchData.C | 19 + .../faMesh/faPatches/faPatch/faPatchData.H | 3 + 4 files changed, 772 insertions(+), 555 deletions(-) diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index 8a9b1650118..16aa1d16843 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -197,6 +197,13 @@ class faMesh static const int quadricsFit_; + // Data Types + + //- Internal class for managing patch/patch bookkeeping + //- during construction + class patchTuple; + + // Private Member Functions //- No copy construct @@ -268,8 +275,9 @@ class faMesh // Helpers - //- Get the polyPatch pairs for the boundary edges (natural order) - List<LabelledItem<edge>> getBoundaryEdgePatchPairs() const; + //- Get list of (proc/patchi/patchEdgei) tuple pairs in an + //- globally consistent ordering + List<Pair<patchTuple>> getBoundaryEdgeConnections() const; //- Create a single patch PtrList<faPatch> createOnePatch @@ -286,14 +294,6 @@ class faMesh const dictionary* defaultPatchDefinition = nullptr ) const; - //- Reorder processor edges using order of the - //- neighbour processorPolyPatch - void reorderProcEdges - ( - faPatchData& patchDef, - const List<LabelledItem<edge>>& bndEdgePatchPairs - ) const; - public: diff --git a/src/finiteArea/faMesh/faMeshPatches.C b/src/finiteArea/faMesh/faMeshPatches.C index 436f3f5678e..b57f817e3ff 100644 --- a/src/finiteArea/faMesh/faMeshPatches.C +++ b/src/finiteArea/faMesh/faMeshPatches.C @@ -35,303 +35,140 @@ License #include "indirectPrimitivePatch.H" #include "edgeHashes.H" #include "LabelledItem.H" - -// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // - -namespace Foam +#include "foamVtkLineWriter.H" +#include "foamVtkIndPatchWriter.H" +#include "foamVtkUIndPatchWriter.H" + +// * * * * * * * * * * * * * * * Local Classes * * * * * * * * * * * * * * // + +// A tuple of (proc, patchi, patchEdgei) in that sort order +// Any finiteArea patches are stored with negated indices, which +// makes them readily identifiable and always sort before normal patches +class Foam::faMesh::patchTuple +: + public FixedList<label, 3> { +public: -// Manage patch pairs with a 'labelled' edge. -// The edge first/second correspond to the owner/neighbour patches. -// The index is a face index on the neighbour patch (FUTURE). - -// Local typedefs - -typedef LabelledItem<edge> patchPairInfo; -typedef List<patchPairInfo> patchPairInfoList; -typedef UIndirectList<patchPairInfo> patchPairInfoUIndList; - + // Constructors -// Handling of dangling coupled edges. -// Tag values to "push" with special -(patchId+2) -struct combineDanglingEdge -{ - const label upperLimit; - - // Set dangling patchId from real patchId - static void setDangling(patchPairInfo& pairing, const label patchId) - { - pairing.first() = pairing.second() = -(patchId + 2); - pairing.setIndex(-1); // Invalidate - } + // Inherit constructors + using FixedList<label, 3>::FixedList; - // Convert dangling patchId to real patchId - static void correct(patchPairInfo& pairing) - { - if (pairing.first() < -1) - { - pairing.first() = -(pairing.first() + 2); - } - if (pairing.second() < -1) + //- Default construct as 'invalid' + patchTuple() { - pairing.second() = -(pairing.second() + 2); + clear(); } - } - //- Construct with upper limit (the number of non-processor patches) - explicit combineDanglingEdge(const label nNonProcessor) - : - upperLimit(nNonProcessor) - {} + // Static Member Functions - // Combine operation: overwrite unused or processor patches with - // 'dangling' patch information only - void operator()(patchPairInfo& x, const patchPairInfo& y) const - { - if (y.first() < -1 && edge::compare(x, y) == 0) + // Globally consistent ordering + // + // 1. sort left/right as lower/higher processor connection + // 2. sort by proc/patch/patch index + static void sort(UList<Pair<patchTuple>>& list) { - if (x.first() == -1 || x.first() >= upperLimit) + for (auto& tuples : list) { - x.first() = y.first(); - } - if (x.second() == -1 || x.second() >= upperLimit) - { - x.second() = y.first(); - x.index() = y.index(); + tuples.sort(); } + Foam::stableSort(list); } - } -}; -// Populate patch pairings according to the boundary edges -void findEdgePatchPairing -( - const polyBoundaryMesh& pbm, - const EdgeMap<label>& edgeToIndex, - patchPairInfoList& patchPairs, - label nMissing = -1 -) -{ - // Count how many slots (both sides) to be filled - if (nMissing < 0) - { - nMissing = 0; - for (const patchPairInfo& pairing : patchPairs) + // Member Functions + + //- Reset to 'invalid' + void clear() { - if (pairing.first() == -1) ++nMissing; - if (pairing.second() == -1) ++nMissing; + procNo(-1); + patchi(labelMax); + patchEdgei(-1); } - } - forAll(pbm, patchID) - { - if (!nMissing) break; // Everything filled - - const polyPatch& pp = pbm[patchID]; - - const bool isProcPatch = isA<processorPolyPatch>(pp); - - // Examine neighbour boundary edges - for (label edgei = pp.nInternalEdges(); edgei < pp.nEdges(); ++edgei) + //- Valid if proc and edge are non-negative + bool valid() const noexcept { - // Lookup global edge of this neighbour, - // find matching owner boundary edge - - const label edgeIndex = - edgeToIndex.lookup(pp.meshEdge(edgei), -1); - - if (edgeIndex != -1) - { - // Add patchId. - // - hash-like so will only insert once - // - also add in the attached patch face (there is only one), - // saved in mesh face numbering for processor patches - - patchPairInfo& pairing = patchPairs[edgeIndex]; - - if (pairing.insert(patchID)) - { - if (isProcPatch) - { - // Save the mesh face - - const label patchFacei = pp.edgeFaces()[edgei][0]; - - pairing.index() = (pp.start() + patchFacei); - } - - --nMissing; - - if (!nMissing) break; // Early exit - } - } + return (procNo() >= 0 && patchEdgei() >= 0); } - } -} - -} // End namespace Foam - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::faMesh::reorderProcEdges -( - faPatchData& patchDef, - const List<LabelledItem<edge>>& bndEdgePatchPairs -) const -{ - if (!patchDef.coupled() || patchDef.edgeLabels_.empty()) - { - return; - } - - const polyBoundaryMesh& pbm = mesh().boundaryMesh(); - - const label procPatchID = patchDef.neighPolyPatchId_; - - const auto* procPatch = isA<processorPolyPatch>(pbm[procPatchID]); - - if (!procPatch) - { - FatalErrorInFunction - << "Internal addressing error. Patch " << procPatchID - << " is not a processor patch" << nl - << abort(FatalError); - } - - // Reorder processor edges using order of neighbour processorPolyPatch - - const label nProcEdges = patchDef.edgeLabels_.size(); - - labelList procFaces(nProcEdges, -1); - - forAll(procFaces, edgei) - { - const label bndEdgei = - (patchDef.edgeLabels_[edgei] - patch().nInternalEdges()); - - procFaces[edgei] = bndEdgePatchPairs[bndEdgei].index(); - } - - // Ascending proc-face numbering - const labelList sortIndices(Foam::sortedOrder(procFaces)); - - if (procFaces.size() && procFaces[sortIndices[0]] < 0) - { - FatalErrorInFunction - << "Internal addressing error. Patch " << procPatchID - << " with negative face" << nl - << abort(FatalError); - } - - const labelList& oldEdgeLabels = patchDef.edgeLabels_; - labelList newEdgeLabels(oldEdgeLabels.size()); + // Processor is the first sort index + label procNo() const { return (*this)[0]; } + void procNo(label val) { (*this)[0] = val; } - // Most of the time, an individual proc-face will only be singly - // attached to the finite-area patch. In rarer case, there could - // multiple connections. For these cases, need to walk the face - // edges - the direction depends on owner vs neighbour side. + // PatchId (-ve for finiteArea patches) is the second sort index + label patchi() const { return (*this)[1]; } + void patchi(label val) { (*this)[1] = val; } - EdgeMap<label> multihit; + // The patch edge index is the third sort index + label patchEdgei() const { return (*this)[2]; } + void patchEdgei(label val) { (*this)[2] = val; } - for (label edgei = 0; edgei < nProcEdges; /*nil*/) - { - const label meshFacei = procFaces[sortIndices[edgei]]; - - // Find all identical faces - label endEdgei = edgei + 1; // one beyond - while - ( - (endEdgei < nProcEdges) - && (meshFacei == procFaces[sortIndices[endEdgei]]) - ) + //- Return the real patchId + label realPatchi() const { - ++endEdgei; + const label id = patchi(); + return (id < 0 ? -(id + 1) : id); } - if (edgei + 1 == endEdgei) + //- Set patchId as finiteArea + void faPatchi(label val) { - // Simplest case - a single connection + patchi(-(val + 1)); + } - newEdgeLabels[edgei] = oldEdgeLabels[sortIndices[edgei]]; + //- Considered to be finiteArea if patchi < 0 + bool is_finiteArea() const noexcept + { + return (patchi() < 0); } - else + + //- Considered to be processor local + bool is_localProc() const noexcept { - multihit.clear(); + return (procNo() == Pstream::myProcNo()); + } +}; - // Map from global edge to patch local edgeId - for (label i = edgei; i < endEdgei; ++i) - { - const label patchEdgei = oldEdgeLabels[sortIndices[i]]; - // The edge in mesh numbering - multihit.insert - ( - patch().meshEdge(patchEdgei), - patchEdgei - ); - } +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // - if (multihit.size() != (endEdgei - edgei)) - { - FatalErrorInFunction - << "Could only hash " << multihit.size() - << " edges from " << (endEdgei - edgei) - << " ... indicates a non-manifold connection" << nl - << multihit << nl - << abort(FatalError); - } +namespace Foam +{ - const face& f = mesh().faces()[meshFacei]; +// Print out edges as point pairs +template<class PatchType> +static void printPatchEdges +( + Ostream& os, + const PatchType& p, + const labelList& edgeIds, + label maxOutput = 10 +) +{ + label nOutput = 0; - forAll(f, fedgei) // Note size() == nEdges() - { - edge e = - ( - patchDef.owner() - ? f.edge(fedgei) // Forward walk - : f.rcEdge(fedgei) // Reverse walk - ); + for (const label patchEdgei : edgeIds) + { + const edge e(p.meshEdge(patchEdgei)); - auto iter = multihit.find(e); - if (iter.found()) - { - newEdgeLabels[edgei++] = iter.val(); - multihit.erase(iter); - if (multihit.empty()) - { - break; - } - } - } + os << " " + << p.points()[e.first()] << ' ' + << p.points()[e.second()] << nl; - if (edgei != endEdgei) - { - FatalErrorInFunction - << "Missed " << (edgei < endEdgei) - << " edges for face: " << meshFacei - << " ... indicates serious geometry issue" << nl - << multihit << nl - << abort(FatalError); - } - if (!multihit.empty()) - { - FatalErrorInFunction - << "Missed edges for face: " << meshFacei - << " ... indicates serious geometry issue" << nl - << multihit << nl - << abort(FatalError); - } + ++nOutput; + if (maxOutput > 0 && nOutput >= maxOutput) + { + os << " ... suppressing further output" << nl; + break; } - edgei = endEdgei; } - - patchDef.edgeLabels_.transfer(newEdgeLabels); } +} // End namespace Foam + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -399,92 +236,272 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createOnePatch } -Foam::List<Foam::LabelledItem<Foam::edge>> -Foam::faMesh::getBoundaryEdgePatchPairs() const +Foam::List<Foam::Pair<Foam::faMesh::patchTuple>> +Foam::faMesh::getBoundaryEdgeConnections() const { const polyBoundaryMesh& pbm = mesh().boundaryMesh(); + const label nNonProcessor = pbm.nNonProcessor(); + const label nInternalEdges = patch().nInternalEdges(); const label nBoundaryEdges = patch().nBoundaryEdges(); + // The output result: + List<Pair<patchTuple>> bndEdgeConnections(nBoundaryEdges); + // Map edges (mesh numbering) back to a boundary index EdgeMap<label> edgeToBoundaryIndex(2*nBoundaryEdges); - // Use labelled 'edge' for accounting of patch pairs - patchPairInfoList bndEdgePatchPairs(nBoundaryEdges); + label nBadEdges(0); + labelHashSet badEdges(2*nBoundaryEdges); - - // Pass 1: - // - setup lookup (edge -> bnd index) - // - add owner patch for each boundary edge { - const SubList<labelList> bndEdgeToPatchFace - ( - patch().edgeFaces(), - patch().nBoundaryEdges(), - patch().nInternalEdges() - ); + // Use labelled 'edge' for local accounting of patch pairs + // It has hash-like insertion properties that are useful. + List<LabelledItem<edge>> patchPairings(nBoundaryEdges); + DebugInFunction + << "Determining required boundary edge connections, " + << "resolving locally attached boundary edges." << endl; + + + // Pass 1: + // - setup lookup (edge -> bnd index) + // - add owner patch for each boundary edge for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) { + const label patchEdgei = (bndEdgei + nInternalEdges); + edgeToBoundaryIndex.insert ( - patch().meshEdge(bndEdgei + nInternalEdges), - edgeToBoundaryIndex.size() + patch().meshEdge(patchEdgei), + bndEdgei ); - // The attached patch face (there is only one): - const label patchFacei = bndEdgeToPatchFace[bndEdgei][0]; + // The attached patch face. Should only be one! + const labelList& edgeFaces = patch().edgeFaces()[patchEdgei]; + + if (edgeFaces.size() != 1) + { + badEdges.insert(patchEdgei); + continue; + } + + const label patchFacei = edgeFaces[0]; + + /// const label meshFacei = faceLabels_[patchFacei]; + /// const label patchId = pbm.whichPatch(meshFacei); + + const label bndFacei = + (faceLabels_[patchFacei] - mesh().nInternalFaces()); + + const label patchId = pbm.patchID()[bndFacei]; - const label meshFacei = faceLabels_[patchFacei]; + // Primary bookkeeping + { + auto& tuple = bndEdgeConnections[bndEdgei].first(); + + tuple.procNo(Pstream::myProcNo()); + tuple.faPatchi(patchId); // Tag as finiteArea patch + tuple.patchEdgei(patchEdgei); + } - const label patchId = pbm.whichPatch(meshFacei); + // Temporary local bookkeeping + { + auto& pairing = patchPairings[bndEdgei]; - bndEdgePatchPairs[bndEdgei].insert(patchId); + pairing.clear(); // invalidate + pairing.setIndex(-1); + pairing.insert(patchId); // 'hash' into the first location + } } - } - // Pass 2: - // - Add in first neighbour patch for the boundary edges - // - examine all possible connecting neighbours + if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0) + { + edgeList dumpEdges(patch().edges(), badEdges.sortedToc()); - findEdgePatchPairing - ( - pbm, - edgeToBoundaryIndex, - bndEdgePatchPairs, - edgeToBoundaryIndex.size() // Number of places to still fill - ); + vtk::lineWriter writer + ( + patch().localPoints(), + dumpEdges, + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.nonManifoldEdges") + ) + ); + + writer.writeGeometry(); + + // CellData + writer.beginCellData(); + writer.writeProcIDs(); + + FatalErrorInFunction + << "Boundary edges not singly connected: " + << nBadEdges << '/' << nBoundaryEdges << nl; + + printPatchEdges + ( + FatalError, + patch(), + badEdges.sortedToc() + ); + + InfoInFunction + << "(debug) wrote " << writer.output().name() << nl; + + FatalError << abort(FatalError); + } + badEdges.clear(); + + // Pass 2: + // Add in first connecting neighbour patch for the boundary edges. + // Need to examine all possibly connecting (non-processor) neighbours, + // but only need to check their boundary edges. + label nMissing = patchPairings.size(); - // Nothing dangling if running in serial - can return already + for (label patchi = 0; patchi < nNonProcessor; ++patchi) + { + if (!nMissing) break; // Early exit + + const polyPatch& pp = pbm[patchi]; + + // Check boundary edges + for + ( + label patchEdgei = pp.nInternalEdges(); + patchEdgei < pp.nEdges(); + ++patchEdgei + ) + { + const label bndEdgei = + edgeToBoundaryIndex.lookup(pp.meshEdge(patchEdgei), -1); + + if (bndEdgei != -1) + { + // Has a matching owner boundary edge + + auto& pairing = patchPairings[bndEdgei]; + + // Add neighbour patchId/patchEdgei + // 'hash' into the second location + // which does not insert the same value twice + if (pairing.insert(patchi)) + { + // The neighbour's patch edge index too + pairing.setIndex(patchEdgei); + + --nMissing; + if (!nMissing) break; // Early exit + } + } + } + } + + + // Combine local bookkeeping into final list + badEdges.clear(); + for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) + { + const auto& pairing = patchPairings[bndEdgei]; + const label nbrPatchi = pairing.second(); + const label nbrPatchEdgei = pairing.index(); + + if (nbrPatchEdgei >= 0) + { + // Add into primary bookkeeping + auto& tuple = bndEdgeConnections[bndEdgei].second(); + + tuple.procNo(Pstream::myProcNo()); + tuple.patchi(nbrPatchi); + tuple.patchEdgei(nbrPatchEdgei); + } + else if (!Pstream::parRun()) + { + badEdges.insert(nInternalEdges + bndEdgei); + } + } + } + + + // ~~~~~~ + // Serial - can return already + // ~~~~~~ if (!Pstream::parRun()) { - return bndEdgePatchPairs; + // Verbose report of missing edges - in serial + + nBadEdges = badEdges.size(); + if (nBadEdges) + { + edgeList dumpEdges(patch().edges(), badEdges.sortedToc()); + + vtk::lineWriter writer + ( + patch().localPoints(), + dumpEdges, + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.invalidEdges") + ) + ); + + writer.writeGeometry(); + + // CellData + writer.beginCellData(); + writer.writeProcIDs(); + + FatalErrorInFunction + << "Boundary edges with missing/invalid neighbours: " + << nBadEdges << '/' << nBoundaryEdges << nl; + + printPatchEdges + ( + FatalError, + patch(), + badEdges.sortedToc() + ); + + InfoInFunction + << "(debug) wrote " << writer.output().name() << nl; + + FatalError << abort(FatalError); + } + + // Globally consistent ordering + patchTuple::sort(bndEdgeConnections); + return bndEdgeConnections; } - // In parallel need to check for "dangling" edges, which are finiteArea - // boundary edges that only exist on one side of a proc boundary. - // Eg, proc boundary coincides with the outer extent of the finiteArea + + // ~~~~~~~~ + // Parallel + // ~~~~~~~~ + + DebugInFunction + << "Creating global coupling data" << endl; const globalMeshData& globalData = mesh().globalData(); const indirectPrimitivePatch& cpp = globalData.coupledPatch(); const mapDistribute& map = globalData.globalEdgeSlavesMap(); + const label nCoupledEdges = cpp.nEdges(); // Construct coupled edge usage with all data - List<unsigned char> coupledEdgesUsed(map.constructSize(), 0u); + List<bool> coupledEdgesUsed(map.constructSize(), false); - forAll(cpp.edges(), coupledEdgei) + // Markup finiteArea boundary edges that are coupled across processors + for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei) { - const auto iter = - edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei)); - - // Used from this side or other other side - coupledEdgesUsed[coupledEdgei] = (iter.found() ? 1 : 2); + coupledEdgesUsed[cppEdgei] = + edgeToBoundaryIndex.found(cpp.meshEdge(cppEdgei)); } - // Save the original (pre-sync) coupling state - const List<unsigned char> coupledEdgesOrig(coupledEdgesUsed); + DebugInFunction + << "Starting sync of boundary edge topology" << endl; globalData.syncData ( @@ -492,180 +509,269 @@ Foam::faMesh::getBoundaryEdgePatchPairs() const globalData.globalEdgeSlaves(), globalData.globalEdgeTransformedSlaves(), // probably not used map, - bitOrEqOp<unsigned char>() + orEqOp<bool>() ); + if (debug) + { + label nAttached = 0; + for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei) + { + if (coupledEdgesUsed[cppEdgei]) + { + ++nAttached; + } + } + + InfoInFunction + << "Approx " + << returnReduce(nAttached, sumOp<label>()) + << " connected boundary edges (total, some duplicates)" << endl; + } - // Check for one-sided edge coupling (coupled value == 3) - // original == 1: - // - coupled to a real finiteArea edge. - // - receive a patch-pair value - // - // original == 2: - // - a "dangled" edge. Information required for other procs. - // - push a patch-pair value + // Combine information - // Map edges (mesh numbering) back to a coupled index. - // These are the edges to 'push' information for. + List<List<patchTuple>> gatheredConnections(map.constructSize()); - EdgeMap<label> edgeToCoupledIndex; + // Map edges (mesh numbering) back to a coupled index in use + EdgeMap<label> edgeToCoupledIndex(2*nCoupledEdges); - label nEdgesPull = 0; + // Pass 1 + // Look for attached boundary edges + // - boundary edges from this side go into gathered connections + // - boundary edges connected from the other side are noted for later - forAll(coupledEdgesUsed, coupledEdgei) + for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei) { - if (coupledEdgesUsed[coupledEdgei] == 3) + if (coupledEdgesUsed[cppEdgei]) { - if (coupledEdgesOrig[coupledEdgei] == 1) + const edge meshEdge(cpp.meshEdge(cppEdgei)); + + const label bndEdgei = + edgeToBoundaryIndex.lookup(meshEdge, -1); + + if (bndEdgei != -1) { - // Coupled side with finiteArea - ++nEdgesPull; + // A boundary finiteEdge edge (known from this side) + + auto& gathered = gatheredConnections[cppEdgei]; + gathered.resize(1); + auto& tuple = gathered.last(); + + tuple = bndEdgeConnections[bndEdgei].first(); } - else if (coupledEdgesOrig[coupledEdgei] == 2) + else { - // Coupled side without finiteArea - edgeToCoupledIndex.insert - ( - cpp.meshEdge(coupledEdgei), - coupledEdgei - ); + // Boundary edge connected from the other side + // - mark for it to be added later + + edgeToCoupledIndex.insert(meshEdge, cppEdgei); } } } - // Nothing to do - can return already - if (returnReduce(edgeToCoupledIndex.empty(), andOp<bool>())) - { - return bndEdgePatchPairs; - } - - // Data locations to pull - labelList patchEdgeLabels(nEdgesPull); - labelList coupledEdgeLabels(nEdgesPull); + // Pass 2 + // - add previously noted boundary edges (connected from other side) + // into gathered connections - // Populate the locations + for (label patchi = 0; patchi < nNonProcessor; ++patchi) { - nEdgesPull = 0; + if (edgeToCoupledIndex.empty()) break; // Early exit + + const polyPatch& pp = pbm[patchi]; - forAll(cpp.edges(), coupledEdgei) + // Check boundary edges + for + ( + label patchEdgei = pp.nInternalEdges(); + patchEdgei < pp.nEdges(); + ++patchEdgei + ) { - if - ( - coupledEdgesUsed[coupledEdgei] == 3 - && coupledEdgesOrig[coupledEdgei] == 1 - ) + const edge meshEdge(pp.meshEdge(patchEdgei)); + + const label cppEdgei = + edgeToCoupledIndex.lookup(meshEdge, -1); + + if (cppEdgei != -1) { - // Pull this edge - const auto iter = - edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei)); + // A known connection - if (iter.found()) - { - patchEdgeLabels[nEdgesPull] = iter.val(); - coupledEdgeLabels[nEdgesPull] = coupledEdgei; - ++nEdgesPull; - } - else - { - // Should be impossible to fail here - FatalErrorInFunction - << "Failed on second lookup of " - << cpp.meshEdge(coupledEdgei) << nl - << abort(FatalError); - } - } - } + auto& gathered = gatheredConnections[cppEdgei]; + gathered.resize(1); + auto& tuple = gathered.last(); - if (nEdgesPull != coupledEdgeLabels.size()) - { - FatalErrorInFunction - << "Failed lookup of some coupled edges" << nl - << abort(FatalError); + tuple.procNo(Pstream::myProcNo()); + tuple.patchi(patchi); + tuple.patchEdgei(patchEdgei); + + // Do not consider again + edgeToCoupledIndex.erase(meshEdge); + + if (edgeToCoupledIndex.empty()) break; // Early exit + } } } - //- Construct edge sync with all data - patchPairInfoList cppEdgeData(map.constructSize()); + DebugInFunction + << "Starting sync of boundary edge information" << endl; - // Fill in for 'push' locations. Only really interested in the owner - // (corresponds to non-proc connection), but grab everything - findEdgePatchPairing + globalMeshData::syncData ( - pbm, - edgeToCoupledIndex, - cppEdgeData, - 2*edgeToCoupledIndex.size() // Accept both sides? + gatheredConnections, + globalData.globalEdgeSlaves(), + globalData.globalEdgeTransformedSlaves(), // probably not used + map, + ListOps::appendEqOp<patchTuple>() ); - const label nNonProcessor = pbm.nNonProcessor(); + DebugInFunction + << "Collating sync information" << endl; - // Adjust patch information to reflect dangling patch neighbour - // Tag with -(value+2) - forAllConstIters(edgeToCoupledIndex, iter) + // Pick out gathered connections and add into primary bookkeeping + for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei) { - const edge& e = iter.key(); - const label coupledEdgei = iter.val(); + const auto& gathered = gatheredConnections[cppEdgei]; - patchPairInfo& pairing = cppEdgeData[coupledEdgei]; - const label ownerPatchId = pairing.first(); + const label bndEdgei = + edgeToBoundaryIndex.lookup(cpp.meshEdge(cppEdgei), -1); - // Some sanity checks - if (ownerPatchId < 0) + if + ( + // A boundary finiteEdge edge (known from this side) + bndEdgei != -1 + + // Gathered a one-to-one connection + && gathered.size() == 2 + ) { - FatalErrorInFunction - << "Error finding dangling edge at " - << cpp.points()[e.first()] << ' ' - << cpp.points()[e.second()] << nl - << abort(FatalError); + const auto& a = gathered[0]; + const auto& b = gathered[1]; + + // Copy second side of connection + auto& connection = bndEdgeConnections[bndEdgei]; + + connection.second() = (connection.first() == b) ? a : b; } - else if (ownerPatchId >= nNonProcessor) + } + + // Check missing/invalid + badEdges.clear(); + for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) + { + const auto& connection = bndEdgeConnections[bndEdgei]; + + if (!connection.second().valid()) { - FatalErrorInFunction - << "Cannot handle edge on processor-processor connection at " - << cpp.points()[e.first()] << ' ' - << cpp.points()[e.second()] << nl - << abort(FatalError); + badEdges.insert(nInternalEdges + bndEdgei); } + } - combineDanglingEdge::setDangling(pairing, ownerPatchId); - // TBD: - // may wish to remember the corresponding proc number, - // if we wish to bridge across 'fan-like' connections. - // - // pairing.setIndex(-(Pstream::myProcNo() + 2)); - } + if (debug & 8) + { + // Boundary edges + { + vtk::lineWriter writer + ( + patch().localPoints(), + patch().boundaryEdges(), + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.boundaryEdges") + ) + ); + writer.writeGeometry(); - // Synchronize edge information + // CellData + writer.beginCellData(); + writer.writeProcIDs(); - const combineDanglingEdge edgeCombineOp(nNonProcessor); + // For each boundary edge - the associate neighbour patch + labelList neighProc(nBoundaryEdges); + labelList neighPatch(nBoundaryEdges); + for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) + { + const auto& connection = bndEdgeConnections[bndEdgei]; - globalMeshData::syncData - ( - cppEdgeData, - globalData.globalEdgeSlaves(), - globalData.globalEdgeTransformedSlaves(), // probably not used - map, - edgeCombineOp - ); + neighProc[bndEdgei] = connection.second().procNo(); + neighPatch[bndEdgei] = connection.second().patchi(); + } + writer.write("neighProc", neighProc); + writer.write("neighPatch", neighPatch); + } - // Combine back from pushed cpp-edge data + // finiteArea + { + vtk::uindirectPatchWriter writer + ( + patch(), + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.faPatch") + ) + ); + + writer.writeGeometry(); - forAll(patchEdgeLabels, i) + // CellData + writer.beginCellData(); + writer.writeProcIDs(); + } + } + + // Verbose report of missing edges + if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0) { - patchPairInfo& pairing = bndEdgePatchPairs[patchEdgeLabels[i]]; - const patchPairInfo& other = cppEdgeData[coupledEdgeLabels[i]]; + edgeList dumpEdges(patch().edges(), badEdges.sortedToc()); + + vtk::lineWriter writer + ( + patch().localPoints(), + dumpEdges, + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.invalidEdges") + ) + ); + + writer.writeGeometry(); + + // CellData + writer.beginCellData(); + writer.writeProcIDs(); + + FatalErrorInFunction + << "Boundary edges with missing/invalid neighbours: " + << nBadEdges << '/' << nBoundaryEdges << nl; + + printPatchEdges + ( + FatalError, + patch(), + badEdges.sortedToc() + ); - edgeCombineOp(pairing, other); + InfoInFunction + << "(debug) wrote " << writer.output().name() << nl; - // Resolve special tagging - combineDanglingEdge::correct(pairing); + FatalError << abort(FatalError); } - return bndEdgePatchPairs; + + // Globally consistent ordering + patchTuple::sort(bndEdgeConnections); + + DebugInFunction + << "Return sorted list of boundary connections" << endl; + + return bndEdgeConnections; } @@ -679,7 +785,7 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList const polyBoundaryMesh& pbm = mesh().boundaryMesh(); // Transcribe into patch definitions - DynamicList<faPatchData> faPatchDefs(bndDict.size() + 4); + DynamicList<faPatchData> faPatchDefs(bndDict.size() + 8); for (const entry& dEntry : bndDict) { if (!dEntry.isDict()) @@ -735,7 +841,7 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList } // Placeholder for any undefined edges - const label undefPatchId = faPatchDefs.size(); + const label undefPatchIndex = faPatchDefs.size(); { faPatchDefs.append(faPatchData()); @@ -752,192 +858,264 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList // ---------------------------------------------------------------------- - const label nInternalEdges = patch().nInternalEdges(); - const label nBoundaryEdges = patch().nBoundaryEdges(); + // Get edge connections (in globally consistent ordering) + List<Pair<patchTuple>> bndEdgeConnections + ( + this->getBoundaryEdgeConnections() + ); - patchPairInfoList bndEdgePatchPairs(this->getBoundaryEdgePatchPairs()); + // Lookup of patchDef for each connection. Initially all -1 (unvisited) + labelList patchDefLookup(bndEdgeConnections.size(), -1); - labelList bndEdgeFaPatchIDs(nBoundaryEdges, -1); + Map<labelHashSet> procConnections; + labelHashSet patchDefsUsed; - for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) - { - const patchPairInfo& patchPair = bndEdgePatchPairs[bndEdgei]; + label nBadEdges(0); + labelHashSet badEdges(2*bndEdgeConnections.size()); - // Find first definition with a matching neighbour and - // possibly with a matching owner. + forAll(bndEdgeConnections, connecti) + { + const Pair<patchTuple>& connection = bndEdgeConnections[connecti]; + const auto& a = connection.first(); + const auto& b = connection.second(); - label bestPatchi = -1; + edge patchPair; - for (label patchi = 0; patchi < faPatchDefs.size(); ++patchi) + if (a.is_finiteArea()) { - const int match = faPatchDefs[patchi].matchPatchPair(patchPair); - if (match == 3) + if (b.is_finiteArea()) { - // Match (owner/neighbour) - done! - bestPatchi = patchi; - break; + // A processor-processor connection + if (a.procNo() == b.procNo()) + { + FatalErrorInFunction + << "Processor-processor addressing error:" << nl + << "Both connections have the same processor: " + << a.procNo() << nl + << abort(FatalError); + } + else if (a.is_localProc()) + { + procConnections(b.procNo()).insert(connecti); + } + else + { + procConnections(a.procNo()).insert(connecti); + } + + continue; } - else if (match == 2 && bestPatchi < 0) + else if (a.is_localProc()) { - // Match (neighbour) - keep looking for exact match - bestPatchi = patchi; + patchPair.first() = a.realPatchi(); + patchPair.second() = b.realPatchi(); } } - - bndEdgeFaPatchIDs[bndEdgei] = bestPatchi; - } + else if (b.is_finiteArea() && b.is_localProc()) + { + patchPair.first() = b.realPatchi(); + patchPair.second() = a.realPatchi(); + } - // Extract which edges map to which patch - // and set edgeLabels for each faPatch - - DynamicList<label> selectEdges(bndEdgeFaPatchIDs.size()); + // Find first definition with a matching neighbour and + // possibly with a matching owner. - for (label patchi = 0; patchi < faPatchDefs.size(); ++patchi) - { - auto& patchDef = faPatchDefs[patchi]; + label bestPatchDefi = -1; - selectEdges.clear(); + const label nPatchDefs = (patchPair.valid() ? faPatchDefs.size() : 0); - forAll(bndEdgeFaPatchIDs, bndEdgei) + for (label patchDefi = 0; patchDefi < nPatchDefs; ++patchDefi) { - if (bndEdgeFaPatchIDs[bndEdgei] == patchi) + const int match = faPatchDefs[patchDefi].matchPatchPair(patchPair); + if (match == 3) { - selectEdges.append(bndEdgei + nInternalEdges); + // Exact match (owner/neighbour) - done! + bestPatchDefi = patchDefi; + break; + } + else if (match == 2 && bestPatchDefi < 0) + { + // Match (neighbour) - keep looking for exact match + bestPatchDefi = patchDefi; } } - patchDef.edgeLabels_ = selectEdges; - } - - // Check for undefined edges - selectEdges.clear(); - - forAll(bndEdgeFaPatchIDs, bndEdgei) - { - if (bndEdgeFaPatchIDs[bndEdgei] == -1) + if (bestPatchDefi < 0) { - selectEdges.append(bndEdgei + nInternalEdges); + bestPatchDefi = undefPatchIndex; // Missed? } - } - - // Save the information - faPatchDefs[undefPatchId].edgeLabels_ = selectEdges; - bool hasUndefined = returnReduce(!selectEdges.empty(), orOp<bool>()); + patchDefLookup[connecti] = bestPatchDefi; + patchDefsUsed.insert(bestPatchDefi); + } - if (hasUndefined) + // Remove undefPatchIndex if not actually needed + if (!returnReduce(patchDefsUsed.found(undefPatchIndex), orOp<bool>())) { - // The initial edges to consider - const labelList& undefinedEdges = - faPatchDefs[undefPatchId].edgeLabels_; + faPatchDefs.remove(undefPatchIndex); + } + else + { + patchDefsUsed.insert(undefPatchIndex); // Parallel consistency - // Check for edges that butt against a processor (or other) patch + // Report locations of undefined edges - labelList edgeNbrPolyPatch(undefinedEdges.size(), -1); - forAll(edgeNbrPolyPatch, edgei) + badEdges.clear(); + forAll(patchDefLookup, connecti) { - const label patchEdgei = undefinedEdges[edgei]; - const label bndEdgei = (patchEdgei - nInternalEdges); - - edgeNbrPolyPatch[edgei] = bndEdgePatchPairs[bndEdgei].second(); - } + if (patchDefLookup[connecti] == undefPatchIndex) + { + const auto& connection = bndEdgeConnections[connecti]; - // Categorize as processor/non-processor associations - labelHashSet procPatchIDs; - labelHashSet nonProcPatchIDs; + const auto& a = connection.first(); + const auto& b = connection.second(); - for (const label polyPatchID : edgeNbrPolyPatch) - { - if (polyPatchID == -1) - { - nonProcPatchIDs.insert(polyPatchID); - } - else if - ( - !nonProcPatchIDs.found(polyPatchID) - && !procPatchIDs.found(polyPatchID) - ) - { - if (isA<processorPolyPatch>(pbm[polyPatchID])) + if (a.is_localProc() && a.is_finiteArea()) { - procPatchIDs.insert(polyPatchID); + badEdges.insert(a.patchEdgei()); } - else + else if (b.is_localProc() && b.is_finiteArea()) { - nonProcPatchIDs.insert(polyPatchID); + badEdges.insert(b.patchEdgei()); } + + Pout<< "Undefined connection: " << flatOutput(connection) + << " connects: " + << pbm[a.realPatchi()].name() << " to " + << pbm[b.realPatchi()].name() << nl; } } - // Select by processor association - for (const label polyPatchID : procPatchIDs.sortedToc()) + if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0) { - selectEdges.clear(); + edgeList dumpEdges(patch().edges(), badEdges.sortedToc()); - forAll(edgeNbrPolyPatch, edgei) - { - if (edgeNbrPolyPatch[edgei] == polyPatchID) - { - selectEdges.append(undefinedEdges[edgei]); - } - } + vtk::lineWriter writer + ( + patch().localPoints(), + dumpEdges, + fileName + ( + mesh().time().globalPath() + / ("faMesh-construct.undefEdges") + ) + ); - faPatchDefs.append(faPatchData()); + writer.writeGeometry(); - auto& patchDef = faPatchDefs.last(); - patchDef.name_ = pbm[polyPatchID].name(); - patchDef.type_ = processorFaPatch::typeName; - patchDef.neighPolyPatchId_ = polyPatchID; // Needed for reorder + // CellData + writer.beginCellData(); + writer.writeProcIDs(); - const auto* ppp = isA<processorPolyPatch>(pbm[polyPatchID]); - if (ppp) - { - patchDef.ownerProcId_ = ppp->myProcNo(); - patchDef.neighProcId_ = ppp->neighbProcNo(); - } + InfoInFunction + << "Had undefined edge connections: " + << nBadEdges << '/' << patch().nBoundaryEdges() << nl + << "(debug) wrote " << writer.output().name() << nl; + } + } + + // Create processor-processor definitions + Map<label> procToDefLookup(2*procConnections.size()); + { + faPatchDefs.reserve(faPatchDefs.size() + procConnections.size()); + + for (const label otherProci : procConnections.sortedToc()) + { + const label patchDefi = faPatchDefs.size(); + procToDefLookup.insert(otherProci, patchDefi); + + // Add entry + faPatchDefs.append(faPatchData()); + auto& patchDef = faPatchDefs.last(); - patchDef.edgeLabels_ = selectEdges; + patchDef.assign_coupled(Pstream::myProcNo(), otherProci); } + } + + + // Extract which edges map to which patch + // and set edgeLabels for each faPatch + DynamicList<label> selectEdges(bndEdgeConnections.size()); - // Check for any remaining undefined edges + for (const label patchDefi : patchDefsUsed.sortedToc()) + { + auto& patchDef = faPatchDefs[patchDefi]; selectEdges.clear(); - // Simply grab any/all (don't worry about which patch) - if (!nonProcPatchIDs.empty()) + // Find the corresponding entries + // and extract the patchEdgeId + + forAll(patchDefLookup, connecti) { - forAll(edgeNbrPolyPatch, edgei) + if (patchDefLookup[connecti] == patchDefi) { - const label polyPatchID = edgeNbrPolyPatch[edgei]; - if (nonProcPatchIDs.found(polyPatchID)) + const auto& a = bndEdgeConnections[connecti].first(); + const auto& b = bndEdgeConnections[connecti].second(); + + if (a.is_localProc() && a.is_finiteArea()) { - selectEdges.append(undefinedEdges[edgei]); + selectEdges.append(a.patchEdgei()); } + else if (b.is_localProc() && b.is_finiteArea()) + { + selectEdges.append(b.patchEdgei()); + } + else + { + FatalErrorInFunction + << "Error in programming logic" << nl + << abort(FatalError); + } + + patchDefLookup[connecti] = -2; // Mark as handled } } - // Complete the information - faPatchDefs[undefPatchId].edgeLabels_ = selectEdges; - - hasUndefined = returnReduce(!selectEdges.empty(), orOp<bool>()); + patchDef.edgeLabels_ = selectEdges; } - // Remove unnecessary entry - if (!hasUndefined) - { - faPatchDefs.remove(undefPatchId); - } - for (auto& patchDef : faPatchDefs) + // Processor patches + for (const label otherProci : procToDefLookup.sortedToc()) { - if (patchDef.coupled()) + const label patchDefi = procToDefLookup[otherProci]; + + auto& patchDef = faPatchDefs[patchDefi]; + selectEdges.clear(); + + // Find the corresponding entries + // and extract the patchEdgeId + + for (const label connecti : procConnections(otherProci).sortedToc()) { - reorderProcEdges(patchDef, bndEdgePatchPairs); - patchDef.neighPolyPatchId_ = -1; // No lookup of neighbour faces + const auto& connection = bndEdgeConnections[connecti]; + const auto& a = connection.first(); + const auto& b = connection.second(); + + if (a.is_localProc()) + { + selectEdges.append(a.patchEdgei()); + } + else if (b.is_localProc()) + { + selectEdges.append(b.patchEdgei()); + } + else + { + FatalErrorInFunction + << "Error in programming logic" << nl + << abort(FatalError); + } + + patchDefLookup[connecti] = -2; // Mark as handled } + + patchDef.edgeLabels_ = selectEdges; } + // Now convert list of definitions to list of patches label nPatches = 0; @@ -962,6 +1140,23 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList ++nPatches; } + + if (debug > 1) + { + label nTotal = 0; + Pout<< "Created new finiteArea patches:" << nl; + for (const faPatch& p : newPatches) + { + Pout<< " size: " << p.size() + << " name:" << p.name() + << " type:" << p.type() << nl; + nTotal += p.size(); + } + + Pout<< "addressed: " << nTotal + << '/' << patch().nBoundaryEdges() << " edges" << endl; + } + return newPatches; } diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C index c83150cea4c..00ad7d679cb 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C @@ -28,6 +28,7 @@ License #include "faPatchData.H" #include "dictionary.H" #include "processorFaPatch.H" +#include "processorPolyPatch.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -108,6 +109,24 @@ void Foam::faPatchData::assign(const faPatch& fap) } +bool Foam::faPatchData::assign_coupled(int ownProci, int neiProci) +{ + clear(); + + if (ownProci == neiProci) + { + return false; + } + + name_ = processorPolyPatch::newName(ownProci, neiProci); + type_ = processorFaPatch::typeName; + ownerProcId_ = ownProci; + neighProcId_ = neiProci; + + return true; +} + + int Foam::faPatchData::matchPatchPair ( const labelPair& patchPair diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H index e66c91723e6..ee90a80522a 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H @@ -98,6 +98,9 @@ public: //- Clear and populate with values from finiteArea patch void assign(const faPatch& fap); + //- Set values consistent with a processor coupling + bool assign_coupled(int ownProci, int neiProci); + //- True if owner/neighbour processor ids are non-equal bool coupled() const noexcept { -- GitLab