single precision changes; corrected intersection with tolerance

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -222,16 +222,16 @@ public:
         ) const;
 
         //- Fast intersection with a ray.
-        //  For a hit, the pointHit.distance() is the line parameter t :
-        //  intersection=p+t*q. Only defined for FULL_RAY or
-        //  HALF_RAY.
+        //  Does face-center decomposition and returns triangle intersection
+        //  point closest to p. See triangle::intersection for details.
         pointHit intersection
         (
             const point& p,
             const vector& q,
             const point& ctr,
             const pointField& meshPoints,
-            const intersection::algorithm alg
+            const intersection::algorithm alg,
+            const scalar tol = 0.0
         ) const;
 
         //- Return nearest point to face

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -136,7 +136,8 @@ Foam::pointHit Foam::face::intersection
     const vector& q,
     const point& ctr,
     const pointField& meshPoints,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     scalar nearestHitDist = VGREAT;
@@ -154,7 +155,7 @@ Foam::pointHit Foam::face::intersection
             meshPoints[f[pI]],
             meshPoints[f[fcIndex(pI)]],
             ctr
-        ).intersection(p, q, alg);
+        ).intersection(p, q, alg, tol);
 
         if (curHit.hit())
         {

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -153,7 +153,7 @@ public:
             inline scalar sweptVol(const triangle& t) const;
 
             //- Return point intersection with a ray.
-            //  For a hit, the distance is signed.  Positive number
+            //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.
             //  In case of miss pointHit is set to nearest point
             //  on triangle and its distance to the distance between
@@ -166,12 +166,14 @@ public:
                 const intersection::direction dir = intersection::VECTOR
             ) const;
 
-            //- Fast intersection with a ray.
-            //  For a hit, the pointHit.distance() is the line parameter t :
-            //  intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or
+            //- Fast intersection with a ray. Distance is normal distance
+            //  to the intersection.
+            //  For a hit, the distance is signed. Positive number
+            //  represents the point in front of triangle. In case of miss
+            //  pointHit position is undefined.
+            //  Only defined for VISIBLE, FULL_RAY or
             //  HALF_RAY. tol increases the virtual size of the triangle
-            //  by a relative factor - can be used to compensate for
-            //  limited precision.
+            //  by a relative factor.
             inline pointHit intersection
             (
                 const point& p,

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -585,7 +585,7 @@ inline bool triangle<Point, PointRef>::classify
 
     bool hit = false;
 
-    if (Foam::mag(u1) < SMALL)
+    if (Foam::mag(u1) < ROOTVSMALL)
     {
         beta = u0/u2;
 

src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C

@@ -73,7 +73,7 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
     const labelList& cFaces = mesh.cells()[this->celli_];
 
     point intersection(vector::zero);
-    scalar trackFraction = VGREAT;
+    scalar distanceSqr = VGREAT;
     label hitFacei = -1;
 
     const vector vec = endPosition-this->position_;
@@ -93,32 +93,21 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
                 intersection::HALF_RAY
             );
 
-            if (inter.hit() && inter.distance() < trackFraction)
+            if (inter.hit())
             {
-                trackFraction = inter.distance();
-                hitFacei = facei;
-                intersection = inter.hitPoint();
+                scalar s = magSqr(inter.hitPoint()-this->position_);
+
+                if (s < distanceSqr)
+                {
+                    distanceSqr = s;
+                    hitFacei = facei;
+                    intersection = inter.hitPoint();
+                }
             }
         }
     }
 
-
-    if (hitFacei != -1)
-    {
-        if (trackFraction > 1.0)
-        {
-            // Nearest intersection beyond endPosition so we hit endPosition.
-            this->position_ = endPosition;
-            this->facei_ = -1;
-            return 1.0;
-        }
-        else
-        {
-            this->position_ = intersection;
-            this->facei_ = hitFacei;
-        }
-    }
-    else
+    if (hitFacei == -1)
     {
         // Did not find any intersection. Fall back to original approximate
         // algorithm
@@ -130,7 +119,24 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
     }
 
 
-    // Normal situation (trackFraction 0..1)
+    scalar trackFraction = Foam::sqrt(distanceSqr/magSqr(vec));
+
+    if (trackFraction >= (1.0-SMALL))
+    {
+        // Nearest intersection beyond endPosition so we hit endPosition.
+        this->position_ = endPosition;
+        this->facei_ = -1;
+        return 1.0;
+    }
+    else
+    {
+        this->position_ = intersection;
+        this->facei_ = hitFacei;
+    }
+
+
+    // Normal situation (trackFraction 0..1). Straight copy
+    // of Particle::trackToFace.
 
     bool internalFace = this->cloud().internalFace(this->facei_);
 

src/meshTools/indexedOctree/indexedOctree.H

@@ -128,6 +128,15 @@ public:
 
 private:
 
