From 160a948abfd10893b7432b59b269f1f8ca58293e Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Mon, 6 May 2024 23:12:38 +0200
Subject: [PATCH 01/10] improve numerical stability of C++ intermediate test
 cases

- add method for calculating ULP distance
- Let the UtilityStructs use this method
- create custom matcher for comparing 2D containers
---
 .../model/GravityModelData.h                  | 49 +++++++++++++------
 src/polyhedralGravity/util/UtilityConstants.h |  1 -
 .../util/UtilityFloatArithmetic.cpp           | 27 ++++++++++
 .../util/UtilityFloatArithmetic.h             | 48 ++++++++++++++++++
 test/model/GoogleTestMatcher.h                | 27 ++++++++++
 test/model/GravityModelBigTest.cpp            |  6 +--
 test/util/UtilityFloatArithmeticTest.cpp      | 31 ++++++++++++
 7 files changed, 169 insertions(+), 20 deletions(-)
 create mode 100644 src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
 create mode 100644 src/polyhedralGravity/util/UtilityFloatArithmetic.h
 create mode 100644 test/model/GoogleTestMatcher.h
 create mode 100644 test/util/UtilityFloatArithmeticTest.cpp

diff --git a/src/polyhedralGravity/model/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h
index 0f61aaa..ecc5d2d 100644
--- a/src/polyhedralGravity/model/GravityModelData.h
+++ b/src/polyhedralGravity/model/GravityModelData.h
@@ -6,7 +6,7 @@
 #include <algorithm>
 #include <tuple>
 #include "polyhedralGravity/util/UtilityContainer.h"
-#include "polyhedralGravity/util/UtilityConstants.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
 namespace polyhedralGravity {
 
@@ -76,20 +76,23 @@ namespace polyhedralGravity {
         double s2;
 
         /**
-         * Checks two Distance structs for equality.
-         * @warning This method compares doubles! So only exact copies will evaluate to true.
+         * Checks two Distance structs for equality with another one by comparing their members values'
+         * ULP distances.
          * @param rhs the other Distance struct
          * @return true if equal
          *
          * @note Just used for testing purpose
          */
         bool operator==(const Distance &rhs) const {
-            return l1 == rhs.l1 && l2 == rhs.l2 && s1 == rhs.s1 && s2 == rhs.s2;
+            return util::floatNear(l1, rhs.l1) &&
+                   util::floatNear(l2, rhs.l2) &&
+                   util::floatNear(s1, rhs.s1) &&
+                   util::floatNear(s2, rhs.s2);
         }
 
         /**
-         * Checks two Distance structs for inequality.
-         * @warning This method compares doubles! So only exact copies will evaluate to false.
+         * Checks two Distance structs for inequality with another one by comparing their members values'
+         * ULP distances.
          * @param rhs the other Distance struct
          * @return false if unequal
          *
@@ -130,20 +133,20 @@ namespace polyhedralGravity {
         double an;
 
         /**
-         * Checks two TranscendentalExpressions for equality.
-         * @warning This method compares doubles! So only exact copies will evaluate to true.
+        * Checks two TranscendentalExpressions for equality with another one by comparing their members values'
+         * ULP distances.
          * @param rhs the other TranscendentalExpressions
          * @return true if equal
          *
          * @note Just used for testing purpose
          */
         bool operator==(const TranscendentalExpression &rhs) const {
-            return ln == rhs.ln && an == rhs.an;
+            return util::floatNear(ln, rhs.ln) && util::floatNear(an, rhs.an);
         }
 
         /**
-         * Checks two TranscendentalExpressions for inequality.
-         * @warning This method compares doubles! So only exact copies will evaluate to false.
+         * Checks two TranscendentalExpressions for inequality with another one by comparing their members values
+         * ULP distances.
          * @param rhs the other TranscendentalExpressions
          * @return false if unequal
          *
@@ -191,20 +194,23 @@ namespace polyhedralGravity {
         double d;
 
         /**
-         * Checking the equality of two this Hessian Plane with another one by comparing their members.
-         * @warning This method compares doubles! So only exact copies will evaluate to true.
+         * Checking the equality of two this Hessian Plane with another one by comparing their members values
+         * ULP distances.
          * @param rhs other HessianPlane
          * @return true if equal
          *
          * @note Just used for testing purpose
          */
         bool operator==(const HessianPlane &rhs) const {
-            return a == rhs.a && b == rhs.b && c == rhs.c && d == rhs.d;
+            return util::floatNear(a, rhs.a) &&
+                   util::floatNear(b, rhs.b) &&
+                   util::floatNear(c, rhs.c) &&
+                   util::floatNear(d, rhs.d);
         }
 
         /**
-         * Checking the inequality of two this Hessian Plane with another one by comparing their members.
-         * @warning This method compares doubles! So only exact copies will evaluate to false.
+         * Checking the inequality of two this Hessian Plane with another one with another one by comparing their members values'
+         * ULP distances.
          * @param rhs other HessianPlane
          * @return true if unequal
          *
@@ -213,6 +219,17 @@ namespace polyhedralGravity {
         bool operator!=(const HessianPlane &rhs) const {
             return !(rhs == *this);
         }
+
+         /**
+         * Pretty output of this struct on the given ostream.
+         * @param os the ostream
+         * @param hessianPlane a HessianPlane
+         * @return os
+         */
+        friend std::ostream &operator<<(std::ostream &os, const HessianPlane &hessianPlane) {
+            os << "a: " << hessianPlane.a << " b: " << hessianPlane.b << " c: " << hessianPlane.c << " d: " << hessianPlane.d;
+            return os;
+        }
     };
 
 }
diff --git a/src/polyhedralGravity/util/UtilityConstants.h b/src/polyhedralGravity/util/UtilityConstants.h
index f9ab381..043c555 100644
--- a/src/polyhedralGravity/util/UtilityConstants.h
+++ b/src/polyhedralGravity/util/UtilityConstants.h
@@ -10,7 +10,6 @@ namespace polyhedralGravity::util {
      */
     constexpr double EPSILON = 1e-14;
 
-
     /**
      * PI with enough precision
      */
diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
new file mode 100644
index 0000000..022edb5
--- /dev/null
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
@@ -0,0 +1,27 @@
+#include "UtilityFloatArithmetic.h"
+
+
+namespace polyhedralGravity::util {
+
+    bool floatNear(double lhs, double rhs, int ulpDistance) {
+        if (lhs == rhs) {
+            return true;
+        }
+        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
+            return false;
+        }
+        return reinterpret_cast<std::int64_t&>(rhs) - reinterpret_cast<std::int64_t&>(lhs) <= ulpDistance;
+    }
+
+    bool floatNear(float lhs, float rhs, int ulpDistance) {
+        if (lhs == rhs) {
+            return true;
+        }
+        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
+            return false;
+        }
+        return reinterpret_cast<std::int32_t&>(rhs) - reinterpret_cast<std::int32_t&>(lhs) <= ulpDistance;
+    }
+
+
+}
\ No newline at end of file
diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.h b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
new file mode 100644
index 0000000..15cec69
--- /dev/null
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
@@ -0,0 +1,48 @@
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+
+namespace polyhedralGravity::util {
+
+    /**
+     * The maximal allowed ULP distance as computed by {@ref polyhedralGravity::util::floatNear}
+     * utilized for FloatingPoint comparisons in this implementation
+     *
+     * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+     */
+    constexpr int MAX_ULP_DISTANCE = 4;
+
+
+    /**
+     * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
+     *
+     * @param lhs The left hand side floating point number to compare.
+     * @param rhs The right hand side floating point number to compare.
+     * @param ulpDistance The maximum acceptable ULP distance between the two floating points
+     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     *
+     * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
+     *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
+     *
+     * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+     */
+    bool floatNear(double lhs, double rhs, int ulpDistance = MAX_ULP_DISTANCE);
+
+    /**
+     * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
+     *
+     * @param lhs The left hand side floating point number to compare.
+     * @param rhs The right hand side floating point number to compare.
+     * @param ulpDistance The maximum acceptable ULP distance between the two floating points
+     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     *
+     * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
+     *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
+     *
+     * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+     */
+    bool floatNear(float lhs, float rhs, int ulpDistance = MAX_ULP_DISTANCE);
+
+}
diff --git a/test/model/GoogleTestMatcher.h b/test/model/GoogleTestMatcher.h
new file mode 100644
index 0000000..db6757e
--- /dev/null
+++ b/test/model/GoogleTestMatcher.h
@@ -0,0 +1,27 @@
+#pragma once
+
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
+
+MATCHER_P2(FloatContainter2D, container, ulpDiff, "Comparing 2D Containers") {
+    if (container.size() != arg.size()) {
+        *result_listener << "The top-level container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
+        return false;
+    }
+    for (size_t idx = 0; idx < std::size(container); ++idx) {
+        if (container[idx].size() != arg[idx].size()) {
+            *result_listener << "The inner container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+            return false;
+        }
+        for (size_t idy = 0; idy < std::size(container[idx]); ++idy) {
+            if (!polyhedralGravity::util::floatNear(container[idx][idy], arg[idx][idy], ulpDiff)) {
+                *result_listener << "The elements at idx = " << idx << " and idy = " << idy << " do not match. Values: " << container[idx][idy] << " != " << arg[idx][idy];
+                return false;
+            }
+        }
+    }
+    return true;
+}
\ No newline at end of file
diff --git a/test/model/GravityModelBigTest.cpp b/test/model/GravityModelBigTest.cpp
index e80285d..0e77204 100644
--- a/test/model/GravityModelBigTest.cpp
+++ b/test/model/GravityModelBigTest.cpp
@@ -10,6 +10,7 @@
 #include "polyhedralGravity/model/Polyhedron.h"
 
 #include "GravityModelVectorUtility.h"
