Skip to content

GitLab

Commit 7bcd87eb authored by graham's avatar graham
Browse files

ENH: triangle and tetrahedron robustness improvements. 

Currently very verbose when a degenerate shape is encountered.
parent aa15f447
Hide whitespace changes
Inline Side-by-side
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
View file @ 7bcd87eb
......@@ -156,29 +156,12 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
if (Foam::mag(denom) < ROOTVSMALL)
{
// Degenerate tet. Use test of the individual triangles.
{
point triCentre = triPointRef(a_, b_, c_).circumCentre();
if (magSqr(d_ - triCentre) < magSqr(a_ - triCentre))
{
return triCentre;
}
}
{
point triCentre = triPointRef(a_, b_, d_).circumCentre();
if (magSqr(c_ - triCentre) < magSqr(a_ - triCentre))
{
return triCentre;
}
}
{
point triCentre = triPointRef(a_, c_, d_).circumCentre();
if (magSqr(b_ - triCentre) < magSqr(a_ - triCentre))
{
return triCentre;
}
}
return triPointRef(b_, c_, d_).circumCentre();
WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
<< "Degenerate tetrahedron:" << nl << *this << nl
<<"Returning centre instead of circumCentre."
<< endl;
return centre();
}
return a_ + 0.5*(a + num/denom);
......@@ -188,16 +171,43 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
template<class Point, class PointRef>
inline scalar tetrahedron<Point, PointRef>::circumRadius() const
{
return Foam::mag(a_ - circumCentre());
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)
{
WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
<< "Degenerate tetrahedron:" << nl << *this << nl
<< "Returning GREAT for circumRadius."
<< endl;
return GREAT;
}
return Foam::mag(0.5*(a + num/denom));
}
template<class Point, class PointRef>
inline scalar tetrahedron<Point, PointRef>::quality() const
{
// Note: 8/(9*sqrt(3)) = 0.5132002393
return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
return
mag()
/(
8.0/(9.0*sqrt(3.0))
*pow3(min(circumRadius(), GREAT))
+ ROOTVSMALL
);
}
......@@ -266,7 +276,8 @@ scalar tetrahedron<Point, PointRef>::barycentric
"const point& pt"
") const"
)
<< "Degenerate tetrahedron - returning 1/4 barycentric coordinates."
<< "Degenerate tetrahedron:" << nl << *this << nl
<< "Returning 1/4 barycentric coordinates."
<< endl;
bary = List<scalar>(4, 0.25);
......
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View file @ 7bcd87eb
......@@ -109,9 +109,9 @@ inline vector triangle<Point, PointRef>::normal() const
template<class Point, class PointRef>
inline Point triangle<Point, PointRef>::circumCentre() const
{
scalar d1 = (c_ - a_)&(b_ - a_);
scalar d2 = -(c_ - b_)&(b_ - a_);
scalar d3 = (c_ - a_)&(c_ - b_);
scalar d1 = (c_ - a_) & (b_ - a_);
scalar d2 = -(c_ - b_) & (b_ - a_);
scalar d3 = (c_ - a_) & (c_ - b_);
scalar c1 = d2*d3;
scalar c2 = d3*d1;
......@@ -119,6 +119,16 @@ inline Point triangle<Point, PointRef>::circumCentre() const
scalar c = c1 + c2 + c3;
if (Foam::mag(c) < ROOTVSMALL)
{
WarningIn("Point triangle<Point, PointRef>::circumCentre() const")
<< "Degenerate triangle:" << nl << *this << nl
<< "Returning centre instead of circumCentre."
<< endl;
return centre();
}
return
(
((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
......@@ -129,14 +139,19 @@ inline Point triangle<Point, PointRef>::circumCentre() const
template<class Point, class PointRef>
inline scalar triangle<Point, PointRef>::circumRadius() const
{
scalar d1 = (c_ - a_) & (b_ - a_);
scalar d2 = - (c_ - b_) & (b_ - a_);
scalar d3 = (c_ - a_) & (c_ - b_);
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)
{
WarningIn("scalar triangle<Point, PointRef>::circumRadius() const")
<< "Degenerate triangle:" << nl << *this << nl
<< "Returning GREAT for circumRadius."
<< endl;
return GREAT;
}
else
......@@ -151,16 +166,7 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
template<class Point, class PointRef>
inline scalar triangle<Point, PointRef>::quality() const
{
// Note: 3*sqr(3)/(4*pi) = 0.4134966716
return
mag()
/ (
constant::mathematical::pi
*Foam::sqr(circumRadius())
*0.4134966716
+ VSMALL
);
return mag()/(Foam::sqr(circumRadius())*3.0*sqrt(3.0)/4.0 + VSMALL);
}
......@@ -264,7 +270,8 @@ scalar triangle<Point, PointRef>::barycentric
"const point& pt"
") const"
)
<< "Degenerate triangle - returning 1/3 barycentric coordinates."
<< "Degenerate triangle:" << nl << *this << nl
<< "Returning 1/3 barycentric coordinates."
<< endl;
bary = List<scalar>(3, 1.0/3.0);
......@@ -490,7 +497,7 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
) const
{
// Adapted from:
// Real-time collision detection, Christer Ericson, 2005, 136-142
// Real-time collision detection, Christer Ericson, 2005, p136-142
// Check if P in vertex region outside A
vector ab = b_ - a_;
......@@ -528,6 +535,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
{
if ((d1 - d3) < ROOTVSMALL)
{
WarningIn
(
"pointHit triangle<Point, PointRef>::nearestPointClassify"
"("
"const point& p,"
"label& nearType,"
"label& nearLabel"
") const"
)
<< "Degenerate triangle:" << nl << *this << nl
<< "d1, d3: " << d1 << ", " << d3 << endl;
// For d1 = d3, a_ and b_ are likely coincident.
nearType = POINT;
nearLabel = 0;
return pointHit(false, a_, Foam::mag(a_ - p), true);
}
// barycentric coordinates (1-v, v, 0)
scalar v = d1/(d1 - d3);
......@@ -556,6 +584,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
{
if ((d2 - d6) < ROOTVSMALL)
{
WarningIn
(
"pointHit triangle<Point, PointRef>::nearestPointClassify"
"("
"const point& p,"
"label& nearType,"
"label& nearLabel"
") const"
)
<< "Degenerate triangle:" << nl << *this << nl
<< "d2, d6: " << d2 << ", " << d6 << endl;
// For d2 = d6, a_ and c_ are likely coincident.
nearType = POINT;
nearLabel = 0;
return pointHit(false, a_, Foam::mag(a_ - p), true);
}
// barycentric coordinates (1-w, 0, w)
scalar w = d2/(d2 - d6);
......@@ -570,6 +619,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
{
if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
{
WarningIn
(
"pointHit triangle<Point, PointRef>::nearestPointClassify"
"("
"const point& p,"
"label& nearType,"
"label& nearLabel"
") const"
)
<< "Degenerate triangle:" << nl << *this << nl
<< "(d4 - d3), (d6 - d5): " << (d4 - d3) << ", " << (d6 - d5)
<< endl;
// For (d4 - d3) = (d6 - d5), b_ and c_ are likely coincident.
nearType = POINT;
nearLabel = 1;
return pointHit(false, b_, Foam::mag(b_ - p), true);
}
// barycentric coordinates (0, 1-w, w)
scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
......@@ -581,6 +652,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
// P inside face region. Compute Q through its barycentric
// coordinates (u, v, w)
if ((va + vb + vc) < ROOTVSMALL)
{
WarningIn
(
"pointHit triangle<Point, PointRef>::nearestPointClassify"
"("
"const point& p,"
"label& nearType,"
"label& nearLabel"
") const"
)
<< "Degenerate triangle:" << nl << *this << nl
<< "va, vb, vc: " << va << ", " << vb << ", " << vc
<< endl;
point nearPt = centre();
nearType = NONE,
nearLabel = -1;
return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
}
scalar denom = 1.0/(va + vb + vc);
scalar v = vb * denom;
scalar w = vc * denom;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment