diff --git a/applications/test/triTet/Test-triTet.C b/applications/test/triTet/Test-triTet.C index 5c850a4377..c85adc83a6 100644 --- a/applications/test/triTet/Test-triTet.C +++ b/applications/test/triTet/Test-triTet.C @@ -7,17 +7,20 @@ using namespace Foam; int main() { - triangle tri + const triangle tri ( vector(0, 0, 0), vector(1, 0, 0), vector(1, 1, 0) ); - Info<< "tri circumCentre = " << tri.circumCentre() << endl; - Info<< "tri circumRadius = " << tri.circumRadius() << endl; + const Tuple2 triCircle = tri.circumCircle(); - tetrahedron tet + Info<< "tri circumCentre = " << triCircle.first() << endl + << "tri circumRadius = " << triCircle.second() << endl + << " tri quality = " << tri.quality() << endl; + + const tetrahedron tet ( vector(1, 0, 0), vector(0, 1, 0), @@ -25,16 +28,11 @@ int main() vector(0.5773502, 0.5773502, 0.5773502) ); - Info<< "tet circumCentre = " << tet.circumCentre() << endl; - Info<< "tet circumRadius = " << tet.circumRadius() << endl; + const Tuple2 tetSphere = tet.circumSphere(); - vector a(Sin); - vector b(Sin); - vector c(Sin); - vector d(Sin); - - Info<< "tet circumRadius = " - << tetrahedron(a, b, c, d).circumRadius() << endl; + Info<< "tet circumCentre = " << tetSphere.first() << endl + << "tet circumRadius = " << tetSphere.second() << endl + << " tet quality = " << tet.quality() << endl; return 0; } diff --git a/applications/utilities/surface/surfaceFeatures/surfaceFeatures.C b/applications/utilities/surface/surfaceFeatures/surfaceFeatures.C index a538fedb0c..15807b4ef0 100644 --- a/applications/utilities/surface/surfaceFeatures/surfaceFeatures.C +++ b/applications/utilities/surface/surfaceFeatures/surfaceFeatures.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2018-2021 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2018-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -616,24 +616,22 @@ namespace Foam forAll(surf, fi) { const triPointRef& tri = surf[fi].tri(surf.points()); - const point& triCentre = tri.circumCentre(); - const scalar radiusSqr = min - ( - sqr(4*tri.circumRadius()), - sqr(searchDistance) - ); + const Tuple2 circle = tri.circumCircle(); + const point& c = circle.first(); + const scalar rSqr = + min(sqr(4*circle.second()), sqr(searchDistance)); pointIndexHitList hitList; - feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList); + feMesh.allNearestFeatureEdges(c, rSqr, hitList); featureProximity[fi] = min ( feMesh.minDisconnectedDist(hitList), featureProximity[fi] ); - feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList); + feMesh.allNearestFeaturePoints(c, rSqr, hitList); featureProximity[fi] = min ( minDist(hitList), diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index d898bd24dc..6f703fa9e6 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H @@ -152,11 +152,8 @@ public: //- Return volume inline scalar mag() const; - //- Return circum-centre - inline Point circumCentre() const; - - //- Return circum-radius - inline scalar circumRadius() const; + //- Return the circum centre and radius + inline Tuple2 circumSphere() const; //- Return quality: Ratio of tetrahedron and circum-sphere // volume, scaled so that a regular tetrahedron has a diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index f68cb207d2..d84dd5d520 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -174,68 +174,42 @@ inline Foam::scalar Foam::tetrahedron::mag() const template -inline Point Foam::tetrahedron::circumCentre() const +inline Foam::Tuple2 +Foam::tetrahedron::circumSphere() const { - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; + const vector a = b_ - a_; + const vector b = c_ - a_; + const vector c = d_ - a_; - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); + const vector ba = b ^ a; + const vector ca = c ^ a; - vector ba = b ^ a; - vector ca = c ^ a; + const scalar lambda = magSqr(c) - (a & c); + const scalar mu = magSqr(b) - (a & b); - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); + const vector num = lambda*ba - mu*ca; + const scalar denom = (c & ba); if (Foam::mag(denom) < rootVSmall) { - // Degenerate tetrahedron, returning centre instead of circumCentre. - - return centre(); + // Degenerate. Return a point far away. + static const scalar sqrt3 = sqrt(scalar(3)); + return Tuple2(point::uniform(great), sqrt3*great); } - - return a_ + 0.5*(a + num/denom); -} - - -template -inline Foam::scalar Foam::tetrahedron::circumRadius() const -{ - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; - - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); - - vector ba = b ^ a; - vector ca = c ^ a; - - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); - - if (Foam::mag(denom) < rootVSmall) + else { - // Degenerate tetrahedron, returning great for circumRadius. - return great; + const vector v = (a + num/denom)/2; + return Tuple2(a_ + v, Foam::mag(v)); } - - return Foam::mag(0.5*(a + num/denom)); } template inline Foam::scalar Foam::tetrahedron::quality() const { - return - mag() - /( - 8.0/(9.0*sqrt(3.0)) - *pow3(min(circumRadius(), great)) - + rootVSmall - ); + const scalar r = circumSphere().second(); + static const scalar sqrt3 = sqrt(scalar(3)); + return mag()/((8.0/27.0)*sqrt3*pow3(min(r, great)) + rootVSmall); } diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 4b81213ea9..f902a18b11 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -140,11 +140,8 @@ public: //- Return unit normal inline vector normal() const; - //- Return circum-centre - inline Point circumCentre() const; - - //- Return circum-radius - inline scalar circumRadius() const; + //- Return the circum centre and radius + inline Tuple2 circumCircle() const; //- Return quality: Ratio of triangle and circum-circle // area, scaled so that an equilateral triangle has a diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 852b1a4fee..03238f62ff 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -116,52 +116,33 @@ inline Foam::vector Foam::triangle::normal() const template -inline Point Foam::triangle::circumCentre() const +inline Foam::Tuple2 +Foam::triangle::circumCircle() const { - scalar d1 = (c_ - a_) & (b_ - a_); - scalar d2 = -(c_ - b_) & (b_ - a_); - scalar d3 = (c_ - a_) & (c_ - b_); + const scalar d1 = (c_ - a_) & (b_ - a_); + const scalar d2 = -(c_ - b_) & (b_ - a_); + const scalar d3 = (c_ - a_) & (c_ - b_); - scalar c1 = d2*d3; - scalar c2 = d3*d1; - scalar c3 = d1*d2; + const scalar c1 = d2*d3; + const scalar c2 = d3*d1; + const scalar c3 = d1*d2; - scalar c = c1 + c2 + c3; + const scalar denom = c1 + c2 + c3; - if (Foam::mag(c) < rootVSmall) + if (Foam::mag(denom) < rootVSmall) { - // Degenerate triangle, returning centre instead of circumCentre. - - return centre(); - } - - return - ( - ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c) - ); -} - - -template -inline Foam::scalar Foam::triangle::circumRadius() const -{ - scalar d1 = (c_ - a_) & (b_ - a_); - scalar d2 = -(c_ - b_) & (b_ - a_); - scalar d3 = (c_ - a_) & (c_ - b_); - - scalar denom = d2*d3 + d3*d1 + d1*d2; - - if (Foam::mag(denom) < vSmall) - { - // Degenerate triangle, returning great for circumRadius. - - return great; + // Degenerate. Return a point far away. + static const scalar sqrt3 = sqrt(scalar(3)); + return Tuple2(point::uniform(great), sqrt3*great); } else { - scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom; - - return 0.5*Foam::sqrt(min(great, max(0, a))); + return + Tuple2 + ( + ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*denom), + sqrt(max(0, (d1 + d2)*(d2 + d3)*(d3 + d1)/(4*denom))) + ); } } @@ -169,15 +150,9 @@ inline Foam::scalar Foam::triangle::circumRadius() const template inline Foam::scalar Foam::triangle::quality() const { - scalar c = circumRadius(); - - if (c < rootVSmall) - { - // zero circumRadius, something has gone wrong. - return small; - } - - return mag()/(Foam::sqr(c)*3.0*sqrt(3.0)/4.0); + const scalar r = circumCircle().second(); + static const scalar sqrt3 = sqrt(scalar(3)); + return mag()/(0.75*sqrt3*sqr(min(r, great)) + rootVSmall); }