Added small utility memcmp_double()

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

opm/core/grid.h

@@ -302,6 +302,7 @@ read_grid(const char *fname);
 
 
 
+
 bool
 grid_equal(const struct UnstructuredGrid * grid1 , const struct UnstructuredGrid * grid2);
 

opm/core/grid/grid.c

@@ -19,6 +19,7 @@
 
 #include "config.h"
 #include <opm/core/grid.h>
+#include <opm/core/utility/memcmp_double.h>
 
 #include <assert.h>
 #include <errno.h>
@@ -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 

opm/core/utility/memcmp_double.c

@@ -0,0 +1,30 @@
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+/**
+   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;
+    }
+}

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

tests/test_ug.cpp

@@ -19,6 +19,7 @@
 #include <algorithm>
 #include <vector>
 #include <opm/core/grid.h>
+#include <opm/core/utility/memcmp_double.h>
 #include <opm/core/grid/cornerpoint_grid.h>  /* compute_geometry */
 #include <opm/core/grid/GridManager.hpp>  /* compute_geometry */
 #include <opm/core/grid/cpgpreprocess/preprocess.h>
@@ -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 );
+    }
+}