diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 3a670f39a4..46a3adaea7 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -156,29 +156,12 @@ inline Point tetrahedron::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::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::circumCentre() const template inline scalar tetrahedron::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::circumCentre() const") + << "Degenerate tetrahedron:" << nl << *this << nl + << "Returning GREAT for circumRadius." + << endl; + + return GREAT; + } + + return Foam::mag(0.5*(a + num/denom)); } template inline scalar tetrahedron::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::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(4, 0.25); diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 7840d8ee7b..b99ba990cd 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -109,9 +109,9 @@ inline vector triangle::normal() const template inline Point triangle::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::circumCentre() const scalar c = c1 + c2 + c3; + if (Foam::mag(c) < ROOTVSMALL) + { + WarningIn("Point triangle::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::circumCentre() const template inline scalar triangle::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::circumRadius() const") + << "Degenerate triangle:" << nl << *this << nl + << "Returning GREAT for circumRadius." + << endl; + return GREAT; } else @@ -151,16 +166,7 @@ inline scalar triangle::circumRadius() const template inline scalar triangle::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::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(3, 1.0/3.0); @@ -490,7 +497,7 @@ pointHit triangle::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::nearestPointClassify if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) { + if ((d1 - d3) < ROOTVSMALL) + { + WarningIn + ( + "pointHit triangle::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::nearestPointClassify if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) { + if ((d2 - d6) < ROOTVSMALL) + { + WarningIn + ( + "pointHit triangle::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::nearestPointClassify if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) { + if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL) + { + WarningIn + ( + "pointHit triangle::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::nearestPointClassify // P inside face region. Compute Q through its barycentric // coordinates (u, v, w) + + if ((va + vb + vc) < ROOTVSMALL) + { + WarningIn + ( + "pointHit triangle::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;