+#include "GoogleTestMatcher.h"
 
 
 /**
@@ -268,7 +269,7 @@ TEST_F(GravityModelBigTest, PlaneUnitNormals) {
 
     auto actualPlaneUnitNormals = polyhedralGravity::GravityModel::calculatePlaneUnitNormals(expectedGij);
 
-    ASSERT_THAT(actualPlaneUnitNormals, ContainerEq(expectedPlaneUnitNormals));
+    ASSERT_THAT(actualPlaneUnitNormals, FloatContainter2D(expectedPlaneUnitNormals, 10));
 }
 
 TEST_F(GravityModelBigTest, SegmentUnitNormals) {
@@ -291,12 +292,11 @@ TEST_F(GravityModelBigTest, PlaneNormalOrientations) {
 
 TEST_F(GravityModelBigTest, HessianPlane) {
     using namespace testing;
-    using namespace polyhedralGravity;
 
     auto actualHessianPlane =
             polyhedralGravity::GravityModel::calculateFacesToHessianPlanes(_computationPoint, _polyhedron);
 
-    ASSERT_EQ(actualHessianPlane, expectedHessianPlanes);
+    ASSERT_THAT(actualHessianPlane, ContainerEq(expectedHessianPlanes));
 }
 
 TEST_F(GravityModelBigTest, PlaneDistances) {
diff --git a/test/util/UtilityFloatArithmeticTest.cpp b/test/util/UtilityFloatArithmeticTest.cpp
new file mode 100644
index 0000000..c155d5e
--- /dev/null
+++ b/test/util/UtilityFloatArithmeticTest.cpp
@@ -0,0 +1,31 @@
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
+
+
+TEST(UtilityFloatArithmeticTest, Test01) {
+    ASSERT_FALSE(polyhedralGravity::util::floatNear(3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::floatNear(-3.0, -4.0));
+    ASSERT_FALSE(polyhedralGravity::util::floatNear(-3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::floatNear(3.0, -4.0));
+    ASSERT_TRUE(polyhedralGravity::util::floatNear(1.0, 1.0));
+}
+
+TEST(UtilityFloatArithmeticTest, Test02) {
+    ASSERT_TRUE(polyhedralGravity::util::floatNear(9.40569e-05, 9.40569e-05));
+    ASSERT_TRUE(polyhedralGravity::util::floatNear(-0.000150712, -0.000150712));
+    ASSERT_TRUE(polyhedralGravity::util::floatNear(0.000135291, 0.000135291));
+    ASSERT_TRUE(polyhedralGravity::util::floatNear(-8.63978e-05, -8.63978e-05));
+}
+
+TEST(UtilityFloatArithmeticTest, Test03) {
+    double x = 300.3;
+    double y = std::nextafter(x, 400.0);
+    double z = std::nextafter(y, 400.0);
+    ASSERT_EQ(1, reinterpret_cast<int64_t&>(y) - reinterpret_cast<int64_t&>(x));
+    ASSERT_EQ(2, reinterpret_cast<int64_t&>(z) - reinterpret_cast<int64_t&>(x));
+}
\ No newline at end of file

From cbc508a8cfb349aa9c17755bcbace8c527384b56 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Tue, 7 May 2024 21:08:29 +0200
Subject: [PATCH 02/10] improves numerical test stability

This improves the stability of the tests and makes them independent of the utilized OS/ Compiler/ Architecture by:
- adding dedicated methods for FloatingPoint comparison using
  - UPL distance
  - Relative Epsilon
- adding custom GoogleTestMatchers using these new methods
- updating the DataStructures Equality methods
---
 .../model/GravityModelData.h                  | 18 ++---
 .../model/GravityModelDetail.cpp              | 24 +++---
 .../model/GravityModelDetail.h                |  1 +
 src/polyhedralGravity/model/Polyhedron.cpp    |  6 +-
 src/polyhedralGravity/model/Polyhedron.h      |  1 +
 src/polyhedralGravity/util/UtilityConstants.h |  8 --
 .../util/UtilityFloatArithmetic.cpp           | 46 ++++++++---
 .../util/UtilityFloatArithmetic.h             | 53 +++++++++----
 test/model/GoogleTestMatcher.h                | 54 +++++++++++--
 test/model/GravityModelBigTest.cpp            | 11 ++-
 test/util/UtilityFloatArithmeticTest.cpp      | 79 ++++++++++++++-----
 11 files changed, 214 insertions(+), 87 deletions(-)

diff --git a/src/polyhedralGravity/model/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h
index ecc5d2d..2d3abb5 100644
--- a/src/polyhedralGravity/model/GravityModelData.h
+++ b/src/polyhedralGravity/model/GravityModelData.h
@@ -84,10 +84,10 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const Distance &rhs) const {
-            return util::floatNear(l1, rhs.l1) &&
-                   util::floatNear(l2, rhs.l2) &&
-                   util::floatNear(s1, rhs.s1) &&
-                   util::floatNear(s2, rhs.s2);
+            return util::almostEqualRelative(l1, rhs.l1) &&
+                   util::almostEqualRelative(l2, rhs.l2) &&
+                   util::almostEqualRelative(s1, rhs.s1) &&
+                   util::almostEqualRelative(s2, rhs.s2);
         }
 
         /**
@@ -141,7 +141,7 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const TranscendentalExpression &rhs) const {
-            return util::floatNear(ln, rhs.ln) && util::floatNear(an, rhs.an);
+            return util::almostEqualRelative(ln, rhs.ln) && util::almostEqualRelative(an, rhs.an);
         }
 
         /**
@@ -202,10 +202,10 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const HessianPlane &rhs) const {
-            return util::floatNear(a, rhs.a) &&
-                   util::floatNear(b, rhs.b) &&
-                   util::floatNear(c, rhs.c) &&
-                   util::floatNear(d, rhs.d);
+            return util::almostEqualRelative(a, rhs.a) &&
+                   util::almostEqualRelative(b, rhs.b) &&
+                   util::almostEqualRelative(c, rhs.c) &&
+                   util::almostEqualRelative(d, rhs.d);
         }
 
         /**
diff --git a/src/polyhedralGravity/model/GravityModelDetail.cpp b/src/polyhedralGravity/model/GravityModelDetail.cpp
index 6c0038b..1f3e70b 100644
--- a/src/polyhedralGravity/model/GravityModelDetail.cpp
+++ b/src/polyhedralGravity/model/GravityModelDetail.cpp
@@ -28,7 +28,7 @@ namespace polyhedralGravity::GravityModel::detail {
         //Calculate N_i * -G_i1 where * is the dot product and then use the inverted sgn
         //We abstain on the double multiplication with -1 in the line above and beyond since two
         //times multiplying with -1 equals no change
-        return sgn(dot(planeUnitNormal, vertex0), util::EPSILON);
+        return sgn(dot(planeUnitNormal, vertex0), util::EPSILON_ZERO_OFFSET);
     }
 
     HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) {
@@ -95,7 +95,7 @@ namespace polyhedralGravity::GravityModel::detail {
                        segmentNormalOrientations.begin(),
                        [&projectionPointOnPlane](const Array3 &segmentUnitNormal, const Array3 &vertex) {
                            using namespace util;
-                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0;
+                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON_ZERO_OFFSET) * -1.0;
                        });
         return segmentNormalOrientations;
     }
@@ -190,8 +190,8 @@ namespace polyhedralGravity::GravityModel::detail {
 
                               //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on
                               // the direction of a specific segment (P coincides with P' and P'')
-                              if (std::abs(distance.s1 - distance.l1) < EPSILON &&
-                                  std::abs(distance.s2 - distance.l2) < EPSILON) {
+                              if (std::abs(distance.s1 - distance.l1) < EPSILON_ZERO_OFFSET &&
+                                  std::abs(distance.s2 - distance.l2) < EPSILON_ZERO_OFFSET) {
                                   //4. Option - Case 2: P is located on the segment from its right side
                                   // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2|
                                   if (distance.s2 < distance.s1) {
@@ -200,7 +200,7 @@ namespace polyhedralGravity::GravityModel::detail {
                                       distance.l1 *= -1.0;
                                       distance.l2 *= -1.0;
                                       return distance;
-                                  } else if (std::abs(distance.s2 - distance.s1) < EPSILON) {
+                                  } else if (std::abs(distance.s2 - distance.s1) < EPSILON_ZERO_OFFSET) {
                                       //4. Option - Case 1: P is located inside the segment (s2 == s1)
                                       // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2|
                                       distance.s1 *= -1.0;
@@ -264,9 +264,9 @@ namespace polyhedralGravity::GravityModel::detail {
                               // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR
                               // the 1D and 3D distances are smaller than some EPSILON
                               // then LN_pq can be set to zero
-                              if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) ||
-                                  (std::abs(distance.s1 + distance.s2) < EPSILON &&
-                                   std::abs(distance.l1 + distance.l2) < EPSILON)) {
+                              if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON_ZERO_OFFSET || r2Norm < EPSILON_ZERO_OFFSET)) ||
+                                  (std::abs(distance.s1 + distance.s2) < EPSILON_ZERO_OFFSET &&
+                                   std::abs(distance.l1 + distance.l2) < EPSILON_ZERO_OFFSET)) {
                                   transcendentalExpressionPerSegment.ln = 0.0;
                               } else {
                                   //Implementation of
@@ -277,7 +277,7 @@ namespace polyhedralGravity::GravityModel::detail {
 
                               //Compute AN_pq according to (15)
                               // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!)
-                              if (planeDistance < EPSILON || segmentDistance < EPSILON) {
+                              if (planeDistance < EPSILON_ZERO_OFFSET || segmentDistance < EPSILON_ZERO_OFFSET) {
                                   transcendentalExpressionPerSegment.an = 0.0;
                               } else {
                                   //Implementation of:
@@ -331,7 +331,7 @@ namespace polyhedralGravity::GravityModel::detail {
             const unsigned int j = thrust::get<2>(tuple);
 
             //segmentNormalOrientation != 0.0
-            if (std::abs(segmentNormalOrientation) > EPSILON) {
+            if (std::abs(segmentNormalOrientation) > EPSILON_ZERO_OFFSET) {
                 return false;
             }
 
@@ -358,14 +358,14 @@ namespace polyhedralGravity::GravityModel::detail {
             j = thrust::get<1>(tuple);
 
             //segmentNormalOrientation != 0.0
-            if (std::abs(segmentNormalOrientation) > EPSILON) {
+            if (std::abs(segmentNormalOrientation) > EPSILON_ZERO_OFFSET) {
                 return false;
             }
 
             r1Norm = projectionPointVertexNorms[(j + 1) % 3];
             r2Norm = projectionPointVertexNorms[j];
             //r1Norm == 0.0 || r2Norm == 0.0
-            return r1Norm < EPSILON || r2Norm < EPSILON;
+            return r1Norm < EPSILON_ZERO_OFFSET || r2Norm < EPSILON_ZERO_OFFSET;
         })) {
             using namespace util;
             //Two segment vectors G_1 and G_2 of this plane
diff --git a/src/polyhedralGravity/model/GravityModelDetail.h b/src/polyhedralGravity/model/GravityModelDetail.h
index 8336afd..ccb3f34 100644
--- a/src/polyhedralGravity/model/GravityModelDetail.h
+++ b/src/polyhedralGravity/model/GravityModelDetail.h
@@ -21,6 +21,7 @@
 #include "polyhedralGravity/util/UtilityConstants.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityThrust.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 #include "polyhedralGravity/output/Logging.h"
 
 namespace polyhedralGravity::GravityModel::detail {
diff --git a/src/polyhedralGravity/model/Polyhedron.cpp b/src/polyhedralGravity/model/Polyhedron.cpp
index a6780da..6f6d5c9 100644
--- a/src/polyhedralGravity/model/Polyhedron.cpp
+++ b/src/polyhedralGravity/model/Polyhedron.cpp
@@ -202,7 +202,7 @@ namespace polyhedralGravity {
         const Array3 rayVector = normal(segmentVector1, segmentVector2);
 
         // The origin of the array has a slight offset in direction of the normal
-        const Array3 rayOrigin = centroid + (rayVector * EPSILON);
+        const Array3 rayOrigin = centroid + (rayVector * EPSILON_ZERO_OFFSET);
 
         // Count every triangular face which is intersected by the ray
         const auto &[begin, end] = this->transformIterator();
@@ -224,7 +224,7 @@ namespace polyhedralGravity {
         const Array3 edge2 = triangle[2] - triangle[0];
         const Array3 h = cross(rayVector, edge2);
         const double a = dot(edge1, h);
-        if (a > -EPSILON && a < EPSILON) {
+        if (a > -EPSILON_ZERO_OFFSET && a < EPSILON_ZERO_OFFSET) {
             return nullptr;
         }
 
@@ -242,7 +242,7 @@ namespace polyhedralGravity {
         }
 
         const double t = f * dot(edge2, q);
-        if (t > EPSILON) {
+        if (t > EPSILON_ZERO_OFFSET) {
             return std::make_unique<Array3>(rayOrigin + rayVector * t);
         } else {
             return nullptr;
diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h
index c80ed30..cc4a461 100644
--- a/src/polyhedralGravity/model/Polyhedron.h
+++ b/src/polyhedralGravity/model/Polyhedron.h
@@ -22,6 +22,7 @@
 #include "thrust/iterator/transform_iterator.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityConstants.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
 namespace polyhedralGravity {
 
diff --git a/src/polyhedralGravity/util/UtilityConstants.h b/src/polyhedralGravity/util/UtilityConstants.h
index 043c555..4292850 100644
--- a/src/polyhedralGravity/util/UtilityConstants.h
+++ b/src/polyhedralGravity/util/UtilityConstants.h
@@ -2,14 +2,6 @@
 
 namespace polyhedralGravity::util {
 
-    /**
-     * The EPSILON used in the polyhedral gravity model.
-     * @related Used to determine if a floating point number is equal to zero as threshold for rounding errors
-     * @related Used for the sgn() function to determine the sign of a double value. Different compilers
-     * produce different results if no EPSILON is applied for the comparison!
-     */
-    constexpr double EPSILON = 1e-14;
-
     /**
      * PI with enough precision
      */
diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
index 022edb5..c1c60ec 100644
--- a/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp
@@ -3,25 +3,49 @@
 
 namespace polyhedralGravity::util {
 
-    bool floatNear(double lhs, double rhs, int ulpDistance) {
+    template<typename FloatType>
+    bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance) {
+        static_assert(std::is_same_v<FloatType, float> || std::is_same_v<FloatType, double>,
+                      "Template argument must be FloatType: Either float or double!");
+
+        // In case the floats are equal in their representation, return true
         if (lhs == rhs) {
             return true;
         }
-        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
+
+        // In case the signs mismatch, return false
+        if (lhs < static_cast<FloatType>(0.0) && rhs > static_cast<FloatType>(0.0) ||
+            lhs > static_cast<FloatType>(0.0) && rhs < static_cast<FloatType>(0.0)) {
             return false;
-        }
-        return reinterpret_cast<std::int64_t&>(rhs) - reinterpret_cast<std::int64_t&>(lhs) <= ulpDistance;
-    }
+            }
 
-    bool floatNear(float lhs, float rhs, int ulpDistance) {
-        if (lhs == rhs) {
-            return true;
+        if constexpr (std::is_same_v<FloatType, float>) {
+            // In case of float, compute ULP distance by interpreting float as 32-bit integer
+            return reinterpret_cast<std::int32_t&>(rhs) - reinterpret_cast<std::int32_t&>(lhs) <= ulpDistance;
         }
-        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
-            return false;
+        else if constexpr (std::is_same_v<FloatType, double>) {
+            // In case of double, compute ULP distance by interpreting double as 64-bit integer
+            return reinterpret_cast<std::int64_t&>(rhs) - reinterpret_cast<std::int64_t&>(lhs) <= ulpDistance;
         }
-        return reinterpret_cast<std::int32_t&>(rhs) - reinterpret_cast<std::int32_t&>(lhs) <= ulpDistance;
+
+        // Due to the static_assert above, this should not happen
+        return false;
+    }
+
+    // Template Instantations for float and double (required since definition is in .cpp file,
+    //      also makes the static assert not strictly necessary)
+    template bool almostEqualUlps<float>(float lhs, float rhs, int ulpDistance);
+    template bool almostEqualUlps<double>(double lhs, double rhs, int ulpDistance);
+
+
+    template<typename FloatType>
+    bool almostEqualRelative(FloatType lhs, FloatType rhs, double epsilon) {
+        const FloatType diff = std::abs(rhs - lhs);
+        const FloatType largerValue = std::max(std::abs(rhs), std::abs(lhs));
+        return diff <= largerValue * epsilon;
     }
 
+    template bool almostEqualRelative<float>(float lhs, float rhs, double epsilon);
+    template bool almostEqualRelative<double>(double lhs, double rhs, double epsilon);
 
 }
\ No newline at end of file
diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.h b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
index 15cec69..8a9f995 100644
--- a/src/polyhedralGravity/util/UtilityFloatArithmetic.h
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
@@ -3,12 +3,32 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
+#include <type_traits>
 
 namespace polyhedralGravity::util {
 
     /**
-     * The maximal allowed ULP distance as computed by {@ref polyhedralGravity::util::floatNear}
-     * utilized for FloatingPoint comparisons in this implementation
+     * The EPSILON used in the polyhedral gravity model to determine a radius around zero/ to use as sligth offset.
+     * @related Used to determine if a floating point number is equal to zero as threshold for rounding errors
+     * @related Used for the sgn() function to determine the sign of a double value. Different compilers
+     * produce different results if no EPSILON is applied for the comparison!
+     */
+    constexpr double EPSILON_ZERO_OFFSET = 1e-14;
+
+    /**
+     * This relative EPSILON is utilized ONLY for testing purposes to compare intermediate values to
+     * Tsoulis' reference implementation Fortran.
+     * It is used in the {@ref polyhedralGravity::util::almostEqualRelative} function.
+     *
+     * @note While in theory no difference at all is observed when compiling this program on Linux using GCC on x86_64,
+     *  the intermediate values change when the program is compiled in different environments.
+     *  Hence, this solves the flakiness of the tests when on different plattforms
+     */
+    constexpr double EPSILON_ALMOST_EQUAL = 1e-10;
+
+    /**
+     * The maximal allowed ULP distance utilized for FloatingPoint comparisons using the
+     * {@ref polyhedralGravity::util::almostEqualUlps} function.
      *
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
@@ -18,31 +38,36 @@ namespace polyhedralGravity::util {
     /**
      * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
      *
+     * @tparam FloatType must be either double or float (ensured by static assertion)
      * @param lhs The left hand side floating point number to compare.
      * @param rhs The right hand side floating point number to compare.
      * @param ulpDistance The maximum acceptable ULP distance between the two floating points
-     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     *      for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
      *
      * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
      *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
-     *
+     * @example The ULP distance between 3.0 and std::nextafter(3.0, INFINITY) would be 1,
+     *      the ULP distance of 3.0 and std::nextafter(std::nextafter(3.0, INFINITY), INFINITY) would be 2, etc.
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
-    bool floatNear(double lhs, double rhs, int ulpDistance = MAX_ULP_DISTANCE);
+    template<typename FloatType>
+    bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance = MAX_ULP_DISTANCE);
 
     /**
-     * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
-     *
-     * @param lhs The left hand side floating point number to compare.
-     * @param rhs The right hand side floating point number to compare.
-     * @param ulpDistance The maximum acceptable ULP distance between the two floating points
-     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     * Function to check if two floating point numbers are relatively equal to each other within a given error range or tolerance.
      *
-     * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
-     *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
+     * @tparam FloatType must be either double or float (ensured by static assertion)
+     * @param lhs The first floating-point number to be compared.
+     * @param rhs The second floating-point number to be compared.
+     * @param epsilon The tolerance for comparison. Two numbers that are less than epsilon apart are considered equal.
+     *                The default value is `EPSILON`.
      *
+     * @return boolean value - Returns `true` if the absolute difference between `lhs` and `rhs` is less than or equal to
+     *                         the relative error factored by the larger of the magnitude of `lhs` and `rhs`. Otherwise, `false`.
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
-    bool floatNear(float lhs, float rhs, int ulpDistance = MAX_ULP_DISTANCE);
+    template<typename FloatType>
+    bool almostEqualRelative(FloatType lhs, FloatType rhs, double epsilon = EPSILON_ALMOST_EQUAL);
+
 
 }
diff --git a/test/model/GoogleTestMatcher.h b/test/model/GoogleTestMatcher.h
index db6757e..2c7b3d0 100644
--- a/test/model/GoogleTestMatcher.h
+++ b/test/model/GoogleTestMatcher.h
@@ -6,22 +6,66 @@
 
 #include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
-MATCHER_P2(FloatContainter2D, container, ulpDiff, "Comparing 2D Containers") {
+MATCHER_P(FloatContainter1D, container, "Comparing 1D Containers") {
     if (container.size() != arg.size()) {
-        *result_listener << "The top-level container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
+        *result_listener << "The container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
+        return false;
+    }
+    for (size_t idx = 0; idx < std::size(container); ++idx) {
+        if (!polyhedralGravity::util::almostEqualRelative(container[idx], arg[idx])) {
+            *result_listener << "The elements at idx = " << idx << " do not match. Values: " << container[idx]<< " != " << arg[idx];
+            return false;
+        }
+    }
+    return true;
+}
+
+// Using the FloatContainter1D would be nice, but this leads to issues with template instatantation
+// Hence, this easy way and a double-for-loop (but at least better fitting messages)
+MATCHER_P(FloatContainter2D, container, "Comparing 2D Containers") {
+    if (container.size() != arg.size()) {
+        *result_listener << "The container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
         return false;
     }
     for (size_t idx = 0; idx < std::size(container); ++idx) {
         if (container[idx].size() != arg[idx].size()) {
-            *result_listener << "The inner container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+            *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
             return false;
         }
         for (size_t idy = 0; idy < std::size(container[idx]); ++idy) {
-            if (!polyhedralGravity::util::floatNear(container[idx][idy], arg[idx][idy], ulpDiff)) {
+            if (!polyhedralGravity::util::almostEqualRelative(container[idx][idy], arg[idx][idy])) {
                 *result_listener << "The elements at idx = " << idx << " and idy = " << idy << " do not match. Values: " << container[idx][idy] << " != " << arg[idx][idy];
                 return false;
             }
         }
     }
     return true;
-}
\ No newline at end of file
+}
+
+// Using the FloatContainter2D would be nice, but this leads to issues with template instatantation
+// Hence, this easy way and a triple-for-loop (but at least better fitting messages)
+MATCHER_P(FloatContainter3D, container, "Comparing 3D Containers") {
+    if (container.size() != arg.size()) {
+        *result_listener << "The container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
+        return false;
+    }
+    for (size_t idx = 0; idx < std::size(container); ++idx) {
+        if (container[idx].size() != arg[idx].size()) {
+            *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+            return false;
+        }
+        for (size_t idy = 0; idy < std::size(container[idx]); ++idy) {
+            if (container[idx][idy].size() != arg[idx][idy].size()) {
+                *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+                return false;
+            }
+            for (size_t idz = 0; idz < std::size(container[idx][idy]); ++idz) {
+                if (!polyhedralGravity::util::almostEqualRelative(container[idx][idy][idz], arg[idx][idy][idz])) {
+                    *result_listener << "The elements at idx = " << idx << " and idy = " << idy << " and idz = " << idz << " do not match. Values: " << container[idx][idy][idz] << " != " << arg[idx][idy][idz];
+                    return false;
+                }
+            }
+        }
+    }
+    return true;
+}
diff --git a/test/model/GravityModelBigTest.cpp b/test/model/GravityModelBigTest.cpp
index 0e77204..9f6b21d 100644
--- a/test/model/GravityModelBigTest.cpp
+++ b/test/model/GravityModelBigTest.cpp
@@ -269,16 +269,15 @@ TEST_F(GravityModelBigTest, PlaneUnitNormals) {
 
     auto actualPlaneUnitNormals = polyhedralGravity::GravityModel::calculatePlaneUnitNormals(expectedGij);
 
-    ASSERT_THAT(actualPlaneUnitNormals, FloatContainter2D(expectedPlaneUnitNormals, 10));
+    ASSERT_THAT(actualPlaneUnitNormals, FloatContainter2D(expectedPlaneUnitNormals));
 }
 
 TEST_F(GravityModelBigTest, SegmentUnitNormals) {
     using namespace testing;
 
-    auto actualSegmentUnitNormals = polyhedralGravity::GravityModel::calculateSegmentUnitNormals(expectedGij,
-                                                                                                 expectedPlaneUnitNormals);
+    auto actualSegmentUnitNormals = polyhedralGravity::GravityModel::calculateSegmentUnitNormals(expectedGij, expectedPlaneUnitNormals);
 
-    ASSERT_THAT(actualSegmentUnitNormals, ContainerEq(expectedSegmentUnitNormals));
+    ASSERT_THAT(actualSegmentUnitNormals, FloatContainter3D(expectedSegmentUnitNormals));
 }
 
 TEST_F(GravityModelBigTest, PlaneNormalOrientations) {
@@ -304,7 +303,7 @@ TEST_F(GravityModelBigTest, PlaneDistances) {
 
     auto actualPlaneDistances = polyhedralGravity::GravityModel::calculatePlaneDistances(expectedHessianPlanes);
 
-    ASSERT_THAT(actualPlaneDistances, ContainerEq(expectedPlaneDistances));
+    ASSERT_THAT(actualPlaneDistances, FloatContainter1D(expectedPlaneDistances));
 }
 
 TEST_F(GravityModelBigTest, OrthogonalProjectionPointsOnPlane) {
@@ -363,7 +362,7 @@ TEST_F(GravityModelBigTest, SegmentDistances) {
             polyhedralGravity::GravityModel::calculateSegmentDistances(expectedOrthogonalProjectionPointsOnPlane,
                                                                        expectedOrthogonalProjectionPointsOnSegment);
 
-    ASSERT_THAT(actualSegmentDistances, ContainerEq(expectedSegmentDistances));
+    ASSERT_THAT(actualSegmentDistances, FloatContainter2D(expectedSegmentDistances));
 }
 
 TEST_F(GravityModelBigTest, DistancesPerSegmentEndpoint) {
diff --git a/test/util/UtilityFloatArithmeticTest.cpp b/test/util/UtilityFloatArithmeticTest.cpp
index c155d5e..10c0fc6 100644
--- a/test/util/UtilityFloatArithmeticTest.cpp
+++ b/test/util/UtilityFloatArithmeticTest.cpp
@@ -7,25 +7,66 @@
 #include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
 
-TEST(UtilityFloatArithmeticTest, Test01) {
-    ASSERT_FALSE(polyhedralGravity::util::floatNear(3.0, 4.0));
-    ASSERT_FALSE(polyhedralGravity::util::floatNear(-3.0, -4.0));
-    ASSERT_FALSE(polyhedralGravity::util::floatNear(-3.0, 4.0));
-    ASSERT_FALSE(polyhedralGravity::util::floatNear(3.0, -4.0));
-    ASSERT_TRUE(polyhedralGravity::util::floatNear(1.0, 1.0));
-}
+TEST(UtilityFloatArithmeticTest, TestAlmostEqualUlps) {
+    // Checking the signess && identity
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(-3.0, -4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(-3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, -4.0));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(1.0, 1.0));
+
+    // Some random values, which are equal
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(9.40569e-05, 9.40569e-05));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(-0.000150712, -0.000150712));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(0.000135291, 0.000135291));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(-8.63978e-05, -8.63978e-05));
+
+    /* FROM HERE IT GETS INTERESTING */
+
+    // Checking the relative difference
+    // The ULP distance is always greater than 4 for these "big" epsilons
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, 3.0 + 1e-9));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, 3.0 + 1e-10));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, 3.0 + 1e-11));
+
+    // Checking the the sensitivity towards the next floating point
+    // Note: The default maximal ULPS distance is set to 4
+    // The ULP distance is one time higher than 4 leading to false, one time lower equal leading to true
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0);
+    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0), 4.0);
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(3.0, fourHops));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, fiveHops));
+
 
