ENH: faMeshDistributor for handling redistribution of finiteArea (#2436)

This commit is contained in:
Mark Olesen
2022-05-12 10:05:12 +02:00
committed by Andrew Heather
parent 15535b8970
commit 4c5c15664e
5 changed files with 1788 additions and 0 deletions

View File

@ -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

View File

@ -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;
}
// ************************************************************************* //

View File

@ -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
// ************************************************************************* //

View File

@ -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
);
}
// ************************************************************************* //

View File

@ -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
// ************************************************************************* //