From 891ac808def2ecfbdf054d3c48f334010cd97273 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 6 Aug 2025 09:07:04 +0200 Subject: [PATCH] FIX: robuster blockMesh, face::centre with SP mode (#3411) - appears to hit single precision overflow with clang-15 in face::center(), cellModel::center() and blockMesh createPoints(). The blockMesh might be particularly sensitive, since the points are frequently defined in millimeters (scaled later), which results in large intermediate summations. Similar to primitiveMesh checks, use double precision for these calculations. ENH: support vector += and -= from compatible types - eg, doubleVector += floatVector is now supported. This streamlines some coding for mixed precision. - To avoid lots of boilerplate, do not yet attempt to support general operations such as `operator+(doubleVector, floatVector)` until they become necessary. --- applications/test/vector/Test-vector.cxx | 56 +++++++++ .../meshes/meshShapes/cellModel/cellModel.C | 70 +++++++---- .../meshes/meshShapes/cellModel/cellModel.H | 6 +- src/OpenFOAM/meshes/meshShapes/face/face.C | 103 +++++++---------- src/OpenFOAM/meshes/meshShapes/face/face.H | 3 + src/OpenFOAM/meshes/meshShapes/face/faceI.H | 19 +++ .../meshes/meshShapes/face/faceTemplates.C | 18 +-- .../primitiveMeshCheck/primitiveMeshCheck.C | 10 +- .../primitiveMeshCheck/primitiveMeshTools.C | 89 ++++++++------ .../primitiveShapes/tetrahedron/tetrahedron.H | 20 ++-- .../tetrahedron/tetrahedronI.H | 15 --- .../primitiveShapes/triangle/triangle.H | 16 +-- .../primitiveShapes/triangle/triangleI.H | 15 --- src/OpenFOAM/primitives/Tensor/Tensor.H | 3 +- src/OpenFOAM/primitives/Tensor/TensorI.H | 4 +- src/OpenFOAM/primitives/Tensor2D/Tensor2D.H | 3 +- src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H | 4 +- src/OpenFOAM/primitives/Vector/Vector.H | 43 +++++-- src/OpenFOAM/primitives/Vector/VectorI.H | 10 +- src/OpenFOAM/primitives/Vector2D/Vector2D.H | 39 ++++++- src/OpenFOAM/primitives/Vector2D/Vector2DI.H | 10 +- .../polyMeshGeometry/polyMeshGeometry.C | 50 +++++--- .../polyTopoChange/addPatchCellLayer.C | 4 +- .../highAspectRatioFvGeometryScheme.C | 8 +- .../stabilised/stabilisedFvGeometryScheme.C | 55 ++++++--- src/mesh/blockMesh/blocks/block/blockCreate.C | 109 ++++++++++-------- src/meshTools/indexedOctree/treeDataFace.C | 2 +- .../indexedOctree/treeDataPrimitivePatch.C | 2 +- .../primitiveMeshGeometry.C | 4 +- 29 files changed, 492 insertions(+), 298 deletions(-) diff --git a/applications/test/vector/Test-vector.cxx b/applications/test/vector/Test-vector.cxx index cc2ad57d36..621479828a 100644 --- a/applications/test/vector/Test-vector.cxx +++ b/applications/test/vector/Test-vector.cxx @@ -33,6 +33,7 @@ Description #include "vectorField.H" #include "boolVector.H" +#include "complexVector.H" #include "labelVector.H" #include "IOstreams.H" #include "FixedList.H" @@ -164,6 +165,30 @@ void testTranscribe(Type& input) } +// Test of adding to vectors, possible of dissimilar types +template +void testAdding(const Target& a, const Source& b) +{ + typedef typename Target::cmptType cmptType1; + typedef typename Source::cmptType cmptType2; + + Info<< " " + << pTraits::typeName << " += " + << pTraits::typeName << " : "; + + if constexpr (std::is_convertible_v) + { + Target res(a); + res += b; + Info<< a << " += " << b << " == " << res << nl; + } + else + { + Info<< "unsupported" << nl; + } +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: @@ -173,6 +198,37 @@ int main(int argc, char *argv[]) Info<<"normalised: " << vector::uniform(VSMALL).normalise() << nl; Info<<"normalised: " << vector::uniform(ROOTVSMALL).normalise() << nl; + { + complexVector cvec(complex(1), complex(0), complex(0)); + doubleVector dvec(1, 0, 0); + floatVector fvec(1, 0, 0); + labelVector lvec(1, 0, 0); + + Info<< nl << "additions:" << nl; + + testAdding(cvec, cvec); + testAdding(cvec, dvec); + testAdding(cvec, fvec); + testAdding(cvec, lvec); + + testAdding(dvec, cvec); + testAdding(dvec, dvec); + testAdding(dvec, fvec); + testAdding(dvec, lvec); + + testAdding(fvec, cvec); + testAdding(fvec, dvec); + testAdding(fvec, fvec); + testAdding(fvec, lvec); + + testAdding(lvec, cvec); + testAdding(lvec, dvec); + testAdding(lvec, fvec); + testAdding(lvec, lvec); + + Info<< nl; + } + { vector vec1(0.5, 0.5, 0.5); vector vec2(0.5, 0.51, -0.5); diff --git a/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C b/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C index da60911c43..e7e63db824 100644 --- a/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C +++ b/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017,2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -29,34 +29,61 @@ License #include "cellModel.H" #include "pyramid.H" +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Average of points +// Note: use double precision to avoid overflows when summing +static inline Foam::doubleVector pointsAverage +( + const UList& points, + const labelUList& pointLabels +) +{ + doubleVector avg(Zero); + + if (const auto n = pointLabels.size(); n) + { + for (const auto pointi : pointLabels) + { + avg += points[pointi]; + } + avg /= n; + } + + return avg; +} + +} // End namespace Foam + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // Foam::vector Foam::cellModel::centre ( - const labelList& pointLabels, + const labelUList& pointLabels, const UList& points ) const { - // Estimate cell centre by averaging the cell points - vector cEst = Zero; - for (const label pointi : pointLabels) - { - cEst += points[pointi]; - } - cEst /= scalar(pointLabels.size()); + // Note: use double precision to avoid overflows when summing + + // Estimated centre by averaging the cell points + const point centrePoint(pointsAverage(points, pointLabels)); // Calculate the centre by breaking the cell into pyramids and // volume-weighted averaging their centres - scalar sumV = 0; - vector sumVc = Zero; + doubleScalar sumV(0); + doubleVector sumVc(Zero); forAll(faces_, facei) { const Foam::face f(pointLabels, faces_[facei]); - const scalar pyrVol = pyramidPointFaceRef(f, cEst).mag(points); + const scalar pyrVol = pyramidPointFaceRef(f, centrePoint).mag(points); if (pyrVol > SMALL) { @@ -67,7 +94,7 @@ Foam::vector Foam::cellModel::centre } sumV -= pyrVol; - sumVc -= pyrVol * pyramidPointFaceRef(f, cEst).centre(points); + sumVc -= pyrVol * pyramidPointFaceRef(f, centrePoint).centre(points); } return sumVc/(sumV + VSMALL); @@ -76,17 +103,14 @@ Foam::vector Foam::cellModel::centre Foam::scalar Foam::cellModel::mag ( - const labelList& pointLabels, + const labelUList& pointLabels, const UList& points ) const { - // Estimate cell centre by averaging the cell points - vector cEst = Zero; - for (const label pointi : pointLabels) - { - cEst += points[pointi]; - } - cEst /= scalar(pointLabels.size()); + // Note: use double precision to avoid overflows when summing + + // Estimated centre by averaging the cell points + const point centrePoint(pointsAverage(points, pointLabels)); // Calculate the magnitude by summing the -mags of the pyramids @@ -94,13 +118,13 @@ Foam::scalar Foam::cellModel::mag // and a pyramid is constructed from an inward pointing face // and the base centre-apex vector - scalar sumV = 0; + scalar sumV(0); forAll(faces_, facei) { const Foam::face f(pointLabels, faces_[facei]); - const scalar pyrVol = pyramidPointFaceRef(f, cEst).mag(points); + const scalar pyrVol = pyramidPointFaceRef(f, centrePoint).mag(points); if (pyrVol > SMALL) { diff --git a/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.H b/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.H index 118589f3a5..8614d39423 100644 --- a/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.H +++ b/src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation - Copyright (C) 2017-2020 OpenCFD Ltd. + Copyright (C) 2017-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -221,14 +221,14 @@ public: //- Centroid of the cell vector centre ( - const labelList& pointLabels, + const labelUList& pointLabels, const UList& points ) const; //- Cell volume scalar mag ( - const labelList& pointLabels, + const labelUList& pointLabels, const UList& points ) const; diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C index 7b238c6c8b..914225687d 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.C +++ b/src/OpenFOAM/meshes/meshShapes/face/face.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2021-2023 OpenCFD Ltd. + Copyright (C) 2021-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -92,7 +92,7 @@ Foam::label Foam::face::mostConcaveAngle angle = constant::mathematical::pi - edgeAngle; } - if (angle > maxAngle) + if (maxAngle < angle) { maxAngle = angle; index = i; @@ -527,34 +527,34 @@ Foam::point Foam::face::centre(const UList& points) const points[operator[](2)] ); } + // else if (nPoints == 0) + // { + // // Should not happen + // return Zero; + // } - point centrePoint = Zero; - for (label pI=0; pIaverage(points)); + + doubleScalar sumA(0); + doubleVector sumAc(Zero); + + for (label pI = 0; pI < nPoints; ++pI) { - centrePoint += points[operator[](pI)]; - } - centrePoint /= nPoints; - - scalar sumA = 0; - vector sumAc = Zero; - - for (label pI=0; pI& points) const } -Foam::vector Foam::face::areaNormal(const UList& p) const +Foam::vector Foam::face::areaNormal(const UList& points) const { const label nPoints = size(); // Calculate the area normal by summing the face triangle area normals. // Changed to deal with small concavity by using a central decomposition - // If the face is a triangle, do a direct calculation to avoid round-off - // error-related problems - + // If the face is a triangle, do a direct calculation if (nPoints == 3) { return triPointRef::areaNormal ( - p[operator[](0)], - p[operator[](1)], - p[operator[](2)] + points[operator[](0)], + points[operator[](1)], + points[operator[](2)] ); } + // else if (nPoints == 0) + // { + // // Should not happen + // return Zero; + // } - label pI; + // Estimated centre by averaging the face points + const point centrePoint(this->average(points)); - point centrePoint = Zero; - for (pI = 0; pI < nPoints; ++pI) + vector n(Zero); + + for (label pI = 0; pI < nPoints; ++pI) { - centrePoint += p[operator[](pI)]; - } - centrePoint /= nPoints; - - vector n = Zero; - - point nextPoint = centrePoint; - - for (pI = 0; pI < nPoints; ++pI) - { - if (pI < nPoints - 1) - { - nextPoint = p[operator[](pI + 1)]; - } - else - { - nextPoint = p[operator[](0)]; - } + const label nextPti(pI == nPoints-1 ? 0 : pI+1); + const auto& thisPoint = points[operator[](pI)]; + const auto& nextPoint = points[operator[](nextPti)]; // Note: for best accuracy, centre point always comes last - // - n += triPointRef::areaNormal - ( - p[operator[](pI)], - nextPoint, - centrePoint - ); + + n += triPointRef::areaNormal(thisPoint, nextPoint, centrePoint); } return n; @@ -689,10 +674,10 @@ Foam::scalar Foam::face::sweptVol // summing their swept volumes. // Changed to deal with small concavity by using a central decomposition - point centreOldPoint = centre(oldPoints); - point centreNewPoint = centre(newPoints); + const point centreOldPoint = centre(oldPoints); + const point centreNewPoint = centre(newPoints); - label nPoints = size(); + const label nPoints = size(); for (label pi=0; pi& points) const; + //- Average of face points: a quick estimate of the face centre. + inline point average(const UList& points) const; + //- Calculate average value at centroid of face template Type average diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceI.H b/src/OpenFOAM/meshes/meshShapes/face/faceI.H index 5f2e3698a6..73392203f8 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceI.H +++ b/src/OpenFOAM/meshes/meshShapes/face/faceI.H @@ -136,6 +136,25 @@ Foam::face::box(const UList& pts) const } +inline Foam::point +Foam::face::average(const UList& points) const +{ + // Note: use double precision to avoid overflows when summing + doubleVector avg(Zero); + + if (const auto n = size(); n) + { + for (const auto pointi : *this) + { + avg += points[pointi]; + } + avg /= n; + } + + return avg; +} + + inline Foam::label Foam::face::nEdges() const noexcept { // for a closed polygon a number of edges is the same as number of points diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C index fdf688e11b..f8ef1a4793 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C @@ -58,8 +58,10 @@ Type Foam::face::average // Calculate the average by breaking the face into triangles and // area-weighted averaging their averages + const label nPoints = size(); + // If the face is a triangle, do a direct calculation - if (size() == 3) + if (nPoints == 3) { return (1.0/3.0) @@ -70,18 +72,14 @@ Type Foam::face::average ); } - label nPoints = size(); + const point centrePoint(this->average(meshPoints)); - point centrePoint = Zero; Type cf = Zero; - for (label pI=0; pI sumClosed(Zero); - solveScalar sumMagClosedBoundary = 0; + solveVector sumClosed(Zero); + solveScalar sumMagClosedBoundary(0); for (label facei = nInternalFaces(); facei < areas.size(); facei++) { if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei]) { - sumClosed += Vector(areas[facei]); + sumClosed += areas[facei]; sumMagClosedBoundary += mag(areas[facei]); } } - reduce(sumClosed, sumOp>()); + reduce(sumClosed, sumOp()); reduce(sumMagClosedBoundary, sumOp()); - Vector openness = sumClosed/(sumMagClosedBoundary + VSMALL); + solveVector openness = sumClosed/(sumMagClosedBoundary + VSMALL); if (cmptMax(cmptMag(openness)) > closedThreshold_) { diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index ad5e0aa63c..751e81c529 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2012-2016 OpenFOAM Foundation - Copyright (C) 2017-2022 OpenCFD Ltd. + Copyright (C) 2017-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -33,6 +33,36 @@ License #include "tetrahedron.H" #include "PrecisionAdaptor.H" +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Average of points +// Note: use double precision to avoid overflows when summing +static inline Foam::doubleVector pointsAverage +( + const UList& points, + const labelUList& pointLabels +) +{ + doubleVector avg(Zero); + + if (const auto n = pointLabels.size(); n) + { + for (const auto pointi : pointLabels) + { + avg += points[pointi]; + } + avg /= n; + } + + return avg; +} + +} // End namespace Foam + + // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // void Foam::primitiveMeshTools::updateFaceCentresAndAreas @@ -48,31 +78,25 @@ void Foam::primitiveMeshTools::updateFaceCentresAndAreas for (const label facei : faceIDs) { - const labelList& f = fs[facei]; + const auto& f = fs[facei]; const label nPoints = f.size(); // If the face is a triangle, do a direct calculation for efficiency // and to avoid round-off error-related problems if (nPoints == 3) { - fCtrs[facei] = triPointRef::centre(p[f[0]], p[f[1]], p[f[2]]); - fAreas[facei] = triPointRef::areaNormal(p[f[0]], p[f[1]], p[f[2]]); + const triPointRef tri(p, f[0], f[1], f[2]); + fCtrs[facei] = tri.centre(); + fAreas[facei] = tri.areaNormal(); } else { - typedef Vector solveVector; + solveVector sumN(Zero); + solveScalar sumA(0); + solveVector sumAc(Zero); - solveVector sumN = Zero; - solveScalar sumA = 0.0; - solveVector sumAc = Zero; - - solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; ++pi) - { - fCentre += solveVector(p[f[pi]]); - } - - fCentre /= nPoints; + // Estimated centre by averaging the face points + const solveVector fCentre(pointsAverage(p, f)); for (label pi = 0; pi < nPoints; ++pi) { @@ -116,12 +140,10 @@ void Foam::primitiveMeshTools::updateCellCentresAndVols scalarField& cellVols_s ) { - typedef Vector solveVector; - PrecisionAdaptor tcellCtrs(cellCtrs_s, false); PrecisionAdaptor tcellVols(cellVols_s, false); - Field& cellCtrs = tcellCtrs.ref(); - Field& cellVols = tcellVols.ref(); + auto& cellCtrs = tcellCtrs.ref(); + auto& cellVols = tcellVols.ref(); // Use current cell centres as estimates for the new cell centres @@ -222,21 +244,18 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas // and to avoid round-off error-related problems if (nPoints == 3) { - fCtrs[facei] = triPointRef::centre(p[f[0]], p[f[1]], p[f[2]]); - fAreas[facei] = triPointRef::areaNormal(p[f[0]], p[f[1]], p[f[2]]); + const triPointRef tri(p, f[0], f[1], f[2]); + fCtrs[facei] = tri.centre(); + fAreas[facei] = tri.areaNormal(); } else { - solveVector sumN = Zero; - solveScalar sumA = Zero; - solveVector sumAc = Zero; + solveVector sumN(Zero); + solveScalar sumA(0); + solveVector sumAc(Zero); - solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; ++pi) - { - fCentre += solveVector(p[f[pi]]); - } - fCentre /= nPoints; + // Estimated centre by averaging the face points + const solveVector fCentre(pointsAverage(p, f)); for (label pi = 0; pi < nPoints; ++pi) { @@ -292,8 +311,8 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols { PrecisionAdaptor tcellCtrs(cellCtrs_s, false); PrecisionAdaptor tcellVols(cellVols_s, false); - Field& cellCtrs = tcellCtrs.ref(); - Field& cellVols = tcellVols.ref(); + auto& cellCtrs = tcellCtrs.ref(); + auto& cellVols = tcellVols.ref(); // Clear the fields for accumulation cellCtrs = Zero; @@ -310,13 +329,13 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols forAll(own, facei) { - cEst[own[facei]] += solveVector(fCtrs[facei]); + cEst[own[facei]] += fCtrs[facei]; ++nCellFaces[own[facei]]; } forAll(nei, facei) { - cEst[nei[facei]] += solveVector(fCtrs[facei]); + cEst[nei[facei]] += fCtrs[facei]; ++nCellFaces[nei[facei]]; } diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index 3a31dd5d21..41f7c87551 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H @@ -182,34 +182,34 @@ public: typedef FixedList tetIntersectionList; - // Classes for use in sliceWithPlane. What to do with decomposition - // of tet. + // Classes for use in sliceWithPlane. + // What to do with decomposition of tetrahedron. //- Dummy - class dummyOp + struct dummyOp { - public: - inline void operator()(const tetPoints&); + void operator()(const tetPoints&) const noexcept + {} }; //- Sum resulting volumes - class sumVolOp + struct sumVolOp { - public: scalar vol_; - inline sumVolOp(); + constexpr sumVolOp() noexcept : vol_(0) {} + + scalar total() const noexcept { return vol_; } inline void operator()(const tetPoints&); }; //- Store resulting tets - class storeOp + struct storeOp { tetIntersectionList& tets_; label& nTets_; - public: inline storeOp(tetIntersectionList&, label&); inline void operator()(const tetPoints&); diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 951dc07193..4fcb0d136b 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -580,21 +580,6 @@ bool Foam::tetrahedron::inside(const point& p) const } -template -inline void Foam::tetrahedron::dummyOp::operator() -( - const tetPoints& -) -{} - - -template -inline Foam::tetrahedron::sumVolOp::sumVolOp() -: - vol_(0.0) -{} - - template inline void Foam::tetrahedron::sumVolOp::operator() ( diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 0f48056c36..134a848e4c 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -211,27 +211,27 @@ public: // What to do with decomposition of triangle. //- Dummy - class dummyOp + struct dummyOp { - public: - inline void operator()(const triPoints&); + void operator()(const triPoints&) const noexcept + {} }; //- Sum resulting areas - class sumAreaOp + struct sumAreaOp { - public: scalar area_; - inline sumAreaOp(); + constexpr sumAreaOp() noexcept : area_(0) {} + + scalar total() const noexcept { return area_; } inline void operator()(const triPoints&); }; //- Store resulting tris - class storeOp + struct storeOp { - public: triIntersectionList& tris_; label& nTris_; diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 1e5a7bd25d..e4a4d2a743 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -1068,21 +1068,6 @@ inline int Foam::triangle::sign } -template -inline void Foam::triangle::dummyOp::operator() -( - const triPoints& -) -{} - - -template -inline Foam::triangle::sumAreaOp::sumAreaOp() -: - area_(0.0) -{} - - template inline void Foam::triangle::sumAreaOp::operator() ( diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H index 601e004416..6a0ec41706 100644 --- a/src/OpenFOAM/primitives/Tensor/Tensor.H +++ b/src/OpenFOAM/primitives/Tensor/Tensor.H @@ -101,7 +101,7 @@ public: // Constructors //- Construct initialized to zero - inline Tensor(const Foam::zero); + inline Tensor(Foam::zero); //- Construct given MatrixSpace of the same rank template @@ -338,6 +338,7 @@ public: //- Deprecated(2018-12) Return vector for given row (0,1) // \deprecated(2018-12) use row() method + FOAM_DEPRECATED_FOR(2018-12, "row()") Vector vectorComponent(const direction cmpt) const { return row(cmpt); diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H index 7c46f93639..b08a1feb35 100644 --- a/src/OpenFOAM/primitives/Tensor/TensorI.H +++ b/src/OpenFOAM/primitives/Tensor/TensorI.H @@ -33,9 +33,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -inline Foam::Tensor::Tensor(const Foam::zero) +inline Foam::Tensor::Tensor(Foam::zero) : - Tensor::msType(Zero) + Tensor::msType(Foam::zero{}) {} diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H index 69f8efd091..4f62657ce4 100644 --- a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H +++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H @@ -96,7 +96,7 @@ public: // Constructors //- Construct initialized to zero - inline Tensor2D(const Foam::zero); + inline Tensor2D(Foam::zero); //- Construct given VectorSpace inline Tensor2D(const VectorSpace, Cmpt, 4>& vs); @@ -241,6 +241,7 @@ public: //- Deprecated(2018-12) Return vector for given row (0,1) // \deprecated(2018-12) use row() method + FOAM_DEPRECATED_FOR(2018-12, "row()") Vector2D vectorComponent(const direction cmpt) const { return row(cmpt); diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H index 1acbb876a8..e478159833 100644 --- a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H +++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H @@ -29,9 +29,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -inline Foam::Tensor2D::Tensor2D(const Foam::zero) +inline Foam::Tensor2D::Tensor2D(Foam::zero) : - Tensor2D::vsType(Zero) + Tensor2D::vsType(Foam::zero{}) {} diff --git a/src/OpenFOAM/primitives/Vector/Vector.H b/src/OpenFOAM/primitives/Vector/Vector.H index 053865a2a8..c1b5b185eb 100644 --- a/src/OpenFOAM/primitives/Vector/Vector.H +++ b/src/OpenFOAM/primitives/Vector/Vector.H @@ -97,7 +97,7 @@ public: // Constructors //- Construct initialized to zero - inline Vector(const Foam::zero); + inline Vector(Foam::zero); //- Copy construct from VectorSpace of the same rank template @@ -115,22 +115,22 @@ public: // Component Access //- Access to the vector x component - const Cmpt& x() const noexcept { return this->v_[X]; } + const Cmpt& x() const noexcept { return this->v_[components::X]; } //- Access to the vector y component - const Cmpt& y() const noexcept { return this->v_[Y]; } + const Cmpt& y() const noexcept { return this->v_[components::Y]; } //- Access to the vector z component - const Cmpt& z() const noexcept { return this->v_[Z]; } + const Cmpt& z() const noexcept { return this->v_[components::Z]; } //- Access to the vector x component - Cmpt& x() noexcept { return this->v_[X]; } + Cmpt& x() noexcept { return this->v_[components::X]; } //- Access to the vector y component - Cmpt& y() noexcept { return this->v_[Y]; } + Cmpt& y() noexcept { return this->v_[components::Y]; } //- Access to the vector z component - Cmpt& z() noexcept { return this->v_[Z]; } + Cmpt& z() noexcept { return this->v_[components::Z]; } // Vector Operations @@ -193,6 +193,35 @@ public: const Vector& a, const Vector& b ); + + + // Member Operators + + //- Inherit VectorSpace += operations + using Vector::vsType::operator+=; + + //- Inherit VectorSpace -= operations + using Vector::vsType::operator-=; + + //- Add compatible vector to this + template + std::enable_if_t, void> + inline operator+=(const Vector& b) + { + this->x() += b.x(); + this->y() += b.y(); + this->z() += b.z(); + } + + //- Subtract compatible vector from this + template + std::enable_if_t, void> + inline operator-=(const Vector& b) + { + this->x() -= b.x(); + this->y() -= b.y(); + this->z() -= b.z(); + } }; diff --git a/src/OpenFOAM/primitives/Vector/VectorI.H b/src/OpenFOAM/primitives/Vector/VectorI.H index 2b5449f680..8a20481754 100644 --- a/src/OpenFOAM/primitives/Vector/VectorI.H +++ b/src/OpenFOAM/primitives/Vector/VectorI.H @@ -29,9 +29,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -inline Foam::Vector::Vector(const Foam::zero) +inline Foam::Vector::Vector(Foam::zero) : - Vector::vsType(Zero) + Vector::vsType(Foam::zero{}) {} @@ -54,9 +54,9 @@ inline Foam::Vector::Vector const Cmpt& vz ) { - this->v_[X] = vx; - this->v_[Y] = vy; - this->v_[Z] = vz; + this->x() = vx; + this->y() = vy; + this->z() = vz; } diff --git a/src/OpenFOAM/primitives/Vector2D/Vector2D.H b/src/OpenFOAM/primitives/Vector2D/Vector2D.H index 01ad98d2f6..1ca1c35a47 100644 --- a/src/OpenFOAM/primitives/Vector2D/Vector2D.H +++ b/src/OpenFOAM/primitives/Vector2D/Vector2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2022 OpenCFD Ltd. + Copyright (C) 2018-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -90,7 +90,7 @@ public: // Constructors //- Construct initialized to zero - inline Vector2D(const Foam::zero); + inline Vector2D(Foam::zero); //- Copy construct from VectorSpace of the same rank inline Vector2D(const VectorSpace, Cmpt, 2>& vs); @@ -107,16 +107,16 @@ public: // Component Access //- Access to the vector x component - const Cmpt& x() const noexcept { return this->v_[X]; } + const Cmpt& x() const noexcept { return this->v_[components::X]; } //- Access to the vector y component - const Cmpt& y() const noexcept { return this->v_[Y]; } + const Cmpt& y() const noexcept { return this->v_[components::Y]; } //- Access to the vector x component - Cmpt& x() noexcept { return this->v_[X]; } + Cmpt& x() noexcept { return this->v_[components::X]; } //- Access to the vector y component - Cmpt& y() noexcept { return this->v_[Y]; } + Cmpt& y() noexcept { return this->v_[components::Y]; } // Vector Operations @@ -170,6 +170,33 @@ public: const Vector2D& a, const Vector2D& b ); + + + // Member Operators + + //- Inherit VectorSpace += operations + using Vector2D::vsType::operator+=; + + //- Inherit VectorSpace -= operations + using Vector2D::vsType::operator-=; + + //- Add compatible 2D vector to this + template + std::enable_if_t, void> + inline operator+=(const Vector2D& b) + { + this->x() += b.x(); + this->y() += b.y(); + } + + //- Subtract compatible 2D vector from this + template + std::enable_if_t, void> + inline operator-=(const Vector2D& b) + { + this->x() -= b.x(); + this->y() -= b.y(); + } }; diff --git a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H b/src/OpenFOAM/primitives/Vector2D/Vector2DI.H index bdd677986f..cf3cb2282d 100644 --- a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H +++ b/src/OpenFOAM/primitives/Vector2D/Vector2DI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2022 OpenCFD Ltd. + Copyright (C) 2018-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -29,9 +29,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -inline Foam::Vector2D::Vector2D(const Foam::zero) +inline Foam::Vector2D::Vector2D(Foam::zero) : - Vector2D::vsType(Zero) + Vector2D::vsType(Foam::zero{}) {} @@ -48,8 +48,8 @@ inline Foam::Vector2D::Vector2D template inline Foam::Vector2D::Vector2D(const Cmpt& vx, const Cmpt& vy) { - this->v_[X] = vx; - this->v_[Y] = vy; + this->x() = vx; + this->y() = vy; } diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index 966f278e07..ecc6b41b57 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -41,6 +41,35 @@ namespace Foam defineTypeNameAndDebug(polyMeshGeometry, 0); } +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Average of points +// Note: use double precision to avoid overflows when summing +static inline Foam::doubleVector pointsAverage +( + const UList& points, + const labelUList& pointLabels +) +{ + doubleVector avg(Zero); + + if (const auto n = pointLabels.size(); n) + { + for (const auto pointi : pointLabels) + { + avg += points[pointi]; + } + avg /= n; + } + + return avg; +} + +} // End namespace Foam + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -61,23 +90,18 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas // and to avoid round-off error-related problems if (nPoints == 3) { - faceCentres_[facei] = - triPointRef::centre(p[f[0]], p[f[1]], p[f[2]]); - faceAreas_[facei] = - triPointRef::areaNormal(p[f[0]], p[f[1]], p[f[2]]); + const triPointRef tri(p, f[0], f[1], f[2]); + faceCentres_[facei] = tri.centre(); + faceAreas_[facei] = tri.areaNormal(); } else { - solveVector sumN = Zero; - solveScalar sumA = Zero; - solveVector sumAc = Zero; + solveVector sumN(Zero); + solveScalar sumA(0); + solveVector sumAc(Zero); - solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; ++pi) - { - fCentre += solveVector(p[f[pi]]); - } - fCentre /= nPoints; + // Estimated centre by averaging the face points + const solveVector fCentre(pointsAverage(p, f)); for (label pi = 0; pi < nPoints; ++pi) { diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C index 848e854ad1..6dc3232522 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C @@ -3001,7 +3001,7 @@ void Foam::addPatchCellLayer::setRefinement // // forAll(own, facei) // { - // cEst[own[facei]] += solveVector(fCtrs[facei]); + // cEst[own[facei]] += fCtrs[facei]; // ++nCellFaces[own[facei]]; // } // @@ -3011,7 +3011,7 @@ void Foam::addPatchCellLayer::setRefinement // { // continue; // } - // cEst[nei[facei]] += solveVector(fCtrs[facei]); + // cEst[nei[facei]] += fCtrs[facei]; // ++nCellFaces[nei[facei]]; // } // diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C index 6e60a2eac3..72c51d5103 100644 --- a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C +++ b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C @@ -246,7 +246,7 @@ void Foam::highAspectRatioFvGeometryScheme::makeAverageCentres forAll(fs, facei) { - const labelList& f = fs[facei]; + const auto& f = fs[facei]; const label nPoints = f.size(); if (nPoints == 3) @@ -255,8 +255,8 @@ void Foam::highAspectRatioFvGeometryScheme::makeAverageCentres } else { - solveScalar sumA = 0.0; - solveVector sumAc = Zero; + solveScalar sumA(0); + solveVector sumAc(Zero); for (label pi = 0; pi < nPoints; pi++) { @@ -282,7 +282,7 @@ void Foam::highAspectRatioFvGeometryScheme::makeAverageCentres sumAc = Zero; for (label pi = 0; pi < nPoints; pi++) { - sumAc += static_cast(p[f[pi]]); + sumAc += p[f[pi]]; } faceCentres[facei] = sumAc/nPoints; } diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C index 186994a50d..c182793577 100644 --- a/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C +++ b/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C @@ -47,6 +47,37 @@ namespace Foam } +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Average of points +// Note: use double precision to avoid overflows when summing +static inline Foam::doubleVector pointsAverage +( + const UList& points, + const labelUList& pointLabels +) +{ + doubleVector avg(Zero); + + if (const auto n = pointLabels.size(); n) + { + for (const auto pointi : pointLabels) + { + avg += points[pointi]; + } + avg /= n; + } + + return avg; +} + +} // End namespace Foam + + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::stabilisedFvGeometryScheme::makeFaceCentresAndAreas @@ -61,32 +92,28 @@ void Foam::stabilisedFvGeometryScheme::makeFaceCentresAndAreas forAll(fs, facei) { - const labelList& f = fs[facei]; - label nPoints = f.size(); + const auto& f = fs[facei]; + const label nPoints = f.size(); // If the face is a triangle, do a direct calculation for efficiency // and to avoid round-off error-related problems if (nPoints == 3) { - fCtrs[facei] = triPointRef::centre(p[f[0]], p[f[1]], p[f[2]]); - fAreas[facei] = triPointRef::areaNormal(p[f[0]], p[f[1]], p[f[2]]); + const triPointRef tri(p, f[0], f[1], f[2]); + fCtrs[facei] = tri.centre(); + fAreas[facei] = tri.areaNormal(); } // For more complex faces, decompose into triangles else { - // Compute an estimate of the centre as the average of the points - solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; pi++) - { - fCentre += solveVector(p[f[pi]]); - } - fCentre /= nPoints; + // Estimated centre by averaging the face points + const solveVector fCentre(pointsAverage(p, f)); // Compute the face area normal and unit normal by summing up the // normals of the triangles formed by connecting each edge to the // point average. - solveVector sumA = Zero; + solveVector sumA(Zero); for (label pi = 0; pi < nPoints; pi++) { const label nextPi(pi == nPoints-1 ? 0 : pi+1); @@ -104,8 +131,8 @@ void Foam::stabilisedFvGeometryScheme::makeFaceCentresAndAreas // the triangle area projected in the direction of the face normal // as the weight, *not* the triangle area magnitude. Only the // former makes the calculation independent of the initial estimate. - solveScalar sumAn = 0.0; - solveVector sumAnc = Zero; + solveScalar sumAn(0); + solveVector sumAnc(Zero); for (label pi = 0; pi < nPoints; pi++) { const label nextPi(pi == nPoints-1 ? 0 : pi+1); diff --git a/src/mesh/blockMesh/blocks/block/blockCreate.C b/src/mesh/blockMesh/blocks/block/blockCreate.C index 175b933118..da28dc3eb6 100644 --- a/src/mesh/blockMesh/blocks/block/blockCreate.C +++ b/src/mesh/blockMesh/blocks/block/blockCreate.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2021 OpenCFD Ltd. + Copyright (C) 2019-2025 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -96,39 +96,39 @@ void Foam::block::createPoints() scalar wx2 = (1 - w1)*w4*(1 - w11) + w1*w5*(1 - w10); scalar wx3 = (1 - w2)*w7*w11 + w2*w6*w10; scalar wx4 = (1 - w3)*(1 - w7)*w8 + w3*(1 - w6)*w9; - - const scalar sumWx = wx1 + wx2 + wx3 + wx4; - wx1 /= sumWx; - wx2 /= sumWx; - wx3 /= sumWx; - wx4 /= sumWx; - + { + const scalar sum(wx1 + wx2 + wx3 + wx4); + wx1 /= sum; + wx2 /= sum; + wx3 /= sum; + wx4 /= sum; + } // y-direction scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11); scalar wy2 = (1 - w5)*w0*(1 - w9) + w5*w1*(1 - w10); scalar wy3 = (1 - w6)*w3*w9 + w6*w2*w10; scalar wy4 = (1 - w7)*(1 - w3)*w8 + w7*(1 - w2)*w11; - - const scalar sumWy = wy1 + wy2 + wy3 + wy4; - wy1 /= sumWy; - wy2 /= sumWy; - wy3 /= sumWy; - wy4 /= sumWy; - + { + const scalar sum(wy1 + wy2 + wy3 + wy4); + wy1 /= sum; + wy2 /= sum; + wy3 /= sum; + wy4 /= sum; + } // z-direction scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7); scalar wz2 = (1 - w9)*w0*(1 - w5) + w9*w3*(1 - w6); scalar wz3 = (1 - w10)*w1*w5 + w10*w2*w6; scalar wz4 = (1 - w11)*(1 - w1)*w4 + w11*(1 - w2)*w7; - - const scalar sumWz = wz1 + wz2 + wz3 + wz4; - wz1 /= sumWz; - wz2 /= sumWz; - wz3 /= sumWz; - wz4 /= sumWz; - + { + const scalar sum(wz1 + wz2 + wz3 + wz4); + wz1 /= sum; + wz2 /= sum; + wz3 /= sum; + wz4 /= sum; + } // Points on straight edges const vector edgex1 = p000 + (p100 - p000)*w0; @@ -147,40 +147,47 @@ void Foam::block::createPoints() const vector edgez4 = p010 + (p011 - p010)*w11; // Add the contributions - points_[vijk] = - ( - wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4 - + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4 - + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4 - )/3; + // Note: use double precision to avoid overflows when summing + doubleVector sumVector(Zero); + + sumVector += wx1*edgex1; + sumVector += wx2*edgex2; + sumVector += wx3*edgex3; + sumVector += wx4*edgex4; + + sumVector += wy1*edgey1; + sumVector += wy2*edgey2; + sumVector += wy3*edgey3; + sumVector += wy4*edgey4; + + sumVector += wz1*edgez1; + sumVector += wz2*edgez2; + sumVector += wz3*edgez3; + sumVector += wz4*edgez4; + sumVector /= 3; // Apply curved-edge correction if block has curved edges if (nCurvedEdges) { // Calculate the correction vectors - const vector corx1 = wx1*(p[0][i] - edgex1); - const vector corx2 = wx2*(p[1][i] - edgex2); - const vector corx3 = wx3*(p[2][i] - edgex3); - const vector corx4 = wx4*(p[3][i] - edgex4); + sumVector += wx1*(p[0][i] - edgex1); // corx1 + sumVector += wx2*(p[1][i] - edgex2); // corx2 + sumVector += wx3*(p[2][i] - edgex3); // corx3 + sumVector += wx4*(p[3][i] - edgex4); // corx4 - const vector cory1 = wy1*(p[4][j] - edgey1); - const vector cory2 = wy2*(p[5][j] - edgey2); - const vector cory3 = wy3*(p[6][j] - edgey3); - const vector cory4 = wy4*(p[7][j] - edgey4); + sumVector += wy1*(p[4][j] - edgey1); // cory1 + sumVector += wy2*(p[5][j] - edgey2); // cory2 + sumVector += wy3*(p[6][j] - edgey3); // cory3 + sumVector += wy4*(p[7][j] - edgey4); // cory4 - const vector corz1 = wz1*(p[8][k] - edgez1); - const vector corz2 = wz2*(p[9][k] - edgez2); - const vector corz3 = wz3*(p[10][k] - edgez3); - const vector corz4 = wz4*(p[11][k] - edgez4); - - points_[vijk] += - ( - corx1 + corx2 + corx3 + corx4 - + cory1 + cory2 + cory3 + cory4 - + corz1 + corz2 + corz3 + corz4 - ); + sumVector += wz1*(p[8][k] - edgez1); // corz1 + sumVector += wz2*(p[9][k] - edgez2); // corz2 + sumVector += wz3*(p[10][k] - edgez3); // corz3 + sumVector += wz4*(p[11][k] - edgez4); // corz4 } + + points_[vijk] = sumVector; } } } @@ -326,9 +333,11 @@ void Foam::block::createPoints() + w11*(1 - w2)*w7 ); - const scalar sumWf = wf4 + wf5; - wf4 /= sumWf; - wf5 /= sumWf; + { + const scalar sum(wf4 + wf5); + wf4 /= sum; + wf5 /= sum; + } points_[vijk] += ( diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C index 7be1099a86..377c387a62 100644 --- a/src/meshTools/indexedOctree/treeDataFace.C +++ b/src/meshTools/indexedOctree/treeDataFace.C @@ -557,7 +557,7 @@ bool Foam::treeDataFace::overlaps if (f.size() == 3) { - const triPointRef tri(points[f[0]], points[f[1]], points[f[2]]); + const triPointRef tri(points, f[0], f[1], f[2]); return searchBox.intersects(tri); } diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C index 7cafbe5dc3..c10d2bd734 100644 --- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C @@ -417,7 +417,7 @@ bool Foam::treeDataPrimitivePatch::overlaps if (f.size() == 3) { - const triPointRef tri(points[f[0]], points[f[1]], points[f[2]]); + const triPointRef tri(points, f[0], f[1], f[2]); return searchBox.intersects(tri); } diff --git a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C index cb6d5bd663..08c9d85d2e 100644 --- a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C +++ b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C @@ -51,8 +51,8 @@ void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas for (label facei : changedFaces) { - const labelList& f = fs[facei]; - label nPoints = f.size(); + const auto& f = fs[facei]; + const label nPoints = f.size(); // If the face is a triangle, do a direct calculation for efficiency // and to avoid round-off error-related problems