diff --git a/src/axom/core/tests/core_utilities.hpp b/src/axom/core/tests/core_utilities.hpp index 3b17d60c85..bb82e7560e 100644 --- a/src/axom/core/tests/core_utilities.hpp +++ b/src/axom/core/tests/core_utilities.hpp @@ -207,3 +207,24 @@ TEST(core_utilities, sort_multiple_values) } } } + +TEST(core_utilities, sign_of) +{ + using axom::utilities::sign_of; + EXPECT_EQ(sign_of(1234), 1); + EXPECT_EQ(sign_of(-1234), -1); + EXPECT_EQ(sign_of(0), 0); + + EXPECT_EQ(sign_of(0.1234), 1); + EXPECT_EQ(sign_of(-0.1234), -1); + EXPECT_EQ(sign_of(0.0), 0); + + EXPECT_EQ(sign_of(1e-9), 1); + EXPECT_EQ(sign_of(-1e-9), -1); + + EXPECT_EQ(sign_of(1e-9, 1e-12), 1); + EXPECT_EQ(sign_of(-1e-9, 1e-12), -1); + + EXPECT_EQ(sign_of(1e-9, 1e-8), 0); + EXPECT_EQ(sign_of(-1e-9, 1e-8), 0); +} diff --git a/src/axom/core/utilities/Utilities.hpp b/src/axom/core/utilities/Utilities.hpp index df17f8d362..6e868a34b4 100644 --- a/src/axom/core/utilities/Utilities.hpp +++ b/src/axom/core/utilities/Utilities.hpp @@ -348,6 +348,18 @@ inline AXOM_HOST_DEVICE bool isNearlyEqualRelative(RealType a, // return abs(a-b) <= max(absThresh, relThresh * maxFabs ); } +/*! + * \brief Sign of a value of any type that supports comparison and + * negation operators. + * \return 0 for v in [-eps, eps], -1 for v < eps and +1 for v > eps + */ +template +inline int sign_of(const T& v, const T& eps = {0}) +{ + assert(eps >= 0); + return v > eps ? 1 : v < -eps ? -1 : 0; +} + /*! * \brief Insertion sort of an array. * \accelerated diff --git a/src/axom/klee/Geometry.cpp b/src/axom/klee/Geometry.cpp index 09c8b03e6e..ff8a0f8c06 100644 --- a/src/axom/klee/Geometry.cpp +++ b/src/axom/klee/Geometry.cpp @@ -154,7 +154,7 @@ void Geometry::populateGeomInfo() m_discreteFunction = axom::Array(2, 2); m_discreteFunction(0, 0) = 0.0; m_discreteFunction(0, 1) = cone.getBaseRadius(); - m_discreteFunction(1, 1) = cone.getLength(); + m_discreteFunction(1, 0) = cone.getLength(); m_discreteFunction(1, 1) = cone.getTopRadius(); m_geomInfo["discreteFunction"].set(m_discreteFunction.data(), m_discreteFunction.size()); m_geomInfo["sorOrigin"].set(cone.getBaseCenter().data(), 3); diff --git a/src/axom/primal/geometry/Cone.hpp b/src/axom/primal/geometry/Cone.hpp index 8eb2e827d4..d7ef33ce74 100644 --- a/src/axom/primal/geometry/Cone.hpp +++ b/src/axom/primal/geometry/Cone.hpp @@ -143,7 +143,8 @@ class Cone } /*! - * \brief Returns the algebraic volume of the cone + * \brief Returns the signed volume of the cone, + * which is negative if the length is negative. * * Volume is only defined when NDIMS == 3. */ diff --git a/src/axom/primal/geometry/CoordinateTransformer.hpp b/src/axom/primal/geometry/CoordinateTransformer.hpp index ef983c28f7..0405e0cd9d 100644 --- a/src/axom/primal/geometry/CoordinateTransformer.hpp +++ b/src/axom/primal/geometry/CoordinateTransformer.hpp @@ -66,6 +66,17 @@ class CoordinateTransformer */ CoordinateTransformer(const numerics::Matrix& matrix) { setMatrix(matrix); } + /*! + * @brief Construct transformer that moves 4 starting points to 4 + * destination points. + * @see setByTerminusPts. + */ + AXOM_HOST_DEVICE CoordinateTransformer(const primal::Point* startPts, + const primal::Point* destPts) + { + setByTerminusPts(startPts, destPts); + } + /*! * @brief Set the matrix, discarding the current transformation. * @param matrix [in] The transformation matrix for homogeneous @@ -100,8 +111,10 @@ class CoordinateTransformer * @param [in] destPts Four destination points. * * The four starting points must define a non-degenerate volume. - * Else the transformation is ill-defined and the transformer - * is set to invalid. + * Else the transformation is ill-defined and the transformer is set + * to invalid. If you think of each 4 points as the vertices of a + * tetrahedron, the transformation moves the starting tet to the + * destination tet. */ AXOM_HOST_DEVICE void setByTerminusPts(const primal::Point* startPts, const primal::Point* destPts) @@ -206,7 +219,6 @@ class CoordinateTransformer */ void applyRotation(const axom::primal::Vector& start, const axom::primal::Vector& end) { - // Note that the rotation matrix is not unique. Vector s = start.unitVector(); Vector e = end.unitVector(); Vector u; // Rotation vector, the cross product of start and end. @@ -214,10 +226,10 @@ class CoordinateTransformer const T sinT = u.norm(); const T cosT = numerics::dot_product(s.data(), e.data(), 3); - // Degenerate: end is parallel to start. - // angle near 0 (identity transform) or pi. if(utilities::isNearlyEqual(sinT, 0.0)) { + // Degenerate: end is parallel to start. + // angle near 0 (identity transform) or pi. if(cosT < 0) { setInvalid(); // Transformation is ill-defined diff --git a/src/axom/primal/geometry/Sphere.hpp b/src/axom/primal/geometry/Sphere.hpp index 9f89e37a32..d8eef6eb74 100644 --- a/src/axom/primal/geometry/Sphere.hpp +++ b/src/axom/primal/geometry/Sphere.hpp @@ -178,6 +178,19 @@ class Sphere AXOM_HOST_DEVICE inline bool intersectsWith(const Sphere& sphere, double TOL = 1.e-9) const; + /*! + * \brief Tests if this sphere completely contains another sphere. + * + * \param [in] other The sphere object to check for containment + * \param [in] margin Amount to add to other's radius before comparing. + * + * Note: a sphere does contain itself. + * + * \return true if this sphere contains the other, false otherwise. + */ + AXOM_HOST_DEVICE + inline bool contains(const Sphere& other, double margin = 0.0) const; + /*! * \brief Prints the Sphere information in the given output stream. * \param [in,out] os the output stream to write to. @@ -233,6 +246,15 @@ AXOM_HOST_DEVICE inline bool Sphere::intersectsWith(const Sphere +AXOM_HOST_DEVICE inline bool Sphere::contains(const Sphere& other, + double margin) const +{ + const T center_sep = VectorType(other.getCenter(), m_center).norm(); + return (m_radius >= center_sep + other.getRadius() + margin); +} + //------------------------------------------------------------------------------ // implementation of free functions //------------------------------------------------------------------------------ diff --git a/src/axom/primal/tests/primal_sphere.cpp b/src/axom/primal/tests/primal_sphere.cpp index 3777b75f80..4313263373 100644 --- a/src/axom/primal/tests/primal_sphere.cpp +++ b/src/axom/primal/tests/primal_sphere.cpp @@ -150,6 +150,34 @@ void check_sphere_intersection() EXPECT_FALSE(S0.intersectsWith(S3)); } +//------------------------------------------------------------------------------ +template +void check_sphere_containment() +{ + using PointType = primal::Point; + using SphereType = primal::Sphere; + + PointType center {3.0, 4.0, 0.0}; + double radius = 5.0; + + SphereType S0(center, radius); + SphereType S1(PointType {1.5, 2.0, 0.0}, radius / 2); + SphereType S2(S1.getCenter(), radius); + + double tol = 1e-12; + + // STEP 0: test fully containing + EXPECT_TRUE(S0.contains(S0, -tol)); + EXPECT_TRUE(S0.contains(S1, -tol)); + + // STEP 1: test barely not containing. + EXPECT_FALSE(S0.contains(S0, tol)); + EXPECT_FALSE(S0.contains(S1, tol)); + + // STEP 2: test partial containment. + EXPECT_FALSE(S0.contains(S2)); +} + //------------------------------------------------------------------------------ template void check_copy_constructor() @@ -244,6 +272,13 @@ TEST(primal_sphere, sphere_sphere_intersection) check_sphere_intersection<3>(); } +//------------------------------------------------------------------------------ +TEST(primal_sphere, sphere_sphere_containment) +{ + check_sphere_containment<2>(); + check_sphere_containment<3>(); +} + //------------------------------------------------------------------------------ int main(int argc, char* argv[]) {