ENH: tetOverlapVolume - improved robustness. Fixes #207

src/meshTools/tetOverlapVolume/tetOverlapVolume.C

@@ -39,12 +39,6 @@ namespace Foam
     defineTypeNameAndDebug(tetOverlapVolume, 0);
 }
 
-// When to consider a tet to be zero volume. We want to avoid doing clipping
-// against negative volume tets. Tet volume can be calculated incorrectly
-// due to truncation errors. The value below works for single and double
-// precision but could probably be tighter for double precision.
-Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/meshTools/tetOverlapVolume/tetOverlapVolume.H

@@ -54,12 +54,6 @@ class polyMesh;
 
 class tetOverlapVolume
 {
-    // Private data
-
-        //- Minimum tet volume to skip test
-        static scalar minTetVolume_;
-
-
     // Private classes
 
         //- tetPoints handling : sum resulting volumes

src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,21 +45,52 @@ void Foam::tetOverlapVolume::tetTetOverlap
     tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
     tetPointRef::dummyOp outside;
 
-    if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
+    // Precompute the tet face areas and exit early if any face area is
+    // too small
+    static FixedList<vector, 4> tetAFaceAreas;
+    static FixedList<scalar, 4> tetAMag2FaceAreas;
+    tetPointRef tetATet = tetA.tet();
+    for (label facei = 0; facei < 4; ++facei)
+    {
+        tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+        tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
+        if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
         {
             return;
         }
+    }
+
+    static FixedList<vector, 4> tetBFaceAreas;
+    static FixedList<scalar, 4> tetBMag2FaceAreas;
+    tetPointRef tetBTet = tetB.tet();
+    for (label facei = 0; facei < 4; ++facei)
+    {
+        tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+        tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
+        if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
+    }
+
+
+    // Face 0
+    {
+        vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
+        plane pl0(tetBTet.tri(0).a(), n, false);
 
-    // face0
-    plane pl0(tetB[1], tetB[3], tetB[2]);
         tetA.tet().sliceWithPlane(pl0, cutInside, outside);
         if (nCutInside == 0)
         {
             return;
         }
+    }
+
+    // Face 1
+    {
+        vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
+        plane pl1(tetBTet.tri(1).a(), n, false);
 
-    // face1
-    plane pl1(tetB[0], tetB[2], tetB[3]);
         nInside = 0;
         for (label i = 0; i < nCutInside; i++)
         {
@@ -70,9 +101,13 @@ void Foam::tetOverlapVolume::tetTetOverlap
         {
             return;
         }
+    }
+
+    // Face 2
+    {
+        vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
+        plane pl2(tetBTet.tri(2).a(), n, false);
 
-    // face2
-    plane pl2(tetB[0], tetB[3], tetB[1]);
         nCutInside = 0;
         for (label i = 0; i < nInside; i++)
         {
@@ -83,14 +118,18 @@ void Foam::tetOverlapVolume::tetTetOverlap
         {
             return;
         }
+    }
 
-    // face3
+    // Face 3
-    plane pl3(tetB[0], tetB[1], tetB[2]);
+    {
+        vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
+        plane pl3(tetBTet.tri(3).a(), n, false);
         for (label i = 0; i < nCutInside; i++)
         {
             const tetPointRef t = cutInsideTets[i].tet();
             t.sliceWithPlane(pl3, insideOp, outside);
         }
+    }
 }