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