diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 3a96a556e7ae77cb176eb7ad94fa486c242c7013..f1fe690f301675d10d99a2d5c520fcf1390dcdd6 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 OpenCFD Ltd.
+ Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -521,13 +521,12 @@ Foam::point Foam::face::centre(const UList<point>& points) const
// If the face is a triangle, do a direct calculation
if (nPoints == 3)
{
- return
- (1.0/3.0)
- *(
- points[operator[](0)]
- + points[operator[](1)]
- + points[operator[](2)]
- );
+ return triPointRef::centre
+ (
+ points[operator[](0)],
+ points[operator[](1)],
+ points[operator[](2)]
+ );
}
@@ -587,12 +586,12 @@ Foam::vector Foam::face::areaNormal(const UList<point>& p) const
if (nPoints == 3)
{
- return triPointRef
+ return triPointRef::areaNormal
(
p[operator[](0)],
p[operator[](1)],
p[operator[](2)]
- ).areaNormal();
+ );
}
label pI;
@@ -621,12 +620,12 @@ Foam::vector Foam::face::areaNormal(const UList<point>& p) const
// Note: for best accuracy, centre point always comes last
//
- n += triPointRef
+ n += triPointRef::areaNormal
(
p[operator[](pI)],
nextPoint,
centrePoint
- ).areaNormal();
+ );
}
return n;
diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
index 18e30c7b770b9c8b90ffe43fc5c689eb680af012..69ae3f889fccae3a6c70947b249ac2b3a79a83b5 100644
--- a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H
@@ -180,31 +180,34 @@ inline Foam::triPointRef Foam::triFace::tri(const UList<point>& points) const
inline Foam::point Foam::triFace::centre(const UList<point>& points) const
{
- return (1.0/3.0)*
+ return triPointRef::centre
(
- points[operator[](0)]
- + points[operator[](1)]
- + points[operator[](2)]
+ points[operator[](0)],
+ points[operator[](1)],
+ points[operator[](2)]
);
}
inline Foam::vector Foam::triFace::areaNormal(const UList<point>& points) const
{
- // As per triPointRef(...).areaNormal()
- return 0.5*
+ return triPointRef::areaNormal
(
- (points[operator[](1)] - points[operator[](0)])
- ^(points[operator[](2)] - points[operator[](0)])
+ points[operator[](0)],
+ points[operator[](1)],
+ points[operator[](2)]
);
}
inline Foam::vector Foam::triFace::unitNormal(const UList<point>& points) const
{
- const vector n(areaNormal(points));
- const scalar s(Foam::mag(n));
- return s < ROOTVSMALL ? Zero : n/s;
+ return triPointRef::unitNormal
+ (
+ points[operator[](0)],
+ points[operator[](1)],
+ points[operator[](2)]
+ );
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index e704fbc8c5257ab4945b787ba18f56d46c7abd1f..5d2d148128603beddc2adb6ddf161c08ab6e6502 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -399,14 +399,14 @@ public:
return UIndirectList<T>(internalValues, faceCells());
}
- //- Slice list to patch
+ //- Slice List to patch, using the number of patch faces
template<class T>
const typename List<T>::subList patchSlice(const UList<T>& l) const
{
return typename List<T>::subList(l, this->size(), start_);
}
- //- Slice Field to patch
+ //- Slice Field to patch, using the number of patch faces
template<class T>
const typename Field<T>::subField patchSlice(const Field<T>& l) const
{
@@ -449,9 +449,9 @@ public:
// Other patch operations
//- Return label of face in patch from global face label
- inline label whichFace(const label l) const
+ label whichFace(const label facei) const noexcept
{
- return l - start_;
+ return facei - start_;
}
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
index 78c40b4926c6a37c490010537a4c21e256574793..ad5e0aa63cd57de8604c3708bfaa9bb9cb5b0cfc 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
@@ -30,6 +30,7 @@ License
#include "primitiveMesh.H"
#include "syncTools.H"
#include "pyramid.H"
+#include "tetrahedron.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@@ -54,8 +55,8 @@ void Foam::primitiveMeshTools::updateFaceCentresAndAreas
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
- fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+ 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]]);
}
else
{
@@ -221,8 +222,8 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
- fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+ 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]]);
}
else
{
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
index 7cc6a63077991213fd205b8ac810fc0a82bbc18e..f9a9ba6293b294a8aa1dabf963ebf5f2182c90a6 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
@@ -168,7 +168,7 @@ public:
//- Operate on a triangle
result operator()(const FixedList<point, 3>& p) const
{
- return result(triPointRef(p[0], p[1], p[2]).areaNormal());
+ return result(triPointRef::areaNormal(p[0], p[1], p[2]));
}
};
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index 7613f380487569b396e356f5766c00d9f3e73f0b..4dc069a1341a040e6d3a995025a59a8ae819b61a 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
@@ -186,28 +186,28 @@ inline Foam::triPointRef Foam::tetrahedron<Point, PointRef>::tri
template<class Point, class PointRef>
inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sa() const
{
- return triangle<Point, PointRef>(b_, c_, d_).areaNormal();
+ return triangle<Point, PointRef>::areaNormal(b_, c_, d_);
}
template<class Point, class PointRef>
inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sb() const
{
- return triangle<Point, PointRef>(a_, d_, c_).areaNormal();
+ return triangle<Point, PointRef>::areaNormal(a_, d_, c_);
}
template<class Point, class PointRef>
inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sc() const
{
- return triangle<Point, PointRef>(a_, b_, d_).areaNormal();
+ return triangle<Point, PointRef>::areaNormal(a_, b_, d_);
}
template<class Point, class PointRef>
inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sd() const
{
- return triangle<Point, PointRef>(a_, c_, b_).areaNormal();
+ return triangle<Point, PointRef>::areaNormal(a_, c_, b_);
}
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 6f1ed0a1aeee7243dbba8b6c602cda297cdac3c5..83d61172e48ab06695208d0979725c11c1e0f99d 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -129,11 +129,26 @@ public:
//- The third vertex
point& c() { return operator[](2); }
+ //- Flip triangle orientation by swapping second and third vertices
+ inline void flip();
+
//- Return as triangle reference
inline triPointRef tri() const;
- //- Flip triangle orientation by swapping second and third vertices
- inline void flip();
+
+ // Properties
+
+ //- Return centre (centroid)
+ inline point centre() const;
+
+ //- The area normal - with magnitude equal to area of triangle
+ inline vector areaNormal() const;
+
+ //- Return unit normal
+ inline vector unitNormal() const;
+
+ //- The magnitude of the triangle area
+ inline scalar mag() const;
//- The bounding box for the tetrahedron
inline treeBoundBox bounds() const;
@@ -268,6 +283,35 @@ public:
const Point& c() const noexcept { return c_; }
+ // Triangle properties (static calculations)
+
+ //- The centre (centroid) of three points
+ inline static Point centre
+ (
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+ );
+
+ //- The area normal for a triangle defined by three points
+ //- (right-hand rule). Magnitude equal to area of the triangle
+ inline static vector areaNormal
+ (
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+ );
+
+ //- The unit normal for a triangle defined by three points
+ //- (right-hand rule).
+ inline static vector unitNormal
+ (
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+ );
+
+
// Properties
//- Return centre (centroid)
@@ -287,7 +331,7 @@ public:
return areaNormal();
}
- //- Return scalar magnitude
+ //- The magnitude of the triangle area
inline scalar mag() const;
//- Return circum-centre
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index 8e20ec1e760daf955c9dce88f1cc765192a6c578..3be9f01eba9d20d2e48b8cd6102e6d33d061cdc0 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -117,33 +117,40 @@ inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-inline Foam::triPointRef Foam::triPoints::tri() const
+template<class Point, class PointRef>
+inline Point Foam::triangle<Point, PointRef>::centre
+(
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+)
{
- return triPointRef(a(), b(), c());
+ return (1.0/3.0)*(p0 + p1 + p2);
}
-inline void Foam::triPoints::flip()
+template<class Point, class PointRef>
+inline Point Foam::triangle<Point, PointRef>::centre() const
{
- // swap pt1 <-> pt2
- point t(b());
- b() = c();
- c() = t;
+ return (1.0/3.0)*(a_ + b_ + c_);
}
-inline Foam::treeBoundBox Foam::triPoints::bounds() const
+inline Foam::point Foam::triPoints::centre() const
{
- treeBoundBox bb;
- bb.add(static_cast<const FixedList<point, 3>&>(*this));
- return bb;
+ return triPointRef::centre(a(), b(), c());
}
template<class Point, class PointRef>
-inline Point Foam::triangle<Point, PointRef>::centre() const
+inline Foam::vector Foam::triangle<Point, PointRef>::areaNormal
+(
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+)
{
- return (1.0/3.0)*(a_ + b_ + c_);
+ return 0.5*((p1 - p0)^(p2 - p0));
}
@@ -154,15 +161,70 @@ inline Foam::vector Foam::triangle<Point, PointRef>::areaNormal() const
}
+inline Foam::vector Foam::triPoints::areaNormal() const
+{
+ return triPointRef::areaNormal(a(), b(), c());
+}
+
+
template<class Point, class PointRef>
-inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal() const
+inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal
+(
+ const Point& p0,
+ const Point& p1,
+ const Point& p2
+)
{
- const vector n(areaNormal());
- const scalar s(Foam::mag(n));
+ const vector n(areaNormal(p0, p1, p2));
+ const scalar s(::Foam::mag(n));
return s < ROOTVSMALL ? Zero : n/s;
}
+template<class Point, class PointRef>
+inline Foam::vector Foam::triangle<Point, PointRef>::unitNormal() const
+{
+ return triPointRef::unitNormal(a_, b_, c_);
+}
+
+
+inline Foam::vector Foam::triPoints::unitNormal() const
+{
+ return triPointRef::unitNormal(a(), b(), c());
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+inline Foam::triPointRef Foam::triPoints::tri() const
+{
+ return triPointRef(a(), b(), c());
+}
+
+
+inline void Foam::triPoints::flip()
+{
+ // swap pt1 <-> pt2
+ point t(b());
+ b() = c();
+ c() = t;
+}
+
+
+inline Foam::treeBoundBox Foam::triPoints::bounds() const
+{
+ treeBoundBox bb;
+ bb.add(static_cast<const FixedList<point, 3>&>(*this));
+ return bb;
+}
+
+
+inline Foam::scalar Foam::triPoints::mag() const
+{
+ return ::Foam::mag(areaNormal());
+}
+
+
template<class Point, class PointRef>
inline Foam::scalar Foam::triangle<Point, PointRef>::mag() const
{
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 89e717596af5d71f9617f5031932a6236d26320b..77e655392ab59c5d46eba0e6f1fcd298dfa016da 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -61,8 +61,10 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
- faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+ 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]]);
}
else
{
@@ -1631,12 +1633,12 @@ bool Foam::polyMeshGeometry::checkFaceTwist
{
vector triArea
(
- triPointRef
+ triPointRef::areaNormal
(
p[f[fpI]],
p[f.nextLabel(fpI)],
fc
- ).areaNormal()
+ )
);
scalar magTri = mag(triArea);
@@ -1732,12 +1734,12 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
forAll(f, fp)
{
- prevN = triPointRef
+ prevN = triPointRef::areaNormal
(
p[f[fp]],
p[f.nextLabel(fp)],
fc
- ).areaNormal();
+ );
scalar magTri = mag(prevN);
@@ -1759,12 +1761,12 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
vector triN
(
- triPointRef
+ triPointRef::areaNormal
(
p[f[fp]],
p[f.nextLabel(fp)],
fc
- ).areaNormal()
+ )
);
scalar magTri = mag(triN);
diff --git a/src/fileFormats/stl/STLtriangleI.H b/src/fileFormats/stl/STLtriangleI.H
index 6ecf80f68eb7fe9aa11b90f9f5d8020172aefcf6..c8ab26c6e60799ccdb6166c888e0b88c5ad002e8 100644
--- a/src/fileFormats/stl/STLtriangleI.H
+++ b/src/fileFormats/stl/STLtriangleI.H
@@ -113,7 +113,7 @@ inline void Foam::STLtriangle::write
)
{
// Calculate the normal ourselves
- const vector norm = triPointRef(p0, p1, p2).unitNormal();
+ const vector norm = triPointRef::unitNormal(p0, p1, p2);
write(os, norm, p0, p1, p2);
}
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index 8430af759710a178d38084330716615231eeb872..a6848fc4944f3e163f40de07317e0a1a12d14f46 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -462,10 +462,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::edgeNormals() const
{
auto tedgeNorm = tmp<vectorField>::New(edgeLengths());
- for (vector& n : tedgeNorm.ref())
- {
- n.normalise();
- }
+ tedgeNorm.ref().normalise();
return tedgeNorm;
}
@@ -473,27 +470,15 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::edgeNormals() const
Foam::tmp<Foam::vectorField> Foam::faPatch::edgeFaceCentres() const
{
- auto tfc = tmp<vectorField>::New(size());
- auto& fc = tfc.ref();
-
- // get reference to global face centres
- const vectorField& gfc =
- boundaryMesh().mesh().areaCentres().internalField();
-
- const labelUList& faceLabels = edgeFaces();
-
- forAll(faceLabels, edgeI)
- {
- fc[edgeI] = gfc[faceLabels[edgeI]];
- }
-
- return tfc;
+ return patchInternalField(boundaryMesh().mesh().areaCentres());
}
Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
{
- return edgeNormals()*(edgeNormals() & (edgeCentres() - edgeFaceCentres()));
+ // Use patch-normal delta for all non-coupled BCs
+ const vectorField nHat(edgeNormals());
+ return nHat*(nHat & (edgeCentres() - edgeFaceCentres()));
}
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
index a1bf87933def5867768582c6c110d8603108165d..7c2cecc140a9b5063c1d60e64213357ef76718cb 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
@@ -308,25 +308,32 @@ public:
//- Patch start in edge list
label start() const;
- //- Patch size is the number of edge labels
+ //- Patch size is the number of edge labels, but can be overloaded
virtual label size() const
{
return labelList::size();
}
//- Return label of edge in patch from global edge label
- inline label whichEdge(const label l) const
+ label whichEdge(const label edgei) const
{
- return l - start();
+ return edgei - start();
}
- //- Slice list to patch
+ //- Slice List to patch, using the virtual patch size
template<class T>
typename List<T>::subList patchSlice(const List<T>& l) const
{
return typename List<T>::subList(l, size(), start());
}
+ //- Slice List to patch, using the number of patch edges
+ template<class T>
+ typename List<T>::subList patchRawSlice(const List<T>& l) const
+ {
+ return typename List<T>::subList(l, nEdges(), start());
+ }
+
//- Write
virtual void write(Ostream&) const;
diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
index 89e08fd3779bf1a6e208c77835b285389e226ba0..1202b80880422b7c916be8802e2fbfad2fe5d67a 100644
--- a/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
+++ b/src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C
@@ -28,6 +28,7 @@ License
#include "highAspectRatioFvGeometryScheme.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
+#include "triangle.H"
#include "syncTools.H"
#include "cellAspectRatio.H"
#include "emptyPolyPatch.H"
@@ -259,7 +260,7 @@ void Foam::highAspectRatioFvGeometryScheme::makeAverageCentres
if (nPoints == 3)
{
- faceCentres[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
+ faceCentres[facei] = triPointRef::centre(p[f[0]], p[f[1]], p[f[2]]);
}
else
{
diff --git a/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C b/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C
index 78d10ceb65ad0dc5649a43b2fa66410656dfdbd9..90c3c4e00f320fefd6712645a42646cb3cfd1a39 100644
--- a/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C
+++ b/src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C
@@ -31,6 +31,7 @@ License
#include "fvMesh.H"
#include "PrecisionAdaptor.H"
#include "primitiveMeshTools.H"
+#include "triangle.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -67,8 +68,8 @@ void Foam::stabilisedFvGeometryScheme::makeFaceCentresAndAreas
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
- fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+ 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]]);
}
// For more complex faces, decompose into triangles
diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
index 4d30c3baaab1673c83025624554b7651c5cc4939..e286fbc5df369d6750e06bf266f8eb9722855f49 100644
--- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
@@ -125,20 +125,7 @@ const Foam::vectorField& Foam::fvPatch::Cf() const
Foam::tmp<Foam::vectorField> Foam::fvPatch::Cn() const
{
- auto tcc = tmp<vectorField>::New(size());
- auto& cc = tcc.ref();
-
- const labelUList& faceCells = this->faceCells();
-
- // get reference to global cell centres
- const vectorField& gcc = boundaryMesh().mesh().cellCentres();
-
- forAll(faceCells, facei)
- {
- cc[facei] = gcc[faceCells[facei]];
- }
-
- return tcc;
+ return patchInternalField(boundaryMesh().mesh().cellCentres());
}
diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H
index cb9569d5747c022ce1c3dd281550f39c6b28efa2..0e575be72dbc7920cb462a35b0270ebd5c9e9a67 100644
--- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H
+++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H
@@ -165,7 +165,7 @@ public:
// Access
//- Return the polyPatch
- const polyPatch& patch() const
+ const polyPatch& patch() const noexcept
{
return polyPatch_;
}
@@ -201,24 +201,39 @@ public:
static wordList constraintTypes();
//- Return the index of this patch in the fvBoundaryMesh
- label index() const
+ label index() const noexcept
{
return polyPatch_.index();
}
//- Return boundaryMesh reference
- const fvBoundaryMesh& boundaryMesh() const
+ const fvBoundaryMesh& boundaryMesh() const noexcept
{
return boundaryMesh_;
}
- //- Slice list to patch
+ //- Slice List to patch, using the virtual patch size
template<class T>
const typename List<T>::subList patchSlice(const List<T>& l) const
{
return typename List<T>::subList(l, size(), start());
}
+ //- Slice List to patch, using the underlying polyPatch information
+ template<class T>
+ const typename List<T>::subList patchRawSlice
+ (
+ const List<T>& l
+ ) const
+ {
+ return typename List<T>::subList
+ (
+ l,
+ polyPatch_.size(),
+ polyPatch_.start()
+ );
+ }
+
//- Return faceCells
virtual const labelUList& faceCells() const;
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
index de522dcdb9f3dd2e60b561e65c80e0cfba619220..6d8cad97f713f2060432527384ce077b9c3a33f7 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C
@@ -67,7 +67,7 @@ correction
const face& f = faces[facei];
- scalar at = triangle<point, const point&>
+ scalar at = triPointRef
(
pi,
points[f[0]],
@@ -84,7 +84,7 @@ correction
for (label pointi=1; pointi<f.size(); pointi++)
{
- at = triangle<point, const point&>
+ at = triPointRef
(
pi,
points[f[pointi]],
@@ -129,7 +129,7 @@ correction
const face& f = faces[facei+fvp.start()];
- scalar at = triangle<point, const point&>
+ scalar at = triPointRef
(
pi,
points[f[0]],
@@ -146,7 +146,7 @@ correction
for (label pointi=1; pointi<f.size(); pointi++)
{
- at = triangle<point, const point&>
+ at = triPointRef
(
pi,
points[f[pointi]],
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 79532b6d8243f7f5c5c480baf43254adc5350109..2414584e953c10d2ea03880ec471547f217a7047 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -315,7 +315,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);
- n = triPointRef(base[0], vertex1[0], vertex2[0]).unitNormal();
+ n = triPointRef::unitNormal(base[0], vertex1[0], vertex2[0]);
// Interpolate the motion of the three face vertices to the current
// coordinates
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
index 8976e57b586dc365a0b6dd4366d35a19fc567a82..d2c358e0c0147d8c50d52905f7cbaff538ca4424 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2018 OpenCFD Ltd.
+ Copyright (C) 2018-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -240,7 +240,7 @@ void Foam::faceAreaIntersect::triangleIntersect
// Cut source triangle with all inward pointing faces of target triangle
// - triangles in workTris1 are inside target triangle
- const scalar srcArea(triArea(src));
+ const scalar srcArea(src.mag());
if (srcArea < ROOTVSMALL)
{
return;
@@ -363,11 +363,11 @@ void Foam::faceAreaIntersect::triangleIntersect
for (label i = 0; i < nWorkTris1; ++i)
{
// Area of intersection
- const scalar currArea = triArea(workTris1[i]);
+ const scalar currArea = workTris1[i].mag();
area += currArea;
// Area-weighted centroid of intersection
- centroid += currArea*triCentroid(workTris1[i]);
+ centroid += currArea*workTris1[i].centre();
if (cacheTriangulation_)
{
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
index fc2d559021e7044c029c5fae8a2e71226e780d9e..7a45177352ebc39ab7fd6487b21cd885a50121c6 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
@@ -131,12 +131,6 @@ private:
const label posI
) const;
- //- Return triangle area
- inline scalar triArea(const triPoints& t) const;
-
- //- Return triangle centre
- inline vector triCentroid(const triPoints& t) const;
-
//- Slice triangle with plane and generate new cut sub-triangles
void triSliceWithPlane
(
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
index 836b87d76200a29f5ec519dfca87d7559d1a50ad..7541b14141cff5feace46313a2f501be49e59d21 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
@@ -107,21 +107,6 @@ inline Foam::point Foam::faceAreaIntersect::planeIntersection
}
-inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
-{
- return mag(0.5*((t[1] - t[0])^(t[2] - t[0])));
-}
-
-
-inline Foam::vector Foam::faceAreaIntersect::triCentroid
-(
- const triPoints& t
-) const
-{
- return (t[0] + t[1] + t[2])/3;
-}
-
-
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar& Foam::faceAreaIntersect::tolerance()
diff --git a/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H b/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H
index 6f73594f60f3996160df5501c80cbb04a91825eb..53a9c56054529483e4b66fe034243f00f27a710c 100644
--- a/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H
+++ b/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H
@@ -93,7 +93,7 @@ inline Foam::vector2D Foam::triangle2D::centre() const
{
const triangle2D& tri = *this;
- return (1/3.0)*(tri[0] + tri[1] + tri[2]);
+ return (1.0/3.0)*(tri[0] + tri[1] + tri[2]);
}
diff --git a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
index a7b0a891139192c9e2fa84a4d6da9ee011bb066a..75ac9dc72b753bd2dadac96be26dc7e2cd65cfe7 100644
--- a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
+++ b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
@@ -58,8 +58,11 @@ void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
- faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+ 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]]);
}
else
{
@@ -1004,12 +1007,12 @@ bool Foam::primitiveMeshGeometry::checkFaceTwist
{
const vector triArea
(
- triPointRef
+ triPointRef::areaNormal
(
p[f[fpI]],
p[f.nextLabel(fpI)],
fc
- ).areaNormal()
+ )
);
const scalar magTri = mag(triArea);
diff --git a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
index b673c238bab331d62bb28069e799ccb16ee379e7..88afcb88850b834ab9c940d934d91f91ae2cacf0 100644
--- a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
+++ b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
@@ -28,7 +28,6 @@ License
#include "faceTriangulation.H"
#include "plane.H"
-
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1e-6;
@@ -160,9 +159,9 @@ bool Foam::faceTriangulation::triangleContainsPoint
const point& pt
)
{
- const scalar area01Pt = triPointRef(p0, p1, pt).areaNormal() & n;
- const scalar area12Pt = triPointRef(p1, p2, pt).areaNormal() & n;
- const scalar area20Pt = triPointRef(p2, p0, pt).areaNormal() & n;
+ const scalar area01Pt = triPointRef::areaNormal(p0, p1, pt) & n;
+ const scalar area12Pt = triPointRef::areaNormal(p1, p2, pt) & n;
+ const scalar area20Pt = triPointRef::areaNormal(p2, p0, pt) & n;
if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
{
diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
index dddd5af17dd540740e21202caf2905aedf255eb4..89db16bd039419a6f3c3ba693505e703eb473354 100644
--- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
@@ -42,7 +42,7 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell
{
// Calculate the normal ourselves, for flexibility and speed
const vector norm =
- triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal();
+ triPointRef::unitNormal(pts[f[0]], pts[f[1]], pts[f[2]]);
// simple triangulation about f[0].
// better triangulation should have been done before
@@ -76,7 +76,7 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell
{
// Calculate the normal ourselves, for flexibility and speed
const vector norm =
- triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal();
+ triPointRef::unitNormal(pts[f[0]], pts[f[1]], pts[f[2]]);
// simple triangulation about f[0].
// better triangulation should have been done before
diff --git a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
index b9c258b8e29048c32bac20ceedccea5d4fd53f89..7bc7868424f4b8a9c1c0b1f4dc60fd0eab021a34 100644
--- a/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
+++ b/src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
@@ -28,6 +28,7 @@ License
\*---------------------------------------------------------------------------*/
#include "cutFace.H"
+#include "triangle.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -274,11 +275,19 @@ void Foam::cutFace::calcSubFaceCentreAndArea
// and to avoid round-off error-related problems
if (nPoints == 3)
{
- subFaceCentre =
- (1.0/3.0)*(subFacePoints[0] + subFacePoints[1] + subFacePoints[2]);
+ subFaceCentre = triPointRef::centre
+ (
+ subFacePoints[0],
+ subFacePoints[1],
+ subFacePoints[2]
+ );
- subFaceArea = 0.5*((subFacePoints[1] - subFacePoints[0]) ^
- (subFacePoints[2] - subFacePoints[0]));
+ subFaceArea = triPointRef::areaNormal
+ (
+ subFacePoints[0],
+ subFacePoints[1],
+ subFacePoints[2]
+ );
}
else if (nPoints > 0)
{