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

FIX: robuster blockMesh, face::centre with SP mode (#3411)

Browse Source 
- 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.
This commit is contained in:
Mark Olesen
2025-08-06 09:07:04 +02:00
parent fe7006acd7
commit 891ac808de
29 changed files with 492 additions and 298 deletions
Download Patch File Download Diff File Expand all files Collapse all files

56
applications/test/vector/Test-vector.cxx
View File

@ -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<class Target, class Source>
void testAdding(const Target& a, const Source& b)
{
typedef typename Target::cmptType cmptType1;
typedef typename Source::cmptType cmptType2;
Info<< " "
<< pTraits<Target>::typeName << " += "
<< pTraits<Source>::typeName << " : ";
if constexpr (std::is_convertible_v<cmptType2, cmptType1>)
{
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);

70
src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C
View File

@ -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<point>& 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<point>& 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<point>& 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)
{

6
src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.H
View File

@ -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<point>& points
) const;
//- Cell volume
scalar mag
(
const labelList& pointLabels,
const labelUList& pointLabels,
const UList<point>& points
) const;

103
src/OpenFOAM/meshes/meshShapes/face/face.C
View File

@ -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<point>& points) const
points[operator[](2)]
);
}
// else if (nPoints == 0)
// {
// // Should not happen
// return Zero;
// }
point centrePoint = Zero;
for (label pI=0; pI<nPoints; ++pI)
// Note: use double precision to avoid overflows when summing
// Estimated centre by averaging the face points
const point centrePoint(this->average(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<nPoints; ++pI)
{
const point& nextPoint = points[operator[]((pI + 1) % nPoints)];
const label nextPti(pI == nPoints-1 ? 0 : pI+1);
const auto& thisPoint = points[operator[](pI)];
const auto& nextPoint = points[operator[](nextPti)];
// Calculate 3*triangle centre
const vector ttc
(
points[operator[](pI)]
+ nextPoint
+ centrePoint
);
const vector ttc(thisPoint + nextPoint + centrePoint);
// Calculate 2*triangle area
const scalar ta = Foam::mag
(
(points[operator[](pI)] - centrePoint)
(thisPoint - centrePoint)
^ (nextPoint - centrePoint)
);
@ -573,58 +573,43 @@ Foam::point Foam::face::centre(const UList<point>& points) const
}
Foam::vector Foam::face::areaNormal(const UList<point>& p) const
Foam::vector Foam::face::areaNormal(const UList<point>& 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<nPoints-1; ++pi)
{

3
src/OpenFOAM/meshes/meshShapes/face/face.H
View File

@ -182,6 +182,9 @@ public:
//- Centre point of face
point centre(const UList<point>& points) const;
//- Average of face points: a quick estimate of the face centre.
inline point average(const UList<point>& points) const;
//- Calculate average value at centroid of face
template<class Type>
Type average

19
src/OpenFOAM/meshes/meshShapes/face/faceI.H
View File

@ -136,6 +136,25 @@ Foam::face::box(const UList<point>& pts) const
}
inline Foam::point
Foam::face::average(const UList<point>& 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

18
src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C
View File

@ -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<nPoints; pI++)
{
centrePoint += meshPoints[operator[](pI)];
cf += fld[operator[](pI)];
}
centrePoint /= nPoints;
cf /= nPoints;
scalar sumA = 0;
@ -89,19 +87,21 @@ Type Foam::face::average
for (label pI=0; pI<nPoints; pI++)
{
const label nextPti(pI == nPoints-1 ? 0 : pI+1);
// Calculate 3*triangle centre field value
Type ttcf =
(
fld[operator[](pI)]
+ fld[operator[]((pI + 1) % nPoints)]
+ fld[operator[](nextPti)]
+ cf
);
// Calculate 2*triangle area
scalar ta = Foam::mag
const scalar ta = Foam::mag
(
(meshPoints[operator[](pI)] - centrePoint)
^ (meshPoints[operator[]((pI + 1) % nPoints)] - centrePoint)
^ (meshPoints[operator[](nextPti)] - centrePoint)
);
sumA += ta;

10
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
View File

@ -57,22 +57,22 @@ bool Foam::primitiveMesh::checkClosedBoundary
// Loop through all boundary faces and sum up the face area vectors.
// For a closed boundary, this should be zero in all vector components
Vector<solveScalar> 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<solveScalar>(areas[facei]);
sumClosed += areas[facei];
sumMagClosedBoundary += mag(areas[facei]);
}
}
reduce(sumClosed, sumOp<Vector<solveScalar>>());
reduce(sumClosed, sumOp<solveVector>());
reduce(sumMagClosedBoundary, sumOp<solveScalar>());
Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
solveVector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
if (cmptMax(cmptMag(openness)) > closedThreshold_)
{

89
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
View File

@ -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<point>& 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<solveScalar> 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<solveScalar> solveVector;
PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s, false);
PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s, false);
Field<solveVector>& cellCtrs = tcellCtrs.ref();
Field<solveScalar>& 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<solveVector, vector> tcellCtrs(cellCtrs_s, false);
PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s, false);
Field<solveVector>& cellCtrs = tcellCtrs.ref();
Field<solveScalar>& 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]];
}

20
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
View File

@ -182,34 +182,34 @@ public:
typedef FixedList<tetPoints, 200> 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&);

15
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
View File