-TEST(UtilityFloatArithmeticTest, Test02) {
-    ASSERT_TRUE(polyhedralGravity::util::floatNear(9.40569e-05, 9.40569e-05));
-    ASSERT_TRUE(polyhedralGravity::util::floatNear(-0.000150712, -0.000150712));
-    ASSERT_TRUE(polyhedralGravity::util::floatNear(0.000135291, 0.000135291));
-    ASSERT_TRUE(polyhedralGravity::util::floatNear(-8.63978e-05, -8.63978e-05));
 }
 
-TEST(UtilityFloatArithmeticTest, Test03) {
-    double x = 300.3;
-    double y = std::nextafter(x, 400.0);
-    double z = std::nextafter(y, 400.0);
-    ASSERT_EQ(1, reinterpret_cast<int64_t&>(y) - reinterpret_cast<int64_t&>(x));
-    ASSERT_EQ(2, reinterpret_cast<int64_t&>(z) - reinterpret_cast<int64_t&>(x));
-}
\ No newline at end of file
+TEST(UtilityFloatArithmeticTest, TestAlmostEqualRelative) {
+    // Checking the signess && identity
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualRelative(3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualRelative(-3.0, -4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualRelative(-3.0, 4.0));
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualRelative(3.0, -4.0));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(1.0, 1.0));
+
+    // Some random values, which are equal
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(9.40569e-05, 9.40569e-05));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(-0.000150712, -0.000150712));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(0.000135291, 0.000135291));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(-8.63978e-05, -8.63978e-05));
+
+    /* FROM HERE IT GETS INTERESTING */
+
+    // Checking the relative difference
+    // Note: 1e-10 is the sensitvity of almostEqualRelative
+    // The method returns true if the relative difference is smaller equal than 1e-10
+    ASSERT_FALSE(polyhedralGravity::util::almostEqualRelative(3.0, 3.0 + 1e-9));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, 3.0 + 1e-10));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, 3.0 + 1e-11));
+
+    // Checking the the sensitivity towards the next floating point
+    // The ULP distnace is really small, this method will always return true
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0);
+    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0), 4.0);
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fourHops));
+    ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fiveHops));
+}

From f1f26d7bcc7adb77dcedc35d4dc53171413957a2 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Tue, 7 May 2024 21:45:53 +0200
Subject: [PATCH 03/10] Update code documentation & Bump version number

Improved the code comments for better clarity and updated some of the examples in the docs. Also corrected minor typos in documentation and comments. The package version number is updated to 3.1.
---
 docs/api/util.rst                                       | 2 +-
 docs/quickstart/examples_cpp.rst                        | 8 ++++----
 setup.py                                                | 2 +-
 src/polyhedralGravity/model/Polyhedron.h                | 8 ++++----
 src/polyhedralGravity/util/UtilityFloatArithmetic.h     | 8 ++++----
 src/polyhedralGravityPython/PolyhedralGravityPython.cpp | 6 +++---
 6 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/docs/api/util.rst b/docs/api/util.rst
index 4c85418..9dc761b 100644
--- a/docs/api/util.rst
+++ b/docs/api/util.rst
@@ -6,7 +6,7 @@ Overview
 
 The namespace :code:`polyhedralGravity::util` contains utility
 for operations on iterable Containers, Constants and syntactic
