From e27a36bec41694408eedb47596c54d326f69b334 Mon Sep 17 00:00:00 2001 From: vbd Date: Fri, 24 Aug 2007 09:12:17 +0000 Subject: [PATCH] staffan : * made some changes testing floating points comparisons --- src/INTERP_KERNEL/Interpolation3D.cxx | 5 +- src/INTERP_KERNEL/TransformedTriangle.cxx | 5 +- src/INTERP_KERNEL/TransformedTriangle.hxx | 6 +++ .../TransformedTriangle_intersect.cxx | 49 +++++++++++++++---- src/INTERP_KERNEL/VectorUtils.hxx | 12 +++-- 5 files changed, 60 insertions(+), 17 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index faf19c113..b3ef55ad7 100644 --- 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 49e7e6013..fd135e74d 100644 --- 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 55e8c5aec..63a95813d 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -3,6 +3,12 @@ #include +#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 35afa8721..49c152e24 100644 --- 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 a85525de0..aa434e3b1 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -7,6 +7,10 @@ #include #include +#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; -- 2.39.2