staffan :

	* made some changes testing floating points comparisons

diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx
index faf19c11385432a84adab88e06e98f5f8ed747e2..b3ef55ad73c82ea3da5463c1d087f701695fae51 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation3D.cxx
+++ b/src/INTERP_KERNEL/Interpolation3D.cxx
@@ -173,8 +173,7 @@ namespace MEDMEM
              {
                const int targetIdx = (*iter)->getIndex();
                const double vol = intersector->intersectCells(srcIdx, targetIdx);
-               if(!epsilonEqual(vol, 0.0, 1.0e-10))
-                 //  if(vol != 0.0)
+               if(!epsilonEqual(vol, 0.0))
                  {
                    volumes->insert(make_pair(targetIdx, vol));
                    LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]);
@@ -301,7 +300,7 @@ namespace MEDMEM
            const int srcIdx = *iter + 1;
            const double vol = intersector->intersectCells(srcIdx, targetIdx);
 
-           if(!epsilonEqual(vol, 0.0, 1.0e-10))
+           if(!epsilonEqual(vol, 0.0))
              {
                matrix[srcIdx - 1].insert(make_pair(targetIdx, vol));
              }

diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx
index 49e7e6013621507e85a5631ef2b5df8816e7e068..fd135e74d89de666d71a591a9454af9e9768882e 100644 (file)
--- a/src/INTERP_KERNEL/TransformedTriangle.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle.cxx
@@ -596,8 +596,11 @@ namespace INTERP_UTILS
     for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
       {
        // ? should have epsilon-equality here?
-       //if(!epsilonEqual(_coords[5*c + coord], 0.0, 1.0e-15))
+#ifdef EPS_TESTING
+       if(!epsilonEqualRelative(_coords[5*c + coord], 0.0, TEST_EPS * _coords[5*c + coord]))
+#else
        if(_coords[5*c + coord] != 0.0)
+#endif
          {
            return false;
          }

diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx
index 55e8c5aecb176d2f2d9b744cde02b985f7538eb5..63a95813dbbeed87ff307c20c9997490786840b3 100644 (file)
--- a/src/INTERP_KERNEL/TransformedTriangle.hxx
+++ b/src/INTERP_KERNEL/TransformedTriangle.hxx
@@ -3,6 +3,12 @@
 
 #include <vector>
 
+#undef EPS_TESTING  // does not give good results
+
+#ifdef EPS_TESTING
+#define TEST_EPS 1.0e-14
+#endif
+
 // Levels : 
 // 1 - overview of algorithm + volume result
 // 2 - algorithm detail

diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
index 35afa87212152c1905f79a66ced2ef2a8089e2dd..49c152e24767ea36a2d5c9e7068e5728edda04ee 100644 (file)
--- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
@@ -6,7 +6,10 @@
 #include "VectorUtils.hxx"
 
 #define SEG_RAY_TABLE 1 // seems correct
-#undef TEST_EPS 0.0//1.0e-14
+
+
+
+
 
 namespace INTERP_UTILS
 {
@@ -200,8 +203,12 @@ namespace INTERP_UTILS
        OZX, XYZ  // ZX
       };
 
-        if(calcStableC(seg,DoubleProduct( edge )) != 0.0)
-    //    if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS))
+#ifdef EPS_TESTING
+    if(!epsilonEqual(calcStableC(seg,DoubleProduct( edge )), 0.0, TEST_EPS))
+
+#else
+      if(calcStableC(seg,DoubleProduct( edge )) != 0.0)
+#endif
       {
        return false;
       } 
@@ -356,8 +363,11 @@ namespace INTERP_UTILS
        const TetraEdge edge = EDGES_FOR_CORNER[3*corner + i];
        const DoubleProduct dp = DoubleProduct( edge );
        const double c = calcStableC(seg, dp);
+#ifdef EPS_TESTING
+       if(!epsilonEqual(c, 0.0, TEST_EPS))
+#else
        if(c != 0.0)
-       //      if(!epsilonEqual(c, 0.0, TEST_EPS))
+#endif
          {
            return false;
          }
@@ -510,7 +520,7 @@ namespace INTERP_UTILS
        assert(pt[i] >= 0.0);
        assert(pt[i] <= 1.0);
       }
