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);
}