diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index 9498c7ade34cfd1992a2a02d5ba30456585b5886..8a9b16501181a82973c33de49a4e98f833713a46 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -269,10 +269,7 @@ class faMesh // Helpers //- Get the polyPatch pairs for the boundary edges (natural order) - List<LabelledItem<edge>> getBoundaryEdgePatchPairs - ( - const labelUList& meshEdges - ) const; + List<LabelledItem<edge>> getBoundaryEdgePatchPairs() const; //- Create a single patch PtrList<faPatch> createOnePatch @@ -294,7 +291,6 @@ class faMesh void reorderProcEdges ( faPatchData& patchDef, - const labelUList& meshEdges, const List<LabelledItem<edge>>& bndEdgePatchPairs ) const; diff --git a/src/finiteArea/faMesh/faMeshPatches.C b/src/finiteArea/faMesh/faMeshPatches.C index 966f58299402e7c706a555c8b52ecc97c119ba4e..889bdd7826c301813a9eb69d7d4b2aa3fb7c023d 100644 --- a/src/finiteArea/faMesh/faMeshPatches.C +++ b/src/finiteArea/faMesh/faMeshPatches.C @@ -52,52 +52,124 @@ typedef List<patchPairInfo> patchPairInfoList; typedef UIndirectList<patchPairInfo> patchPairInfoUIndList; -// Synchronize edge patch pairs. -// - only propagate real (non-processor) patch ids - -struct syncEdgePatchPairs +// Handling of dangling coupled edges. +// Tag values to "push" with special -(patchId+2) +struct combineDanglingEdge { const label upperLimit; - explicit syncEdgePatchPairs(const label nNonProcessor) + // Set dangling patchId from real patchId + static void setDangling(patchPairInfo& pairing, const label patchId) + { + pairing.first() = pairing.second() = -(patchId + 2); + pairing.setIndex(-1); // Invalidate + } + + // Convert dangling patchId to real patchId + static void correct(patchPairInfo& pairing) + { + if (pairing.first() < -1) + { + pairing.first() = -(pairing.first() + 2); + } + if (pairing.second() < -1) + { + pairing.second() = -(pairing.second() + 2); + } + } + + //- Construct with upper limit (the number of non-processor patches) + explicit combineDanglingEdge(const label nNonProcessor) : upperLimit(nNonProcessor) {} - void insert(edge& e, const label i) const + + // Combine operation: overwrite unused or processor patches with + // 'dangling' patch information only + void operator()(patchPairInfo& x, const patchPairInfo& y) const { - // This could probably be simpler - if (i >= 0 && i < upperLimit && !e.found(i)) + if (y.first() < -1 && edge::compare(x, y) == 0) { - if (e.first() == -1) - { - e.first() = i; - } - else if (e.second() == -1) + if (x.first() == -1 || x.first() >= upperLimit) { - e.second() = i; + x.first() = y.first(); } - else if (upperLimit < e.first()) + if (x.second() == -1 || x.second() >= upperLimit) { - e.first() = i; - } - else if (upperLimit < e.second()) - { - e.second() = i; + x.second() = y.first(); + x.index() = y.index(); } } } +}; - void operator()(edge& x, const edge& y) const + +// 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) { - if (edge::compare(x, y) == 0) + nMissing = 0; + for (const patchPairInfo& pairing : patchPairs) { - insert(x, y.first()); - insert(x, y.second()); + if (pairing.first() == -1) ++nMissing; + if (pairing.second() == -1) ++nMissing; } } -}; + 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) + { + // 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 + } + } + } + } +} } // End namespace Foam @@ -107,7 +179,6 @@ struct syncEdgePatchPairs void Foam::faMesh::reorderProcEdges ( faPatchData& patchDef, - const labelUList& meshEdges, const List<LabelledItem<edge>>& bndEdgePatchPairs ) const { @@ -117,38 +188,45 @@ void Foam::faMesh::reorderProcEdges } const polyBoundaryMesh& pbm = mesh().boundaryMesh(); - const labelListList& edgeFaces = mesh().edgeFaces(); + + 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 procPatchID = patchDef.neighPolyPatchId_; const label nProcEdges = patchDef.edgeLabels_.size(); labelList procFaces(nProcEdges, -1); forAll(procFaces, edgei) { - const label localEdgei = patchDef.edgeLabels_[edgei]; - const label meshEdgei = meshEdges[localEdgei]; + const label bndEdgei = + (patchDef.edgeLabels_[edgei] - patch().nInternalEdges()); - for (const label meshFacei : edgeFaces[meshEdgei]) - { - if - ( - !faceLabels_.found(meshFacei) - && (procPatchID == pbm.whichPatch(meshFacei)) - ) - { - // The edge's proc-face - procFaces[edgei] = meshFacei; - break; - } - } + 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()); @@ -161,14 +239,14 @@ void Foam::faMesh::reorderProcEdges for (label edgei = 0; edgei < nProcEdges; /*nil*/) { - const label procFacei = procFaces[sortIndices[edgei]]; + const label meshFacei = procFaces[sortIndices[edgei]]; // Find all identical faces label endEdgei = edgei + 1; // one beyond while ( (endEdgei < nProcEdges) - && (procFacei == procFaces[sortIndices[endEdgei]]) + && (meshFacei == procFaces[sortIndices[endEdgei]]) ) { ++endEdgei; @@ -184,16 +262,16 @@ void Foam::faMesh::reorderProcEdges { multihit.clear(); - // Map from global edge to local edgeId + // Map from global edge to patch local edgeId for (label i = edgei; i < endEdgei; ++i) { - label localEdgei = oldEdgeLabels[sortIndices[i]]; - label meshEdgei = meshEdges[localEdgei]; + const label patchEdgei = oldEdgeLabels[sortIndices[i]]; + // The edge in mesh numbering multihit.insert ( - mesh().edges()[meshEdgei], - localEdgei + patch().meshEdge(patchEdgei), + patchEdgei ); } @@ -207,7 +285,7 @@ void Foam::faMesh::reorderProcEdges << abort(FatalError); } - const face& f = mesh().faces()[procFacei]; + const face& f = mesh().faces()[meshFacei]; forAll(f, fedgei) // Note size() == nEdges() { @@ -234,7 +312,7 @@ void Foam::faMesh::reorderProcEdges { FatalErrorInFunction << "Missed " << (edgei < endEdgei) - << " edges for face: " << procFacei + << " edges for face: " << meshFacei << " ... indicates serious geometry issue" << nl << multihit << nl << abort(FatalError); @@ -242,7 +320,7 @@ void Foam::faMesh::reorderProcEdges if (!multihit.empty()) { FatalErrorInFunction - << "Missed edges for face: " << procFacei + << "Missed edges for face: " << meshFacei << " ... indicates serious geometry issue" << nl << multihit << nl << abort(FatalError); @@ -322,13 +400,9 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createOnePatch Foam::List<Foam::LabelledItem<Foam::edge>> -Foam::faMesh::getBoundaryEdgePatchPairs -( - const labelUList& meshEdges -) const +Foam::faMesh::getBoundaryEdgePatchPairs() const { const polyBoundaryMesh& pbm = mesh().boundaryMesh(); - const labelListList& edgeFaces = mesh().edgeFaces(); const label nInternalEdges = patch().nInternalEdges(); const label nBoundaryEdges = patch().nBoundaryEdges(); @@ -336,155 +410,260 @@ Foam::faMesh::getBoundaryEdgePatchPairs // Map edges (mesh numbering) back to a boundary index EdgeMap<label> edgeToBoundaryIndex(2*nBoundaryEdges); - - // Use 'edge' for accounting + // Use labelled 'edge' for accounting of patch pairs patchPairInfoList bndEdgePatchPairs(nBoundaryEdges); - // Get pair of polyPatches for each boundary edge - for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) + + // Pass 1: + // - setup lookup (edge -> bnd index) + // - add owner patch for each boundary edge { - edgeToBoundaryIndex.insert + const SubList<labelList> bndEdgeToPatchFace ( - patch().meshEdge(bndEdgei + nInternalEdges), - edgeToBoundaryIndex.size() + patch().edgeFaces(), + patch().nBoundaryEdges(), + patch().nInternalEdges() ); - const label patchEdgei = (bndEdgei + nInternalEdges); - const label meshEdgei = meshEdges[patchEdgei]; + for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei) + { + edgeToBoundaryIndex.insert + ( + patch().meshEdge(bndEdgei + nInternalEdges), + edgeToBoundaryIndex.size() + ); - patchPairInfo& patchPair = bndEdgePatchPairs[bndEdgei]; + // The attached patch face (there is only one): + const label patchFacei = bndEdgeToPatchFace[bndEdgei][0]; + + const label meshFacei = faceLabels_[patchFacei]; - for (const label meshFacei : edgeFaces[meshEdgei]) - { const label patchId = pbm.whichPatch(meshFacei); - // Note: negative labels never insert - patchPair.insert(patchId); + bndEdgePatchPairs[bndEdgei].insert(patchId); } } + // Pass 2: + // - Add in first neighbour patch for the boundary edges + // - examine all possible connecting neighbours - // Synchronize edge information - we want the 'global' patch connectivity + findEdgePatchPairing + ( + pbm, + edgeToBoundaryIndex, + bndEdgePatchPairs, + edgeToBoundaryIndex.size() // Number of places to still fill + ); - // Looks like PatchTools::matchEdges - const indirectPrimitivePatch& cpp = mesh().globalData().coupledPatch(); + // Nothing dangling if running in serial - can return already + if (!Pstream::parRun()) + { + return bndEdgePatchPairs; + } - labelList patchEdgeLabels(nBoundaryEdges); - labelList coupledEdgeLabels(nBoundaryEdges); + // 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 - { - label nMatches = 0; - forAll(cpp.edges(), coupledEdgei) - { - const edge coupledMeshEdge(cpp.meshEdge(coupledEdgei)); + const globalMeshData& globalData = mesh().globalData(); + const indirectPrimitivePatch& cpp = globalData.coupledPatch(); + const mapDistribute& map = globalData.globalEdgeSlavesMap(); - const auto iter = edgeToBoundaryIndex.cfind(coupledMeshEdge); + // Construct coupled edge usage with all data + List<unsigned char> coupledEdgesUsed(map.constructSize(), 0u); - if (iter.found()) - { - patchEdgeLabels[nMatches] = iter.val(); - coupledEdgeLabels[nMatches] = coupledEdgei; - ++nMatches; - } - } + forAll(cpp.edges(), coupledEdgei) + { + const auto iter = + edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei)); - patchEdgeLabels.resize(nMatches); - coupledEdgeLabels.resize(nMatches); + // Used from this side or other other side + coupledEdgesUsed[coupledEdgei] = (iter.found() ? 1 : 2); } + // Save the original (pre-sync) coupling state + const List<unsigned char> coupledEdgesOrig(coupledEdgesUsed); - const globalMeshData& globalData = mesh().globalData(); - const mapDistribute& map = globalData.globalEdgeSlavesMap(); + globalData.syncData + ( + coupledEdgesUsed, + globalData.globalEdgeSlaves(), + globalData.globalEdgeTransformedSlaves(), // probably not used + map, + bitOrEqOp<unsigned char>() + ); - //- Construct with all data - patchPairInfoList cppEdgeData(map.constructSize()); - // Convert patch-edge data into cpp-edge data - patchPairInfoUIndList(cppEdgeData, coupledEdgeLabels) = - patchPairInfoUIndList(bndEdgePatchPairs, patchEdgeLabels); + // 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 + // Map edges (mesh numbering) back to a coupled index. + // These are the edges to 'push' information for. - // Also 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 end of the finiteArea + EdgeMap<label> edgeToCoupledIndex; - { - boolList edgeInUse(map.constructSize(), false); + label nEdgesPull = 0; - for (const label coupledEdgei : coupledEdgeLabels) + forAll(coupledEdgesUsed, coupledEdgei) + { + if (coupledEdgesUsed[coupledEdgei] == 3) { - edgeInUse[coupledEdgei] = true; + if (coupledEdgesOrig[coupledEdgei] == 1) + { + // Coupled side with finiteArea + ++nEdgesPull; + } + else if (coupledEdgesOrig[coupledEdgei] == 2) + { + // Coupled side without finiteArea + edgeToCoupledIndex.insert + ( + cpp.meshEdge(coupledEdgei), + coupledEdgei + ); + } } + } - // Retain pre-synchronized state for later xor - const boolList nonSyncEdgeInUse - ( - SubList<bool>(edgeInUse, cpp.nEdges()) - ); + // Nothing to do - can return already + if (returnReduce(edgeToCoupledIndex.empty(), andOp<bool>())) + { + return bndEdgePatchPairs; + } - globalData.syncData - ( - edgeInUse, - globalData.globalEdgeSlaves(), - globalData.globalEdgeTransformedSlaves(), - map, - orEqOp<bool>() - ); + // Data locations to pull + labelList patchEdgeLabels(nEdgesPull); + labelList coupledEdgeLabels(nEdgesPull); - // Check for anything coupled from the other side, - // but not originally from this. - // - // Process these dangling edges, obtain the attached pair - // of polyPatches for each. - // - // These edges have no finiteArea correspondence on this processor, - // but the information is obviously needed for other processors + // Populate the locations + { + nEdgesPull = 0; - forAll(nonSyncEdgeInUse, coupledEdgei) + forAll(cpp.edges(), coupledEdgei) { - // Coupled, but originating from elsewhere - if (edgeInUse[coupledEdgei] && !nonSyncEdgeInUse[coupledEdgei]) + if + ( + coupledEdgesUsed[coupledEdgei] == 3 + && coupledEdgesOrig[coupledEdgei] == 1 + ) { - // Already default initialized - patchPairInfo& patchPair = cppEdgeData[coupledEdgei]; - - const label meshEdgei = - cpp.meshEdge - ( - coupledEdgei, - mesh().edges(), - mesh().pointEdges() - ); - - for (const label meshFacei : edgeFaces[meshEdgei]) - { - const label patchId = pbm.whichPatch(meshFacei); + // Pull this edge + const auto iter = + edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei)); - // Note: negative labels never insert - patchPair.insert(patchId); + 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); } } } + + if (nEdgesPull != coupledEdgeLabels.size()) + { + FatalErrorInFunction + << "Failed lookup of some coupled edges" << nl + << abort(FatalError); + } } - // Convert patch-edge data into cpp-edge data - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + //- Construct edge sync with all data + patchPairInfoList cppEdgeData(map.constructSize()); - globalData.syncData + // Fill in for 'push' locations. Only really interested in the owner + // (corresponds to non-proc connection), but grab everything + findEdgePatchPairing + ( + pbm, + edgeToCoupledIndex, + cppEdgeData, + 2*edgeToCoupledIndex.size() // Accept both sides? + ); + + + const label nNonProcessor = pbm.nNonProcessor(); + + // Adjust patch information to reflect dangling patch neighbour + // Tag with -(value+2) + forAllConstIters(edgeToCoupledIndex, iter) + { + const edge& e = iter.key(); + const label coupledEdgei = iter.val(); + + patchPairInfo& pairing = cppEdgeData[coupledEdgei]; + const label ownerPatchId = pairing.first(); + + // Some sanity checks + if (ownerPatchId < 0) + { + FatalErrorInFunction + << "Error finding dangling edge at " + << cpp.points()[e.first()] << ' ' + << cpp.points()[e.second()] << nl + << abort(FatalError); + } + else if (ownerPatchId >= nNonProcessor) + { + FatalErrorInFunction + << "Cannot handle edge on processor-processor connection at " + << cpp.points()[e.first()] << ' ' + << cpp.points()[e.second()] << nl + << abort(FatalError); + } + + 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)); + } + + + // Synchronize edge information + + const combineDanglingEdge edgeCombineOp(nNonProcessor); + + globalMeshData::syncData ( cppEdgeData, globalData.globalEdgeSlaves(), - globalData.globalEdgeTransformedSlaves(), + globalData.globalEdgeTransformedSlaves(), // probably not used map, - syncEdgePatchPairs(pbm.nNonProcessor()) + edgeCombineOp ); - // Back from cpp-edge to patch-edge data - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - patchPairInfoUIndList(bndEdgePatchPairs, patchEdgeLabels) = - patchPairInfoUIndList(cppEdgeData, coupledEdgeLabels); + // Combine back from pushed cpp-edge data + + forAll(patchEdgeLabels, i) + { + patchPairInfo& pairing = bndEdgePatchPairs[patchEdgeLabels[i]]; + const patchPairInfo& other = cppEdgeData[coupledEdgeLabels[i]]; + + edgeCombineOp(pairing, other); + + // Resolve special tagging + combineDanglingEdge::correct(pairing); + } return bndEdgePatchPairs; } @@ -499,13 +678,6 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList { const polyBoundaryMesh& pbm = mesh().boundaryMesh(); - const labelListList& edgeFaces = mesh().edgeFaces(); - const labelList meshEdges - ( - patch().meshEdges(mesh().edges(), mesh().pointEdges()) - ); - - // Transcribe into patch definitions DynamicList<faPatchData> faPatchDefs(bndDict.size() + 4); for (const entry& dEntry : bndDict) @@ -576,10 +748,7 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList const label nInternalEdges = patch().nInternalEdges(); const label nBoundaryEdges = patch().nBoundaryEdges(); - patchPairInfoList bndEdgePatchPairs - ( - getBoundaryEdgePatchPairs(meshEdges) - ); + patchPairInfoList bndEdgePatchPairs(this->getBoundaryEdgePatchPairs()); labelList bndEdgeFaPatchIDs(nBoundaryEdges, -1); @@ -653,23 +822,10 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList labelList edgeNbrPolyPatch(undefinedEdges.size(), -1); forAll(edgeNbrPolyPatch, edgei) { - const label localEdgei = undefinedEdges[edgei]; - const label meshEdgei = meshEdges[localEdgei]; - label polyPatchID; + const label patchEdgei = undefinedEdges[edgei]; + const label bndEdgei = (patchEdgei - nInternalEdges); - for (const label meshFacei : edgeFaces[meshEdgei]) - { - if - ( - !faceLabels_.found(meshFacei) - && (polyPatchID = pbm.whichPatch(meshFacei)) != -1 - ) - { - // Found the edge's off-patch neighbour face - edgeNbrPolyPatch[edgei] = polyPatchID; - break; - } - } + edgeNbrPolyPatch[edgei] = bndEdgePatchPairs[bndEdgei].second(); } // Categorize as processor/non-processor associations @@ -762,7 +918,7 @@ Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList { if (patchDef.coupled()) { - reorderProcEdges(patchDef, meshEdges, bndEdgePatchPairs); + reorderProcEdges(patchDef, bndEdgePatchPairs); patchDef.neighPolyPatchId_ = -1; // No lookup of neighbour faces } }