diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C index 67f18252ff..67f687a577 100644 --- a/src/finiteArea/faMesh/faMesh.C +++ b/src/finiteArea/faMesh/faMesh.C @@ -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]; } diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index 7a4a17c8e3..bd391977ac 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -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 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 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 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 diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C index 7fd15f4db4..2fbfd468d9 100644 --- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C +++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C @@ -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::faMesh::calcRawEdgeNormals(int order) const +{ + // Return edge normals with flat boundary addressing + auto tedgeNormals = tmp::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 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(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 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(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(fap)) { // Correct processor patch points - const auto& procPatch = refCast(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(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()); + Pstream::combineReduce(gpNormals, plusEqOp()); // Extract local data forAll(addr, i) @@ -1304,14 +1386,11 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const } - if (returnReduce(nbrBoundaryAdjust.any(), orOp())) + 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::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; + } } }