Merge pull request #39 from esa/stable-cpp-tests

Improve Test Stability & Minor fixes to documentation

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
 -------------

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)

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.
 

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="[email protected]",
     description="Package to compute full gravity tensor of a given constant density polyhedron for arbitrary points "

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 ensuring that the members are
+         * almost equal.
          * @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::almostEqualRelative(l1, rhs.l1) &&
+                   util::almostEqualRelative(l2, rhs.l2) &&
+                   util::almostEqualRelative(s1, rhs.s1) &&
+                   util::almostEqualRelative(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 ensuring that the members are
+         * not almost equal.
          * @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 ensuring that the members are
+         * almost equal.
          * @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::almostEqualRelative(ln, rhs.ln) && util::almostEqualRelative(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 ensuring that the members are
+         * not almost equal.
          * @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 ensuring that the members are
+         * almost equal.
          * @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::almostEqualRelative(a, rhs.a) &&
+                   util::almostEqualRelative(b, rhs.b) &&
+                   util::almostEqualRelative(c, rhs.c) &&
+                   util::almostEqualRelative(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 by ensuring that the members are
+         * not almost equal.
          * @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;
+        }
     };
 
 }

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

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 {

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;

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 {
 
@@ -150,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,
@@ -169,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,
@@ -187,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,
@@ -203,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,

src/polyhedralGravity/util/UtilityConstants.h

@@ -2,15 +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
      */