From 09b8e8888711d14b45de003e53f4258e1aa5b593 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 22 Sep 2011 11:30:24 +0100
Subject: [PATCH] Added stabilisation in division by zero; optimising compiler

---
 .../algorithms/octree/octree/treeBoundBox.C   | 46 +++++++++++++++----
 1 file changed, 38 insertions(+), 8 deletions(-)

diff --git a/src/OpenFOAM/algorithms/octree/octree/treeBoundBox.C b/src/OpenFOAM/algorithms/octree/octree/treeBoundBox.C
index 1229403dd..2675d3401 100644
--- a/src/OpenFOAM/algorithms/octree/octree/treeBoundBox.C
+++ b/src/OpenFOAM/algorithms/octree/octree/treeBoundBox.C
@@ -318,8 +318,21 @@ bool Foam::treeBoundBox::intersects
 ) const
 {
     const direction endBits = posBits(end);
+
     pt = start;
 
+    // Note
+    // Optimising compilers (Intel 12) reorganises the code in
+    // a way which executes all divisions before doing the if-check.
+    // This causes a floating point exception of division by zero.
+    // Furthermore, division by VSMALL (first number larger than zero) may
+    // in some cases cause floating point exception overflow, and stabilisation
+    // is required.  Added by HJ, 21/Sep/2011
+
+    scalar s;
+
+    static const scalar smallNumber = SMALL;
+
     while (true)
     {
         direction ptBits = posBits(pt);
@@ -328,6 +341,7 @@ bool Foam::treeBoundBox::intersects
         {
             // pt inside bb
             ptOnFaces = faceBits(pt);
+
             return true;
         }
 
@@ -335,6 +349,7 @@ bool Foam::treeBoundBox::intersects
         {
             // pt and end in same block outside of bb
             ptOnFaces = faceBits(pt);
+
             return false;
         }
 
@@ -343,7 +358,9 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=min, n=-1,0,0
             if (Foam::mag(overallVec.x()) > VSMALL)
             {
-                scalar s = (min().x() - overallStart.x())/overallVec.x();
+                s = (min().x() - overallStart.x())/
+                    stabilise(overallVec.x(), smallNumber);
+
                 pt.x() = min().x();
                 pt.y() = overallStart.y() + overallVec.y()*s;
                 pt.z() = overallStart.z() + overallVec.z()*s;
@@ -360,7 +377,9 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=max, n=1,0,0
             if (Foam::mag(overallVec.x()) > VSMALL)
             {
-                scalar s = (max().x() - overallStart.x())/overallVec.x();
+                s = (max().x() - overallStart.x())/
+                    stabilise(overallVec.x(), smallNumber);
+
                 pt.x() = max().x();
                 pt.y() = overallStart.y() + overallVec.y()*s;
                 pt.z() = overallStart.z() + overallVec.z()*s;
@@ -375,7 +394,9 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=min, n=0,-1,0
             if (Foam::mag(overallVec.y()) > VSMALL)
             {
-                scalar s = (min().y() - overallStart.y())/overallVec.y();
+                s = (min().y() - overallStart.y())/
+                    stabilise(overallVec.y(), smallNumber);
+
                 pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = min().y();
                 pt.z() = overallStart.z() + overallVec.z()*s;
@@ -390,7 +411,9 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=max, n=0,1,0
             if (Foam::mag(overallVec.y()) > VSMALL)
             {
-                scalar s = (max().y() - overallStart.y())/overallVec.y();
+                s = (max().y() - overallStart.y())/
+                    stabilise(overallVec.y(), smallNumber);
+
                 pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = max().y();
                 pt.z() = overallStart.z() + overallVec.z()*s;
@@ -405,7 +428,8 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=min, n=0,0,-1
             if (Foam::mag(overallVec.z()) > VSMALL)
             {
-                scalar s = (min().z() - overallStart.z())/overallVec.z();
+                s = (min().z() - overallStart.z())/
+                    stabilise(overallVec.z(), smallNumber);
                 pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = overallStart.y() + overallVec.y()*s;
                 pt.z() = min().z();
@@ -420,7 +444,9 @@ bool Foam::treeBoundBox::intersects
             // Intersect with plane V=max, n=0,0,1
             if (Foam::mag(overallVec.z()) > VSMALL)
             {
-                scalar s = (max().z() - overallStart.z())/overallVec.z();
+                s = (max().z() - overallStart.z())/
+                    stabilise(overallVec.z(), smallNumber);
+
                 pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = overallStart.y() + overallVec.y()*s;
                 pt.z() = max().z();
@@ -442,7 +468,8 @@ bool Foam::treeBoundBox::intersects
 ) const
 {
     direction ptBits;
-    return intersects(start, end-start, start, end, pt, ptBits);
+
+    return intersects(start, end - start, start, end, pt, ptBits);
 }
 
 
@@ -458,7 +485,7 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
     //
     // Compare all components against min and max of bb
     //
-    for (direction cmpt=0; cmpt<3; cmpt++)
+    for (direction cmpt = 0; cmpt < 3; cmpt++)
     {
         if (pt[cmpt] < min()[cmpt])
         {
@@ -496,6 +523,7 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
 {
     direction faceBits = 0;
+
     if (pt.x() == min().x())
     {
         faceBits |= LEFTBIT;
@@ -522,6 +550,7 @@ Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
     {
         faceBits |= FRONTBIT;
     }
+
     return faceBits;
 }
 
@@ -557,6 +586,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= FRONTBIT;
     }
+
     return posBits;
 }