+    // Static data
+
+        //- Relative peturbation tolerance. Determines when point is
+        //  considered to be close to face/edge of bb of node.
+        //  The tolerance is relative to the bounding box of the smallest
+        //  node.
+        static scalar perturbTol_;
+
+
     // Private data
 
         //- Underlying shapes for geometric queries.
@@ -144,8 +153,9 @@ private:
 
     // Private Member Functions
 
-        //- Like above but now bb is implicitly provided as parent bb + mid
-        //  + octant
+        //- Helper: does bb intersect a sphere around sample? Or is any
+        //  corner point of bb closer than nearestDistSqr to sample.
+        //  (bb is implicitly provided as parent bb + octant)
         static bool overlaps
         (
             const treeBoundBox& parentBb,
@@ -215,6 +225,60 @@ private:
                 point& nearestPoint
             ) const;
 
+            //- Return bbox of octant
+            treeBoundBox subBbox
+            (
+                const label parentNodeI,
+                const direction octant
+            ) const;
+
+            //- Helper: take a point on/close to face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take a point on face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const direction,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take point on face(s) of bb and push it away from
+            //  edges of face.
+            static point pushPointIntoFace
+            (
+                const treeBoundBox& bb,
+                const vector& dir,          // end-start
+                const point& pt
+            );
+
+            ////- Push point on multiple faces away from any corner/edge.
+            //static void checkMultipleFaces
+            //(
+            //    const treeBoundBox& bb,
+            //    const vector& dir,          // end-start
+            //    pointIndexHit& faceHitInfo,
+            //    direction& faceID
+            //);
+
+            //- Walk to parent of node+octant.
+            bool walkToParent
+            (
+                const label nodeI,
+                const direction octant,
+
+                label& parentNodeI,
+                label& parentOctant
+            ) const;
+
             //- Walk tree to neighbouring node. Return false if edge of tree
             //  hit.
             bool walkToNeighbour
@@ -225,8 +289,8 @@ private:
                 direction& octant
             ) const;
 
-            //- Categorize point on bounding box.
-            static direction getFace(const treeBoundBox&, const point&);
+            //- Debug: return verbose the bounding box faces
+            static word faceString(const direction faceID);
 
             //- Traverse a node. If intersects a triangle return first
             // intersection point.
@@ -235,9 +299,11 @@ private:
             void traverseNode
             (
                 const bool findAny,
+                const point& treeStart,
+                const vector& treeVec,
+
                 const point& start,
                 const point& end,
-                const vector& dir,
                 const label nodeI,
                 const direction octantI,
 
@@ -252,7 +318,8 @@ private:
                 const point& treeStart,
                 const point& treeEnd,
                 const label startNodeI,
-                const direction startOctantI
+                const direction startOctantI,
+                const bool verbose = false
             ) const;
 
             //- Find any or nearest intersection of line between start and end.
@@ -310,6 +377,10 @@ private:
 
 public:
 
+        //- Get the perturbation tolerance
+        static scalar& perturbTol();
+
+
     // Constructors
 
         //- Construct null
@@ -505,6 +576,7 @@ public:
                 const point& sample
             );
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or

src/meshTools/indexedOctree/treeDataFace.C

@@ -528,7 +528,7 @@ bool Foam::treeDataFace::intersects
         intersection::HALF_RAY
     );
 
-    if (inter.hit() && inter.distance() <= 1)
+    if (inter.hit() && magSqr(inter.hitPoint()-start) <= magSqr(dir))
     {
         // Note: no extra test on whether intersection is in front of us
         // since using half_ray

src/meshTools/indexedOctree/treeDataTriSurface.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "treeDataTriSurface.H"
@@ -226,7 +224,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
 
     if (!pHit.hit())
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "treeBb:" << treeBb
             << " sample:" << sample
             << " pHit:" << pHit
@@ -238,7 +236,8 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
         surface_,
         sample,
         pHit.index(),
-        pHit.hitPoint()
+        pHit.hitPoint(),
+        indexedOctree<treeDataTriSurface>::perturbTol()
     );
 
     if (t == triSurfaceTools::UNKNOWN)
@@ -255,12 +254,13 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
     }
     else
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "problem" << abort(FatalError);
         return indexedOctree<treeDataTriSurface>::UNKNOWN;
     }
 }
 
+
 // Check if any point on triangle is inside cubeBb.
 bool Foam::treeDataTriSurface::overlaps
 (
@@ -305,6 +305,7 @@ bool Foam::treeDataTriSurface::overlaps
     // Now we have the difficult case: all points are outside but connecting
     // edges might go through cube. Use fast intersection of bounding box.
 
+    //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
     return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
 }
 
@@ -445,9 +446,17 @@ bool Foam::treeDataTriSurface::intersects
 
     const vector dir(end - start);
 
-    pointHit inter = tri.intersection(start, dir, intersection::HALF_RAY);
+    // Use relative tolerance (from octree) to determine intersection.
 
