mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: use fallback value if calculated Le() is degenerate (#2592)
- the Le vector is calculated from (edgeVec ^ edgeNorm) and should be oriented in direction (faceCentre -> edgeCentre). If, however, the edgeNorm value is bad for any reason, the cross-product falls apart and Le vector is calculated as a zero vector! For these cases, revert to using (faceCentre -> edgeCentre) as a better approximation than a zero vector. In the future, will very likely switch calculating the edge normals directly from the attached faces, instead of from the attached points as is currently done, which should improve robustness. ENH: expose fa:geometryOrder as a registered OptimisationSwitch ENN: reuse polyMesh data (eg, faceCentres) if possible in faMesh STYLE: add code lambdas and static functions to isolate logic
This commit is contained in:
Mark Olesen
@ -40,12 +40,24 @@ License
|
||||
#include "processorFaPatch.H"
|
||||
#include "wedgeFaPatch.H"
|
||||
#include "faPatchData.H"
|
||||
#include "registerSwitch.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(faMesh, 0);
|
||||
|
||||
int faMesh::geometryOrder_
|
||||
(
|
||||
debug::optimisationSwitch("fa:geometryOrder", 2)
|
||||
);
|
||||
registerOptSwitch
|
||||
(
|
||||
"fa:geometryOrder",
|
||||
int,
|
||||
faMesh::geometryOrder_
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
@ -53,8 +65,6 @@ const Foam::word Foam::faMesh::prefix("finite-area");
|
||||
|
||||
Foam::word Foam::faMesh::meshSubDir = "faMesh";
|
||||
|
||||
int Foam::faMesh::geometryOrder_ = 1; // 1: Standard treatment
|
||||
|
||||
const int Foam::faMesh::quadricsFit_ = 0; // Tuning (experimental)
|
||||
|
||||
|
||||
@ -228,7 +238,7 @@ void Foam::faMesh::clearGeomNotAreas() const
|
||||
deleteDemandDrivenData(patchStartsPtr_);
|
||||
deleteDemandDrivenData(LePtr_);
|
||||
deleteDemandDrivenData(magLePtr_);
|
||||
deleteDemandDrivenData(centresPtr_);
|
||||
deleteDemandDrivenData(faceCentresPtr_);
|
||||
deleteDemandDrivenData(edgeCentresPtr_);
|
||||
deleteDemandDrivenData(faceAreaNormalsPtr_);
|
||||
deleteDemandDrivenData(edgeAreaNormalsPtr_);
|
||||
@ -373,7 +383,7 @@ Foam::faMesh::faMesh
|
||||
patchStartsPtr_(nullptr),
|
||||
LePtr_(nullptr),
|
||||
magLePtr_(nullptr),
|
||||
centresPtr_(nullptr),
|
||||
faceCentresPtr_(nullptr),
|
||||
edgeCentresPtr_(nullptr),
|
||||
faceAreaNormalsPtr_(nullptr),
|
||||
edgeAreaNormalsPtr_(nullptr),
|
||||
@ -479,7 +489,7 @@ Foam::faMesh::faMesh
|
||||
patchStartsPtr_(nullptr),
|
||||
LePtr_(nullptr),
|
||||
magLePtr_(nullptr),
|
||||
centresPtr_(nullptr),
|
||||
faceCentresPtr_(nullptr),
|
||||
edgeCentresPtr_(nullptr),
|
||||
faceAreaNormalsPtr_(nullptr),
|
||||
edgeAreaNormalsPtr_(nullptr),
|
||||
@ -560,7 +570,7 @@ Foam::faMesh::faMesh
|
||||
patchStartsPtr_(nullptr),
|
||||
LePtr_(nullptr),
|
||||
magLePtr_(nullptr),
|
||||
centresPtr_(nullptr),
|
||||
faceCentresPtr_(nullptr),
|
||||
edgeCentresPtr_(nullptr),
|
||||
faceAreaNormalsPtr_(nullptr),
|
||||
edgeAreaNormalsPtr_(nullptr),
|
||||
@ -780,12 +790,12 @@ const Foam::edgeScalarField& Foam::faMesh::magLe() const
|
||||
|
||||
const Foam::areaVectorField& Foam::faMesh::areaCentres() const
|
||||
{
|
||||
if (!centresPtr_)
|
||||
if (!faceCentresPtr_)
|
||||
{
|
||||
calcAreaCentres();
|
||||
calcFaceCentres();
|
||||
}
|
||||
|
||||
return *centresPtr_;
|
||||
return *faceCentresPtr_;
|
||||
}
|
||||
|
||||
|
||||
@ -999,14 +1009,11 @@ bool Foam::faMesh::movePoints()
|
||||
|
||||
bool Foam::faMesh::correctPatchPointNormals(const label patchID) const
|
||||
{
|
||||
if ((patchID < 0) || (patchID >= boundary().size()))
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "patchID is not in the valid range"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
if (correctPatchPointNormalsPtr_)
|
||||
if
|
||||
(
|
||||
bool(correctPatchPointNormalsPtr_)
|
||||
&& patchID >= 0 && patchID < boundary().size()
|
||||
)
|
||||
{
|
||||
return (*correctPatchPointNormalsPtr_)[patchID];
|
||||
}
|
||||
|
||||
@ -287,7 +287,7 @@ class faMesh
|
||||
mutable edgeScalarField* magLePtr_;
|
||||
|
||||
//- Face centres
|
||||
mutable areaVectorField* centresPtr_;
|
||||
mutable areaVectorField* faceCentresPtr_;
|
||||
|
||||
//- Edge centres
|
||||
mutable edgeVectorField* edgeCentresPtr_;
|
||||
@ -328,9 +328,6 @@ class faMesh
|
||||
|
||||
// Static Private Data
|
||||
|
||||
//- Geometry treatment (0: primitive, 1: standard)
|
||||
static int geometryOrder_;
|
||||
|
||||
//- Quadrics fit for pointAreaNormals (experimental)
|
||||
static const int quadricsFit_;
|
||||
|
||||
@ -378,6 +375,11 @@ class faMesh
|
||||
//- are related to the areaMesh
|
||||
void calcWhichPatchFaces() const;
|
||||
|
||||
//- Calculate the edge normals (direct calculation)
|
||||
//- with flat boundary addressing
|
||||
// Possible communication
|
||||
tmp<vectorField> calcRawEdgeNormals(int calcType) const;
|
||||
|
||||
//- Calculate edge lengths
|
||||
// Triggers communication via calcEdgeAreaNormals
|
||||
void calcLe() const;
|
||||
@ -388,7 +390,7 @@ class faMesh
|
||||
|
||||
//- Calculate face centres
|
||||
// No communication
|
||||
void calcAreaCentres() const;
|
||||
void calcFaceCentres() const;
|
||||
|
||||
//- Calculate edge centres
|
||||
// No communication
|
||||
@ -505,6 +507,12 @@ public:
|
||||
typedef faBoundaryMesh BoundaryMesh;
|
||||
|
||||
|
||||
// Tuning switches
|
||||
|
||||
//- Geometry treatment
|
||||
static int geometryOrder_;
|
||||
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("faMesh");
|
||||
|
||||
@ -568,8 +576,8 @@ public:
|
||||
|
||||
// Static Functions
|
||||
|
||||
//- Return the current geometry treatment (0: primitive, 1: standard)
|
||||
// A zero level is with restricted neighbour information
|
||||
//- Return the current geometry treatment
|
||||
// A zero level or negative is with restricted neighbour information
|
||||
static int geometryOrder() noexcept
|
||||
{
|
||||
return geometryOrder_;
|
||||
@ -577,7 +585,7 @@ public:
|
||||
|
||||
//- Set the preferred geometry treatment
|
||||
// \return the previous value
|
||||
static int geometryOrder(const int order) noexcept
|
||||
static int geometryOrder(int order) noexcept
|
||||
{
|
||||
int old(geometryOrder_);
|
||||
geometryOrder_ = order;
|
||||
@ -770,13 +778,13 @@ public:
|
||||
//- Face centres of boundary halo neighbours
|
||||
const pointField& haloFaceCentres() const;
|
||||
|
||||
//- Face normals of boundary halo neighbours
|
||||
//- Face unit-normals of boundary halo neighbours
|
||||
const vectorField& haloFaceNormals() const;
|
||||
|
||||
//- Face centres of boundary halo neighbours for specified patch
|
||||
tmp<pointField> haloFaceCentres(const label patchi) const;
|
||||
|
||||
//- Face normals of boundary halo neighbours for specified patch
|
||||
//- Face unit-normals of boundary halo neighbours for specified patch
|
||||
tmp<vectorField> haloFaceNormals(const label patchi) const;
|
||||
|
||||
|
||||
@ -786,7 +794,7 @@ public:
|
||||
labelList faceCells() const;
|
||||
|
||||
|
||||
// Storage management
|
||||
// Storage Management
|
||||
|
||||
//- Remove all files from mesh instance
|
||||
void removeFiles(const fileName& instanceDir) const;
|
||||
@ -795,6 +803,46 @@ public:
|
||||
void removeFiles() const;
|
||||
|
||||
|
||||
//- Has face areas: S()
|
||||
bool hasFaceAreas() const noexcept { return bool(SPtr_); }
|
||||
|
||||
//- Has face centres: areaCentres()
|
||||
bool hasAreaCentres() const noexcept { return bool(faceCentresPtr_); }
|
||||
|
||||
//- Has edge centres: edgeCentres()
|
||||
bool hasAdgeCentres() const noexcept { return bool(edgeCentresPtr_); }
|
||||
|
||||
//- Has edge length vectors: Le()
|
||||
bool hasLe() const noexcept { return bool(LePtr_); }
|
||||
|
||||
//- Has edge length magnitudes: magLe()
|
||||
bool hasMagLe() const noexcept { return bool(magLePtr_); }
|
||||
|
||||
//- Has face area normals: faceAreaNormals()
|
||||
bool hasFaceAreaNormals() const noexcept
|
||||
{
|
||||
return bool(faceAreaNormalsPtr_);
|
||||
}
|
||||
|
||||
//- Has edge area normals: edgeAreaNormals()
|
||||
bool hasEdgeAreaNormals() const noexcept
|
||||
{
|
||||
return bool(edgeAreaNormalsPtr_);
|
||||
}
|
||||
|
||||
//- Has point area normals: pointAreaNormals()
|
||||
bool hasPointAreaNormals() const noexcept
|
||||
{
|
||||
return bool(pointAreaNormalsPtr_);
|
||||
}
|
||||
|
||||
//- Has face curvatures: faceCurvatures()
|
||||
bool hasFaceCurvatures() const noexcept
|
||||
{
|
||||
return bool(faceCurvaturesPtr_);
|
||||
}
|
||||
|
||||
|
||||
// Mesh motion and morphing
|
||||
|
||||
//- Is mesh moving
|
||||
|
||||
@ -104,34 +104,6 @@ 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).
|
||||
//
|
||||
@ -221,6 +193,43 @@ static inline vector calcEdgeNormalFromFace
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// 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
|
||||
);
|
||||
}
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
@ -269,7 +278,7 @@ void Foam::faMesh::calcPatchStarts() const
|
||||
if (patchStartsPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "patchStartsPtr_ already allocated"
|
||||
<< "patchStarts already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
@ -317,6 +326,126 @@ void Foam::faMesh::calcWhichPatchFaces() const
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
//- Calculate the 'Le' vector from faceCentre to edge centre
|
||||
// using the edge normal to correct for curvature
|
||||
//
|
||||
// Normalise and rescaled to the edge length
|
||||
static inline vector calcLeVector
|
||||
(
|
||||
const point& faceCentre,
|
||||
const linePointRef& edgeLine,
|
||||
const vector& edgeNormal // (unit or area normal)
|
||||
)
|
||||
{
|
||||
const vector centreToEdge(edgeLine.centre() - faceCentre);
|
||||
|
||||
vector leVector(edgeLine.vec() ^ edgeNormal);
|
||||
|
||||
scalar s(mag(leVector));
|
||||
|
||||
if (s < ROOTVSMALL)
|
||||
{
|
||||
// The calculated edgeNormal somehow degenerate and thus a
|
||||
// bad cross-product?
|
||||
// Revert to basic centre -> edge
|
||||
|
||||
leVector = centreToEdge;
|
||||
leVector.removeCollinear(edgeLine.unitVec());
|
||||
s = mag(leVector);
|
||||
|
||||
if (s < ROOTVSMALL)
|
||||
{
|
||||
// Unlikely that this should happen
|
||||
return Zero;
|
||||
}
|
||||
|
||||
leVector *= edgeLine.mag()/s;
|
||||
}
|
||||
else
|
||||
{
|
||||
// The additional orientation is probably unnecessary
|
||||
leVector *= edgeLine.mag()/s * sign(leVector & centreToEdge);
|
||||
}
|
||||
|
||||
return leVector;
|
||||
}
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::vectorField> Foam::faMesh::calcRawEdgeNormals(int order) const
|
||||
{
|
||||
// Return edge normals with flat boundary addressing
|
||||
auto tedgeNormals = tmp<vectorField>::New(nEdges_);
|
||||
auto& edgeNormals = tedgeNormals.ref();
|
||||
|
||||
// Need face centres
|
||||
const areaVectorField& fCentres = areaCentres();
|
||||
|
||||
// Also need local points
|
||||
const pointField& localPoints = points();
|
||||
|
||||
|
||||
{
|
||||
// Simple (primitive) edge normal calculations.
|
||||
// These are primarly designed to avoid any communication
|
||||
// but are thus necessarily inconsistent across processor boundaries!
|
||||
|
||||
WarningInFunction
|
||||
<< "Using geometryOrder < 1 : "
|
||||
"simplified edge area-normals, without processor connectivity"
|
||||
<< endl;
|
||||
|
||||
// Internal (edge normals) - contributions from owner/neighbour
|
||||
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
|
||||
{
|
||||
const linePointRef edgeLine(edges_[edgei].line(localPoints));
|
||||
|
||||
edgeNormals[edgei] =
|
||||
calcEdgeNormalFromFace
|
||||
(
|
||||
edgeLine,
|
||||
fCentres[edgeOwner()[edgei]],
|
||||
fCentres[edgeNeighbour()[edgei]]
|
||||
);
|
||||
}
|
||||
|
||||
// Boundary (edge normals) - like about but only has owner
|
||||
for (label edgei = nInternalEdges_; edgei < nEdges_; ++edgei)
|
||||
{
|
||||
const linePointRef edgeLine(edges_[edgei].line(localPoints));
|
||||
|
||||
edgeNormals[edgei] =
|
||||
calcEdgeNormalFromFace
|
||||
(
|
||||
edgeLine,
|
||||
fCentres[edgeOwner()[edgei]]
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Remove collinear components and normalise
|
||||
|
||||
forAll(edgeNormals, edgei)
|
||||
{
|
||||
const linePointRef edgeLine(edges_[edgei].line(localPoints));
|
||||
|
||||
edgeNormals[edgei].removeCollinear(edgeLine.unitVec());
|
||||
edgeNormals[edgei].normalise();
|
||||
}
|
||||
|
||||
return tedgeNormals;
|
||||
}
|
||||
|
||||
|
||||
void Foam::faMesh::calcLe() const
|
||||
{
|
||||
DebugInFunction
|
||||
@ -345,187 +474,97 @@ void Foam::faMesh::calcLe() const
|
||||
|
||||
edgeVectorField& Le = *LePtr_;
|
||||
|
||||
// Need face centres
|
||||
const areaVectorField& fCentres = areaCentres();
|
||||
|
||||
// Also need local points
|
||||
const pointField& localPoints = points();
|
||||
|
||||
|
||||
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!
|
||||
// The edge normals with flat boundary addressing
|
||||
// (which _may_ use communication)
|
||||
vectorField edgeNormals
|
||||
(
|
||||
calcRawEdgeNormals(faMesh::geometryOrder())
|
||||
);
|
||||
|
||||
// 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;
|
||||
// Calculate the Le vectors.
|
||||
// Can do inplace (overwrite with the edgeNormals)
|
||||
|
||||
UIndirectList<vector> fCentres(mesh().faceCentres(), faceLabels());
|
||||
|
||||
// Flat addressing
|
||||
vectorField edgeNormals(nEdges_);
|
||||
|
||||
// Internal (edge normals)
|
||||
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
|
||||
vectorField& leVectors = edgeNormals;
|
||||
forAll(leVectors, 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]]
|
||||
);
|
||||
leVectors[edgei] = calcLeVector
|
||||
(
|
||||
fCentres[edgeOwner()[edgei]],
|
||||
edges_[edgei].line(localPoints),
|
||||
edgeNormals[edgei]
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Now use these edge normals for calculating Le
|
||||
// Copy internal field
|
||||
Le.primitiveFieldRef() =
|
||||
vectorField::subList(leVectors, nInternalEdges_);
|
||||
|
||||
// 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)
|
||||
// Transcribe boundary field
|
||||
auto& bfld = Le.boundaryFieldRef();
|
||||
|
||||
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));
|
||||
}
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
bfld[patchi] = fap.patchRawSlice(leVectors);
|
||||
}
|
||||
|
||||
// Done
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
// Longer forms.
|
||||
// Using edgeAreaNormals, which uses pointAreaNormals (communication!)
|
||||
|
||||
const edgeVectorField& eCentres = edgeCentres();
|
||||
const areaVectorField& fCentres = areaCentres();
|
||||
const edgeVectorField& edgeNormals = edgeAreaNormals();
|
||||
|
||||
// Internal (edge vector)
|
||||
|
||||
else
|
||||
{
|
||||
vectorField& internalFld = Le.ref();
|
||||
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
|
||||
// Using edgeAreaNormals,
|
||||
// which _may_ use pointAreaNormals (communication!)
|
||||
|
||||
const edgeVectorField& edgeNormals = edgeAreaNormals();
|
||||
|
||||
// Internal (edge vector)
|
||||
{
|
||||
vector& leVector = internalFld[edgei];
|
||||
|
||||
vector edgeVec = edges_[edgei].vec(localPoints);
|
||||
const scalar magEdge(mag(edgeVec));
|
||||
|
||||
if (magEdge < ROOTVSMALL)
|
||||
vectorField& fld = Le.primitiveFieldRef();
|
||||
for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
|
||||
{
|
||||
// Too small
|
||||
leVector = Zero;
|
||||
continue;
|
||||
fld[edgei] = calcLeVector
|
||||
(
|
||||
fCentres[edgeOwner()[edgei]],
|
||||
edges_[edgei].line(localPoints),
|
||||
edgeNormals[edgei]
|
||||
);
|
||||
}
|
||||
|
||||
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)
|
||||
{
|
||||
const labelUList& bndEdgeFaces = boundary()[patchi].edgeFaces();
|
||||
|
||||
const edgeList::subList bndEdges =
|
||||
boundary()[patchi].patchSlice(edges_);
|
||||
|
||||
const vectorField& bndEdgeNormals =
|
||||
edgeNormals.boundaryField()[patchi];
|
||||
|
||||
vectorField& patchLe = Le.boundaryFieldRef()[patchi];
|
||||
const vectorField& patchECentres = eCentres.boundaryField()[patchi];
|
||||
|
||||
forAll(patchLe, bndEdgei)
|
||||
// Boundary (edge vector)
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
vector& leVector = patchLe[bndEdgei];
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
vectorField& pfld = Le.boundaryFieldRef()[patchi];
|
||||
|
||||
vector edgeVec = bndEdges[bndEdgei].vec(localPoints);
|
||||
const scalar magEdge(mag(edgeVec));
|
||||
const vectorField& bndEdgeNormals =
|
||||
edgeNormals.boundaryField()[patchi];
|
||||
|
||||
if (magEdge < ROOTVSMALL)
|
||||
label edgei = fap.start();
|
||||
|
||||
forAll(pfld, patchEdgei)
|
||||
{
|
||||
// Too small
|
||||
leVector = Zero;
|
||||
continue;
|
||||
pfld[patchEdgei] = calcLeVector
|
||||
(
|
||||
fCentres[edgeOwner()[edgei]],
|
||||
edges_[edgei].line(localPoints),
|
||||
bndEdgeNormals[patchEdgei]
|
||||
);
|
||||
|
||||
++edgei;
|
||||
}
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -539,7 +578,7 @@ void Foam::faMesh::calcMagLe() const
|
||||
if (magLePtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "magLePtr_ already allocated"
|
||||
<< "magLe() already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
@ -590,19 +629,19 @@ void Foam::faMesh::calcMagLe() const
|
||||
}
|
||||
|
||||
|
||||
void Foam::faMesh::calcAreaCentres() const
|
||||
void Foam::faMesh::calcFaceCentres() const
|
||||
{
|
||||
DebugInFunction
|
||||
<< "Calculating face centres" << endl;
|
||||
|
||||
if (centresPtr_)
|
||||
if (faceCentresPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "centresPtr_ already allocated"
|
||||
<< "areaCentres already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
centresPtr_ =
|
||||
faceCentresPtr_ =
|
||||
new areaVectorField
|
||||
(
|
||||
IOobject
|
||||
@ -616,16 +655,32 @@ void Foam::faMesh::calcAreaCentres() const
|
||||
dimLength
|
||||
);
|
||||
|
||||
areaVectorField& centres = *centresPtr_;
|
||||
areaVectorField& centres = *faceCentresPtr_;
|
||||
|
||||
// Need local points
|
||||
const pointField& localPoints = points();
|
||||
const faceList& localFaces = faces();
|
||||
|
||||
|
||||
// Internal (face centres)
|
||||
// Could also obtain from volume mesh faceCentres()
|
||||
forAll(localFaces, facei)
|
||||
{
|
||||
centres.ref()[facei] = localFaces[facei].centre(localPoints);
|
||||
if (mesh().hasFaceCentres())
|
||||
{
|
||||
// The volume mesh has faceCentres, can reuse them
|
||||
|
||||
centres.primitiveFieldRef()
|
||||
= UIndirectList<vector>(mesh().faceCentres(), faceLabels());
|
||||
}
|
||||
else
|
||||
{
|
||||
// Calculate manually
|
||||
auto iter = centres.primitiveFieldRef().begin();
|
||||
|
||||
for (const face& f : faces())
|
||||
{
|
||||
*iter = f.centre(localPoints);
|
||||
++iter;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Boundary (edge centres)
|
||||
@ -654,7 +709,7 @@ void Foam::faMesh::calcEdgeCentres() const
|
||||
if (edgeCentresPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "edgeCentresPtr_ already allocated"
|
||||
<< "edgeCentres already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
@ -676,6 +731,7 @@ void Foam::faMesh::calcEdgeCentres() const
|
||||
|
||||
const pointField& localPoints = points();
|
||||
|
||||
|
||||
// Internal (edge centres)
|
||||
{
|
||||
auto iter = centres.primitiveFieldRef().begin();
|
||||
@ -713,7 +769,7 @@ void Foam::faMesh::calcS() const
|
||||
if (SPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "SPtr_ already allocated"
|
||||
<< "S() already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
@ -730,15 +786,35 @@ void Foam::faMesh::calcS() const
|
||||
*this,
|
||||
dimArea
|
||||
);
|
||||
auto& S = *SPtr_;
|
||||
auto& areas = *SPtr_;
|
||||
|
||||
const pointField& localPoints = points();
|
||||
const faceList& localFaces = faces();
|
||||
|
||||
// Could also obtain from volume mesh faceAreas()
|
||||
forAll(S, facei)
|
||||
// No access to fvMesh::magSf(), only polyMesh::faceAreas()
|
||||
if (mesh().hasFaceAreas())
|
||||
{
|
||||
S[facei] = localFaces[facei].mag(localPoints);
|
||||
// The volume mesh has faceAreas, can reuse them
|
||||
UIndirectList<vector> meshFaceAreas(mesh().faceAreas(), faceLabels());
|
||||
|
||||
auto& fld = areas.field();
|
||||
|
||||
forAll(fld, facei)
|
||||
{
|
||||
fld[facei] = Foam::mag(meshFaceAreas[facei]);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Calculate manually
|
||||
|
||||
const pointField& localPoints = points();
|
||||
|
||||
auto iter = areas.field().begin();
|
||||
|
||||
for (const face& f : faces())
|
||||
{
|
||||
*iter = f.mag(localPoints);
|
||||
++iter;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -751,7 +827,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
|
||||
if (faceAreaNormalsPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "faceAreaNormalsPtr_ already allocated"
|
||||
<< "faceAreaNormals already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
@ -772,23 +848,43 @@ void Foam::faMesh::calcFaceAreaNormals() const
|
||||
areaVectorField& faceNormals = *faceAreaNormalsPtr_;
|
||||
|
||||
const pointField& localPoints = points();
|
||||
const faceList& localFaces = faces();
|
||||
|
||||
// Internal (faces)
|
||||
// Could also obtain from volume mesh Sf() + normalise
|
||||
vectorField& nInternal = faceNormals.ref();
|
||||
forAll(localFaces, faceI)
|
||||
// Internal
|
||||
{
|
||||
nInternal[faceI] = localFaces[faceI].unitNormal(localPoints);
|
||||
auto& fld = faceNormals.primitiveFieldRef();
|
||||
|
||||
if (mesh().hasFaceAreas())
|
||||
{
|
||||
// The volume mesh has faceAreas, can reuse them
|
||||
fld = UIndirectList<vector>(mesh().faceAreas(), faceLabels());
|
||||
}
|
||||
else
|
||||
{
|
||||
// Calculate manually
|
||||
|
||||
auto iter = fld.begin();
|
||||
|
||||
for (const face& f : faces())
|
||||
{
|
||||
*iter = f.areaNormal(localPoints);
|
||||
++iter;
|
||||
}
|
||||
}
|
||||
|
||||
// Make unit normals
|
||||
fld.normalise();
|
||||
}
|
||||
|
||||
|
||||
// Boundary - copy from edges
|
||||
|
||||
const auto& edgeNormalsBoundary = edgeAreaNormals().boundaryField();
|
||||
|
||||
forAll(boundary(), patchI)
|
||||
{
|
||||
faceNormals.boundaryFieldRef()[patchI] = edgeNormalsBoundary[patchI];
|
||||
const auto& edgeNormalsBoundary = edgeAreaNormals().boundaryField();
|
||||
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
faceNormals.boundaryFieldRef()[patchi]
|
||||
= edgeNormalsBoundary[patchi];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1131,9 +1227,9 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
const label nVerts(f.size());
|
||||
|
||||
point centrePoint(Zero);
|
||||
for (label i = 0; i < nVerts; ++i)
|
||||
for (const label fp : f)
|
||||
{
|
||||
centrePoint += points[f[i]];
|
||||
centrePoint += points[fp];
|
||||
}
|
||||
centrePoint /= nVerts;
|
||||
|
||||
@ -1153,9 +1249,8 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
}
|
||||
|
||||
|
||||
// Handle the boundary edges
|
||||
|
||||
bitSet nbrBoundaryAdjust(boundary().size(), true);
|
||||
// Boundary edge corrections
|
||||
bitSet nbrBoundaryAdjust;
|
||||
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
@ -1193,14 +1288,10 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
&result[wedgePatch.axisPoint()]
|
||||
);
|
||||
}
|
||||
|
||||
// Handled
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
else if (Pstream::parRun() && isA<processorFaPatch>(fap))
|
||||
{
|
||||
// Correct processor patch points
|
||||
|
||||
const auto& procPatch = refCast<const processorFaPatch>(fap);
|
||||
|
||||
const labelList& patchPoints = procPatch.pointLabels();
|
||||
@ -1244,25 +1335,16 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
result[patchPoints[pti]] +=
|
||||
patchPointNormals[nbrPatchPoints[pti]];
|
||||
}
|
||||
|
||||
// Handled
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
else if (fap.coupled())
|
||||
{
|
||||
// Coupled - no further action for neighbour side
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
// TBD:
|
||||
/// else if (isA<emptyFaPatch>(fap))
|
||||
/// {
|
||||
/// // Ignore this boundary
|
||||
/// nbrBoundaryAdjust.unset(patchi);
|
||||
/// }
|
||||
else if (!correctPatchPointNormals(patchi))
|
||||
else if (correctPatchPointNormals(patchi) && !fap.coupled())
|
||||
{
|
||||
// No corrections
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
// Neighbour correction requested
|
||||
nbrBoundaryAdjust.set(patchi);
|
||||
}
|
||||
}
|
||||
|
||||
@ -1289,7 +1371,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
gpNormals[addr[i]] += spNormals[i];
|
||||
}
|
||||
|
||||
Pstream::combineAllGather(gpNormals, plusEqOp<vectorField>());
|
||||
Pstream::combineReduce(gpNormals, plusEqOp<vectorField>());
|
||||
|
||||
// Extract local data
|
||||
forAll(addr, i)
|
||||
@ -1304,14 +1386,11 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
}
|
||||
|
||||
|
||||
if (returnReduce(nbrBoundaryAdjust.any(), orOp<bool>()))
|
||||
if (returnReduceOr(nbrBoundaryAdjust.any()))
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
PoutInFunction
|
||||
<< "Apply " << nbrBoundaryAdjust.count()
|
||||
<< " boundary neighbour corrections" << nl;
|
||||
}
|
||||
DebugInFunction
|
||||
<< "Apply " << nbrBoundaryAdjust.count()
|
||||
<< " boundary neighbour corrections" << nl;
|
||||
|
||||
// Apply boundary points correction
|
||||
|
||||
@ -1324,7 +1403,6 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
for (const label patchi : nbrBoundaryAdjust)
|
||||
{
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
const labelList& edgeLabels = fap.edgeLabels();
|
||||
|
||||
if (fap.ngbPolyPatchIndex() < 0)
|
||||
{
|
||||
@ -1334,7 +1412,8 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
for (const label edgei : edgeLabels)
|
||||
// NB: haloFaceNormals uses primitivePatch edge indexing
|
||||
for (const label edgei : fap.edgeLabels())
|
||||
{
|
||||
const edge& e = patch().edges()[edgei];
|
||||
|
||||
@ -1359,9 +1438,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
forAllConstIters(fpNormals, iter)
|
||||
{
|
||||
const label pointi = iter.key();
|
||||
vector fpnorm = normalised(iter.val());
|
||||
|
||||
result[pointi].removeCollinear(fpnorm);
|
||||
result[pointi].removeCollinear(normalised(iter.val()));
|
||||
}
|
||||
}
|
||||
|
||||
@ -1957,39 +2034,48 @@ Foam::tmp<Foam::edgeScalarField> Foam::faMesh::edgeLengthCorrection() const
|
||||
|
||||
const vectorField& pointNormals = pointAreaNormals();
|
||||
|
||||
forAll(correction.internalField(), edgeI)
|
||||
const auto angleCorrection =
|
||||
[](const vector& a, const vector& b) -> scalar
|
||||
{
|
||||
return Foam::cos(0.5*Foam::asin(Foam::mag(a ^ b)));
|
||||
};
|
||||
|
||||
|
||||
// Internal
|
||||
{
|
||||
scalar sinAlpha = mag
|
||||
(
|
||||
pointNormals[edges()[edgeI].start()]^
|
||||
pointNormals[edges()[edgeI].end()]
|
||||
);
|
||||
auto& fld = correction.primitiveFieldRef();
|
||||
|
||||
scalar alpha = asin(sinAlpha);
|
||||
|
||||
correction.ref()[edgeI] = cos(0.5*alpha);
|
||||
forAll(fld, edgei)
|
||||
{
|
||||
fld[edgei] = angleCorrection
|
||||
(
|
||||
pointNormals[edges_[edgei].start()],
|
||||
pointNormals[edges_[edgei].end()]
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
forAll(boundary(), patchI)
|
||||
// Boundary
|
||||
{
|
||||
const edgeList::subList patchEdges
|
||||
(
|
||||
boundary()[patchI].patchSlice(edges())
|
||||
);
|
||||
auto& bfld = correction.boundaryFieldRef();
|
||||
|
||||
forAll(patchEdges, edgeI)
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
scalar sinAlpha =
|
||||
mag
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
scalarField& pfld = bfld[patchi];
|
||||
|
||||
label edgei = fap.start();
|
||||
|
||||
forAll(pfld, patchEdgei)
|
||||
{
|
||||
pfld[patchEdgei] = angleCorrection
|
||||
(
|
||||
pointNormals[patchEdges[edgeI].start()]
|
||||
^ pointNormals[patchEdges[edgeI].end()]
|
||||
pointNormals[edges_[edgei].start()],
|
||||
pointNormals[edges_[edgei].end()]
|
||||
);
|
||||
|
||||
scalar alpha = asin(sinAlpha);
|
||||
|
||||
correction.boundaryFieldRef()[patchI][edgeI] = cos(0.5*alpha);
|
||||
++edgei;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user