-sugar for using the thrird party dependency :code:`thrust`.
+sugar for using the third party dependency :code:`thrust`.
 
 Documentation
 -------------
diff --git a/docs/quickstart/examples_cpp.rst b/docs/quickstart/examples_cpp.rst
index 40b9135..7cc49fc 100644
--- a/docs/quickstart/examples_cpp.rst
+++ b/docs/quickstart/examples_cpp.rst
@@ -79,7 +79,7 @@ and we check the if the plane unit normals are actually outwards pointing
 
         // Returns either a single of vector of results
         // Here, we only have one point. Thus we get a single result
-        Polyhedron polyhedron{vertices, faces, density, PolyhedronIntegrity::VERIFY};
+        Polyhedron polyhedron{vertices, faces, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY};
         const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, point, false);
 
 
@@ -131,12 +131,12 @@ The result will always be fine.
 
         // Reading the configuration from a yaml file
         std::shared_ptr<ConfigSource> config = std::make_shared<YAMLConfigReader>("config.yaml");
-        Polyhedron poly = config->getDataSource()->getPolyhedron();
+        auto polyhedralSource = config->getDataSource()->getPolyhedralSource();
         double density = config->getDensity();
         std::array<double, 3> point = config->getPointsOfInterest()[0];
 
-        Polyhedron polyhedron{files, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL};
-        GravityResult result = GravityModel::evaluate(poly, density, point);
+        Polyhedron polyhedron{polyhedralSource, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL};
+        const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, point);
 
 
 GravityEvaluable (with caching)