-    if (inter.hit() && inter.distance() <= 1)
+    pointHit inter = tri.intersection
+    (
+        start,
+        dir,
+        intersection::HALF_RAY,
+        indexedOctree<treeDataTriSurface>::perturbTol()
+    );
+
+    if (inter.hit() && magSqr(inter.hitPoint()-start) <= magSqr(dir))
     {
         // Note: no extra test on whether intersection is in front of us
         // since using half_ray.
@@ -459,6 +468,16 @@ bool Foam::treeDataTriSurface::intersects
     {
         return false;
     }
+
+
+    //- Using exact intersections
+    //pointHit info = f.tri(points).intersectionExact(start, end);
+    //
+    //if (info.hit())
+    //{
+    //    intersectionPoint = info.hitPoint();
+    //}
+    //return info.hit();
 }
 
 

src/meshTools/octree/treeBoundBox.C

@@ -309,12 +309,14 @@ bool Foam::treeBoundBox::overlaps
 //     This makes sure that posBits tests 'inside'
 bool Foam::treeBoundBox::intersects
 (
+    const point& overallStart,
+    const vector& overallVec,
     const point& start,
     const point& end,
-    point& pt
+    point& pt,
+    direction& ptOnFaces
 ) const
 {
-    const vector vec(end - start);
     const direction endBits = posBits(end);
     pt = start;
 
@@ -325,80 +327,106 @@ bool Foam::treeBoundBox::intersects
         if (ptBits == 0)
         {
             // pt inside bb
+            ptOnFaces = faceBits(pt);
             return true;
         }
 
         if ((ptBits & endBits) != 0)
         {
             // pt and end in same block outside of bb
+            ptOnFaces = faceBits(pt);
             return false;
         }
 
         if (ptBits & LEFTBIT)
         {
             // Intersect with plane V=min, n=-1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
+            {
+                scalar s = (min().x() - overallStart.x())/overallVec.x();
+                pt.x() = min().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (min().x() - pt.x())/vec.x();
+                // Vector not in x-direction. But still intersecting bb planes.
+                // So must be close - just snap to bb.
                 pt.x() = min().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-        if (ptBits & RIGHTBIT)
+        else if (ptBits & RIGHTBIT)
         {
             // Intersect with plane V=max, n=1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
+            {
+                scalar s = (max().x() - overallStart.x())/overallVec.x();
+                pt.x() = max().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().x() - pt.x())/vec.x();
                 pt.x() = max().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BOTTOMBIT)
+        else if (ptBits & BOTTOMBIT)
         {
             // Intersect with plane V=min, n=0,-1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
             {
-                scalar s = (min().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
+                scalar s = (min().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = min().y();
-                pt.z() = pt.z() + vec.z()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
+            {
+                pt.x() = min().y();
             }
         }
-        if (ptBits & TOPBIT)
+        else if (ptBits & TOPBIT)
         {
             // Intersect with plane V=max, n=0,1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
+            {
+                scalar s = (max().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = max().y();
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
                 pt.y() = max().y();
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BACKBIT)
+        else if (ptBits & BACKBIT)
         {
             // Intersect with plane V=min, n=0,0,-1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (min().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = min().z();
+            }
+            else
             {
-                scalar s = (min().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = min().z();
             }
         }
-        if (ptBits & FRONTBIT)
+        else if (ptBits & FRONTBIT)
         {
             // Intersect with plane V=max, n=0,0,1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (max().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = max().z();
+            }
+            else
             {
-                scalar s = (max().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = max().z();
             }
         }
@@ -406,6 +434,18 @@ bool Foam::treeBoundBox::intersects
 }
 
 
+bool Foam::treeBoundBox::intersects
+(
+    const point& start,
+    const point& end,
+    point& pt
+) const
+{
+    direction ptBits;
+    return intersects(start, end-start, start, end, pt, ptBits);
+}
+
+
 // this.bb fully contains bb
 bool Foam::treeBoundBox::contains(const treeBoundBox& bb) const
 {
@@ -452,6 +492,40 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 }
 
 
+// Code position of pt on bounding box faces
+Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
+{
+    direction faceBits = 0;
+    if (pt.x() == min().x())
+    {
+        faceBits |= LEFTBIT;
+    }
+    else if (pt.x() == max().x())
+    {
+        faceBits |= RIGHTBIT;
+    }
+
+    if (pt.y() == min().y())
+    {
+        faceBits |= BOTTOMBIT;
+    }
+    else if (pt.y() == max().y())
+    {
+        faceBits |= TOPBIT;
+    }
+
+    if (pt.z() == min().z())
+    {
+        faceBits |= BACKBIT;
+    }
+    else if (pt.z() == max().z())
+    {
+        faceBits |= FRONTBIT;
+    }
+    return faceBits;
+}
+
+
 // Code position of point relative to box
 Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
 {
@@ -461,7 +535,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {