diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 754bef2c..18dd617f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -130,6 +130,7 @@ list (APPEND MAIN_SOURCE_FILES opm/core/utility/VelocityInterpolation.cpp opm/core/utility/WachspressCoord.cpp opm/core/utility/miscUtilities.cpp + opm/core/utility/memcmp_double.c opm/core/utility/miscUtilitiesBlackoil.cpp opm/core/utility/NullStream.cpp opm/core/utility/parameters/Parameter.cpp @@ -425,6 +426,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/utility/Exceptions.hpp opm/core/utility/Factory.hpp opm/core/utility/MonotCubicInterpolator.hpp + opm/core/utility/memcmp_double.h opm/core/utility/NonuniformTableLinear.hpp opm/core/utility/NullStream.hpp opm/core/utility/RegionMapping.hpp diff --git a/opm/core/grid.h b/opm/core/grid.h index 36568045..59d50b23 100644 --- a/opm/core/grid.h +++ b/opm/core/grid.h @@ -302,7 +302,8 @@ read_grid(const char *fname); -bool + +bool grid_equal(const struct UnstructuredGrid * grid1 , const struct UnstructuredGrid * grid2); #ifdef __cplusplus diff --git a/opm/core/grid/grid.c b/opm/core/grid/grid.c index d07017f9..02a956db 100644 --- a/opm/core/grid/grid.c +++ b/opm/core/grid/grid.c @@ -19,6 +19,7 @@ #include "config.h" #include +#include #include #include @@ -542,31 +543,6 @@ read_grid(const char *fname) return G; } -/** - Ahhh - the joys of comparing floating point numbers .... - - http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm -*/ -static int memcmp_double(const double * p1 , const double *p2 , size_t num_elements) { - if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0) - return 0; - else { - const double abs_epsilon = 1e-8; - const double rel_epsilon = 1e-5; - - for (size_t index = 0; index < num_elements; index++) { - double diff = fabs(p1[index] - p2[index]); - if (diff > abs_epsilon) { - double sum = fabs(p1[index]) + fabs(p2[index]); - - if (diff > sum * rel_epsilon) - return 1; - } - } - - return 0; - } -} bool diff --git a/opm/core/utility/memcmp_double.c b/opm/core/utility/memcmp_double.c new file mode 100644 index 00000000..448895ab --- /dev/null +++ b/opm/core/utility/memcmp_double.c @@ -0,0 +1,30 @@ +#include +#include +#include + +/** + Ahhh - the joys of comparing floating point numbers .... + + http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +*/ + +int memcmp_double(const double * p1 , const double *p2 , size_t num_elements) { + if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0) + return 0; + else { + const double abs_epsilon = 1e-8; + const double rel_epsilon = 1e-5; + + for (size_t index = 0; index < num_elements; index++) { + double diff = fabs(p1[index] - p2[index]); + if (diff > abs_epsilon) { + double sum = fabs(p1[index]) + fabs(p2[index]); + + if (diff > sum * rel_epsilon) + return 1; + } + } + + return 0; + } +} diff --git a/opm/core/utility/memcmp_double.h b/opm/core/utility/memcmp_double.h new file mode 100644 index 00000000..86e140ed --- /dev/null +++ b/opm/core/utility/memcmp_double.h @@ -0,0 +1,22 @@ +/* + Will compare all the elements in the two double pointers p1 and + p2. If all elements are 'sufficiently' the function will return 0, + otherwise it will return 1. +*/ + +#ifndef MEMCMP_DOUBLE_H +#define MEMCMP_DOUBLE_H + +#ifdef __cplusplus +extern "C" { +#endif + + +int memcmp_double(const double * p1 , const double *p2 , size_t num_elements); + + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/tests/test_ug.cpp b/tests/test_ug.cpp index f374cab7..73a2f154 100644 --- a/tests/test_ug.cpp +++ b/tests/test_ug.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include /* compute_geometry */ #include /* compute_geometry */ #include @@ -164,3 +165,66 @@ BOOST_AUTO_TEST_CASE(TOPS_Fully_Specified) { BOOST_CHECK( grid_equal( cgrid1 , cgrid2 )); } + + + +BOOST_AUTO_TEST_CASE(compare_double) { + const double abs_epsilon = 1e-8; + const double rel_epsilon = 1e-5; + + double v1,v2; + /* Should be equal: */ + { + v1 = 0.0; + v2 = 0.0; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + v1 = 1e-12; + v2 = v1 + 0.5*abs_epsilon; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + v1 = 7.0; + v2 = 7.0; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + v1 = -7.0; + v2 = -7.0; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + v1 = 0; + v2 = 0.5 * abs_epsilon; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + + v1 = 1e7; + v2 = 1e7 + 2*abs_epsilon; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + + v1 = 1e7; + v2 = 1e7*(1 + 2*rel_epsilon); + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 0 ); + } + + /* Should be different: */ + { + v1 = 0; + v2 = 1.5 * abs_epsilon; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 1 ); + + v1 = 1e-8; + v2 = v1 + 1.5*abs_epsilon; + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 1 ); + + v1 = 1; + v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 1 ); + + v1 = 10; + v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 1 ); + + v1 = 1e7; + v2 = 1e7*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK_EQUAL( memcmp_double( &v1 , &v2 , 1 ) , 1 ); + } +}