diff --git a/setup.py b/setup.py
index aaa5f30..70d486d 100644
--- a/setup.py
+++ b/setup.py
@@ -175,7 +175,7 @@ def build_extension(self, ext):
 # --------------------------------------------------------------------------------
 setup(
     name="polyhedral_gravity",
-    version="3.0",
+    version="3.1",
     author="Jonas Schuhmacher",
     author_email="jonas.schuhmacher@tum.de",
     description="Package to compute full gravity tensor of a given constant density polyhedron for arbitrary points "
diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h
index cc4a461..ef74b7a 100644
--- a/src/polyhedralGravity/model/Polyhedron.h
+++ b/src/polyhedralGravity/model/Polyhedron.h
@@ -151,7 +151,7 @@ namespace polyhedralGravity {
          *
          * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
          * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
-         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         * @throws std::invalid_argument dpending on the {@link integrity} flag
          */
         Polyhedron(
                 const std::vector<Array3> &vertices,
@@ -170,7 +170,7 @@ namespace polyhedralGravity {
          *
          * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
          * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
-         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         * @throws std::invalid_argument dpending on the {@link integrity} flag
          */
         Polyhedron(
                 const PolyhedralSource &polyhedralSource,
@@ -188,7 +188,7 @@ namespace polyhedralGravity {
          *
          * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
          * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
-         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         * @throws std::invalid_argument dpending on the {@link integrity} flag
          */
         Polyhedron(const PolyhedralFiles &polyhedralFiles, double density,
                    const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
@@ -204,7 +204,7 @@ namespace polyhedralGravity {
          *
          * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
          * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
-         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         * @throws std::invalid_argument dpending on the {@link integrity} flag
          */
         Polyhedron(const std::variant<PolyhedralSource, PolyhedralFiles> &polyhedralSource, double density,
                    const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.h b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
index 8a9f995..60507c7 100644
--- a/src/polyhedralGravity/util/UtilityFloatArithmetic.h
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
@@ -18,7 +18,7 @@ namespace polyhedralGravity::util {
     /**
      * This relative EPSILON is utilized ONLY for testing purposes to compare intermediate values to
      * Tsoulis' reference implementation Fortran.
-     * It is used in the {@ref polyhedralGravity::util::almostEqualRelative} function.
+     * It is used in the {@link polyhedralGravity::util::almostEqualRelative} function.
      *
      * @note While in theory no difference at all is observed when compiling this program on Linux using GCC on x86_64,
      *  the intermediate values change when the program is compiled in different environments.
@@ -28,7 +28,7 @@ namespace polyhedralGravity::util {
 
     /**
      * The maximal allowed ULP distance utilized for FloatingPoint comparisons using the
-     * {@ref polyhedralGravity::util::almostEqualUlps} function.
+     * {@link polyhedralGravity::util::almostEqualUlps} function.
      *
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
@@ -42,7 +42,7 @@ namespace polyhedralGravity::util {
      * @param lhs The left hand side floating point number to compare.
      * @param rhs The right hand side floating point number to compare.
      * @param ulpDistance The maximum acceptable ULP distance between the two floating points
-     *      for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     *      for which they would be considered near each other. This is optional and by default, it will be {@link MAX_ULP_DISTANCE}.
      *
      * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
      *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
@@ -60,7 +60,7 @@ namespace polyhedralGravity::util {
      * @param lhs The first floating-point number to be compared.
      * @param rhs The second floating-point number to be compared.
      * @param epsilon The tolerance for comparison. Two numbers that are less than epsilon apart are considered equal.
-     *                The default value is `EPSILON`.
+     *                The default value is {@link EPSILON_ALMOST_EQUAL}.
      *
      * @return boolean value - Returns `true` if the absolute difference between `lhs` and `rhs` is less than or equal to
      *                         the relative error factored by the larger of the magnitude of `lhs` and `rhs`. Otherwise, `false`.
diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
index 75a785c..e6f4d79 100644
--- a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
+++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
@@ -152,8 +152,8 @@ PYBIND11_MODULE(polyhedral_gravity, m) {
 
                                         * :code:`AUTOMATIC` (Default): Prints to stdout and throws ValueError if normal_orientation is wrong/ inconsisten
                                         * :code:`VERIFY`: Like :code:`AUTOMATIC`, but does not print to stdout
-                                        * :code:`DISABLE`: Recommened, when you know the mesh to avoid to pay :math:`O(n^2)` runtime. Disables ALL checks
-                                        * :code`HEAL`: Automatically fixes the normal_orientation and vertex ordering to the correct values
+                                        * :code:`DISABLE`: Recommened, when you are familiar with the mesh to avoid :math:`O(n^2)` runtime cost. Disables ALL checks
+                                        * :code:`HEAL`: Automatically fixes the normal_orientation and vertex ordering to the correct values
 
             Raises:
                 ValueError: If the faces array does not contain a reference to vertex 0 indicating an index start at 1
@@ -182,7 +182,7 @@ PYBIND11_MODULE(polyhedral_gravity, m) {
                 This utility is mainly for diagnostics and debugging purposes. If the polyhedron is constrcuted with `integrity_check`
                 set to :code:`AUTOMATIC` or :code:`VERIFY`, the construction fails anyways.
                 If set to :code:`HEAL`, this method should return an empty set (but maybe a different ordering than initially specified)
-                Only if set to code:`DISABLE`, then this method might actually return a set with faulty indices.
+                Only if set to :code:`DISABLE`, then this method might actually return a set with faulty indices.
                 Hence, if you want to know your mesh error. Construct the polyhedron with :code:`integrity_check=DISABLE` and call this method.
             )mydelimiter")
             .def("__getitem__", &Polyhedron::getResolvedFace, R"mydelimiter(

From 17d8fb424b6284e2a09a9d54f6bca54de574e2a6 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Tue, 7 May 2024 21:48:47 +0200
Subject: [PATCH 04/10] typos

---
 test/model/GoogleTestMatcher.h           | 4 ++--
 test/util/UtilityFloatArithmeticTest.cpp | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/model/GoogleTestMatcher.h b/test/model/GoogleTestMatcher.h
index 2c7b3d0..cd7a294 100644
--- a/test/model/GoogleTestMatcher.h
+++ b/test/model/GoogleTestMatcher.h
@@ -21,7 +21,7 @@ MATCHER_P(FloatContainter1D, container, "Comparing 1D Containers") {
 }
 
 // Using the FloatContainter1D would be nice, but this leads to issues with template instatantation
-// Hence, this easy way and a double-for-loop (but at least better fitting messages)
+// Hence, this is the easy way and a double-for-loop (but at least better fitting messages)
 MATCHER_P(FloatContainter2D, container, "Comparing 2D Containers") {
     if (container.size() != arg.size()) {
         *result_listener << "The container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
@@ -43,7 +43,7 @@ MATCHER_P(FloatContainter2D, container, "Comparing 2D Containers") {
 }
 
 // Using the FloatContainter2D would be nice, but this leads to issues with template instatantation
-// Hence, this easy way and a triple-for-loop (but at least better fitting messages)
+// Hence, this is the easy way and a triple-for-loop (but at least better fitting messages)
 MATCHER_P(FloatContainter3D, container, "Comparing 3D Containers") {
     if (container.size() != arg.size()) {
         *result_listener << "The container sizes do not match. Sizes: " << container.size() << " != " << arg.size();
diff --git a/test/util/UtilityFloatArithmeticTest.cpp b/test/util/UtilityFloatArithmeticTest.cpp
index 10c0fc6..4dc28eb 100644
--- a/test/util/UtilityFloatArithmeticTest.cpp
+++ b/test/util/UtilityFloatArithmeticTest.cpp
@@ -64,7 +64,7 @@ TEST(UtilityFloatArithmeticTest, TestAlmostEqualRelative) {
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, 3.0 + 1e-11));
 
     // Checking the the sensitivity towards the next floating point
-    // The ULP distnace is really small, this method will always return true
+    // A ULP distance of 4 or 5 is still really small than our relative sensitivity of 1e-10
     const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0);
     const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0), 4.0);
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fourHops));

From 6c5f6c12a54792e0ef319357f93fd194eb9bfe39 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 15:14:17 +0200
Subject: [PATCH 05/10] update the UML component diagram to match version 3.x

---
 paper/figures/PolyhedralGravityModel.png          |   3 ---
 ...nt_Diagram_Polyhedral_Gravity_Model.drawio.pdf | Bin 0 -> 130 bytes
 paper/paper.md                                    |   2 +-
 3 files changed, 1 insertion(+), 4 deletions(-)
 delete mode 100644 paper/figures/PolyhedralGravityModel.png
 create mode 100644 paper/figures/UML_Component_Diagram_Polyhedral_Gravity_Model.drawio.pdf

diff --git a/paper/figures/PolyhedralGravityModel.png b/paper/figures/PolyhedralGravityModel.png
deleted file mode 100644
index 80685f2..0000000
--- a/paper/figures/PolyhedralGravityModel.png
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:6399ec6bcb2f1183eb2907f7a73d0ede81749e042eb9cd565f9d631f916e363a
-size 238730
diff --git a/paper/figures/UML_Component_Diagram_Polyhedral_Gravity_Model.drawio.pdf b/paper/figures/UML_Component_Diagram_Polyhedral_Gravity_Model.drawio.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..8c1a6a33ceefa5ac9711cef08da1a483fe2444e8
GIT binary patch
literal 130
zcmWN?!4bkB5CFhGRnUL|E^t7)0qziHR5F5jSiSCNFMiKHUb3xq&O@ntU$;k{+yC~-
zTNzI^PcG^*Vsw(TC2)8c$|gAH2x|h~p<>V}NAHoM$%V7^5j7>W<O~)Fnka&00!_Q2
MDrWz(gn$mk59_!mod5s;

literal 0
HcmV?d00001

diff --git a/paper/paper.md b/paper/paper.md
index 235c717..ca646b6 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -67,7 +67,7 @@ On a mathematical level, the implemented model follows the line integral approac
 
 Implementation-wise, it makes use of the inherent parallelization opportunity of the approach as it iterates over the mesh elements. This parallelization is achieved via *thrust*, which allows utilizing *OpenMP* and *Intel TBB*. On a finer scale, individual, costly operations have been investigated, and, e.g., the \texttt{arctan} operations have been vectorized using *xsimd*. On the application side, the user can choose between the functional interface for evaluating the full gravity tensor or the object-oriented \texttt{GravityEvaluable}, providing the same functionality while implementing a caching mechanism to avoid recomputing mesh properties that can be shared between multipoint evaluation, such as the face normals.
 
-![UML Component Diagram of the implementation. External dependencies are depicted in blue. Internal components are colored in grey.\label{fig:implementation}](figures/PolyhedralGravityModel.png)
+![UML Component Diagram of the implementation. External dependencies are depicted in gray. Components of the polyhedral gravity model are colored in blue and green.\label{fig:implementation}](figures/UML_Component_Diagram_Polyhedral_Gravity_Model.drawio.pdf)
 
 Extensive tests using GoogleTest for the C++ side and pytest for the Python interface are employed via GitHub Actions to ensure the (continued) correctness of the implementation.
 

From 25cd0d01c037c2f7b7ab2e4c55687039bb8a206b Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 15:41:48 +0200
Subject: [PATCH 06/10] update comments

---
 .../model/GravityModelData.h                  | 24 +++++++++----------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/polyhedralGravity/model/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h
index 2d3abb5..65dee5c 100644
--- a/src/polyhedralGravity/model/GravityModelData.h
+++ b/src/polyhedralGravity/model/GravityModelData.h
@@ -76,8 +76,8 @@ namespace polyhedralGravity {
         double s2;
 
         /**
-         * Checks two Distance structs for equality with another one by comparing their members values'
-         * ULP distances.
+         * Checks two Distance structs for equality with another one by ensuring that the members are
+         * almost equal.
          * @param rhs the other Distance struct
          * @return true if equal
          *
@@ -91,8 +91,8 @@ namespace polyhedralGravity {
         }
 
         /**
-         * Checks two Distance structs for inequality with another one by comparing their members values'
-         * ULP distances.
+         * Checks two Distance structs for inequality with another one by ensuring that the members are
+         * not almost equal.
          * @param rhs the other Distance struct
          * @return false if unequal
          *
@@ -133,8 +133,8 @@ namespace polyhedralGravity {
         double an;
 
         /**
-        * Checks two TranscendentalExpressions for equality with another one by comparing their members values'
-         * ULP distances.
+         * Checks two TranscendentalExpressions for equality with another one by ensuring that the members are
+         * almost equal.
          * @param rhs the other TranscendentalExpressions
          * @return true if equal
          *
@@ -145,8 +145,8 @@ namespace polyhedralGravity {
         }
 
         /**
-         * Checks two TranscendentalExpressions for inequality with another one by comparing their members values
-         * ULP distances.
+         * Checks two TranscendentalExpressions for inequality with another one by ensuring that the members are
+         * not almost equal.
          * @param rhs the other TranscendentalExpressions
          * @return false if unequal
          *
@@ -194,8 +194,8 @@ namespace polyhedralGravity {
         double d;
 
         /**
-         * Checking the equality of two this Hessian Plane with another one by comparing their members values
-         * ULP distances.
+         * Checking the equality of two this Hessian Plane with another one by ensuring that the members are
+         * almost equal.
          * @param rhs other HessianPlane
          * @return true if equal
          *
@@ -209,8 +209,8 @@ namespace polyhedralGravity {
         }
 
         /**
-         * Checking the inequality of two this Hessian Plane with another one with another one by comparing their members values'
-         * ULP distances.
+         * Checking the inequality of two this Hessian Plane with another one by ensuring that the members are
+         * not almost equal.
          * @param rhs other HessianPlane
          * @return true if unequal
          *

From f4ce1675c00773330fe9c3ce6a8ed4912efc6be1 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 17:08:32 +0200
Subject: [PATCH 07/10] Update
 src/polyhedralGravity/util/UtilityFloatArithmetic.h
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Pablo Gómez <contact@pablo-gomez.net>
---
 src/polyhedralGravity/util/UtilityFloatArithmetic.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.h b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
index 60507c7..589f022 100644
--- a/src/polyhedralGravity/util/UtilityFloatArithmetic.h
+++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.h
@@ -8,7 +8,7 @@
 namespace polyhedralGravity::util {
 
     /**
-     * The EPSILON used in the polyhedral gravity model to determine a radius around zero/ to use as sligth offset.
+     * The EPSILON used in the polyhedral gravity model to determine a radius around zero/ to use as slight offset.
      * @related Used to determine if a floating point number is equal to zero as threshold for rounding errors
      * @related Used for the sgn() function to determine the sign of a double value. Different compilers
      * produce different results if no EPSILON is applied for the comparison!

From 5b7d79745d2cd99b613b6d714cb5a9087e4450ee Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 17:13:55 +0200
Subject: [PATCH 08/10] use infinity

---
 test/util/UtilityFloatArithmeticTest.cpp | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/test/util/UtilityFloatArithmeticTest.cpp b/test/util/UtilityFloatArithmeticTest.cpp
index 4dc28eb..25853a7 100644
--- a/test/util/UtilityFloatArithmeticTest.cpp
+++ b/test/util/UtilityFloatArithmeticTest.cpp
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <limits>
 
 #include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
@@ -32,8 +33,9 @@ TEST(UtilityFloatArithmeticTest, TestAlmostEqualUlps) {
     // Checking the the sensitivity towards the next floating point
     // Note: The default maximal ULPS distance is set to 4
     // The ULP distance is one time higher than 4 leading to false, one time lower equal leading to true
-    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0);
-    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0), 4.0);
+    const auto inf = std::numeric_limits<double>::infinity();
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf);
+    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf), inf);
     ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(3.0, fourHops));
     ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, fiveHops));
 
@@ -65,8 +67,9 @@ TEST(UtilityFloatArithmeticTest, TestAlmostEqualRelative) {
 
     // Checking the the sensitivity towards the next floating point
     // A ULP distance of 4 or 5 is still really small than our relative sensitivity of 1e-10
-    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0);
-    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, 4.0), 4.0), 4.0), 4.0), 4.0);
+    const auto inf = std::numeric_limits<double>::infinity();
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf);
+    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf), inf);
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fourHops));
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fiveHops));
 }

From 675a1e49b8f34e635aede0c8bd45b830f9f3e442 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 17:23:14 +0200
Subject: [PATCH 09/10] add statement of need

---
 test/model/GoogleTestMatcher.h | 34 +++++++++++++++++++++++++++-------
 1 file changed, 27 insertions(+), 7 deletions(-)

diff --git a/test/model/GoogleTestMatcher.h b/test/model/GoogleTestMatcher.h
index cd7a294..0f6cb23 100644
--- a/test/model/GoogleTestMatcher.h
+++ b/test/model/GoogleTestMatcher.h
@@ -2,9 +2,17 @@
 
 #include "gtest/gtest.h"
 #include "gmock/gmock.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
+/*
+ * This file provides several matchers to be used with ASSERT_TAHT(..) statements for multi-dimensional containers
+ * with flotaing points as content.
+ * GoogleTest provides matchers like DoubleNear(..), ContainerEq(..) and Pointwise(..). These work great when working
+ * with 1D containers. However, there is a lack of these matchers for multi-dimensional containers.
+ * This files provides Matcher for comparing nested containers with floating points where the comparison is
+ * conducted with an epsilon to counter the effect of potential rounding erros due to floating point arithmetic.
+ */
 
-#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
 MATCHER_P(FloatContainter1D, container, "Comparing 1D Containers") {
     if (container.size() != arg.size()) {
@@ -13,7 +21,9 @@ MATCHER_P(FloatContainter1D, container, "Comparing 1D Containers") {
     }
     for (size_t idx = 0; idx < std::size(container); ++idx) {
         if (!polyhedralGravity::util::almostEqualRelative(container[idx], arg[idx])) {
-            *result_listener << "The elements at idx = " << idx << " do not match. Values: " << container[idx]<< " != " << arg[idx];
+            *result_listener
+                    << "The elements at idx = " << idx << " do not match. Values: "
+                    << container[idx] << " != " << arg[idx];
             return false;
         }
     }
@@ -29,12 +39,16 @@ MATCHER_P(FloatContainter2D, container, "Comparing 2D Containers") {
     }
     for (size_t idx = 0; idx < std::size(container); ++idx) {
         if (container[idx].size() != arg[idx].size()) {
-            *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+            *result_listener
+                    << "The container sizes at idx = " << idx << " do not match. Sizes: "
+                    << container[idx].size() << " != " << arg[idx].size();
             return false;
         }
         for (size_t idy = 0; idy < std::size(container[idx]); ++idy) {
             if (!polyhedralGravity::util::almostEqualRelative(container[idx][idy], arg[idx][idy])) {
-                *result_listener << "The elements at idx = " << idx << " and idy = " << idy << " do not match. Values: " << container[idx][idy] << " != " << arg[idx][idy];
+                *result_listener
+                        << "The elements at idx = " << idx << " and idy = " << idy << " do not match. Values: "
+                        << container[idx][idy] << " != " << arg[idx][idy];
                 return false;
             }
         }
@@ -51,17 +65,23 @@ MATCHER_P(FloatContainter3D, container, "Comparing 3D Containers") {
     }
     for (size_t idx = 0; idx < std::size(container); ++idx) {
         if (container[idx].size() != arg[idx].size()) {
-            *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+            *result_listener
+                    << "The container sizes at idx = " << idx << " do not match. Sizes: "
+                    << container[idx].size() << " != " << arg[idx].size();
             return false;
         }
         for (size_t idy = 0; idy < std::size(container[idx]); ++idy) {
             if (container[idx][idy].size() != arg[idx][idy].size()) {
-                *result_listener << "The container sizes at idx = " << idx << " do not match. Sizes: " << container[idx].size() << " != " << arg[idx].size();
+                *result_listener
+                        << "The container sizes at idx = " << idx << " do not match. Sizes: "
+                        << container[idx].size() << " != " << arg[idx].size();
                 return false;
             }
             for (size_t idz = 0; idz < std::size(container[idx][idy]); ++idz) {
                 if (!polyhedralGravity::util::almostEqualRelative(container[idx][idy][idz], arg[idx][idy][idz])) {
-                    *result_listener << "The elements at idx = " << idx << " and idy = " << idy << " and idz = " << idz << " do not match. Values: " << container[idx][idy][idz] << " != " << arg[idx][idy][idz];
+                    *result_listener
+                            << "The elements at idx = " << idx << " and idy = " << idy << " and idz = " << idz << " do not match. Values: "
+                            << container[idx][idy][idz] << " != " << arg[idx][idy][idz];
                     return false;
                 }
             }

From 69126932a1617ae23d064d765046328b4398c2d8 Mon Sep 17 00:00:00 2001
From: Jonas Schuhmacher <jonas.schuhmacher@tum.de>
Date: Wed, 8 May 2024 17:27:27 +0200
Subject: [PATCH 10/10] refactor constant

---
 test/util/UtilityFloatArithmeticTest.cpp | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/test/util/UtilityFloatArithmeticTest.cpp b/test/util/UtilityFloatArithmeticTest.cpp
index 25853a7..cdf08e5 100644
--- a/test/util/UtilityFloatArithmeticTest.cpp
+++ b/test/util/UtilityFloatArithmeticTest.cpp
@@ -7,6 +7,8 @@
 
 #include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
+/** Contains the positive infinity values for doubles */
+constexpr static double INF = std::numeric_limits<double>::infinity();
 
 TEST(UtilityFloatArithmeticTest, TestAlmostEqualUlps) {
     // Checking the signess && identity
@@ -33,9 +35,9 @@ TEST(UtilityFloatArithmeticTest, TestAlmostEqualUlps) {
     // Checking the the sensitivity towards the next floating point
     // Note: The default maximal ULPS distance is set to 4
     // The ULP distance is one time higher than 4 leading to false, one time lower equal leading to true
-    const auto inf = std::numeric_limits<double>::infinity();
-    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf);
-    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf), inf);
+    // Can be constexpr starting with C++23
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, INF), INF), INF), INF);
+    const double fiveHops = std::nextafter(fourHops, INF);
     ASSERT_TRUE(polyhedralGravity::util::almostEqualUlps(3.0, fourHops));
     ASSERT_FALSE(polyhedralGravity::util::almostEqualUlps(3.0, fiveHops));
 
@@ -67,9 +69,9 @@ TEST(UtilityFloatArithmeticTest, TestAlmostEqualRelative) {
 
     // Checking the the sensitivity towards the next floating point
     // A ULP distance of 4 or 5 is still really small than our relative sensitivity of 1e-10
-    const auto inf = std::numeric_limits<double>::infinity();
-    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf);
-    const double fiveHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, inf), inf), inf), inf), inf);
+    // Can be constexpr starting with C++23
+    const double fourHops = std::nextafter(std::nextafter(std::nextafter(std::nextafter(3.0, INF), INF), INF), INF);
+    const double fiveHops = std::nextafter(fourHops, INF);
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fourHops));
     ASSERT_TRUE(polyhedralGravity::util::almostEqualRelative(3.0, fiveHops));
 }