Changes supporting GPU shaping branch

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);
+}

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 <typename T>
+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

src/axom/klee/Geometry.cpp

@@ -154,7 +154,7 @@ void Geometry::populateGeomInfo()
     m_discreteFunction = axom::Array<double, 2>(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);

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.
    */

src/axom/primal/geometry/CoordinateTransformer.hpp

@@ -66,6 +66,17 @@ class CoordinateTransformer
    */
   CoordinateTransformer(const numerics::Matrix<T>& matrix) { setMatrix(matrix); }
 
+  /*!
+   * @brief Construct transformer that moves 4 starting points to 4
+   * destination points.
+   * @see setByTerminusPts.
+   */
+  AXOM_HOST_DEVICE CoordinateTransformer(const primal::Point<T, 3>* startPts,
+                                         const primal::Point<T, 3>* 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<T, 3>* startPts,
                                          const primal::Point<T, 3>* destPts)
@@ -206,18 +219,17 @@ class CoordinateTransformer
    */
   void applyRotation(const axom::primal::Vector<T, 3>& start, const axom::primal::Vector<T, 3>& end)
   {
-    // Note that the rotation matrix is not unique.
     Vector<T, 3> s = start.unitVector();
     Vector<T, 3> e = end.unitVector();
     Vector<T, 3> u;  // Rotation vector, the cross product of start and end.
     numerics::cross_product(s.data(), e.data(), u.data());
     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

src/axom/primal/geometry/Sphere.hpp

@@ -178,6 +178,19 @@ class Sphere
   AXOM_HOST_DEVICE
   inline bool intersectsWith(const Sphere<T, NDIMS>& 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<T, NDIMS>& 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<T, NDIMS>::intersectsWith(const Sphere<T, ND
           utilities::isNearlyEqual(distance_squared, sum_of_radii_2, TOL));
 }
 
+//------------------------------------------------------------------------------
+template <typename T, int NDIMS>
+AXOM_HOST_DEVICE inline bool Sphere<T, NDIMS>::contains(const Sphere<T, NDIMS>& 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
 //------------------------------------------------------------------------------

src/axom/primal/tests/primal_sphere.cpp

@@ -150,6 +150,34 @@ void check_sphere_intersection()
   EXPECT_FALSE(S0.intersectsWith(S3));
 }
 
+//------------------------------------------------------------------------------
+template <int NDIMS>
+void check_sphere_containment()
+{
+  using PointType = primal::Point<double, NDIMS>;
+  using SphereType = primal::Sphere<double, NDIMS>;
+
+  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 <int NDIMS>
 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[])
 {