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 config = std::make_shared("config.yaml"); - Polyhedron poly = config->getDataSource()->getPolyhedron(); + auto polyhedralSource = config->getDataSource()->getPolyhedralSource(); double density = config->getDensity(); std::array 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/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 0000000..8c1a6a3 Binary files /dev/null and b/paper/figures/UML_Component_Diagram_Polyhedral_Gravity_Model.drawio.pdf differ 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. 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/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h index 0f61aaa..65dee5c 100644 --- a/src/polyhedralGravity/model/GravityModelData.h +++ b/src/polyhedralGravity/model/GravityModelData.h @@ -6,7 +6,7 @@ #include #include #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; + } }; } 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(rayOrigin + rayVector * t); } else { return nullptr; diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h index c80ed30..ef74b7a 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 { @@ -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 &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, double density, const NormalOrientation &orientation = NormalOrientation::OUTWARDS, diff --git a/src/polyhedralGravity/util/UtilityConstants.h b/src/polyhedralGravity/util/UtilityConstants.h index f9ab381..4292850 100644 --- a/src/polyhedralGravity/util/UtilityConstants.h +++ b/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 */ diff --git a/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp new file mode 100644 index 0000000..c1c60ec --- /dev/null +++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.cpp @@ -0,0 +1,51 @@ +#include "UtilityFloatArithmetic.h" + + +namespace polyhedralGravity::util { + + template + bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance) { + static_assert(std::is_same_v || std::is_same_v, + "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; + } + + // In case the signs mismatch, return false + if (lhs < static_cast(0.0) && rhs > static_cast(0.0) || + lhs > static_cast(0.0) && rhs < static_cast(0.0)) { + return false; + } + + if constexpr (std::is_same_v) { + // In case of float, compute ULP distance by interpreting float as 32-bit integer + return reinterpret_cast(rhs) - reinterpret_cast(lhs) <= ulpDistance; + } + else if constexpr (std::is_same_v) { + // In case of double, compute ULP distance by interpreting double as 64-bit integer + return reinterpret_cast(rhs) - reinterpret_cast(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 lhs, float rhs, int ulpDistance); + template bool almostEqualUlps(double lhs, double rhs, int ulpDistance); + + + template + 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 lhs, float rhs, double epsilon); + template bool almostEqualRelative(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 new file mode 100644 index 0000000..589f022 --- /dev/null +++ b/src/polyhedralGravity/util/UtilityFloatArithmetic.h @@ -0,0 +1,73 @@ +#pragma once + +#include +#include +#include +#include + +namespace polyhedralGravity::util { + + /** + * 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! + */ + 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 {@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. + * 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 + * {@link polyhedralGravity::util::almostEqualUlps} function. + * + * @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. + * + * @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 {@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. + * @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/ + */ + template + bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance = MAX_ULP_DISTANCE); + + /** + * Function to check if two floating point numbers are relatively equal to each other within a given error range or tolerance. + * + * @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 {@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`. + * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ + */ + template + bool almostEqualRelative(FloatType lhs, FloatType rhs, double epsilon = EPSILON_ALMOST_EQUAL); + + +} 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( diff --git a/test/model/GoogleTestMatcher.h b/test/model/GoogleTestMatcher.h new file mode 100644 index 0000000..0f6cb23 --- /dev/null +++ b/test/model/GoogleTestMatcher.h @@ -0,0 +1,91 @@ +#pragma once + +#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. + */ + + +MATCHER_P(FloatContainter1D, container, "Comparing 1D 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 (!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 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(); + 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 (!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; +} + +// Using the FloatContainter2D would be nice, but this leads to issues with template instatantation +// 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(); + 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 e80285d..9f6b21d 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,16 +269,15 @@ TEST_F(GravityModelBigTest, PlaneUnitNormals) { auto actualPlaneUnitNormals = polyhedralGravity::GravityModel::calculatePlaneUnitNormals(expectedGij); - ASSERT_THAT(actualPlaneUnitNormals, ContainerEq(expectedPlaneUnitNormals)); + 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) { @@ -291,12 +291,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) { @@ -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 new file mode 100644 index 0000000..cdf08e5 --- /dev/null +++ b/test/util/UtilityFloatArithmeticTest.cpp @@ -0,0 +1,77 @@ +#include "gtest/gtest.h" +#include "gmock/gmock.h" + +#include +#include +#include + +#include "polyhedralGravity/util/UtilityFloatArithmetic.h" + +/** Contains the positive infinity values for doubles */ +constexpr static double INF = std::numeric_limits::infinity(); + +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 + // 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)); + + +} + +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 + // A ULP distance of 4 or 5 is still really small than our relative sensitivity of 1e-10 + // 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)); +}