-    assert(epsilonEqual(pt[0] + pt[1] + pt[2] - 1.0, 0.0, 1.0e-9));
+    assert(epsilonEqualRelative(pt[0] + pt[1] + pt[2], 1.0));
   }
     
   /**
@@ -565,7 +575,11 @@ namespace INTERP_UTILS
     
     //? epsilon-equality here?
     // cond. 1
+#ifdef EPS_TESTING
+    if(epsilonEqual(cVal0, 0.0, TEST_EPS))
+#else
     if(cVal0 != 0.0) 
+#endif
       {
        LOG(4, "SR fails at cond 1 cVal0 = "  << cVal0 );
        return false;
@@ -697,8 +711,11 @@ namespace INTERP_UTILS
 
     // if two or more c-values are zero we disallow x-edge intersection
     // Grandy, p.446
-        const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
-    //?    const int numZeros = (epsilonEqual(cPQ,0.0) ? 1 : 0) + (epsilonEqual(cQR,0.0) ? 1 : 0) + (epsilonEqual(cRP, 0.0) ? 1 : 0);
+#ifdef EPS_TESTING
+    const int numZeros = (epsilonEqual(cPQ,0.0) ? 1 : 0) + (epsilonEqual(cQR,0.0) ? 1 : 0) + (epsilonEqual(cRP, 0.0) ? 1 : 0);
+#else
+    const int numZeros = (cPQ == 0.0 ? 1 : 0) + (cQR == 0.0 ? 1 : 0) + (cRP == 0.0 ? 1 : 0);
+#endif
     
     if(numZeros >= 2 ) 
       {
@@ -738,7 +755,11 @@ namespace INTERP_UTILS
 
     //? should equality with zero use epsilon?
     LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 );
+#ifdef EPS_TESTING
+    return (t1*t2 <= 0.0) && (!epsilonEqualRelative(t1, t2, TEST_EPS * std::max(t1, t2)));
+#else
     return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
+#endif
   }
 
   /**
@@ -763,12 +784,13 @@ namespace INTERP_UTILS
     const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]);
     const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]);
 
-    if(false)
-    //if(epsilonEqual(c1, 0.0, TEST_EPS) || epsilonEqual(c2, 0.0, TEST_EPS) || epsilonEqual(c3, 0.0, TEST_EPS))
+#ifdef EPS_TSTING
+    if(epsilonEqual(c1, 0.0, TEST_EPS) || epsilonEqual(c2, 0.0, TEST_EPS) || epsilonEqual(c3, 0.0, TEST_EPS))
       {
        return false;
       }
     else
+#endif
       {
        return (c1*c3 > 0.0) && (c2*c3 > 0.0);
       }
@@ -792,8 +814,11 @@ namespace INTERP_UTILS
 
     //? should we use epsilon-equality here in second test?
     LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
-    //return (coord1*coord2 <= 0.0) && epsilonEqual(coord1,coord2);
+#ifdef EPS_TESTING
+    return (coord1*coord2 <= 0.0) && epsilonEqualRelative(coord1,coord2, TEST_EPS * std::max(coord1, coord2));
+#else
     return (coord1*coord2 <= 0.0) && (coord1 != coord2);
+#endif
   }
 
   bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const
@@ -803,7 +828,11 @@ namespace INTERP_UTILS
     const double coord2 = _coords[5*( (seg + 1) % 3) + 4];
     //? should we use epsilon-equality here in second test?
     LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
+#ifdef EPS_TESTING
+    return (coord1*coord2 <= 0.0) && epsilonEqualRelative(coord1,coord2, TEST_EPS * std::max(coord1, coord2));
+#else
     return (coord1*coord2 <= 0.0) && (coord1 != coord2);
+#endif
   }
 
   /**

diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx
index a85525de0656951695cd5691aa017f53641eef4d..aa434e3b1fac524dccbf5ef9f602c84a2411fbd8 100644 (file)
--- a/src/INTERP_KERNEL/VectorUtils.hxx
+++ b/src/INTERP_KERNEL/VectorUtils.hxx
@@ -7,6 +7,10 @@
 #include <numeric>
 #include <map>
 
+#define VOL_PREC 1.0e-6
+#define DEFAULT_REL_TOL 1.0e-6
+#define DEFAULT_ABS_TOL 1.0e-12
+
 namespace INTERP_UTILS
 {
  
@@ -53,15 +57,17 @@ namespace INTERP_UTILS
 
   }
 
-  inline bool epsilonEqual(const double x, const double y, const double errTol = 1.0e-12)
+  /// Should be used for comparisons to zero
+  inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
   {
     return std::abs(x - y) <= errTol;
   }
 
-  inline bool epsilonEqualRelative(const double x, const double y, const double relTol = 1.0e-6, const double absTol = 1.0e-17)
+  // Should be used for comparisons between numbers that could be large
+  inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
   {
     // necessary for comparing values close to zero
-    // in order to avoid division
+    // in order to avoid division by very small numbers
     if(std::abs(x - y) < absTol)
       {
        return true;