From 9fa37ba068e7f4a6f0d436bd2a5ed2d2aef9a694 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 23 Sep 2022 16:25:12 +0200
Subject: [PATCH] ENH: add static centre(), {area,unit}Normal() methods to
triangle
- commonly used calculations
ENH: add faPatch::patchRawSlice method
- slices using the nEdges() instead of the virtual size(),
which provides similar functionality as finite-volume has with
its distinction between polyPatch vs fvPatch patchSlice
- use patchInternal for obtaining faPatch, fvPatch information
---
src/OpenFOAM/meshes/meshShapes/face/face.C | 23 +++--
.../meshes/meshShapes/triFace/triFaceI.H | 25 ++---
.../polyPatches/polyPatch/polyPatch.H | 8 +-
.../primitiveMeshCheck/primitiveMeshTools.C | 9 +-
src/OpenFOAM/meshes/primitiveShapes/cut/cut.H | 2 +-
.../tetrahedron/tetrahedronI.H | 8 +-
.../primitiveShapes/triangle/triangle.H | 50 +++++++++-
.../primitiveShapes/triangle/triangleI.H | 94 +++++++++++++++----
.../polyMeshGeometry/polyMeshGeometry.C | 18 ++--
src/fileFormats/stl/STLtriangleI.H | 2 +-
.../faMesh/faPatches/faPatch/faPatch.C | 25 +----
.../faMesh/faPatches/faPatch/faPatch.H | 15 ++-
.../highAspectRatioFvGeometryScheme.C | 3 +-
.../stabilised/stabilisedFvGeometryScheme.C | 5 +-
.../fvMesh/fvPatches/fvPatch/fvPatch.C | 15 +--
.../fvMesh/fvPatches/fvPatch/fvPatch.H | 23 ++++-
.../schemes/pointLinear/pointLinear.C | 8 +-
src/lagrangian/basic/particle/particleI.H | 2 +-
.../faceAreaIntersect/faceAreaIntersect.C | 8 +-
.../faceAreaIntersect/faceAreaIntersect.H | 6 --
.../faceAreaIntersect/faceAreaIntersectI.H | 15 ---
.../AMIInterpolation/triangle2D/triangle2DI.H | 2 +-
.../primitiveMeshGeometry.C | 11 ++-
.../faceTriangulation/faceTriangulation.C | 7 +-
.../surfaceFormats/stl/STLsurfaceFormat.C | 4 +-
.../geometricVoF/cellCuts/cutFace/cutFace.C | 17 +++-
26 files changed, 251 insertions(+), 154 deletions(-)
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 3a96a556e7a..f1fe690f301 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 18e30c7b770..69ae3f889fc 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 e704fbc8c52..5d2d1481286 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 78c40b4926c..ad5e0aa63cd 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 7cc6a630779..f9a9ba6293b 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 7613f380487..4dc069a1341 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 6f1ed0a1aee..83d61172e48 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 8e20ec1e760..3be9f01eba9 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 89e717596af..77e655392ab 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 6ecf80f68eb..c8ab26c6e60 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 8430af75971..a6848fc4944 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 a1bf87933de..7c2cecc140a 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 89e08fd3779..1202b808804 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 78d10ceb65a..90c3c4e00f3 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 4d30c3baaab..e286fbc5df3 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 cb9569d5747..0e575be72db 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 de522dcdb9f..6d8cad97f71 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 79532b6d824..2414584e953 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 8976e57b586..d2c358e0c01 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 fc2d559021e..7a45177352e 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 836b87d7620..7541b14141c 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 6f73594f60f..53a9c560545 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 a7b0a891139..75ac9dc72b7 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 b673c238bab..88afcb88850 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 dddd5af17dd..89db16bd039 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 b9c258b8e29..7bc7868424f 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)
{
--
GitLab