From 7872ed89ef587eb2fb5b6e076eb0c0a326f70908 Mon Sep 17 00:00:00 2001 From: Mark Olesen <Mark.Olesen@esi-group.com> Date: Thu, 12 May 2022 10:05:12 +0200 Subject: [PATCH] ENH: faMeshDistributor for handling redistribution of finiteArea (#2436) --- src/finiteArea/Make/files | 3 + .../distributed/faMeshDistributor.C | 200 ++++ .../distributed/faMeshDistributor.H | 241 +++++ .../distributed/faMeshDistributorNew.C | 916 ++++++++++++++++++ .../distributed/faMeshDistributorTemplates.C | 428 ++++++++ 5 files changed, 1788 insertions(+) create mode 100644 src/finiteArea/distributed/faMeshDistributor.C create mode 100644 src/finiteArea/distributed/faMeshDistributor.H create mode 100644 src/finiteArea/distributed/faMeshDistributorNew.C create mode 100644 src/finiteArea/distributed/faMeshDistributorTemplates.C diff --git a/src/finiteArea/Make/files b/src/finiteArea/Make/files index 7db25ed3f3a..3ab2e2063c2 100644 --- a/src/finiteArea/Make/files +++ b/src/finiteArea/Make/files @@ -24,6 +24,9 @@ $(faPatches)/constraint/wedge/wedgeFaPatch.C $(faPatches)/constraint/cyclic/cyclicFaPatch.C $(faPatches)/constraint/symmetry/symmetryFaPatch.C +distributed/faMeshDistributor.C +distributed/faMeshDistributorNew.C + ensight = output/ensight $(ensight)/ensightFaMesh.C diff --git a/src/finiteArea/distributed/faMeshDistributor.C b/src/finiteArea/distributed/faMeshDistributor.C new file mode 100644 index 00000000000..8094d290e5d --- /dev/null +++ b/src/finiteArea/distributed/faMeshDistributor.C @@ -0,0 +1,200 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "faMeshDistributor.H" +#include "BitOps.H" +#include "fileOperation.H" +#include "areaFields.H" +#include "edgeFields.H" +#include "faMeshTools.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +int Foam::faMeshDistributor::verbose_ = 0; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::faMeshDistributor::createPatchMaps() const +{ + const faBoundaryMesh& oldPatches = srcMesh_.boundary(); + const faBoundaryMesh& newPatches = tgtMesh_.boundary(); + + patchEdgeMaps_.clear(); + patchEdgeMaps_.resize(oldPatches.size()); + + // area: edgeMap (volume: faceMap) + const auto& faEdgeMap = distMap_.faceMap(); + + // The logical edge ranges per patch [target mesh] + List<labelRange> ranges = newPatches.patchRanges(); + + forAll(oldPatches, patchi) + { + if (!isA<processorFaPatch>(oldPatches[patchi])) + { + // Map non-processor only + + // Copy full map + patchEdgeMaps_.set + ( + patchi, + new mapDistributeBase(faEdgeMap) + ); + + // Retain patch elements + patchEdgeMaps_[patchi].compactRemoteData + ( + bitSet(ranges[patchi]), + UPstream::msgType(), + true // Also renumber/resize the compact maps + ); + } + } +} + + +void Foam::faMeshDistributor::createInternalEdgeMap() const +{ + // area: edgeMap (volume: faceMap) + const auto& faEdgeMap = distMap_.faceMap(); + + // Copy full map + internalEdgeMap_.reset(new mapDistributeBase(faEdgeMap)); + + // Retain internal edges + internalEdgeMap_->compactRemoteData + ( + bitSet(tgtMesh_.nInternalEdges(), true), + UPstream::msgType(), + true // Also renumber/resize the compact maps + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::faMeshDistributor::faMeshDistributor +( + const faMesh& srcMesh, + const faMesh& tgtMesh, + const mapDistributePolyMesh& distMap, + const bool isWriteProc +) +: + srcMesh_(srcMesh), + tgtMesh_(tgtMesh), + distMap_(distMap), + internalEdgeMap_(), + patchEdgeMaps_(), + isWriteProc_(isWriteProc) +{ + #ifdef FULLDEBUG + { + Pout<< "Create from nFaces:" << srcMesh.faceLabels().size() + << " to:" << tgtMesh.faceLabels().size() << endl; + + vectorField oldFaceCentres(srcMesh_.areaCentres()); + vectorField newFaceCentres(tgtMesh_.areaCentres()); + + // volume: cells, area: faces + distMap_.distributeCellData(oldFaceCentres); + vectorField diff(newFaceCentres - oldFaceCentres); + + Pout<< "diff faces: " << diff << endl; + + vectorField oldEdgeCentres + ( + faMeshTools::flattenEdgeField(srcMesh_.edgeCentres()) + ); + vectorField newEdgeCentres + ( + faMeshTools::flattenEdgeField(tgtMesh_.edgeCentres()) + ); + + Pout<< "distributed edges: " << oldEdgeCentres.size() << " from " + << srcMesh.nEdges() << " to " << tgtMesh.nEdges() << endl; + + // volume: faces, area: edges + distMap_.distributeFaceData(oldEdgeCentres); + + diff = (newEdgeCentres - oldEdgeCentres); + + Pout<< "diff edges: " << diff << endl; + + Info<< "Patch edge maps" << endl; + forAll(patchEdgeMaps_, patchi) + { + if (patchEdgeMaps_.set(patchi)) + { + Pout<< "patch " << patchi << " : " + << patchEdgeMaps_[patchi].info() << endl; + } + } + + Info<< nl << "Detailed patch maps" << endl; + + forAll(patchEdgeMaps_, patchi) + { + if (patchEdgeMaps_.set(patchi)) + { + Info<< "patch " << patchi << " : " + << patchEdgeMaps_[patchi] << endl; + } + } + } + #endif +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::label Foam::faMeshDistributor::distributeAllFields +( + const IOobjectList& objects, + const wordRes& selected +) const +{ + label nTotal = 0; + + nTotal += distributeAreaFields<scalar>(objects, selected); + nTotal += distributeAreaFields<vector>(objects, selected); + nTotal += distributeAreaFields<symmTensor>(objects, selected); + nTotal += distributeAreaFields<sphericalTensor>(objects, selected); + nTotal += distributeAreaFields<tensor>(objects, selected); + + nTotal += distributeEdgeFields<scalar>(objects, selected); + nTotal += distributeEdgeFields<vector>(objects, selected); + nTotal += distributeEdgeFields<symmTensor>(objects, selected); + nTotal += distributeEdgeFields<sphericalTensor>(objects, selected); + nTotal += distributeEdgeFields<tensor>(objects, selected); + + return nTotal; +} + + +// ************************************************************************* // diff --git a/src/finiteArea/distributed/faMeshDistributor.H b/src/finiteArea/distributed/faMeshDistributor.H new file mode 100644 index 00000000000..47afc8c7c2d --- /dev/null +++ b/src/finiteArea/distributed/faMeshDistributor.H @@ -0,0 +1,241 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::faMeshDistributor + +Description + Holds a reference to the original mesh (the baseMesh) + and optionally to a subset of that mesh (the subMesh) + with mapping lists for points, faces, and cells. + +SourceFiles + faMeshDistributor.C + faMeshDistributorNew.C + faMeshDistributorTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_faMeshDistributor_H +#define Foam_faMeshDistributor_H + +#include "faMesh.H" +#include "mapDistributePolyMesh.H" +#include "areaFieldsFwd.H" +#include "edgeFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class IOobjectList; + +/*---------------------------------------------------------------------------*\ + Class faMeshDistributor Declaration +\*---------------------------------------------------------------------------*/ + +class faMeshDistributor +{ + // Private Data + + //- The source mesh reference + const faMesh& srcMesh_; + + //- The destination mesh reference + const faMesh& tgtMesh_; + + //- Distribution map reference (faMesh) + const mapDistributePolyMesh& distMap_; + + //- Internal edge mapper + mutable std::unique_ptr<mapDistributeBase> internalEdgeMap_; + + //- Patch edge mappers + mutable PtrList<mapDistributeBase> patchEdgeMaps_; + + //- Do I need to write (eg, master only for reconstruct) + bool isWriteProc_; + + + // Private Member Functions + + //- Construct internal edge mapping + void createInternalEdgeMap() const; + + //- Construct per-patch edge mapping + void createPatchMaps() const; + + static mapDistributePolyMesh createReconstructMap + ( + const faMesh& mesh, + const autoPtr<faMesh>& baseMeshPtr, + const labelUList& faceProcAddr, + const labelUList& edgeProcAddr, + const labelUList& pointProcAddr, + const labelUList& boundaryProcAddr + ); + + //- No copy construct + faMeshDistributor(const faMeshDistributor&) = delete; + + //- No copy assignment + void operator=(const faMeshDistributor&) = delete; + + +public: + + //- Output verbosity when writing + static int verbose_; + + + // Constructors + + //- Construct from components + faMeshDistributor + ( + const faMesh& srcMesh, + const faMesh& tgtMesh, + const mapDistributePolyMesh& faDistMap, + const bool isWriteProc = false + ); + + + // Static Methods + + //- Distribute mesh according to the given (volume) mesh distribution. + // Uses 'tgtPolyMesh' for the new mesh + static mapDistributePolyMesh distribute + ( + const faMesh& oldMesh, + const mapDistributePolyMesh& distMap, //! From polyMesh + const polyMesh& tgtPolyMesh, + autoPtr<faMesh>& newMeshPtr + ); + + //- Distribute mesh according to the given (volume) mesh distribution. + // Re-uses polyMesh from oldMesh for the new mesh + static mapDistributePolyMesh distribute + ( + const faMesh& oldMesh, + const mapDistributePolyMesh& distMap, //! From polyMesh + autoPtr<faMesh>& newMeshPtr + ); + + + // Member Functions + + //- Get status of write enabled (on this proc) + bool isWriteProc() const noexcept + { + return isWriteProc_; + } + + //- Change status of write enabled (on this proc) + bool isWriteProc(const bool on) noexcept + { + bool old(isWriteProc_); + isWriteProc_ = on; + return old; + } + + + // Field Mapping + + //- Read, distribute and write all/selected point field types + //- (scalar, vector, ... types) + label distributeAllFields + ( + const IOobjectList& objects, + const wordRes& selectedFields = wordRes() + ) const; + + //- Distribute area field + template<class Type> + tmp<GeometricField<Type, faPatchField, areaMesh>> + distributeField + ( + const GeometricField<Type, faPatchField, areaMesh>& fld + ) const; + + //- Distribute edge field + template<class Type> + tmp<GeometricField<Type, faePatchField, edgeMesh>> + distributeField + ( + const GeometricField<Type, faePatchField, edgeMesh>& fld + ) const; + + //- Read and distribute area field + template<class Type> + tmp<GeometricField<Type, faPatchField, areaMesh>> + distributeAreaField + ( + const IOobject& fieldObject + ) const; + + //- Read and distribute edge field + template<class Type> + tmp<GeometricField<Type, faePatchField, edgeMesh>> + distributeEdgeField + ( + const IOobject& fieldObject + ) const; + + //- Read, distribute and write all/selected area fields + template<class Type> + label distributeAreaFields + ( + const IOobjectList& objects, + const wordRes& selectedFields = wordRes() + ) const; + + //- Read, distribute and write all/selected area fields + template<class Type> + label distributeEdgeFields + ( + const IOobjectList& objects, + const wordRes& selectedFields = wordRes() + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +#ifdef NoRepository +# include "faMeshDistributorTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/distributed/faMeshDistributorNew.C b/src/finiteArea/distributed/faMeshDistributorNew.C new file mode 100644 index 00000000000..00e433a0759 --- /dev/null +++ b/src/finiteArea/distributed/faMeshDistributorNew.C @@ -0,0 +1,916 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "faMeshDistributor.H" +#include "globalIndex.H" +#include "BitOps.H" +#include "ListOps.H" +#include "mapDistributePolyMesh.H" +#include "processorFaPatch.H" +#include "labelPairHashes.H" + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::mapDistributePolyMesh +Foam::faMeshDistributor::distribute +( + const faMesh& oldMesh, + const mapDistributePolyMesh& distMap, + const polyMesh& tgtPolyMesh, + autoPtr<faMesh>& newMeshPtr +) +{ + newMeshPtr.reset(nullptr); + + const uindirectPrimitivePatch& oldAreaPatch = oldMesh.patch(); + + // Original (patch) sizes before any changes + const label nOldPoints = oldAreaPatch.nPoints(); + const label nOldEdges = oldAreaPatch.nEdges(); + const label nOldFaces = oldAreaPatch.nFaces(); + + + // ------------------------ + // Step 1: The face mapping + // ------------------------ + // Relatively straightforward. + // - a subset of face selections on the polyMesh faceMap + + mapDistributeBase faFaceMap(distMap.faceMap()); + + { + // Retain faces needed for the faMesh - preserve original order + + // Note can use compactLocalData without problem since + // finiteArea is only defined on real boundary faces so there + // is no danger of sending an internal or processor face. + + labelList oldToNewSub; + labelList oldToNewConstruct; + + faFaceMap.compactLocalData + ( + oldMesh.faceLabels(), + oldToNewSub, + oldToNewConstruct, + distMap.nOldFaces(), + UPstream::msgType() + ); + + + // The receiving side: + // Mapped faces (>= 0) are the polyMesh face labels that define + // the faMesh. Since the compact order is implicitly sorted in + // ascending order, it tends to form contiguous ranges on the + // polyPatches, which serves our purpose nicely. + + labelList newFaceLabels + ( + ListOps::findIndices(oldToNewConstruct, labelRange::ge0()) + ); + + // Set up to read-if-present + IOobject io(tgtPolyMesh); + io.readOpt(IOobject::READ_IF_PRESENT); + + newMeshPtr.reset + ( + new faMesh(tgtPolyMesh, std::move(newFaceLabels), io) + ); + } + + // The face map is now complete. + + // The new faMesh and the corresponding primitive patch + auto& newMesh = newMeshPtr(); + const uindirectPrimitivePatch& newAreaPatch = newMesh.patch(); + + + // ------------------------ + // Step 2: The edge mapping + // ------------------------ + // Use globally unique addressing to identify both sides of the edges. + + // A local bookkeeping struct of (face0 face1 edge0 edge1) + // selected for a stable and unique sort + + struct faceEdgeTuple : public FixedList<label, 4> + { + // Inherit constructors + using FixedList<label, 4>::FixedList; + + // Default construct as 'invalid' + faceEdgeTuple() : FixedList<label, 4>(-1) {} + + // Is empty if face indices are negative + bool empty() const + { + return (face0() < 0 && face1() < 0); + } + + // Global face numbers are the first sort index + label face0() const { return (*this)[0]; } + label face1() const { return (*this)[1]; } + + // Additional sorting based on edges + label edge0() const { return (*this)[2]; } + label edge1() const { return (*this)[3]; } + + label getFace(int i) const { return (*this)[(i ? 1 : 0)]; } + label getEdge(int i) const { return (*this)[(i ? 3 : 2)]; } + + labelPair getFaces() const { return labelPair(face0(), face1()); } + labelPair getPair0() const { return labelPair(face0(), edge0()); } + // labelPair getPair1() const { return labelPair(face1(), edge1()); } + + // Canonically sorted order + void canonical() + { + if (face1() >= 0 && face0() >= face1()) + { + std::swap((*this)[0], (*this)[1]); + std::swap((*this)[2], (*this)[3]); + } + } + + // Slot a face-edge pair into a free location in the tuple + void add(const labelPair& faceEdge) + { + if ((*this)[0] < 0) // face0 + { + (*this)[0] = faceEdge.first(); // face + (*this)[2] = faceEdge.second(); // edge + } + else if ((*this)[1] < 0) // face1 + { + (*this)[1] = faceEdge.first(); // face + (*this)[3] = faceEdge.second(); // edge + } + } + + // Combine operation + void combine(const faceEdgeTuple& y) + { + auto& x = *this; + + if (y.empty() || x == y) + { + // Nothing to do + } + else if (x.empty()) + { + x = y; + } + else // Both non-empty, but non-identical (should not occur) + { + FatalErrorInFunction + << "Unexpected edge matching: " + << x << " vs. " << y << endl + << exit(FatalError); + } + } + + // Combine operation + void operator()(faceEdgeTuple& x, const faceEdgeTuple& y) const + { + x.combine(y); + } + }; + + + // ------------------------ + // (Edge mapping) + // ------------------------ + // 1. + // Create a list of destination edges for each face, + // appended by a unique face identifier. + // Using global face numbering from the *target* mesh. + + const globalIndex uniqFaceIndex(newAreaPatch.nFaces()); + + labelListList dstFaceEdgeIds(newAreaPatch.nFaces()); + + forAll(newAreaPatch.faceEdges(), facei) + { + const labelList& fEdges = newAreaPatch.faceEdges()[facei]; + labelList& dstEdges = dstFaceEdgeIds[facei]; + + dstEdges.resize(fEdges.size() + 1); + + // Leading entries are the destination (patch) edge indices + labelSubList(dstEdges, fEdges.size()) = fEdges; + + // Last entry is the face id + dstEdges.last() = uniqFaceIndex.toGlobal(facei); + } + + // Send back to the original mesh locations + faFaceMap.reverseDistribute(nOldFaces, dstFaceEdgeIds); + + + // 2. + // Walk all original faces and their edges to generate a edge lookup + // table with the destination face/edge information. + // Eg ((globFacei, localEdgei), (globFacei, localEdgei)) + + // NB: currently no provision for face flips, which would potentially + // reverse the local edge order. + + List<faceEdgeTuple> halfEdgeLookup(nOldEdges, faceEdgeTuple()); + + forAll(oldAreaPatch.faceEdges(), facei) + { + const labelList& fEdges = oldAreaPatch.faceEdges()[facei]; + + // The corresponding destination edges (+ uniqFacei). + const labelList& dstFaceEdges = dstFaceEdgeIds[facei]; + + // Last entry has the unique destination face id + const label uniqFacei = dstFaceEdges.last(); + + forAll(fEdges, faceEdgei) + { + const label srcEdgei = fEdges[faceEdgei]; + const label dstEdgei = dstFaceEdges[faceEdgei]; + + halfEdgeLookup[srcEdgei].add(labelPair(uniqFacei, dstEdgei)); + } + } + + // 3. + // Add patch indices - encoded as '-(patchId + 2)' for + // non-processor boundaries (will be needed later). + + // Also record which procs get patches from here [proc] -> [patches] + // Use for building patchMap + Map<labelHashSet> patchMapInfo; + + label nNonProcessor(0); + { + const faBoundaryMesh& oldBndMesh = oldMesh.boundary(); + + forAll(oldBndMesh, patchi) + { + const faPatch& fap = oldBndMesh[patchi]; + if (isA<processorFaPatch>(fap)) + { + break; + } + ++nNonProcessor; + + // Encoded as -(patchId + 2) + const labelPair encodedPatchId(-(patchi+2), -(patchi+2)); + + for (const label srcEdgei : fap.edgeLabels()) + { + faceEdgeTuple& facePairings = halfEdgeLookup[srcEdgei]; + + facePairings.add(encodedPatchId); + + label dstFacei = facePairings.face0(); + + label proci = uniqFaceIndex.whichProcID(dstFacei); + + patchMapInfo(proci).insert(patchi); + } + } + + Pstream::broadcast(nNonProcessor); + } + + // Processor-processor connections + if (Pstream::parRun()) + { + const label startOfRequests = Pstream::nRequests(); + + const faBoundaryMesh& oldBndMesh = oldMesh.boundary(); + + // Setup sends + for (const faPatch& fap : oldBndMesh) + { + const auto* fapp = isA<processorFaPatch>(fap); + + if (fapp) + { + // Send (dstFacei dstEdgei) per patch edge. + // Since we are walking a boundary edge, the first location + // from the previously established lookup provides + // our information. + + List<labelPair> edgeTuples(fap.edgeLabels().size()); + + label edgeTuplei = 0; + + for (const label edgei : fap.edgeLabels()) + { + edgeTuples[edgeTuplei] = halfEdgeLookup[edgei].getPair0(); + ++edgeTuplei; + } + + fapp->send<labelPair> + ( + Pstream::commsTypes::nonBlocking, + edgeTuples + ); + } + } + + // Wait for all to finish + Pstream::waitRequests(startOfRequests); + + // Receive values + for (const faPatch& fap : oldBndMesh) + { + const auto* fapp = isA<processorFaPatch>(fap); + + if (fapp) + { + // Receive (dstFacei dstEdgei) per patch edge. + // - slot into the second location of the lookup + + List<labelPair> edgeTuples(fap.edgeLabels().size()); + + fapp->receive<labelPair> + ( + Pstream::commsTypes::nonBlocking, + edgeTuples + ); + + label edgeTuplei = 0; + + for (const label srcEdgei : fap.edgeLabels()) + { + halfEdgeLookup[srcEdgei].add(edgeTuples[edgeTuplei]); + ++edgeTuplei; + } + } + } + } + + // Globally consistent order + for (auto& tuple : halfEdgeLookup) + { + tuple.canonical(); + } + + + // Now ready to actually construct the edgeMap + + mapDistributeBase faEdgeMap + ( + newAreaPatch.nEdges(), // constructSize + labelListList(), // subMap + labelListList(), // constructMap + false, // subHasFlip + false, // constructHasFlip + faFaceMap.comm() + ); + + { + // Pass 1. + // Count the number of edges to be sent to each proc + + labelList nProcEdges(Pstream::nProcs(faFaceMap.comm()), Zero); + + forAll(halfEdgeLookup, srcEdgei) + { + const faceEdgeTuple& facePairings = halfEdgeLookup[srcEdgei]; + + labelPair dstProcs(-1, -1); + + for (int sidei = 0; sidei < 2; ++sidei) + { + const label dstFacei = facePairings.getFace(sidei); + //const label dstEdgei = facePairings.getEdge(sidei); + + if (dstFacei < 0) + { + continue; + } + + label proci = uniqFaceIndex.whichProcID(dstFacei); + dstProcs[sidei] = proci; + + if (proci != -1 && dstProcs[0] != dstProcs[1]) + { + ++nProcEdges[proci]; + } + } + } + + auto& edgeSubMap = faEdgeMap.subMap(); + auto& edgeCnstrMap = faEdgeMap.constructMap(); + + edgeSubMap.resize(nProcEdges.size()); + edgeCnstrMap.resize(nProcEdges.size()); + + labelListList remoteEdges(nProcEdges.size()); + + forAll(nProcEdges, proci) + { + edgeSubMap[proci].resize(nProcEdges[proci], -1); + remoteEdges[proci].resize(nProcEdges[proci], -1); + } + + + // Pass 2. + // Fill in the maps + + nProcEdges = Zero; // Reset counter + + forAll(halfEdgeLookup, srcEdgei) + { + const faceEdgeTuple& facePairings = halfEdgeLookup[srcEdgei]; + + labelPair dstProcs(-1, -1); + + for (int sidei = 0; sidei < 2; ++sidei) + { + const label dstFacei = facePairings.getFace(sidei); + const label dstEdgei = facePairings.getEdge(sidei); + + if (dstFacei < 0) + { + continue; + } + + label proci = uniqFaceIndex.whichProcID(dstFacei); + dstProcs[sidei] = proci; + + if (proci != -1 && dstProcs[0] != dstProcs[1]) + { + edgeSubMap[proci][nProcEdges[proci]] = srcEdgei; + remoteEdges[proci][nProcEdges[proci]] = dstEdgei; + ++nProcEdges[proci]; + } + } + } + + // The remoteEdges are what we know locally about what will be + // received, but not what is actually received. + // So need an all-to-all exchange + + Pstream::exchange<labelList, label> + ( + remoteEdges, + edgeCnstrMap, + UPstream::msgType(), + faEdgeMap.comm() + ); + } + + // The edge map is now complete [in PrimitivePatch edge order] + + + if (oldMesh.hasInternalEdgeLabels()) + { + // If there are gaps in the edge numbering or the + // internal edge labels are out of sequence would + // have to use compactLocalData etc before sending + // But just issue an error for now + + FatalErrorInFunction + << "Originating faMesh has gaps in the edge addressing" + << " this is currently unsupported" + << abort(FatalError); + } + + + // ------------------------ + // Patch edge labels + // ------------------------ + + faPatchList newFaPatches; + + // Distribute the edge lookups. + // Needs full version for the combine operation + + mapDistributeBase::distribute<faceEdgeTuple, faceEdgeTuple, identityOp> + ( + Pstream::commsTypes::nonBlocking, + List<labelPair>::null(), + faEdgeMap.constructSize(), + + faEdgeMap.subMap(), + faEdgeMap.subHasFlip(), + + faEdgeMap.constructMap(), + faEdgeMap.constructHasFlip(), + + halfEdgeLookup, + faceEdgeTuple(), // nullValue + faceEdgeTuple(), // CombineOp + identityOp(), // NegateOp + UPstream::msgType(), + faEdgeMap.comm() + ); + + { + // Pass 1. + // Count the number of edges for each patch type + + Map<label> nNonProcEdges(2*nNonProcessor); + LabelPairMap<label> nProcEdges(2*nNonProcessor); + + forAll(halfEdgeLookup, edgei) + { + const auto& tuple = halfEdgeLookup[edgei]; + + labelPair target(tuple.getFaces()); + + if (target[1] < -1) + { + // Neighbour face was patchId encoded value + label patchi = mag(target[1])-2; + ++nNonProcEdges(patchi); + } + else if (target[0] >= 0 && target[1] >= 0) + { + // From global face to proc id + target[0] = uniqFaceIndex.whichProcID(target[0]); + target[1] = uniqFaceIndex.whichProcID(target[1]); + + // A connection between different procs but involving + // myProc? + if + ( + target[0] != target[1] + && + ( + target[0] == Pstream::myProcNo() + || target[1] == Pstream::myProcNo() + ) + ) + { + ++nProcEdges(target); + } + } + } + + label nPatches = (nNonProcessor + nProcEdges.size()); + + newFaPatches.resize(nPatches); + + labelList nEdgeLabels(nPatches, Zero); + labelListList newEdgeLabels(nPatches); + + LabelPairMap<label> procPairToPatchId(2*nProcEdges.size()); + + // Map processor-pair to index in patches + { + nPatches = nNonProcessor; + + for (const labelPair& twoProcs : nProcEdges.sortedToc()) + { + procPairToPatchId.set(twoProcs, nPatches); + ++nPatches; + } + } + + // Presize edgeLabels arrays + { + nPatches = 0; + + for (label patchi = 0; patchi < nNonProcessor; ++patchi) + { + label nLabels = nNonProcEdges.lookup(nPatches, 0); + + nEdgeLabels[nPatches] = nLabels; + + newEdgeLabels[nPatches].resize(nLabels); + + ++nPatches; + } + + // Processor patches + for (const labelPair& twoProcs : nProcEdges.sortedToc()) + { + label nLabels = nProcEdges.lookup(twoProcs, 0); + + nEdgeLabels[nPatches] = nLabels; + + newEdgeLabels[nPatches].resize(nLabels); + + ++nPatches; + } + } + + nEdgeLabels = Zero; + + // Populate edgeLabels arrays - walk in canonically sorted + // order to ensure that both sides of processor edges + // correspond. + + const labelList order(Foam::sortedOrder(halfEdgeLookup)); + + for (const label edgei : order) + { + const auto& tuple = halfEdgeLookup[edgei]; + + labelPair target(tuple.getFaces()); + + label patchi = -1; + + if (target[1] < -1) + { + // Neighbour face was patchId encoded value + patchi = mag(target[1])-2; + } + else if (target[0] >= 0 && target[1] >= 0) + { + // From global face to proc id + target[0] = uniqFaceIndex.whichProcID(target[0]); + target[1] = uniqFaceIndex.whichProcID(target[1]); + + // A connection between different procs but involving + // myProc? + if + ( + target[0] != target[1] + && + ( + target[0] == Pstream::myProcNo() + || target[1] == Pstream::myProcNo() + ) + ) + { + patchi = procPairToPatchId.lookup(target, -1); + } + } + + if (patchi >= 0) + { + const label idx = nEdgeLabels[patchi]++; + newEdgeLabels[patchi][idx] = edgei; + } + } + + + // Clone all non-processor patches + + nPatches = 0; + for (label patchi = 0; patchi < nNonProcessor; ++patchi) + { + newFaPatches.set + ( + nPatches, + oldMesh.boundary()[patchi].clone + ( + newMesh.boundary(), + newEdgeLabels[nPatches], // edgeLabels + nPatches, + oldMesh.boundary()[patchi].ngbPolyPatchIndex() + ) + ); + ++nPatches; + } + + // Any processor patches + for (const labelPair& twoProcs : nProcEdges.sortedToc()) + { + label nbrProcNo = + ( + twoProcs[1] == Pstream::myProcNo() + ? twoProcs[0] : twoProcs[1] + ); + + newFaPatches.set + ( + nPatches, + new processorFaPatch + ( + newEdgeLabels[nPatches], // edgeLabels + nPatches, + newMesh.boundary(), + -1, // nbrPolyPatchi + Pstream::myProcNo(), + nbrProcNo + ) + ); + ++nPatches; + } + } + + + newMesh.addFaPatches(newFaPatches); + newMesh.init(true); + + + // At this stage we now have a complete mapping overview in + // terms of the PrimitivePatch edge ordering which now must be + // changed into faMesh edge ordering. + + { + labelList oldToNewSub; + labelList oldToNewConstruct; + + // Assume we use all patch edges for the faMesh too. + oldToNewSub.resize(oldAreaPatch.nEdges(), -1); + oldToNewConstruct.resize(newAreaPatch.nEdges(), -1); + + // Remap old edges + label nEdges = 0; + + // Internal edgeLabels + for + ( + label edgei = 0; + edgei < oldAreaPatch.nInternalEdges(); + ++edgei + ) + { + oldToNewSub[edgei] = nEdges++; + } + + // Boundary edgeLabels + for (const faPatch& fap : oldMesh.boundary()) + { + for (const label edgei : fap.edgeLabels()) + { + oldToNewSub[edgei] = nEdges++; + } + } + + // New edges + nEdges = 0; + + // Internal edgeLabels + for + ( + label edgei = 0; + edgei < newAreaPatch.nInternalEdges(); + ++edgei + ) + { + oldToNewConstruct[edgei] = nEdges++; + } + + // Boundary edgeLabels + for (const faPatch& fap : newMesh.boundary()) + { + for (const label edgei : fap.edgeLabels()) + { + oldToNewConstruct[edgei] = nEdges++; + } + } + + mapDistributeBase::renumberMap + ( + faEdgeMap.subMap(), + oldToNewSub, + faEdgeMap.subHasFlip() + ); + + mapDistributeBase::renumberMap + ( + faEdgeMap.constructMap(), + oldToNewConstruct, + faEdgeMap.constructHasFlip() + ); + } + + + // ------------------------ + // Patch mapping + // ------------------------ + + mapDistributeBase faPatchMap + ( + newMesh.boundary().size(), // constructSize + labelListList(), // subMap + labelListList(), // constructMap + false, // subHasFlip + false, // constructHasFlip + faFaceMap.comm() + ); + + // For patch maps, would normally transcribe from patchMapInfo + // gathered earlier. However, existing practice (APR-2022) for + // faMesh decomposition is to map all non-processor patches + + { + // Map all non-processor patches + const label nProcs = UPstream::nProcs(faPatchMap.comm()); + + faPatchMap.subMap().resize(nProcs, identity(nNonProcessor)); + faPatchMap.constructMap().resize(nProcs, identity(nNonProcessor)); + } + + /// { + /// // Transcribe from patchMapInfo gathered earlier. + /// // - transform Map of labelHashSet to labelListList + /// + /// labelListList sendToRemote(Pstream::nProcs(faPatchMap.comm())); + /// + /// forAllConstIters(patchMapInfo, iter) + /// { + /// const label proci = iter.key(); + /// sendToRemote[proci] = iter.val().sortedToc(); + /// } + /// + /// auto& patchSubMap = faPatchMap.subMap(); + /// auto& patchCnstrMap = faPatchMap.constructMap(); + /// + /// patchSubMap = sendToRemote; + /// patchCnstrMap.resize(patchSubMap.size()); + /// + /// // Change sendToRemote into recv-from-remote by using + /// // all-to-all exchange + /// + /// Pstream::exchange<labelList, label> + /// ( + /// sendToRemote, + /// patchCnstrMap, + /// UPstream::msgType(), + /// faPatchMap.comm() + /// ); + /// } + + + // ------------------------ + // Point mapping + // ------------------------ + + mapDistributeBase faPointMap(distMap.pointMap()); + + { + // Retain meshPoints needed for the faMesh - preserve original order + // Need both sides (local/remote) for correct compaction maps + // without dangling points. + + labelList oldToNewSub; + labelList oldToNewConstruct; + + faPointMap.compactData + ( + oldAreaPatch.meshPoints(), + newAreaPatch.meshPoints(), + oldToNewSub, + oldToNewConstruct, + distMap.nOldPoints(), + UPstream::msgType() + ); + } + + + return mapDistributePolyMesh + ( + // Mesh before changes + nOldPoints, + nOldEdges, // area: nOldEdges (volume: nOldFaces) + nOldFaces, // area: nOldFaces (volume: nOldCells) + + labelList(oldMesh.boundary().patchStarts()), + labelList(), // oldPatchNMeshPoints [unused] + + mapDistribute(std::move(faPointMap)), + mapDistribute(std::move(faEdgeMap)), // area: edgeMap (volume: faceMap) + mapDistribute(std::move(faFaceMap)), // area: faceMap (volume: cellMap) + mapDistribute(std::move(faPatchMap)) + ); +} + + +Foam::mapDistributePolyMesh +Foam::faMeshDistributor::distribute +( + const faMesh& oldMesh, + const mapDistributePolyMesh& distMap, + autoPtr<faMesh>& newMeshPtr +) +{ + return faMeshDistributor::distribute + ( + oldMesh, + distMap, + oldMesh.mesh(), // polyMesh + newMeshPtr + ); +} + + +// ************************************************************************* // diff --git a/src/finiteArea/distributed/faMeshDistributorTemplates.C b/src/finiteArea/distributed/faMeshDistributorTemplates.C new file mode 100644 index 00000000000..6f556f3dab5 --- /dev/null +++ b/src/finiteArea/distributed/faMeshDistributorTemplates.C @@ -0,0 +1,428 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "Time.H" +#include "emptyFaPatchField.H" +#include "emptyFaePatchField.H" +#include "IOobjectList.H" +#include "polyMesh.H" +#include "polyPatch.H" +#include "processorFaPatch.H" +#include "mapDistribute.H" +#include "mapDistributePolyMesh.H" +#include "areaFields.H" +#include "edgeFields.H" + +#include "distributedFieldMapper.H" +#include "distributedFaPatchFieldMapper.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>> +Foam::faMeshDistributor::distributeField +( + const GeometricField<Type, faPatchField, areaMesh>& fld +) const +{ + typedef typename + GeometricField<Type, faPatchField, areaMesh>::Patch + PatchFieldType; + + if (tgtMesh_.boundary().size() && patchEdgeMaps_.empty()) + { + createPatchMaps(); + } + + // Create internalField by remote mapping + + const distributedFieldMapper mapper + ( + labelUList::null(), + distMap_.cellMap() // area: faceMap (volume: cellMap) + ); + + DimensionedField<Type, areaMesh> internalField + ( + IOobject + ( + fld.name(), + tgtMesh_.time().timeName(), + fld.local(), + tgtMesh_.thisDb(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + tgtMesh_, + fld.dimensions(), + Field<Type>(fld.internalField(), mapper) + ); + + internalField.oriented() = fld.oriented(); + + + // Create patchFields by remote mapping + + PtrList<PatchFieldType> newPatchFields(tgtMesh_.boundary().size()); + + const auto& bfld = fld.boundaryField(); + + forAll(bfld, patchi) + { + if (patchEdgeMaps_.set(patchi)) + { + // Clone local patch field + + const distributedFaPatchFieldMapper mapper + ( + labelUList::null(), + patchEdgeMaps_[patchi] + ); + + // Map into local copy + newPatchFields.set + ( + patchi, + PatchFieldType::New + ( + bfld[patchi], + tgtMesh_.boundary()[patchi], + DimensionedField<Type, areaMesh>::null(), + mapper + ) + ); + } + } + + // Add empty patchFields on remaining patches (this also handles + // e.g. processorPatchFields or any other constraint type patches) + forAll(newPatchFields, patchi) + { + if (!newPatchFields.set(patchi)) + { + newPatchFields.set + ( + patchi, + PatchFieldType::New + ( + emptyFaPatchField<Type>::typeName, + tgtMesh_.boundary()[patchi], + DimensionedField<Type, areaMesh>::null() + ) + ); + } + } + + + auto tresult = tmp<GeometricField<Type, faPatchField, areaMesh>>::New + ( + std::move(internalField), + newPatchFields + ); + auto& result = tresult.ref(); + + result.boundaryFieldRef().template evaluateCoupled<processorFaPatch>(); + + return tresult; +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>> +Foam::faMeshDistributor::distributeField +( + const GeometricField<Type, faePatchField, edgeMesh>& fld +) const +{ + typedef typename + GeometricField<Type, faePatchField, edgeMesh>::Patch + PatchFieldType; + + if (!internalEdgeMap_) + { + createInternalEdgeMap(); + } + + + // Create internalField by remote mapping + + const distributedFieldMapper mapper + ( + labelUList::null(), + *(internalEdgeMap_) + ); + + DimensionedField<Type, edgeMesh> internalField + ( + IOobject + ( + fld.name(), + tgtMesh_.time().timeName(), + fld.local(), + tgtMesh_.thisDb(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + tgtMesh_, + fld.dimensions(), + Field<Type>(fld.internalField(), mapper) + ); + + internalField.oriented() = fld.oriented(); + + + // Create patchFields by remote mapping + + PtrList<PatchFieldType> newPatchFields(tgtMesh_.boundary().size()); + + const auto& bfld = fld.boundaryField(); + + forAll(bfld, patchi) + { + if (patchEdgeMaps_.set(patchi)) + { + // Clone local patch field + + const distributedFaPatchFieldMapper mapper + ( + labelUList::null(), + patchEdgeMaps_[patchi] + ); + + // Map into local copy + newPatchFields.set + ( + patchi, + PatchFieldType::New + ( + bfld[patchi], + tgtMesh_.boundary()[patchi], + DimensionedField<Type, edgeMesh>::null(), + mapper + ) + ); + } + } + + // Add empty patchFields on remaining patches (this also handles + // e.g. processorPatchFields or any other constraint type patches) + forAll(newPatchFields, patchi) + { + if (!newPatchFields.set(patchi)) + { + newPatchFields.set + ( + patchi, + PatchFieldType::New + ( + emptyFaePatchField<Type>::typeName, + tgtMesh_.boundary()[patchi], + DimensionedField<Type, edgeMesh>::null() + ) + ); + } + } + + + return tmp<GeometricField<Type, faePatchField, edgeMesh>>::New + ( + std::move(internalField), + newPatchFields + ); +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>> +Foam::faMeshDistributor::distributeAreaField +( + const IOobject& fieldObject +) const +{ + // Read field + GeometricField<Type, faPatchField, areaMesh> fld + ( + fieldObject, + srcMesh_ + ); + + // Redistribute + return distributeField(fld); +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>> +Foam::faMeshDistributor::distributeEdgeField +( + const IOobject& fieldObject +) const +{ + // Read field + GeometricField<Type, faePatchField, edgeMesh> fld + ( + fieldObject, + srcMesh_ + ); + + // Redistribute + return distributeField(fld); +} + + +template<class Type> +Foam::label Foam::faMeshDistributor::distributeAreaFields +( + const IOobjectList& objects, + const wordRes& selectedFields +) const +{ + typedef GeometricField<Type, faPatchField, areaMesh> fieldType; + + label nFields = 0; + + for + ( + const IOobject& io : + ( + selectedFields.empty() + ? objects.sorted<fieldType>() + : objects.sorted<fieldType>(selectedFields) + ) + ) + { + if (verbose_) + { + if (!nFields) + { + Info<< " Reconstructing " + << fieldType::typeName << "s\n" << nl; + } + Info<< " " << io.name() << nl; + } + ++nFields; + + tmp<fieldType> tfld(distributeAreaField<Type>(io)); + if (isWriteProc_) + { + tfld().write(); + } + } + + if (nFields && verbose_) Info<< endl; + return nFields; +} + + +template<class Type> +Foam::label Foam::faMeshDistributor::distributeEdgeFields +( + const IOobjectList& objects, + const wordRes& selectedFields +) const +{ + typedef GeometricField<Type, faePatchField, edgeMesh> fieldType; + + label nFields = 0; + + for + ( + const IOobject& io : + ( + selectedFields.empty() + ? objects.sorted<fieldType>() + : objects.sorted<fieldType>(selectedFields) + ) + ) + { + if (verbose_) + { + if (!nFields) + { + Info<< " Reconstructing " + << fieldType::typeName << "s\n" << nl; + } + Info<< " " << io.name() << nl; + } + ++nFields; + + tmp<fieldType> tfld(distributeEdgeField<Type>(io)); + if (isWriteProc_) + { + tfld().write(); + } + } + + if (nFields && verbose_) Info<< endl; + return nFields; +} + + + +#if 0 +template<class Type> +void Foam::faMeshDistributor::redistributeAndWrite +( + PtrList<GeometricField<Type, faPatchField, areaMesh>>& flds +) const +{ + for (auto& fld : flds) + { + Pout<< "process: " << fld.name() << endl; + + tmp<GeometricField<Type, faPatchField, areaMesh>> tfld = + this->distributeField(fld); + + if (isWriteProc_) + { + tfld().write(); + } + } +} + + +template<class Type> +void Foam::faMeshDistributor::redistributeAndWrite +( + PtrList<GeometricField<Type, faePatchField, edgeMesh>>& flds +) const +{ + for (auto& fld : flds) + { + tmp<GeometricField<Type, faePatchField, edgeMesh>> tfld = + this->distributeField(fld); + + if (isWriteProc_) + { + tfld().write(); + } + } +} +#endif + + +// ************************************************************************* // -- GitLab