@ -580,21 +580,6 @@ bool Foam::tetrahedron<Point, PointRef>::inside(const point& p) const
}
template<class Point, class PointRef>
inline void Foam::tetrahedron<Point, PointRef>::dummyOp::operator()
(
const tetPoints&
)
{}
template<class Point, class PointRef>
inline Foam::tetrahedron<Point, PointRef>::sumVolOp::sumVolOp()
:
vol_(0.0)
{}
template<class Point, class PointRef>
inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
(

16
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
View File

@ -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_;

15
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View File

@ -1068,21 +1068,6 @@ inline int Foam::triangle<Point, PointRef>::sign
}
template<class Point, class PointRef>
inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
(
const triPoints&
)
{}
template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::sumAreaOp::sumAreaOp()
:
area_(0.0)
{}
template<class Point, class PointRef>
inline void Foam::triangle<Point, PointRef>::sumAreaOp::operator()
(

3
src/OpenFOAM/primitives/Tensor/Tensor.H
View File

@ -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<class Cmpt2>
@ -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<Cmpt> vectorComponent(const direction cmpt) const
{
return row(cmpt);

4
src/OpenFOAM/primitives/Tensor/TensorI.H
View File

@ -33,9 +33,9 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::Tensor<Cmpt>::Tensor(const Foam::zero)
inline Foam::Tensor<Cmpt>::Tensor(Foam::zero)
:
Tensor::msType(Zero)
Tensor::msType(Foam::zero{})
{}

3
src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
View File

@ -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<Tensor2D<Cmpt>, 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<Cmpt> vectorComponent(const direction cmpt) const
{
return row(cmpt);

4
src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
View File

@ -29,9 +29,9 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::Tensor2D<Cmpt>::Tensor2D(const Foam::zero)
inline Foam::Tensor2D<Cmpt>::Tensor2D(Foam::zero)
:
Tensor2D::vsType(Zero)
Tensor2D::vsType(Foam::zero{})
{}

43
src/OpenFOAM/primitives/Vector/Vector.H
View File

@ -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<class Cmpt2>
@ -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<Cmpt>& a,
const Vector<Cmpt>& b
);
// Member Operators
//- Inherit VectorSpace += operations
using Vector::vsType::operator+=;
//- Inherit VectorSpace -= operations
using Vector::vsType::operator-=;
//- Add compatible vector to this
template<class Cmpt2>
std::enable_if_t<std::is_convertible_v<Cmpt2, Cmpt>, void>
inline operator+=(const Vector<Cmpt2>& b)
{
this->x() += b.x();
this->y() += b.y();
this->z() += b.z();
}
//- Subtract compatible vector from this
template<class Cmpt2>
std::enable_if_t<std::is_convertible_v<Cmpt2, Cmpt>, void>
inline operator-=(const Vector<Cmpt2>& b)
{
this->x() -= b.x();
this->y() -= b.y();
this->z() -= b.z();
}
};

10
src/OpenFOAM/primitives/Vector/VectorI.H
View File

@ -29,9 +29,9 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline Foam::Vector<Cmpt>::Vector(const Foam::zero)
inline Foam::Vector<Cmpt>::Vector(Foam::zero)
:
Vector::vsType(Zero)
Vector::vsType(Foam::zero{})
{}
@ -54,9 +54,9 @@ inline Foam::Vector<Cmpt>::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;
}

39
src/OpenFOAM/primitives/Vector2D/Vector2D.H
View File

@ -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<Vector2D<Cmpt>, 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<Cmpt>& a,
const Vector2D<Cmpt>& b
);
// Member Operators
//- Inherit VectorSpace += operations
using Vector2D::vsType::operator+=;
//- Inherit VectorSpace -= operations
using Vector2D::vsType::operator-=;
//- Add compatible 2D vector to this
template<class Cmpt2>
std::enable_if_t<std::is_convertible_v<Cmpt2, Cmpt>, void>
inline operator+=(const Vector2D<Cmpt2>& b)
{
this->x() += b.x();
this->y() += b.y();
}
//- Subtract compatible 2D vector from this
template<class Cmpt2>
std::enable_if_t<std::is_convertible_v<Cmpt2, Cmpt>, void>
inline operator-=(const Vector2D<Cmpt2>& b)
{
this->x() -= b.x();
this->y() -= b.y();
}
};

10
src/OpenFOAM/primitives/Vector2D/Vector2DI.H
View File

@ -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<class Cmpt>
inline Foam::Vector2D<Cmpt>::Vector2D(const Foam::zero)
inline Foam::Vector2D<Cmpt>::Vector2D(Foam::zero)
:
Vector2D::vsType(Zero)
Vector2D::vsType(Foam::zero{})
{}
@ -48,8 +48,8 @@ inline Foam::Vector2D<Cmpt>::Vector2D
template<class Cmpt>
inline Foam::Vector2D<Cmpt>::Vector2D(const Cmpt& vx, const Cmpt& vy)
{
this->v_[X] = vx;
this->v_[Y] = vy;
this->x() = vx;
this->y() = vy;
}

50
src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
View File

@ -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<point>& 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)
{

4
src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
View File

@ -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]];
// }
//

8
src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
View File

@ -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<solveVector>(p[f[pi]]);
sumAc += p[f[pi]];
}
faceCentres[facei] = sumAc/nPoints;
}

55
src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C
View File

@ -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<point>& 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);

109
src/mesh/blockMesh/blocks/block/blockCreate.C
View File

@ -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] +=
(

2
src/meshTools/indexedOctree/treeDataFace.C
View File

@ -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);
}

2
src/meshTools/indexedOctree/treeDataPrimitivePatch.C
View File

@ -417,7 +417,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::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);
}

4
src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
View File

@ -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