[MOAB-dev] r3545 - MOAB/trunk

[bmsmith6 at wisc.edu]

Author: bmsmith
Date: 2010-02-15 11:44:27 -0600 (Mon, 15 Feb 2010)
New Revision: 3545

Modified:
   MOAB/trunk/MBOrientedBox.cpp
   MOAB/trunk/MBOrientedBox.hpp
Log:
Bug fix for MBOrientedBox::intersect_ray. 
  A squared distance (x^2) is larger or smaller than x depending if x<1 or x>1.
Changed tolerances to be more consistent in MBOrientedBox::intersect_ray.

This fix decreased the number of lost particles for ITER Module13 model by ~5x.


Modified: MOAB/trunk/MBOrientedBox.cpp
===================================================================
--- MOAB/trunk/MBOrientedBox.cpp	2010-02-14 04:56:42 UTC (rev 3544)
+++ MOAB/trunk/MBOrientedBox.cpp	2010-02-15 17:44:27 UTC (rev 3545)
@@ -477,16 +477,24 @@
   const MBCartVect cx = center - b;
   double dist_s = cx % m;
   double dist_sq = cx % cx - (dist_s*dist_s);
-  double max_diagsq = outer_radius_squared();
+  double max_diagsq = outer_radius_squared(reps);
   
     // if greater than the longest diagonal, we don't hit
-  if (dist_sq > max_diagsq+reps)
+  if (dist_sq > max_diagsq)
     return false;
-  if (len && dist_s - max_diagsq > *len)
-    return false;
-  
+
+  // If the closest possible hit is farther than len, we don't want the hit.
+  // Problem: the previous method was wrong because max_diagsq will be greater
+  // than max_diag if max_diag>1 but less than max_diag if max_diag<1.
+  // Be careful with absolute value, squaring distances, and subtracting squared
+  // distances.
+  if (len) {
+    const double temp = fabs(dist_s) - *len;
+    if(0.0<temp && temp*temp>max_diagsq) return false;
+  } 
+
     // if smaller than shortest diagonal, we do hit
-  if (dist_sq < inner_radius_squared() - reps && dist_s >= 0.0)
+  if (dist_sq < inner_radius_squared(reps) && dist_s >= 0.0)
     return true;
     
     // get transpose of axes

Modified: MOAB/trunk/MBOrientedBox.hpp
===================================================================
--- MOAB/trunk/MBOrientedBox.hpp	2010-02-14 04:56:42 UTC (rev 3544)
+++ MOAB/trunk/MBOrientedBox.hpp	2010-02-15 17:44:27 UTC (rev 3545)
@@ -53,8 +53,8 @@
 
   inline double inner_radius() const; //!< radius of inscribed sphere
   inline double outer_radius() const; //!< radius of circumscribed sphere
-  inline double outer_radius_squared() const; //!< square of radius of circumsphere
-  inline double inner_radius_squared() const; //!< square of radius if inscribed sphere
+  inline double outer_radius_squared(const double reps) const; //!< square of (radius+at least epsilon) of circumsphere
+  inline double inner_radius_squared(const double reps) const; //!< square of (radius-epsilon) of inscribed sphere
   inline double volume() const;               //!< volume of box
   inline MBCartVect dimensions() const;       //!< number of dimensions for which box is not flat
   inline double area() const;                 //!< largest side area
@@ -176,24 +176,30 @@
 #endif
 }
 
-double MBOrientedBox::outer_radius_squared() const
+// Add at least epsilon to the radius, before squaring it.
+double MBOrientedBox::outer_radius_squared(const double reps) const
 {
 #if MB_ORIENTED_BOX_OUTER_RADIUS
-  return radius * radius;
+  return (radius+reps)*(radius+reps);
 #elif MB_ORIENTED_BOX_UNIT_VECTORS
-  return length % length;
+  MBCartVect tmp(length[0]+reps,length[1]+reps,length[2]+reps);
+  return tmp % tmp;
 #else
-  const MBCartVect half_diag = axis[0] + axis[1] + axis[2];
+  MBCartVect half_diag = axis[0] + axis[1] + axis[2];
+  half_diag += MBCartVect(reps,reps,reps);
   return half_diag % half_diag;
 #endif
 }
 
-double MBOrientedBox::inner_radius_squared() const
+// Subtract epsilon from the length of the shortest axis, before squaring it.
+double MBOrientedBox::inner_radius_squared(const double reps) const
 {
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-  return length[0] * length[0];
+  return (length[0]-reps) * (length[0]-reps);
 #else
-  return axis[0] % axis[0];
+  MBCartVect tmp = axis[0];
+  tmp -= MBCartVect(reps,reps,reps);
+  return (tmp % tmp);
 #endif
 }