zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: add faMeshTools for new/load mesh, proc addressing etc

Browse Source 
ENH: simple faMeshSubset (zero-sized meshes only)

ENH: additional access methods for faMesh, primitive geometry mode

- wrapped walking of boundary edgeLabels as list of list
  (similar to edgeFaces).

- primitive finiteArea geometry mode with reduced communication:
  primarily interesting for decomposition/redistribution (#2436)

ENH: extra vtk debug outputs for checkFaMesh

- report per-processor sizes in the mesh summary
This commit is contained in:
Mark Olesen
2022-05-12 09:44:35 +02:00
parent 68b692fdfc
commit 06b353f8cd
19 changed files with 2600 additions and 379 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2021 OpenCFD Ltd.
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -46,6 +46,8 @@ Original Authors
#include "edgeFields.H"
#include "processorFaPatch.H"
#include "foamVtkIndPatchWriter.H"
#include "foamVtkLineWriter.H"
#include "faMeshTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -66,11 +68,27 @@ int main(int argc, char *argv[])
"Write mesh as a vtp (vtk) file for display or debugging"
);
argList::addOption
(
"geometryOrder",
"N",
"Test different geometry order - experimental!!",
true // Advanced option
);
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedPolyMesh.H"
int geometryOrder(1);
if (args.readIfPresent("geometryOrder", geometryOrder))
{
Info<< "Setting faMesh::geometryOrder = " << geometryOrder << nl
<< "(experimental)" << nl << endl;
faMesh::geometryOrder(geometryOrder);
}
// Create
faMesh aMesh(mesh);

94
applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -16,7 +16,7 @@ Description
\*---------------------------------------------------------------------------*/
{
// finiteArea
// finiteArea - faces
vtk::uindirectPatchWriter writer
(
aMesh.patch(),
@ -29,9 +29,23 @@ Description
writer.writeGeometry();
globalIndex procAddr(aMesh.nFaces());
labelList cellIDs;
if (Pstream::master())
{
cellIDs.resize(procAddr.totalSize());
for (const labelRange& range : procAddr.ranges())
{
auto slice = cellIDs.slice(range);
slice = identity(range);
}
}
// CellData
writer.beginCellData(3);
writer.beginCellData(4);
writer.writeProcIDs();
writer.write("cellID", cellIDs);
writer.write("area", aMesh.S().field());
writer.write("normal", aMesh.faceAreaNormals());
@ -43,5 +57,79 @@ Description
<< "Wrote faMesh in vtk format: " << writer.output().name() << nl;
}
{
// finiteArea - edges
vtk::lineWriter writer
(
aMesh.points(),
aMesh.edges(),
// vtk::formatType::INLINE_ASCII,
fileName
(
aMesh.mesh().time().globalPath() / "finiteArea-edges"
)
);
writer.writeGeometry();
// CellData
writer.beginCellData(4);
writer.writeProcIDs();
{
// Use primitive patch order
Field<scalar> fld
(
faMeshTools::flattenEdgeField(aMesh.magLe(), true)
);
writer.write("magLe", fld);
}
// PointData
writer.beginPointData(1);
writer.write("normal", aMesh.pointAreaNormals());
Info<< nl
<< "Wrote faMesh in vtk format: " << writer.output().name() << nl;
}
{
// finiteArea - edgeCentres
// (no other convenient way to display vectors on the edges)
vtk::lineWriter writer
(
aMesh.edgeCentres(),
edgeList::null(),
// vtk::formatType::INLINE_ASCII,
fileName
(
aMesh.mesh().time().globalPath() / "finiteArea-edgesCentres"
)
);
writer.writeGeometry();
// PointData
writer.beginPointData(4);
{
// Use primitive patch order
Field<vector> fld
(
faMeshTools::flattenEdgeField(aMesh.Le(), true)
);
writer.write("Le", fld);
}
{
// Use primitive patch order
Field<vector> fld
(
faMeshTools::flattenEdgeField(aMesh.edgeAreaNormals(), true)
);
writer.write("normal", fld);
}
Info<< nl
<< "Wrote faMesh in vtk format: " << writer.output().name() << nl;
}
// ************************************************************************* //

99
applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -20,24 +20,7 @@ Description
const label nNonProcessor = patches.nNonProcessor();
const label nPatches = patches.size();
Info<< "----------------" << nl
<< "Mesh Information" << nl
<< "----------------" << nl
<< " " << "boundingBox: " << boundBox(aMesh.points()) << nl;
Info<< " Number of points: "
<< returnReduce(aMesh.nPoints(), sumOp<label>()) << nl
<< " Number of faces: "
<< returnReduce(aMesh.nFaces(), sumOp<label>()) << nl;
Info<< " Number of edges: "
<< returnReduce(aMesh.nEdges(), sumOp<label>()) << nl
<< " Number of internal edges: "
<< returnReduce(aMesh.nInternalEdges(), sumOp<label>()) << nl;
label nProcEdges = 0;
label nLocalProcEdges = 0;
if (Pstream::parRun())
{
for (const faPatch& fap : patches)
@ -46,22 +29,84 @@ Description
if (cpp)
{
nProcEdges += fap.nEdges();
nLocalProcEdges += fap.nEdges();
}
}
}
const label nBoundEdges = aMesh.nBoundaryEdges() - nProcEdges;
const labelList nFaces
(
UPstream::listGatherValues<label>(aMesh.nFaces())
);
const labelList nPoints
(
UPstream::listGatherValues<label>(aMesh.nPoints())
);
const labelList nEdges
(
UPstream::listGatherValues<label>(aMesh.nEdges())
);
const labelList nIntEdges
(
UPstream::listGatherValues<label>(aMesh.nInternalEdges())
);
Info<< " Number of boundary edges: "
<< returnReduce(nBoundEdges - nProcEdges, sumOp<label>()) << nl;
// The "real" (non-processor) boundary edges
const labelList nBndEdges
(
UPstream::listGatherValues<label>
(
aMesh.nBoundaryEdges() - nLocalProcEdges
)
);
const labelList nProcEdges
(
UPstream::listGatherValues<label>(nLocalProcEdges)
);
if (Pstream::parRun())
// Format output as
// Number of faces: ...
// per-proc: (...)
const auto reporter =
[&](const char* tag, const labelList& list)
{
Info<< " Number of " << tag << ": " << sum(list) << nl;
if (Pstream::parRun())
{
int padding = static_cast<int>
(
// strlen(" Number of ") - strlen("per-proc")
(12 - 8)
+ strlen(tag)
);
do { Info<< ' '; } while (--padding > 0);
Info<< "per-proc: " << flatOutput(list) << nl;
}
};
Info<< "----------------" << nl
<< "Mesh Information" << nl
<< "----------------" << nl
<< " " << "boundingBox: " << boundBox(aMesh.points()) << nl;
if (Pstream::master())
{
Info<< " Number of processor edges: "
<< returnReduce(nProcEdges, sumOp<label>()) << nl;
}
reporter("faces", nFaces);
reporter("points", nPoints);
reporter("edges", nEdges);
reporter("internal edges", nIntEdges);
reporter("boundary edges", nBndEdges);
if (Pstream::parRun())
{
reporter("processor edges", nProcEdges);
}
}
Info<< "----------------" << nl
<< "Patches" << nl

4
src/finiteArea/Make/files
View File

@ -9,6 +9,10 @@ faMesh/faMeshBoundaryHalo.C
faMesh/faBoundaryMesh/faBoundaryMesh.C
faMesh/faBoundaryMesh/faBoundaryMeshEntries.C
faMesh/faMeshSubset/faMeshSubset.C
faMesh/faMeshTools/faMeshTools.C
faMesh/faMeshTools/faMeshToolsProcAddr.C
faPatches = faMesh/faPatches
$(faPatches)/faPatch/faPatch.C
$(faPatches)/faPatch/faPatchData.C

16
src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
View File

@ -235,6 +235,22 @@ void Foam::faBoundaryMesh::calcGeometry()
}
Foam::UPtrList<const Foam::labelUList>
Foam::faBoundaryMesh::edgeLabels() const
{
const faPatchList& patches = *this;
UPtrList<const labelUList> list(patches.size());
forAll(list, patchi)
{
list.set(patchi, &patches[patchi].edgeLabels());
}
return list;
}
Foam::UPtrList<const Foam::labelUList>
Foam::faBoundaryMesh::edgeFaces() const
{

3
src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
View File

@ -132,6 +132,9 @@ public:
return mesh_;
}
//- Return a list of edgeLabels for each patch
UPtrList<const labelUList> edgeLabels() const;
//- Return a list of edgeFaces for each patch
UPtrList<const labelUList> edgeFaces() const;

8
src/finiteArea/faMesh/faMesh.C
View File

@ -337,10 +337,12 @@ Foam::faMesh::faMesh
*this
),
comm_(Pstream::worldComm),
curTimeIndex_(time().timeIndex()),
patchPtr_(nullptr),
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
curTimeIndex_(time().timeIndex()),
SPtr_(nullptr),
S0Ptr_(nullptr),
S00Ptr_(nullptr),
@ -439,10 +441,12 @@ Foam::faMesh::faMesh
label(0)
),
comm_(Pstream::worldComm),
curTimeIndex_(time().timeIndex()),
patchPtr_(nullptr),
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
curTimeIndex_(time().timeIndex()),
SPtr_(nullptr),
S0Ptr_(nullptr),
S00Ptr_(nullptr),

164
src/finiteArea/faMesh/faMesh.H
View File

