From 6cd78ff36855d451512a55ffb0d66e5a5ea3aeed Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 8 Dec 2010 16:00:59 +0000 Subject: [PATCH] ENH: tetrahedron : handle degenerate tets --- .../tetrahedron/tetrahedronI.H | 35 +++++++++++++++++-- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 0939c060e1..3a670f39a4 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -21,12 +21,11 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Description - \*---------------------------------------------------------------------------*/ #include "triangle.H" #include "IOstreams.H" +#include "triPointRef.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -152,7 +151,37 @@ inline Point tetrahedron::circumCentre() const vector ba = b ^ a; vector ca = c ^ a; - return a_ + 0.5*(a + (lambda*ba - mu*ca)/((c & ba) + ROOTVSMALL)); + vector num = lambda*ba - mu*ca; + scalar denom = (c & ba); + + 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(); + } + + return a_ + 0.5*(a + num/denom); }