@ -209,37 +209,43 @@ class faMesh
faBoundaryMesh boundary_;
// Primitive mesh data
// Primitive mesh data
//- Edges, addressing into local point list
edgeList edges_;
//- Edges, addressing into local point list
edgeList edges_;
//- Edge owner
labelList edgeOwner_;
//- Edge owner
labelList edgeOwner_;
//- Edge neighbour
labelList edgeNeighbour_;
//- Edge neighbour
labelList edgeNeighbour_;
// Primitive size data
// Primitive size data
//- Total number of points
mutable label nPoints_;
//- Total number of points
mutable label nPoints_;
//- Total number of edges
mutable label nEdges_;
//- Total number of edges
mutable label nEdges_;
//- Number of internal edges
mutable label nInternalEdges_;
//- Number of internal edges
mutable label nInternalEdges_;
//- Number of faces
mutable label nFaces_;
//- Number of faces
mutable label nFaces_;
// Communication support
// Communication support, updating
//- Communicator used for parallel communication
label comm_;
//- Communicator used for parallel communication
label comm_;
//- Current time index for motion
// Note. The whole mechanism will be replaced once the
// dimensionedField is created and the dimensionedField
// will take care of the old-time levels.
mutable label curTimeIndex_;
// Demand-driven data
@ -253,11 +259,8 @@ class faMesh
//- Ldu addressing data
mutable faMeshLduAddressing* lduPtr_;
//- Current time index for motion
// Note. The whole mechanism will be replaced once the
// dimensionedField is created and the dimensionedField
// will take care of the old-time levels.
mutable label curTimeIndex_;
// Geometric Data
//- Face areas
mutable DimensionedField<scalar, areaMesh>* SPtr_;
@ -545,6 +548,22 @@ public:
// Static Functions
//- Return the current geometry treatment (0: primitive, 1: standard)
// A zero level is with restricted neighbour information
static int geometryOrder() noexcept
{
return geometryOrder_;
}
//- Set the preferred geometry treatment
// \return the previous value
static int geometryOrder(const int order) noexcept
{
int old(geometryOrder_);
geometryOrder_ = order;
return old;
}
//- Read construction from polyMesh if all files are available
static autoPtr<faMesh> TryNew(const polyMesh& pMesh);
@ -571,73 +590,79 @@ public:
bool init(const bool doInit);
// Database
// Database
//- Return access to polyMesh
inline const polyMesh& mesh() const;
//- Return access to polyMesh
inline const polyMesh& mesh() const;
//- Interface to referenced polyMesh (similar to GeoMesh)
const polyMesh& operator()() const { return mesh(); }
//- Interface to referenced polyMesh (similar to GeoMesh)
const polyMesh& operator()() const { return mesh(); }
//- Return the local mesh directory (dbDir()/meshSubDir)
fileName meshDir() const;
//- Return the local mesh directory (dbDir()/meshSubDir)
fileName meshDir() const;
//- Return reference to time
const Time& time() const;
//- Return reference to time
const Time& time() const;
//- Return the current instance directory for points
// Used in the construction of geometric mesh data dependent
// on points
const fileName& pointsInstance() const;
//- Return the current instance directory for points
// Used in the construction of geometric mesh data dependent
// on points
const fileName& pointsInstance() const;
//- Return the current instance directory for faces
const fileName& facesInstance() const;
//- Return the current instance directory for faces
const fileName& facesInstance() const;
// Communication support
// Communication support
//- Return communicator used for parallel communication
inline label comm() const noexcept;
//- Return communicator used for parallel communication
inline label comm() const noexcept;
//- Return communicator used for parallel communication
inline label& comm() noexcept;
//- Return communicator used for parallel communication
inline label& comm() noexcept;
// Mesh size parameters
// Access: Mesh size parameters
//- Number of local mesh points
inline label nPoints() const noexcept;
//- Number of local mesh points
inline label nPoints() const noexcept;
//- Number of local mesh edges
inline label nEdges() const noexcept;
//- Number of local mesh edges
inline label nEdges() const noexcept;
//- Number of internal faces
inline label nInternalEdges() const noexcept;
//- Number of internal faces
inline label nInternalEdges() const noexcept;
//- Number of boundary edges (== nEdges - nInternalEdges)
inline label nBoundaryEdges() const noexcept;
//- Number of boundary edges (== nEdges - nInternalEdges)
inline label nBoundaryEdges() const noexcept;
//- Number of patch faces
inline label nFaces() const noexcept;
//- Number of patch faces
inline label nFaces() const noexcept;
// Primitive mesh data
// Access: Primitive mesh data
//- Return local patch points
inline const pointField& points() const;
//- Return local points
inline const pointField& points() const;
//- Return local patch edges with reordered boundary
inline const edgeList& edges() const noexcept;
//- Return local edges with reordered boundary
inline const edgeList& edges() const noexcept;
//- Return local patch faces
inline const faceList& faces() const;
//- Sub-list of local internal edges
inline const edgeList::subList internalEdges() const;
//- Edge owner addressing
inline const labelList& edgeOwner() const noexcept;
//- Return local faces
inline const faceList& faces() const;
//- Edge neighbour addressing
inline const labelList& edgeNeighbour() const noexcept;
//- Edge owner addressing
inline const labelList& edgeOwner() const noexcept;
//- Edge neighbour addressing
inline const labelList& edgeNeighbour() const noexcept;
//- True if the internalEdges use an ordering that does not
//- correspond 1-to-1 with the patch internalEdges.
inline bool hasInternalEdgeLabels() const noexcept;
// Access
@ -649,7 +674,7 @@ public:
virtual const objectRegistry& thisDb() const;
//- Name function is needed to disambiguate those inherited
// from base classes
//- from base classes
const word& name() const
{
return thisDb().name();
@ -688,7 +713,7 @@ public:
}
//- True if given edge label is internal to the mesh
bool isInternalEdge(const label edgeIndex) const
bool isInternalEdge(const label edgeIndex) const noexcept
{
return (edgeIndex < nInternalEdges_);
}
@ -817,6 +842,9 @@ public:
boolList& correctPatchPointNormals() const;
// Storage management
//- Write mesh
virtual bool write(const bool valid = true) const;

708
src/finiteArea/faMesh/faMeshDemandDrivenData.C
View File

@ -37,6 +37,7 @@ License
#include "scalarMatrices.H"
#include "processorFaPatchFields.H"
#include "emptyFaPatchFields.H"
#include "triPointRef.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
@ -102,6 +103,121 @@ static inline vector areaInvDistSqrWeightedNormalDualEdge
);
}
// Calculate transform tensor with reference vector (unitAxis1)
// and direction vector (axis2).
//
// This is nearly identical to the meshTools axesRotation
// with an E3_E1 transformation with the following exceptions:
//
// - axis1 (e3 == unitAxis1): is already normalized (unit vector)
// - axis2 (e1 == dirn): no difference
// - transformation is row-vectors, not column-vectors
static inline tensor rotation_e3e1
(
const vector& unitAxis1,
vector dirn
)
{
dirn.removeCollinear(unitAxis1);
dirn.normalise();
// Set row vectors
return tensor
(
dirn,
(unitAxis1^dirn),
unitAxis1
);
}
// Simple area-weighted normal calculation for the specified edge vector
// and its owner/neighbour face centres (internal edges).
//
// Uses four triangles since adjacent faces may be non-planar
// and/or the edge centre is skewed from the face centres.
/*---------------------------------------------------------------------------*\
* * |
* /|\ | Triangles (0,1) are on the owner side.
* / | \ |
* / | \ | Triangles (2,3) are on the neighbour side.
* / 1|3 \ |
* own *----|----* nbr | Boundary edges are the same, but without a neighbour
* \ 0|2 / |
* \ | / |
* \ | / |
* \|/ |
* * |
\*---------------------------------------------------------------------------*/
static inline vector calcEdgeNormalFromFace
(
const linePointRef& edgeVec,
const point& ownCentre,
const point& neiCentre
)
{
const scalar magEdge(edgeVec.mag());
if (magEdge < ROOTVSMALL)
{
return Zero;
}
const vector edgeCtr = edgeVec.centre();
vector edgeNorm
(
// owner
triPointRef(edgeCtr, ownCentre, edgeVec.first()).areaNormal()
+ triPointRef(edgeCtr, edgeVec.last(), ownCentre).areaNormal()
// neighbour
+ triPointRef(edgeCtr, edgeVec.first(), neiCentre).areaNormal()
+ triPointRef(edgeCtr, neiCentre, edgeVec.last()).areaNormal()
);
// Requires a unit-vector (already tested its mag)
edgeNorm.removeCollinear(edgeVec.vec()/magEdge);
edgeNorm.normalise();
// Scale with the original edge length
return magEdge*edgeNorm;
}
// As above, but for boundary edgess (no neighbour)
static inline vector calcEdgeNormalFromFace
(
const linePointRef& edgeVec,
const point& ownCentre
)
{
const scalar magEdge(edgeVec.mag());
if (magEdge < ROOTVSMALL)
{
return Zero;
}
const vector edgeCtr = edgeVec.centre();
vector edgeNorm
(
// owner
triPointRef(edgeCtr, ownCentre, edgeVec.first()).areaNormal()
+ triPointRef(edgeCtr, edgeVec.last(), ownCentre).areaNormal()
);
// Requires a unit-vector (already tested its mag)
edgeNorm.removeCollinear(edgeVec.vec()/magEdge);
edgeNorm.normalise();
// Scale with the original edge length
return magEdge*edgeNorm;
}
} // End namespace Foam
@ -189,71 +305,187 @@ void Foam::faMesh::calcLe() const
edgeVectorField& Le = *LePtr_;
const pointField& localPoints = points();
const pointField& pPoints = points();
const edgeList& pEdges = edges();
if (faMesh::geometryOrder() < 1)
{
// Simple (primitive) edge normal calculations.
// These are primarly designed to avoid any communication
// but are thus necessarily inconsistent across processor boundaries!
// Reasonable to assume that the volume mesh already has faceCentres
// eg, used magSf somewhere.
// Can use these instead of triggering our calcAreaCentres().
WarningInFunction
<< "Using geometryOrder < 1 : "
"simplified edge area-normals for Le() calculation"
<< endl;
UIndirectList<vector> fCentres(mesh().faceCentres(), faceLabels());
// Flat addressing
vectorField edgeNormals(nEdges_);
// Internal (edge normals)
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
{
edgeNormals[edgei] =
calcEdgeNormalFromFace
(
edges_[edgei].line(localPoints),
fCentres[owner()[edgei]],
fCentres[neighbour()[edgei]]
);
}
// Boundary (edge normals). Like above, but only has owner
for (label edgei = nInternalEdges_; edgei < nEdges_; ++edgei)
{
edgeNormals[edgei] =
calcEdgeNormalFromFace
(
edges_[edgei].line(localPoints),
fCentres[owner()[edgei]]
);
}
// Now use these edge normals for calculating Le
// Internal (edge vector)
{
vectorField& internalFld = Le.ref();
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
{
vector& leVector = internalFld[edgei];
vector edgeVec = edges_[edgei].vec(localPoints);
const scalar magEdge(mag(edgeVec));
if (magEdge < ROOTVSMALL)
{
// Too small
leVector = Zero;
continue;
}
const vector edgeCtr = edges_[edgei].centre(localPoints);
const vector& edgeNorm = edgeNormals[edgei];
const vector& ownCentre = fCentres[owner()[edgei]];
leVector = magEdge*normalised(edgeVec ^ edgeNorm);
leVector *= sign(leVector & (edgeCtr - ownCentre));
}
}
// Boundary (edge vector)
forAll(boundary(), patchi)
{
const labelUList& bndEdgeFaces = boundary()[patchi].edgeFaces();
const edgeList::subList bndEdges =
boundary()[patchi].patchSlice(edges_);
vectorField& patchLe = Le.boundaryFieldRef()[patchi];
forAll(patchLe, bndEdgei)
{
vector& leVector = patchLe[bndEdgei];
const label meshEdgei(boundary()[patchi].start() + bndEdgei);
vector edgeVec = bndEdges[bndEdgei].vec(localPoints);
const scalar magEdge(mag(edgeVec));
if (magEdge < ROOTVSMALL)
{
// Too small
leVector = Zero;
continue;
}
const vector edgeCtr = bndEdges[bndEdgei].centre(localPoints);
const vector& edgeNorm = edgeNormals[meshEdgei];
const vector& ownCentre = fCentres[bndEdgeFaces[bndEdgei]];
leVector = magEdge*normalised(edgeVec ^ edgeNorm);
leVector *= sign(leVector & (edgeCtr - ownCentre));
}
}
// Done
return;
}
// Longer forms.
// Using edgeAreaNormals, which uses pointAreaNormals (communication!)
const edgeVectorField& eCentres = edgeCentres();
const areaVectorField& fCentres = areaCentres();
const edgeVectorField& edgeNormals = edgeAreaNormals();
vectorField& leInternal = Le.ref();
const vectorField& edgeNormalsInternal = edgeNormals.internalField();
const vectorField& fCentresInternal = fCentres.internalField();
const vectorField& eCentresInternal = eCentres.internalField();
const scalarField& magLeInternal = magLe().internalField();
// Internal (edge vector)
forAll(leInternal, edgeI)
{
leInternal[edgeI] =
pEdges[edgeI].vec(pPoints) ^ edgeNormalsInternal[edgeI];
vectorField& internalFld = Le.ref();
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
{
vector& leVector = internalFld[edgei];
leInternal[edgeI] *=
- sign
(
leInternal[edgeI] &
(
fCentresInternal[owner()[edgeI]]
- eCentresInternal[edgeI]
)
);
vector edgeVec = edges_[edgei].vec(localPoints);
const scalar magEdge(mag(edgeVec));
leInternal[edgeI] *=
magLeInternal[edgeI]/mag(leInternal[edgeI]);
if (magEdge < ROOTVSMALL)
{
// Too small
leVector = Zero;
continue;
}
const vector& edgeCtr = eCentres[edgei];
const vector& edgeNorm = edgeNormals[edgei];
const vector& ownCentre = fCentres[owner()[edgei]];
leVector = magEdge*normalised(edgeVec ^ edgeNorm);
leVector *= sign(leVector & (edgeCtr - ownCentre));
}
}
forAll(boundary(), patchI)
forAll(boundary(), patchi)
{
const labelUList& bndEdgeFaces = boundary()[patchI].edgeFaces();
const labelUList& bndEdgeFaces = boundary()[patchi].edgeFaces();
const edgeList::subList bndEdges =
boundary()[patchI].patchSlice(pEdges);
boundary()[patchi].patchSlice(edges_);
const vectorField& bndEdgeNormals =
edgeNormals.boundaryField()[patchI];
edgeNormals.boundaryField()[patchi];
vectorField& patchLe = Le.boundaryFieldRef()[patchI];
const vectorField& patchECentres = eCentres.boundaryField()[patchI];
vectorField& patchLe = Le.boundaryFieldRef()[patchi];
const vectorField& patchECentres = eCentres.boundaryField()[patchi];
forAll(patchLe, edgeI)
forAll(patchLe, bndEdgei)
{
patchLe[edgeI] =
bndEdges[edgeI].vec(pPoints) ^ bndEdgeNormals[edgeI];
vector& leVector = patchLe[bndEdgei];
patchLe[edgeI] *=
- sign
(
patchLe[edgeI]&
(
fCentresInternal[bndEdgeFaces[edgeI]]
- patchECentres[edgeI]
)
);
vector edgeVec = bndEdges[bndEdgei].vec(localPoints);
const scalar magEdge(mag(edgeVec));
patchLe[edgeI] *=
magLe().boundaryField()[patchI][edgeI]
/mag(patchLe[edgeI]);
if (magEdge < ROOTVSMALL)
{
// Too small
leVector = Zero;
continue;
}
const vector& edgeCtr = patchECentres[bndEdgei];
const vector& edgeNorm = bndEdgeNormals[bndEdgei];
const vector& ownCentre = fCentres[bndEdgeFaces[bndEdgei]];
leVector = magEdge*normalised(edgeVec ^ edgeNorm);
leVector *= sign(leVector & (edgeCtr - ownCentre));
}
}
}
@ -289,23 +521,30 @@ void Foam::faMesh::calcMagLe() const
const pointField& localPoints = points();
label edgei = 0;
for (const edge& e : patch().internalEdges())
// Internal (edge length)
{
magLe.ref()[edgei] = e.mag(localPoints);
++edgei;
auto iter = magLe.primitiveFieldRef().begin();
for (const edge& e : internalEdges())
{
*iter = e.mag(localPoints);
++iter;
}
}
forAll(boundary(), patchI)
// Boundary (edge length)
{
const edgeList::subList patchEdges =
boundary()[patchI].patchSlice(edges());
auto& bfld = magLe.boundaryFieldRef();
forAll(patchEdges, edgeI)
forAll(boundary(), patchi)
{
magLe.boundaryFieldRef()[patchI][edgeI] =
patchEdges[edgeI].mag(localPoints);
auto iter = bfld[patchi].begin();
for (const edge& e : boundary()[patchi].patchSlice(edges_))
{
*iter = e.mag(localPoints);
++iter;
}
}
}
}
@ -343,21 +582,25 @@ void Foam::faMesh::calcAreaCentres() const
const faceList& localFaces = faces();
// Internal (face centres)
forAll(localFaces, faceI)
// Could also obtain from volume mesh faceCentres()
forAll(localFaces, facei)
{
centres.ref()[faceI] = localFaces[faceI].centre(localPoints);
centres.ref()[facei] = localFaces[facei].centre(localPoints);
}
// Boundary (edge centres)
forAll(boundary(), patchI)
{
const edgeList::subList patchEdges =
boundary()[patchI].patchSlice(edges());
auto& bfld = centres.boundaryFieldRef();
forAll(patchEdges, edgeI)
forAll(boundary(), patchi)
{
centres.boundaryFieldRef()[patchI][edgeI] =
patchEdges[edgeI].centre(localPoints);
auto iter = bfld[patchi].begin();
for (const edge& e : boundary()[patchi].patchSlice(edges_))
{
*iter = e.centre(localPoints);
++iter;
}
}
}
}
@ -389,28 +632,34 @@ void Foam::faMesh::calcEdgeCentres() const
dimLength
);
edgeVectorField& edgeCentres = *edgeCentresPtr_;
edgeVectorField& centres = *edgeCentresPtr_;
const pointField& localPoints = points();
// Internal edges
label edgei = 0;
for (const edge& e : patch().internalEdges())
// Internal (edge centres)
{
edgeCentres.ref()[edgei] = e.centre(localPoints);
++edgei;
auto iter = centres.primitiveFieldRef().begin();
for (const edge& e : internalEdges())
{
*iter = e.centre(localPoints);
++iter;
}
}
// Boundary edges
forAll(boundary(), patchI)
// Boundary (edge centres)
{
const edgeList::subList patchEdges =
boundary()[patchI].patchSlice(edges());
auto& bfld = centres.boundaryFieldRef();
forAll(patchEdges, edgeI)
forAll(boundary(), patchi)
{
edgeCentres.boundaryFieldRef()[patchI][edgeI] =
patchEdges[edgeI].centre(localPoints);
auto iter = bfld[patchi].begin();
for (const edge& e : boundary()[patchi].patchSlice(edges_))
{
*iter = e.centre(localPoints);
++iter;
}
}
}
}
@ -446,9 +695,10 @@ void Foam::faMesh::calcS() const
const pointField& localPoints = points();
const faceList& localFaces = faces();
forAll(S, faceI)
// Could also obtain from volume mesh faceAreas()
forAll(S, facei)
{
S[faceI] = localFaces[faceI].mag(localPoints);
S[facei] = localFaces[facei].mag(localPoints);
}
}
@ -485,6 +735,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
const faceList& localFaces = faces();
// Internal (faces)
// Could also obtain from volume mesh Sf() + normalise
vectorField& nInternal = faceNormals.ref();
forAll(localFaces, faceI)
{
@ -493,8 +744,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
// Boundary - copy from edges
const edgeVectorField::Boundary& edgeNormalsBoundary =
edgeAreaNormals().boundaryField();
const auto& edgeNormalsBoundary = edgeAreaNormals().boundaryField();
forAll(boundary(), patchI)
{
@ -528,46 +778,50 @@ void Foam::faMesh::calcEdgeAreaNormals() const
*this,
dimless
);
edgeVectorField& edgeAreaNormals = *edgeAreaNormalsPtr_;
// Point area normals
// Starting from point area normals
const vectorField& pointNormals = pointAreaNormals();
// Internal edges
forAll(edgeAreaNormals.internalField(), edgei)
{
const edge& e = edges()[edgei];
const edge& e = edges_[edgei];
const vector edgeVec = e.unitVec(points());
vector& n = edgeAreaNormals.ref()[edgei];
vector& edgeNorm = edgeAreaNormals.ref()[edgei];
n = (pointNormals[e.first()] + pointNormals[e.second()]);
// Average of both ends
edgeNorm = (pointNormals[e.first()] + pointNormals[e.second()]);
n -= edgeVec*(edgeVec & n);
n.normalise();
edgeNorm.removeCollinear(edgeVec);
edgeNorm.normalise();
}
// Boundary
auto& bfld = edgeAreaNormals.boundaryFieldRef();
forAll(boundary(), patchi)
{
auto& pfld = bfld[patchi];
const edgeList::subList patchEdges =
boundary()[patchi].patchSlice(edges());
boundary()[patchi].patchSlice(edges_);
vectorField& edgeNorms = edgeAreaNormals.boundaryFieldRef()[patchi];
forAll(patchEdges, edgei)
forAll(patchEdges, bndEdgei)
{
const edge& e = patchEdges[edgei];
const edge& e = patchEdges[bndEdgei];
const vector edgeVec = e.unitVec(points());
vector& n = edgeNorms[edgei];
vector& edgeNorm = pfld[bndEdgei];
n = (pointNormals[e.first()] + pointNormals[e.second()]);
// Average of both ends
edgeNorm = (pointNormals[e.first()] + pointNormals[e.second()]);
n -= edgeVec*(edgeVec & n);
n.normalise();
edgeNorm.removeCollinear(edgeVec);
edgeNorm.normalise();
}
}
}
@ -624,9 +878,14 @@ void Foam::faMesh::calcEdgeTransformTensors() const
<< abort(FatalError);
}
edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges());
FieldField<Field, tensor>& edgeTransformTensors =
*edgeTransformTensorsPtr_;
edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges_);
auto& edgeTransformTensors = *edgeTransformTensorsPtr_;
// Initialize all transformation tensors
for (label edgei = 0; edgei < nEdges_; ++edgei)
{
edgeTransformTensors.set(edgei, new Field<tensor>(3, tensor::I));
}
const areaVectorField& Nf = faceAreaNormals();
const areaVectorField& Cf = areaCentres();
@ -635,193 +894,129 @@ void Foam::faMesh::calcEdgeTransformTensors() const
const edgeVectorField& Ce = edgeCentres();
// Internal edges transformation tensors
for (label edgeI=0; edgeI<nInternalEdges(); ++edgeI)
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
{
edgeTransformTensors.set(edgeI, new Field<tensor>(3, I));
const label ownFacei = owner()[edgei];
const label neiFacei = neighbour()[edgei];
auto& tensors = edgeTransformTensors[edgei];
vector E = Ce.internalField()[edgeI];
vector edgeCtr = Ce.internalField()[edgei];
if (skew())
{
E -= skewCorrectionVectors().internalField()[edgeI];
edgeCtr -= skewCorrectionVectors().internalField()[edgei];
}
// Edge transformation tensor
vector il = E - Cf.internalField()[owner()[edgeI]];
il -= Ne.internalField()[edgeI]
*(Ne.internalField()[edgeI]&il);
il /= mag(il);
vector kl = Ne.internalField()[edgeI];
vector jl = kl^il;
edgeTransformTensors[edgeI][0] =
tensor
tensors[0] =
rotation_e3e1
(
il.x(), il.y(), il.z(),
jl.x(), jl.y(), jl.z(),
kl.x(), kl.y(), kl.z()
Ne.internalField()[edgei],
(edgeCtr - Cf.internalField()[ownFacei])
);
// Owner transformation tensor
il = E - Cf.internalField()[owner()[edgeI]];
il -= Nf.internalField()[owner()[edgeI]]
*(Nf.internalField()[owner()[edgeI]]&il);
il /= mag(il);
kl = Nf.internalField()[owner()[edgeI]];
jl = kl^il;
edgeTransformTensors[edgeI][1] =
tensor
tensors[1] =
rotation_e3e1
(
il.x(), il.y(), il.z(),
jl.x(), jl.y(), jl.z(),
kl.x(), kl.y(), kl.z()
Nf.internalField()[ownFacei],
(edgeCtr - Cf.internalField()[ownFacei])
);
// Neighbour transformation tensor
il = Cf.internalField()[neighbour()[edgeI]] - E;
il -= Nf.internalField()[neighbour()[edgeI]]
*(Nf.internalField()[neighbour()[edgeI]]&il);
il /= mag(il);
kl = Nf.internalField()[neighbour()[edgeI]];
jl = kl^il;
edgeTransformTensors[edgeI][2] =
tensor
tensors[2] =
rotation_e3e1
(
il.x(), il.y(), il.z(),
jl.x(), jl.y(), jl.z(),
kl.x(), kl.y(), kl.z()
Nf.internalField()[neiFacei],
(Cf.internalField()[neiFacei] - edgeCtr)
);
}
// Boundary edges transformation tensors
forAll(boundary(), patchI)
forAll(boundary(), patchi)
{
if (boundary()[patchI].coupled())
const labelUList& edgeFaces = boundary()[patchi].edgeFaces();
if (boundary()[patchi].coupled())
{
const labelUList& edgeFaces =
boundary()[patchI].edgeFaces();
vectorField nbrCf(Cf.boundaryField()[patchi].patchNeighbourField());
vectorField nbrNf(Nf.boundaryField()[patchi].patchNeighbourField());
vectorField ngbCf(Cf.boundaryField()[patchI].patchNeighbourField());
vectorField ngbNf(Nf.boundaryField()[patchI].patchNeighbourField());
forAll(edgeFaces, edgeI)
forAll(edgeFaces, bndEdgei)
{
edgeTransformTensors.set
(
boundary()[patchI].start() + edgeI,
new Field<tensor>(3, I)
);
const label ownFacei = edgeFaces[bndEdgei];
const label meshEdgei(boundary()[patchi].start() + bndEdgei);
vector E = Ce.boundaryField()[patchI][edgeI];
auto& tensors = edgeTransformTensors[meshEdgei];
vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
if (skew())
{
E -= skewCorrectionVectors()
.boundaryField()[patchI][edgeI];
edgeCtr -= skewCorrectionVectors()
.boundaryField()[patchi][bndEdgei];
}
// Edge transformation tensor
vector il = E - Cf.internalField()[edgeFaces[edgeI]];
il -= Ne.boundaryField()[patchI][edgeI]
*(Ne.boundaryField()[patchI][edgeI]&il);
il /= mag(il);
vector kl = Ne.boundaryField()[patchI][edgeI];
vector jl = kl^il;
edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
tensor(il, jl, kl);
tensors[0] =
rotation_e3e1
(
Ne.boundaryField()[patchi][bndEdgei],
(edgeCtr - Cf.internalField()[ownFacei])
);
// Owner transformation tensor
il = E - Cf.internalField()[edgeFaces[edgeI]];
il -= Nf.internalField()[edgeFaces[edgeI]]
*(Nf.internalField()[edgeFaces[edgeI]]&il);
il /= mag(il);
kl = Nf.internalField()[edgeFaces[edgeI]];
jl = kl^il;
edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
tensor(il, jl, kl);
tensors[1] =
rotation_e3e1
(
Nf.internalField()[ownFacei],
(edgeCtr - Cf.internalField()[ownFacei])
);
// Neighbour transformation tensor
il = ngbCf[edgeI] - E;
il -= ngbNf[edgeI]*(ngbNf[edgeI]&il);
il /= mag(il);
kl = ngbNf[edgeI];
jl = kl^il;
edgeTransformTensors[boundary()[patchI].start() + edgeI][2] =
tensor(il, jl, kl);
tensors[2] =
rotation_e3e1
(
nbrNf[bndEdgei],
(nbrCf[bndEdgei] - edgeCtr)
);
}
}
else
{
const labelUList& edgeFaces = boundary()[patchI].edgeFaces();
forAll(edgeFaces, edgeI)
forAll(edgeFaces, bndEdgei)
{
edgeTransformTensors.set
(
boundary()[patchI].start() + edgeI,
new Field<tensor>(3, I)
);
const label ownFacei = edgeFaces[bndEdgei];
const label meshEdgei(boundary()[patchi].start() + bndEdgei);
vector E = Ce.boundaryField()[patchI][edgeI];
auto& tensors = edgeTransformTensors[meshEdgei];
vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
if (skew())
{
E -= skewCorrectionVectors()
.boundaryField()[patchI][edgeI];
edgeCtr -= skewCorrectionVectors()
.boundaryField()[patchi][bndEdgei];
}
// Edge transformation tensor
vector il = E - Cf.internalField()[edgeFaces[edgeI]];
il -= Ne.boundaryField()[patchI][edgeI]
*(Ne.boundaryField()[patchI][edgeI]&il);
il /= mag(il);
vector kl = Ne.boundaryField()[patchI][edgeI];
vector jl = kl^il;
edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
tensor(il, jl, kl);
tensors[0] =
rotation_e3e1
(
Ne.boundaryField()[patchi][bndEdgei],
(edgeCtr - Cf.internalField()[ownFacei])
);
// Owner transformation tensor
il = E - Cf.internalField()[edgeFaces[edgeI]];
tensors[1] =
rotation_e3e1
(
Nf.internalField()[ownFacei],
(edgeCtr - Cf.internalField()[ownFacei])
);
il -= Nf.internalField()[edgeFaces[edgeI]]
*(Nf.internalField()[edgeFaces[edgeI]]&il);
il /= mag(il);
kl = Nf.internalField()[edgeFaces[edgeI]];
jl = kl^il;
edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
tensor(il, jl, kl);
// Neighbour transformation tensor
tensors[2] = tensor::I;
}
}
}
@ -882,7 +1077,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
DebugInFunction
<< "Calculating pointAreaNormals : face dual method" << endl;
result.resize(nPoints());
result.resize_nocopy(nPoints());
result = Zero;
const pointField& points = patch().localPoints();
@ -945,7 +1140,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
for (const label pti : patchPoints)
{
result[pti] -= N*(N&result[pti]);
result[pti].removeCollinear(N);
}
// Axis point correction
@ -1124,17 +1319,13 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
forAllConstIters(fpNormals, iter)
{
const label pointi = iter.key();
vector fpnorm = iter.val();
vector fpnorm = normalised(iter.val());
fpnorm.normalise();
result[pointi] -= fpnorm*(fpnorm & result[pointi]);
result[pointi].removeCollinear(fpnorm);
}
}
for (vector& n : result)
{
n.normalise();
}
result.normalise();
}
@ -1469,7 +1660,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
const vector& origin = points[curPoint];
const vector axis = normalised(result[curPoint]);
vector dir(allPoints[0] - points[curPoint]);
dir -= axis*(axis&dir);
dir.removeCollinear(axis);
dir.normalise();
const coordinateSystem cs("cs", origin, axis, dir);
@ -1633,7 +1824,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
const vector& origin = points[curPoint];
const vector axis = normalised(result[curPoint]);
vector dir(allPoints[0] - points[curPoint]);
dir -= axis*(axis&dir);
dir.removeCollinear(axis);
dir.normalise();
coordinateSystem cs("cs", origin, axis, dir);
@ -1710,22 +1901,19 @@ Foam::tmp<Foam::edgeScalarField> Foam::faMesh::edgeLengthCorrection() const
DebugInFunction
<< "Calculating edge length correction" << endl;
tmp<edgeScalarField> tcorrection
auto tcorrection = tmp<edgeScalarField>::New
(
new edgeScalarField
IOobject
(
IOobject
(
"edgeLengthCorrection",
mesh().pointsInstance(),
meshSubDir,
mesh()
),
*this,
dimless
)
"edgeLengthCorrection",
mesh().pointsInstance(),
meshSubDir,
mesh()
),
*this,
dimless
);
edgeScalarField& correction = tcorrection.ref();
auto& correction = tcorrection.ref();
const vectorField& pointNormals = pointAreaNormals();

16
src/finiteArea/faMesh/faMeshI.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2021 OpenCFD Ltd.
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -73,7 +73,7 @@ inline Foam::label Foam::faMesh::nInternalEdges() const noexcept
inline Foam::label Foam::faMesh::nBoundaryEdges() const noexcept
{
return nEdges_ - nInternalEdges_;
return (nEdges_ - nInternalEdges_);
}
@ -95,6 +95,12 @@ inline const Foam::edgeList& Foam::faMesh::edges() const noexcept
}
inline const Foam::edgeList::subList Foam::faMesh::internalEdges() const
{
return edgeList::subList(edges_, nInternalEdges_);
}
inline const Foam::faceList& Foam::faMesh::faces() const
{
return patch().localFaces();
@ -139,6 +145,12 @@ inline Foam::uindirectPrimitivePatch& Foam::faMesh::patch()
}
inline bool Foam::faMesh::hasInternalEdgeLabels() const noexcept
{
return false;
}
inline const Foam::List<Foam::labelPair>&
Foam::faMesh::boundaryConnections() const
{

152
src/finiteArea/faMesh/faMeshSubset/faMeshSubset.C Normal file
View File

@ -0,0 +1,152 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "faMeshSubset.H"
#include "boolList.H"
#include "BitOps.H"
#include "Pstream.H"
#include "emptyFaPatch.H"
#include "cyclicFaPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::word Foam::faMeshSubset::exposedPatchName("oldInternalEdges");
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::faMeshSubset::checkHasSubMesh() const
{
if (!subMeshPtr_)
{
FatalErrorInFunction
<< "Mesh is not subsetted!" << nl
<< abort(FatalError);
return false;
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::faMeshSubset::faMeshSubset(const faMesh& baseMesh)
:
baseMesh_(baseMesh),
subMeshPtr_(nullptr),
edgeFlipMapPtr_(nullptr),
pointMap_(),
faceMap_(),
cellMap_(),
patchMap_()
{}
Foam::faMeshSubset::faMeshSubset(const faMesh& baseMesh, const Foam::zero)
:
faMeshSubset(baseMesh)
{
reset(Foam::zero{});
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::faMeshSubset::clear()
{
subMeshPtr_.reset(nullptr);
edgeFlipMapPtr_.reset(nullptr);
pointMap_.clear();
faceMap_.clear();
cellMap_.clear();
patchMap_.clear();
}
void Foam::faMeshSubset::reset()
{
clear();
}
void Foam::faMeshSubset::reset(const Foam::zero)
{
clear();
// Create zero-sized subMesh
subMeshPtr_.reset
(
new faMesh
(
baseMesh_.mesh(), // The polyMesh
// IOobject
// (
// baseMesh_.name(),
// baseMesh_.time().timeName(),
// baseMesh_.time(),
// IOobject::READ_IF_PRESENT, // Read fa* if present
// IOobject::NO_WRITE
// ),
Foam::zero{} // zero-sized
)
);
auto& newSubMesh = subMeshPtr_();
// Clone non-processor patches
{
const faBoundaryMesh& oldBoundary = baseMesh_.boundary();
const faBoundaryMesh& newBoundary = newSubMesh.boundary();
faPatchList newPatches(oldBoundary.nNonProcessor());
patchMap_ = identity(newPatches.size());
forAll(newPatches, patchi)
{
newPatches.set
(
patchi,
oldBoundary[patchi].clone
(
newBoundary,
labelList(), // edgeLabels
patchi,
oldBoundary[patchi].ngbPolyPatchIndex()
)
);
}
newSubMesh.addFaPatches(newPatches);
}
}
// ************************************************************************* //

222
src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H Normal file
View File

@ -0,0 +1,222 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::faMeshSubset
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.
Caution
Currently not really functional for subsetting beyond handling
a simple zero-sized mesh.
SourceFiles
faMeshSubset.C
faMeshSubsetI.H
faMeshSubsetTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_faMeshSubset_H
#define Foam_faMeshSubset_H
#include "areaFaMesh.H"
#include "edgeFaMesh.H"
#include "GeometricField.H"
#include "bitSet.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward Declarations
class mapDistributePolyMesh;
/*---------------------------------------------------------------------------*\
Class faMeshSubset Declaration
\*---------------------------------------------------------------------------*/
class faMeshSubset
{
// Private Data
//- The base mesh to subset from
const faMesh& baseMesh_;
//- Demand-driven subset mesh (pointer)
autoPtr<faMesh> subMeshPtr_;
//- Optional edge mapping array with flip encoded (-1/+1)
mutable autoPtr<labelList> edgeFlipMapPtr_;
//- Point mapping array
labelList pointMap_;
//- Face mapping array
labelList faceMap_;
//- Cell mapping array
labelList cellMap_;
//- Patch mapping array
labelList patchMap_;
// Private Member Functions
//- Calculate face flip map
void calcEdgeFlipMap() const;
protected:
// Protected Member Functions
//- FatalError if subset has not been performed
bool checkHasSubMesh() const;
//- No copy construct
faMeshSubset(const faMeshSubset&) = delete;
//- No copy assignment
void operator=(const faMeshSubset&) = delete;
public:
// Static Data Members
//- Name for exposed internal edges (default: oldInternalEdges)
static word exposedPatchName;
// Constructors
//- Construct using the entire mesh (no subset)
explicit faMeshSubset(const faMesh& baseMesh);
//- Construct a zero-sized subset mesh, non-processor patches only
faMeshSubset(const faMesh& baseMesh, const Foam::zero);
// Member Functions
// Access
//- Original mesh
inline const faMesh& baseMesh() const noexcept;
//- Return baseMesh or subMesh, depending on the current state.
inline const faMesh& mesh() const noexcept;
//- Have subMesh?
inline bool hasSubMesh() const noexcept;
//- Return reference to subset mesh
inline const faMesh& subMesh() const;
//- Return reference to subset mesh
inline faMesh& subMesh();
//- Return point map
inline const labelList& pointMap() const;
//- Return face map
inline const labelList& faceMap() const;
//- Return edge map with sign to encode flipped edges
inline const labelList& edgeFlipMap() const;
//- Return cell map
inline const labelList& cellMap() const;
//- Return patch map
inline const labelList& patchMap() const;
// Edit
//- Reset subMesh and all maps
void clear();
//- Reset subMesh and all maps. Same as clear()
void reset();
//- Reset to a zero-sized subset mesh, non-processor patches only
void reset(const Foam::zero);
// Field Mapping (static functions)
//- Map area field.
// Optionally allow unmapped faces not to produce a warning
template<class Type>
static tmp<GeometricField<Type, faPatchField, areaMesh>>
interpolate
(
const GeometricField<Type, faPatchField, areaMesh>&,
const faMesh& sMesh,
const bool allowUnmapped = false
);
// Field Mapping
//- Map area field.
// Optionally allow unmapped faces not to produce a warning
template<class Type>
tmp<GeometricField<Type, faPatchField, areaMesh>>
interpolate
(
const GeometricField<Type, faPatchField, areaMesh>&,
const bool allowUnmapped = false
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "faMeshSubsetI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "faMeshSubsetTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

107
src/finiteArea/faMesh/faMeshSubset/faMeshSubsetI.H Normal file
View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::faMesh& Foam::faMeshSubset::baseMesh() const noexcept
{
return baseMesh_;
}
inline const Foam::faMesh& Foam::faMeshSubset::mesh() const noexcept
{
return subMeshPtr_ ? *subMeshPtr_ : baseMesh_;
}
inline bool Foam::faMeshSubset::hasSubMesh() const noexcept
{
return bool(subMeshPtr_);
}
inline const Foam::faMesh& Foam::faMeshSubset::subMesh() const
{
checkHasSubMesh();
return *subMeshPtr_;
}
inline Foam::faMesh& Foam::faMeshSubset::subMesh()
{
checkHasSubMesh();
return *subMeshPtr_;
}
inline const Foam::labelList& Foam::faMeshSubset::pointMap() const
{
checkHasSubMesh();
return pointMap_;
}
inline const Foam::labelList& Foam::faMeshSubset::faceMap() const
{
checkHasSubMesh();
return faceMap_;
}
inline const Foam::labelList& Foam::faMeshSubset::edgeFlipMap() const
{
if (!edgeFlipMapPtr_)
{
// calcEdgeFlipMap();
}
return *edgeFlipMapPtr_;
}
inline const Foam::labelList& Foam::faMeshSubset::cellMap() const
{
checkHasSubMesh();
return cellMap_;
}
inline const Foam::labelList& Foam::faMeshSubset::patchMap() const
{
checkHasSubMesh();
return patchMap_;
}
// ************************************************************************* //

169
src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C Normal file
View File

@ -0,0 +1,169 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "faMeshSubset.H"
#include "areaFaMesh.H"
#include "edgeFaMesh.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "emptyFaPatchFields.H"
#include "directFaPatchFieldMapper.H"
#include "flipOp.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
template<class Type>
Foam::tmp
<
Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>
>
Foam::faMeshSubset::interpolate
(
const GeometricField<Type, faPatchField, areaMesh>& vf,
const faMesh& sMesh,
const bool allowUnmapped
)
{
// 1. Create the complete field with dummy patch fields
PtrList<faPatchField<Type>> patchFields(sMesh.boundary().size());
forAll(patchFields, patchi)
{
patchFields.set
(
patchi,
faPatchField<Type>::New
(
calculatedFaPatchField<Type>::typeName,
sMesh.boundary()[patchi],
DimensionedField<Type, areaMesh>::null()
)
);
}
auto tresult = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"subset"+vf.name(),
sMesh.time().timeName(),
sMesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sMesh,
vf.dimensions(),
Field<Type>(),
// Field<Type>(vf.primitiveField(), cellMap),
patchFields
);
auto& result = tresult.ref();
result.oriented() = vf.oriented();
// 2. Change the faPatchFields to the correct type using a mapper
// constructor (with reference to the now correct internal field)
auto& bf = result.boundaryFieldRef();
forAll(bf, patchi)
{
// Construct addressing
const faPatch& subPatch = sMesh.boundary()[patchi];
labelList directAddressing;
directFaPatchFieldMapper mapper(directAddressing);
// allowUnmapped : special mode for if we do not want to be
// warned for unmapped faces (e.g. from faMeshDistribute).
const bool hasUnmapped = mapper.hasUnmapped();
if (allowUnmapped)
{
mapper.hasUnmapped() = false;
}
bf.set
(
patchi,
faPatchField<Type>::New
(
vf.boundaryField()[patchi],
subPatch,
result(),
mapper
)
);
if (allowUnmapped && hasUnmapped)
{
// Set unmapped values to zeroGradient. This is the default
// action for unmapped faPatchFields. Note that this bypasses
// any special logic for handling unmapped faPatchFields but
// since this is only used inside faMeshDistribute ...
tmp<Field<Type>> tfld(bf[patchi].patchInternalField());
const Field<Type>& fld = tfld();
Field<Type> value(bf[patchi]);
forAll(directAddressing, i)
{
if (directAddressing[i] == -1)
{
value[i] = fld[i];
}
}
bf[patchi].faPatchField<Type>::operator=(value);
}
}
return tresult;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp
<
Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>
>
Foam::faMeshSubset::interpolate
(
const GeometricField<Type, faPatchField, areaMesh>& vf,
const bool allowUnmapped
) const
{
if (subMeshPtr_)
{
return interpolate(vf, *subMeshPtr_);
}
return vf;
}
// ************************************************************************* //

470
src/finiteArea/faMesh/faMeshTools/faMeshTools.C Normal file
View File

@ -0,0 +1,470 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2012-2016 OpenFOAM Foundation
Copyright (C) 2015-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 "faMeshTools.H"
#include "faBoundaryMeshEntries.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "polyMesh.H"
#include "processorFaPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::faMeshTools::unregisterMesh(const faMesh& mesh)
{
auto& obr = const_cast<objectRegistry&>(mesh.thisDb());
// Checkout by name (casting ambiguity)
obr.checkOut(faMesh::typeName);
obr.checkOut("faBoundaryMesh");
obr.checkOut("faSchemes");
obr.checkOut("faSolution");
}
void Foam::faMeshTools::forceDemandDriven(faMesh& mesh)
{
(void)mesh.globalData();
(void)mesh.Le();
(void)mesh.magLe();
(void)mesh.areaCentres();
(void)mesh.edgeCentres();
(void)mesh.faceAreaNormals();
(void)mesh.edgeAreaNormals();
(void)mesh.pointAreaNormals();
(void)mesh.faceCurvatures();
(void)mesh.edgeTransformTensors();
}
Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
(
const IOobject& io,
const polyMesh& pMesh,
const bool masterOnlyReading,
const bool verbose
)
{
// Region name
// ~~~~~~~~~~~
const fileName meshSubDir
(
(pMesh.name() == polyMesh::defaultRegion ? word::null : pMesh.name())
/ faMesh::meshSubDir
);
fileName facesInstance;
// Patch types
// ~~~~~~~~~~~
// Read and scatter master patches (without reading master mesh!)
PtrList<entry> patchEntries;
if (Pstream::master())
{
const bool oldParRun = Pstream::parRun(false);
facesInstance = io.time().findInstance
(
meshSubDir,
"faceLabels",
IOobject::MUST_READ
);
patchEntries = faBoundaryMeshEntries
(
IOobject
(
"faBoundary",
facesInstance,
meshSubDir,
io.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
Pstream::parRun(oldParRun);
}
// Broadcast information to all
Pstream::broadcasts
(
UPstream::worldComm,
patchEntries,
facesInstance
);
// Dummy meshes
// ~~~~~~~~~~~~
// Set up to read-if-present. Note: does not search for mesh so set
// instance explicitly
IOobject meshIO(io);
meshIO.instance() = facesInstance;
meshIO.readOpt(IOobject::READ_IF_PRESENT);
// For mesh components (faceLabels, ...)
IOobject cmptIO(meshIO, "faceLabels", meshSubDir);
cmptIO.readOpt(IOobject::MUST_READ);
cmptIO.writeOpt(IOobject::NO_WRITE);
cmptIO.registerObject(false);
// Check who has a mesh
const fileName meshDir = io.time().path()/facesInstance/meshSubDir;
bool haveMesh = isDir(meshDir);
if (masterOnlyReading && !Pstream::master())
{
haveMesh = false;
meshIO.readOpt(IOobject::NO_READ);
}
if (!haveMesh)
{
cmptIO.readOpt(IOobject::NO_READ);
}
// Read mesh
// ~~~~~~~~~
// Now all processors use supplied points,faces etc
// Note: solution, schemes are also using the supplied IOobject so
// on slave will be NO_READ, on master READ_IF_PRESENT. This will
// conflict with e.g. timeStampMaster reading so switch off.
const auto oldCheckType = IOobject::fileModificationChecking;
IOobject::fileModificationChecking = IOobject::timeStamp;
// faceLabels
cmptIO.rename("faceLabels");
labelIOList faceLabels(cmptIO);
auto meshPtr = autoPtr<faMesh>::New
(
pMesh,
std::move(faceLabels),
meshIO
);
auto& mesh = *meshPtr;
IOobject::fileModificationChecking = oldCheckType;
// Some processors without patches? - add patches
if (returnReduce(mesh.boundary().empty(), orOp<bool>()))
{
// Use patchEntries, which were read on master and broadcast
faPatchList patches(patchEntries.size());
label nPatches = 0;
const bool isEmptyMesh = (mesh.faceLabels().empty());
forAll(patchEntries, patchi)
{
const entry& e = patchEntries[patchi];
const word type(e.dict().get<word>("type"));
const word& name = e.keyword();
if
(
type == processorFaPatch::typeName
)
{
// Stop at the first processor patch.
// - logic will not work with inter-mixed proc-patches anyhow
break;
}
else
{
dictionary patchDict(e.dict());
if (isEmptyMesh)
{
patchDict.set("edgeLabels", labelList());
}
patches.set
(
patchi,
faPatch::New
(
name,
patchDict,
nPatches++,
mesh.boundary()
)
);
}
}
patches.resize(nPatches);
mesh.addFaPatches(patches, false); // No parallel comms
}
// Recreate basic geometry, globalMeshData etc.
mesh.init(false);
(void)mesh.globalData();
return meshPtr;
}
Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
(
const IOobject& io,
const polyMesh& pMesh,
const bool decompose,
const bool verbose
)
{
// Region name
// ~~~~~~~~~~~
const fileName meshSubDir
(
(pMesh.name() == polyMesh::defaultRegion ? word::null : pMesh.name())
/ faMesh::meshSubDir
);
// Patch types
// ~~~~~~~~~~~
// Read and scatter master patches (without reading master mesh!)
PtrList<entry> patchEntries;
if (Pstream::master())
{
const bool oldParRun = Pstream::parRun(false);
patchEntries = faBoundaryMeshEntries
(
IOobject
(
"faBoundary",
io.instance(),
meshSubDir,
io.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
Pstream::parRun(oldParRun);
}
// Broadcast: send patches to all
Pstream::broadcast(patchEntries); // == worldComm;
/// Info<< patchEntries << nl;
// Dummy meshes
// ~~~~~~~~~~~~
// Check who has or needs a mesh.
// For 'decompose', only need mesh on master.
// Otherwise check for presence of the "faceLabels" file
bool haveMesh =
(
decompose
? Pstream::master()
: fileHandler().isFile
(
fileHandler().filePath
(
io.time().path()/io.instance()/meshSubDir/"faceLabels"
)
)
);
if (!haveMesh)
{
const bool oldParRun = Pstream::parRun(false);
// Create dummy mesh - on procs that don't already have a mesh
faMesh dummyMesh
(
pMesh,
labelList(),
IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE)
);
// Add patches
faPatchList patches(patchEntries.size());
label nPatches = 0;
forAll(patchEntries, patchi)
{
const entry& e = patchEntries[patchi];
const word type(e.dict().get<word>("type"));
const word& name = e.keyword();
if
(
type == processorFaPatch::typeName
)
{
// Stop at the first processor patch.
// - logic will not work with inter-mixed proc-patches anyhow
break;
}
else
{
dictionary patchDict(e.dict());
patchDict.set("edgeLabels", labelList());
patches.set
(
patchi,
faPatch::New
(
name,
patchDict,
nPatches++,
dummyMesh.boundary()
)
);
}
}
patches.resize(nPatches);
dummyMesh.addFaPatches(patches, false); // No parallel comms
// Bad hack, but the underlying polyMesh is NO_WRITE
// so it does not create the faMesh subDir for us...
Foam::mkDir(dummyMesh.boundary().path());
//Pout<< "Writing dummy mesh to " << dummyMesh.boundary().path()
// << endl;
dummyMesh.write();
Pstream::parRun(oldParRun); // Restore parallel state
}
// Read mesh
// ~~~~~~~~~
// Now all processors have a (possibly zero size) mesh so read in
// parallel
/// Pout<< "Reading area mesh from " << io.objectRelPath() << endl;
auto meshPtr = autoPtr<faMesh>::New(pMesh, false);
faMesh& mesh = *meshPtr;
// Sync patches
// ~~~~~~~~~~~~
if (!Pstream::master() && haveMesh)
{
// Check master names against mine
const faBoundaryMesh& patches = mesh.boundary();
forAll(patchEntries, patchi)
{
const entry& e = patchEntries[patchi];
const word type(e.dict().get<word>("type"));
const word& name = e.keyword();
if
(
type == processorFaPatch::typeName
)
{
break;
}
if (patchi >= patches.size())
{
FatalErrorInFunction
<< "Non-processor patches not synchronised."
<< endl
<< "Processor " << Pstream::myProcNo()
<< " has only " << patches.size()
<< " patches, master has "
<< patchi
<< exit(FatalError);
}
if
(
type != patches[patchi].type()
|| name != patches[patchi].name()
)
{
FatalErrorInFunction
<< "Non-processor patches not synchronised."
<< endl
<< "Master patch " << patchi
<< " name:" << type
<< " type:" << type << endl
<< "Processor " << Pstream::myProcNo()
<< " patch " << patchi
<< " has name:" << patches[patchi].name()
<< " type:" << patches[patchi].type()
<< exit(FatalError);
}
}
}
// Recreate basic geometry, globalMeshData etc.
mesh.init(false);
(void)mesh.globalData();
/// #if 0
/// faMeshTools::forceDemandDriven(mesh);
/// faMeshTools::unregisterMesh(mesh);
/// #endif
// Do some checks.
// Check if the boundary definition is unique
// and processor patches are correct
mesh.boundary().checkDefinition(verbose);
mesh.boundary().checkParallelSync(verbose);
return meshPtr;
}
// ************************************************************************* //

154
src/finiteArea/faMesh/faMeshTools/faMeshTools.H Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::faMeshTools
Description
A collection of tools for operating on an faMesh.
SourceFiles
faMeshTools.C
faMeshToolsProcAddr.C
faMeshToolsTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_faMeshTools_H
#define Foam_faMeshTools_H
#include "faMesh.H"
#include "areaFieldsFwd.H"
#include "edgeFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward Declarations
class mapDistributePolyMesh;
class writeHandler;
class polyMesh;
class IOobject;
/*---------------------------------------------------------------------------*\
Class faMeshTools Declaration
\*---------------------------------------------------------------------------*/
class faMeshTools
{
public:
//- Unregister the faMesh from its associated polyMesh
//- to prevent triggering on polyMesh changes etc.
static void unregisterMesh(const faMesh& mesh);
//- Force creation of everything that might vaguely be used by patches.
// This is fairly horrible.
static void forceDemandDriven(faMesh& mesh);
//- Read mesh or create dummy mesh (0 faces, >0 patches).
// Works in two modes according to masterOnlyReading:
// true : create a dummy mesh for all procs
// false: checks locally for mesh directories and only creates dummy mesh
// if not present
static autoPtr<faMesh> newMesh
(
const IOobject& io,
const polyMesh& pMesh,
const bool masterOnlyReading,
const bool verbose = false
);
// Read mesh if available. Otherwise create empty mesh with same non-proc
// patches as proc0 mesh. Requires:
// - all processors to have all patches (and in same order).
// - io.instance() set to facesInstance
static autoPtr<faMesh> loadOrCreateMesh
(
const IOobject& io,
const polyMesh& pMesh,
const bool decompose,
const bool verbose = false
);
//- Read decompose/reconstruct addressing
static mapDistributePolyMesh readProcAddressing
(
const faMesh& mesh,
const autoPtr<faMesh>& baseMeshPtr
);
//- Write decompose/reconstruct addressing
//
// \param mesh
// \param faDistMap
// \param decompose running in decompose vs reconstruct mode
// \param writeHandler file handler
// \param procMesh (optional) processor mesh in reconstruct mode
//
// \note Since the faMesh holds a reference to a polyMesh,
// in reconstruct mode it will refer to the base mesh, but
// we need a means to proc addressing into the processor locations.
// This is the purpose of the additional procMesh reference
static void writeProcAddressing
(
const faMesh& mesh,
const mapDistributePolyMesh& faDistMap,
const bool decompose,
autoPtr<fileOperation>&& writeHandler,
const faMesh* procMesh = nullptr
);
//- Flatten an edge field into linear addressing
// Optionally use primitive patch edge ordering
template<class Type>
static tmp<Field<Type>> flattenEdgeField
(
const GeometricField<Type, faePatchField, edgeMesh>& fld,
const bool primitiveOrdering = false
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "faMeshToolsTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

468
src/finiteArea/faMesh/faMeshTools/faMeshToolsProcAddr.C Normal file
View File

@ -0,0 +1,468 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "faMeshTools.H"
#include "BitOps.H"
#include "fileOperation.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "IOmapDistributePolyMesh.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Create a reconstruct map.
static mapDistributePolyMesh createReconstructMap
(
const faMesh& mesh,
const autoPtr<faMesh>& baseMeshPtr,
const labelUList& faceProcAddr,
const labelUList& edgeProcAddr,
const labelUList& pointProcAddr,
const labelUList& boundaryProcAddr
)
{
const label nOldPoints = mesh.nPoints();
const label nOldFaces = mesh.nFaces();
const label nOldEdges = mesh.nEdges();
///Pout<< "old sizes"
/// << " points:" << nOldPoints
/// << " faces:" << nOldFaces
/// << " edges:" << nOldEdges << nl;
const faBoundaryMesh& oldBndMesh = mesh.boundary();
labelList oldPatchStarts(oldBndMesh.patchStarts());
// Patches: purge -1 entries
labelList patchProcAddr
(
IndirectList<label>::subset_if
(
boundaryProcAddr,
labelRange::ge0()
)
);
labelListList faceSubMap(Pstream::nProcs());
faceSubMap[Pstream::masterNo()] = identity(nOldFaces);
labelListList edgeSubMap(Pstream::nProcs());
edgeSubMap[Pstream::masterNo()] = identity(nOldEdges);
labelListList pointSubMap(Pstream::nProcs());
pointSubMap[Pstream::masterNo()] = identity(nOldPoints);
labelListList patchSubMap(Pstream::nProcs());
patchSubMap[Pstream::masterNo()] = patchProcAddr;
// Gather addressing on the master
labelListList faceAddressing(Pstream::nProcs());
faceAddressing[Pstream::myProcNo()] = faceProcAddr;
Pstream::gatherList(faceAddressing);
labelListList edgeAddressing(Pstream::nProcs());
edgeAddressing[Pstream::myProcNo()] = edgeProcAddr;
Pstream::gatherList(edgeAddressing);
labelListList pointAddressing(Pstream::nProcs());
pointAddressing[Pstream::myProcNo()] = pointProcAddr;
Pstream::gatherList(pointAddressing);
labelListList patchAddressing(Pstream::nProcs());
patchAddressing[Pstream::myProcNo()] = patchProcAddr;
Pstream::gatherList(patchAddressing);
// NB: can only have a reconstruct on master!
if (Pstream::master() && baseMeshPtr && baseMeshPtr->nFaces())
{
const faMesh& baseMesh = *baseMeshPtr;
const label nNewPoints = baseMesh.nPoints();
const label nNewFaces = baseMesh.nFaces();
const label nNewEdges = baseMesh.nEdges();
const label nNewPatches = baseMesh.boundary().size();
/// Pout<< "new sizes"
/// << " points:" << nNewPoints
/// << " faces:" << nNewFaces
/// << " edges:" << nNewEdges << nl;
mapDistribute faFaceMap
(
nNewFaces,
std::move(faceSubMap),
std::move(faceAddressing),
false, // subHasFlip
false // constructHasFlip
);
mapDistribute faEdgeMap
(
nNewEdges,
std::move(edgeSubMap),
std::move(edgeAddressing),
false, // subHasFlip
false // constructHasFlip
);
mapDistribute faPointMap
(
nNewPoints,
std::move(pointSubMap),
std::move(pointAddressing)
);
mapDistribute faPatchMap
(
nNewPatches,
std::move(patchSubMap),
std::move(patchAddressing)
);
return mapDistributePolyMesh
(
// Mesh before changes
nOldPoints,
nOldEdges, // area: nOldEdges (volume: nOldFaces)
nOldFaces, // area: nOldFaces (volume: nOldCells)
std::move(oldPatchStarts),
labelList(), // oldPatchNMeshPoints [unused]
mapDistribute(std::move(faPointMap)),
mapDistribute(std::move(faEdgeMap)), // edgeMap (volume: faceMap)
mapDistribute(std::move(faFaceMap)), // faceMap (volume: cellMap)
mapDistribute(std::move(faPatchMap))
);
}
else
{
mapDistribute faFaceMap
(
0, // nNewFaces
std::move(faceSubMap),
labelListList(Pstream::nProcs()), // constructMap
false, // subHasFlip
false // constructHasFlip
);
mapDistribute faEdgeMap
(
0, // nNewEdges
std::move(edgeSubMap),
labelListList(Pstream::nProcs()), // constructMap
false, // subHasFlip
false // constructHasFlip
);
mapDistribute faPointMap
(
0, // nNewPoints
std::move(pointSubMap),
labelListList(Pstream::nProcs()) // constructMap
);
mapDistribute faPatchMap
(
0, // nNewPatches
std::move(patchSubMap),
labelListList(Pstream::nProcs()) // constructMap
);
return mapDistributePolyMesh
(
// Mesh before changes
nOldPoints,
nOldEdges, // area: nOldEdges (volume: nOldFaces)
nOldFaces, // area: nOldFaces (volume: nOldCells)
std::move(oldPatchStarts),
labelList(), // oldPatchNMeshPoints [unused]
mapDistribute(std::move(faPointMap)),
mapDistribute(std::move(faEdgeMap)), // edgeMap (volume: faceMap)
mapDistribute(std::move(faFaceMap)), // faceMap (volume: cellMap)
mapDistribute(std::move(faPatchMap))
);
}
}
} // End namespace Foam
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::mapDistributePolyMesh
Foam::faMeshTools::readProcAddressing
(
const faMesh& mesh,
const autoPtr<faMesh>& baseMeshPtr
)
{
IOobject ioAddr
(
"procAddressing",
mesh.facesInstance(),
faMesh::meshSubDir,
mesh.thisDb(),
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false // no register
);
//if (ioAddr.typeHeaderOk<labelIOList>(true))
//{
// Pout<< "Reading addressing from " << io.name() << " at "
// << mesh.facesInstance() << nl << endl;
// distMap.reset(new IOmapDistributePolyMesh(io));
//}
//else
{
Info<< "Reading (face|edge|face|point|boundary)ProcAddressing from "
<< mesh.facesInstance().c_str() << '/'
<< faMesh::meshSubDir << nl << endl;
ioAddr.rename("faceProcAddressing");
labelIOList faceProcAddressing(ioAddr, Zero);
ioAddr.rename("edgeProcAddressing");
labelIOList edgeProcAddressing(ioAddr, Zero);
ioAddr.rename("pointProcAddressing");
labelIOList pointProcAddressing(ioAddr, Zero);
ioAddr.rename("boundaryProcAddressing");
labelIOList boundaryProcAddressing(ioAddr, Zero);
if
(
mesh.nFaces() != faceProcAddressing.size()
|| mesh.nEdges() != edgeProcAddressing.size()
|| mesh.nPoints() != pointProcAddressing.size()
|| mesh.boundary().size() != boundaryProcAddressing.size()
)
{
FatalErrorInFunction
<< "Read addressing inconsistent with mesh sizes" << nl
<< "faces:" << mesh.nFaces()
<< " addressing:" << faceProcAddressing.objectRelPath()
<< " size:" << faceProcAddressing.size() << nl
<< "edges:" << mesh.nEdges()
<< " addressing:" << edgeProcAddressing.objectRelPath()
<< " size:" << edgeProcAddressing.size() << nl
<< "points:" << mesh.nPoints()
<< " addressing:" << pointProcAddressing.objectRelPath()
<< " size:" << pointProcAddressing.size()
<< "patches:" << mesh.boundary().size()
<< " addressing:" << boundaryProcAddressing.objectRelPath()
<< " size:" << boundaryProcAddressing.size()
<< exit(FatalError);
}
return createReconstructMap
(
mesh,
baseMeshPtr,
faceProcAddressing,
edgeProcAddressing,
pointProcAddressing,
boundaryProcAddressing
);
}
}
void Foam::faMeshTools::writeProcAddressing
(
const faMesh& mesh,
const mapDistributePolyMesh& map,
const bool decompose,
autoPtr<fileOperation>&& writeHandler,
const faMesh* procMesh
)
{
Info<< "Writing ("
<< (decompose ? "decompose" : "reconstruct")
<< ") procAddressing files to "
<< mesh.facesInstance().c_str() << '/'
<< faMesh::meshSubDir << endl;
IOobject ioAddr
(
"procAddressing",
mesh.facesInstance(),
faMesh::meshSubDir,
(procMesh && !decompose ? procMesh->thisDb() : mesh.thisDb()),
IOobject::NO_READ,
IOobject::NO_WRITE,
false // no register
);
// faceProcAddressing (faMesh)
ioAddr.rename("faceProcAddressing");
labelIOList faceMap(ioAddr, Zero);
// edgeProcAddressing (faMesh)
ioAddr.rename("edgeProcAddressing");
labelIOList edgeMap(ioAddr, Zero);
// pointProcAddressing (faMesh)
ioAddr.rename("pointProcAddressing");
labelIOList pointMap(ioAddr, Zero);
// boundaryProcAddressing (faMesh)
ioAddr.rename("boundaryProcAddressing");
labelIOList patchMap(ioAddr, Zero);
if (decompose)
{
// Decompose
// - forward map: [undecomposed] -> [decomposed]
// area:faces (volume:cells)
faceMap = identity(map.nOldCells());
map.cellMap().distribute(faceMap);
// area:edges (volume:faces)
edgeMap = identity(map.nOldFaces());
map.faceMap().distribute(edgeMap);
pointMap = identity(map.nOldPoints());
map.distributePointData(pointMap);
patchMap = identity(map.patchMap().constructSize());
map.patchMap().mapDistributeBase::distribute
(
Pstream::commsTypes::nonBlocking,
label(-1), // nullValue for new patches...
patchMap,
flipOp() // negate op
);
}
else // reconstruct
{
// Reconstruct
// - reverse map: [undecomposed] <- [decomposed]
// area:faces (volume:cells)
faceMap = identity(mesh.nFaces());
map.cellMap().reverseDistribute(map.nOldCells(), faceMap);
// area:edges (volume:faces)
edgeMap = identity(mesh.patch().nEdges());
map.faceMap().reverseDistribute(map.nOldFaces(), edgeMap);
pointMap = identity(mesh.nPoints());
map.pointMap().reverseDistribute(map.nOldPoints(), pointMap);
patchMap = identity(mesh.boundary().size());
map.patchMap().mapDistributeBase::reverseDistribute
(
Pstream::commsTypes::nonBlocking,
map.oldPatchSizes().size(),
label(-1), // nullValue for unmapped patches...
patchMap
);
}
autoPtr<fileOperation> defaultHandler;
if (writeHandler)
{
defaultHandler = fileHandler(std::move(writeHandler));
}
// If we want procAddressing, need to manually write it ourselves
// since it was not registered anywhere
IOmapDistributePolyMeshRef procAddrMap
(
IOobject
(
"procAddressing",
mesh.facesInstance(),
faMesh::meshSubDir,
mesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false // no register
),
map
);
if (decompose)
{
// Write into proc directories
procAddrMap.write();
}
else
{
// Reconstruct: "procAddressing" only meaningful for rank 0
// and written into base (serial) location (if at all).
if (Pstream::master())
{
const bool oldParRun = Pstream::parRun(false);
procAddrMap.write();
Pstream::parRun(oldParRun);
}
}
const bool faceOk = faceMap.write();
const bool edgeOk = edgeMap.write();
const bool pointOk = pointMap.write();
const bool patchOk = patchMap.write();
if (defaultHandler)
{
writeHandler = fileHandler(std::move(defaultHandler));
}
if (!edgeOk || !faceOk || !pointOk || !patchOk)
{
WarningInFunction
<< "Failed to write some of "
<< faceMap.objectRelPath() << ", "
<< edgeMap.objectRelPath() << ", "
<< pointMap.objectRelPath() << ", "
<< patchMap.objectRelPath() << endl;
}
}
// ************************************************************************* //

80
src/finiteArea/faMesh/faMeshTools/faMeshToolsTemplates.C Normal file
View File

@ -0,0 +1,80 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "edgeFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::faMeshTools::flattenEdgeField
(
const GeometricField<Type, faePatchField, edgeMesh>& fld,
const bool primitiveOrdering
)
{
const faMesh& mesh = fld.mesh();
auto tresult = tmp<Field<Type>>::New(mesh.nEdges(), Zero);
auto& result = tresult.ref();
// Internal field
result.slice(0, fld.size()) = fld;
if (primitiveOrdering)
{
// Boundary field in primitive patch order
forAll(fld.boundaryField(), patchi)
{
UIndirectList<Type>
(
result,
mesh.boundary()[patchi].edgeLabels()
) = fld.boundaryField()[patchi];
}
}
else
{
// Boundary field in sub-list (slice) order
label start = fld.size();
forAll(fld.boundaryField(), patchi)
{
const label len = mesh.boundary()[patchi].size();
result.slice(start, len) = fld.boundaryField()[patchi];
start += len;
}
}
return tresult;
}
// ************************************************************************* //

25
src/finiteArea/faMesh/faMeshUpdate.C
View File

@ -37,8 +37,6 @@ License
void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
{
DebugInFunction << "Updating mesh" << endl;
// if (!mpm.morphing())
// {
// // No topo change
@ -48,18 +46,14 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
// Create fa mesh mapper, using the old mesh
const faMeshMapper mapper(*this, mpm);
// Rebuild mesh
// Cast away const for interface reasons. HJ, 12/Aug/2011
faMesh& m = const_cast<faMesh&>(*this);
// ~~~~~~~~~~~~
// Clear existing mesh data
clearOut();
// Set new labels
m.faceLabels_ = mapper.areaMap().newFaceLabels();
faceLabels_ = mapper.areaMap().newFaceLabels();
const uindirectPrimitivePatch& bp = patch();
@ -121,12 +115,12 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
}
// Set new edges for all patches
forAll(m.boundary_, patchI)
forAll(boundary_, patchI)
{
m.boundary_[patchI].resetEdges(patchEdges[patchI]);
boundary_[patchI].resetEdges(patchEdges[patchI]);
}
m.setPrimitiveMeshData();
setPrimitiveMeshData();
// Create global mesh data
if (Pstream::parRun())
@ -135,10 +129,10 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
}
// Calculate topology for the patches (processor-processor comms etc.)
m.boundary_.updateMesh();
boundary_.updateMesh();
// Calculate the geometry for the patches (transformation tensors etc.)
m.boundary_.calcGeometry();
boundary_.calcGeometry();
// Map fields
@ -156,7 +150,7 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
{
// Map all the areaFields in the objectRegistry
// Map areaFields in the objectRegistry
MapGeometricFields<scalar, faPatchField, faMeshMapper, areaMesh>(mapper);
MapGeometricFields<vector, faPatchField, faMeshMapper, areaMesh>(mapper);
MapGeometricFields<sphericalTensor, faPatchField, faMeshMapper, areaMesh>
@ -165,7 +159,7 @@ void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
(mapper);
MapGeometricFields<tensor, faPatchField, faMeshMapper, areaMesh>(mapper);
// Map all the edgeFields in the objectRegistry
// Map edgeFields in the objectRegistry
MapGeometricFields<scalar, faePatchField, faMeshMapper, edgeMesh>(mapper);
MapGeometricFields<vector, faePatchField, faMeshMapper, edgeMesh>(mapper);
MapGeometricFields<sphericalTensor, faePatchField, faMeshMapper, edgeMesh>
@ -227,7 +221,6 @@ void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
}
}
}
}