Improve clipping code performance

src/axom/core/utilities/Utilities.hpp

@@ -348,6 +348,16 @@ 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.
+ */
+template <typename T>
+inline int sign_of(const T& v, const T& 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 algebraic 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,16 @@ class CoordinateTransformer
    */
   CoordinateTransformer(const numerics::Matrix<T>& matrix) { setMatrix(matrix); }
 
+  /*!
+   * @brief Contruct transformer that moves 4 starting points to 4
+   * destination points.
+   */
+  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

src/axom/primal/geometry/Sphere.hpp

@@ -178,6 +178,18 @@ 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 that this sphere must contain the other sphere by.
+   *   Positive means that the other sphere is "more inside".
+   *
+   * \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 +245,14 @@ 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>& sphere, double TOL) const
+{
+  const T center_sep = VectorType(sphere.getCenter(), m_center).norm();
+  return (m_radius > center_sep + sphere.getRadius() + TOL);
+}
+
 //------------------------------------------------------------------------------
 //  implementation of free functions
 //------------------------------------------------------------------------------

src/axom/primal/operators/detail/clip_impl.hpp

@@ -233,7 +233,7 @@ AXOM_HOST_DEVICE void poly_clip_vertices(Polyhedron<T, NDIMS>& poly,
           // A probable AMD GPU 6.3 compiler bug caused the first assert to fail.
           // The second assert is equivalent because we expect addVertex call to
           // consume the next index and increment the number of vertices.
-          // However, when the compilers no longar cause false failures, we should
+          // However, when the compilers no longer cause false failures, we should
           // restore the first assert.
           // SLIC_ASSERT(newVertexIndex == expectedVertexIndex);
           SLIC_ASSERT(newVertexIndex == poly.numVertices() - 1);

src/axom/primal/tests/primal_sphere.cpp

@@ -150,6 +150,31 @@ 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 {0.0, 0.0, 0.0};
+  double tol = 1e-12;
+
+  // STEP 0: test fully containing
+  SphereType S0;
+  EXPECT_TRUE(S0.contains(S0, -tol));
+
+  // STEP 1: test barely not containing.
+  center[0] = tol;
+  SphereType S1(center);
+  EXPECT_FALSE(S0.contains(S1));
+
+  // STEP 2: test partial containment.
+  center[0] = 0.5;
+  SphereType S3(center);
+  EXPECT_FALSE(S0.contains(S3));
+}
+
 //------------------------------------------------------------------------------
 template <int NDIMS>
 void check_copy_constructor()
@@ -244,6 +269,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[])
 {

src/axom/quest/CMakeLists.txt

@@ -60,6 +60,7 @@ set( quest_headers
     interface/signed_distance.hpp
 
     util/mesh_helpers.hpp
+    util/make_clipper_strategy.hpp
     MeshViewUtil.hpp
     )
 
@@ -90,6 +91,7 @@ set( quest_sources
     interface/signed_distance.cpp
 
     util/mesh_helpers.cpp
+    util/make_clipper_strategy.cpp
    )
 
 
@@ -100,6 +102,7 @@ set( quest_depends_on
     )
 
 blt_list_append(TO quest_depends_on IF AXOM_ENABLE_SIDRE ELEMENTS sidre)
+blt_list_append(TO quest_depends_on IF AXOM_ENABLE_BUMP ELEMENTS bump)
 blt_list_append(TO quest_depends_on IF SPARSEHASH_FOUND ELEMENTS sparsehash)
 blt_list_append(TO quest_depends_on IF MFEM_FOUND ELEMENTS mfem)
 
@@ -150,11 +153,24 @@ if(AXOM_ENABLE_KLEE AND AXOM_ENABLE_SIDRE)
         if(RAJA_FOUND)
             list(APPEND quest_headers MeshClipper.hpp
                                       MeshClipperStrategy.hpp
+                                      detail/clipping/Plane3DClipper.hpp
                                       detail/clipping/TetClipper.hpp
+                                      detail/clipping/TetMeshClipper.hpp
+                                      detail/clipping/HexClipper.hpp
+                                      detail/clipping/SphereClipper.hpp
+                                      detail/clipping/FSorClipper.hpp
+                                      detail/clipping/SorClipper.hpp
                                       detail/clipping/MeshClipperImpl.hpp)
             list(APPEND quest_sources MeshClipper.cpp
+                                      MeshClipperStrategy.cpp
+                                      detail/clipping/Plane3DClipper.cpp
                                       detail/clipping/TetClipper.cpp
-                                      MeshClipperStrategy.cpp)
+                                      detail/clipping/TetMeshClipper.cpp
+                                      detail/clipping/HexClipper.cpp
+                                      detail/clipping/SphereClipper.cpp
+                                      detail/clipping/FSorClipper.cpp
+                                      detail/clipping/SorClipper.cpp
+                                      detail/clipping/MeshClipperImpl.hpp)
         endif()
     endif()
 

src/axom/quest/Discretize.cpp

@@ -108,6 +108,8 @@ OctType new_inscribed_oct(const SphereType& sphere, OctType& o, int s, int t, in
  *
  * This routine allocates an array pointed to by \a out.  The caller is responsible
  * to free the array.
+ *
+ * TODO: If possible, port to GPU (rewrite to be data-parallel).
  */
 bool discretize(const SphereType& sphere, int levels, axom::Array<OctType>& out, int& octcount)
 {
@@ -125,7 +127,7 @@ bool discretize(const SphereType& sphere, int levels, axom::Array<OctType>& out,
 
   octcount = count_sphere_octahedra(levels);
 
-  out = axom::Array<OctType>(octcount, octcount);
+  out.resize(octcount);
 
   // index points to an octahedron of the last generation.  We'll generate
   // new octahedra based on out[index].

src/axom/quest/Discretize.hpp

@@ -50,7 +50,8 @@ bool discretize(const SphereType& s, int levels, axom::Array<OctType>& out, int&
 /*!
  * \brief Given a 2D polyline revolved around the positive X-axis, allocate
  *   and return a list of Octahedra approximating the shape.
- * \param [in] polyline The polyline to revolve around the X-axis
+ * \param [in] polyline The polyline to revolve around the X-axis.
+ *   Data should be in a host-accessible memory space.
  * \param [in] len The number of points in \a polyline
  * \param [in] levels The number of refinements to perform, in addition to
  *   a central level-zero octahedron in each segment

src/axom/quest/InOutOctree.hpp

@@ -17,7 +17,6 @@
 #include "axom/slic.hpp"
 #include "axom/slam.hpp"
 #include "axom/primal.hpp"
-#include "axom/mint.hpp"
 #include "axom/spin.hpp"
 
 #include "detail/inout/BlockData.hpp"

src/axom/quest/IntersectionShaper.hpp

@@ -825,6 +825,7 @@ class IntersectionShaper : public Shaper
 
       // Generate the Octahedra
       // (octahedra m_octs will be on device)
+      m_octs = axom::Array<OctahedronType>(0, 0, axom::execution_space<ExecSpace>::allocatorID());
       const bool disc_status =
         axom::quest::discretize<ExecSpace>(polyline, polyline_size, m_level, m_octs, m_octcount);
 
@@ -1967,13 +1968,15 @@ class IntersectionShaper : public Shaper
     if(m_bpGrp)
     {
       auto fieldsGrp = m_bpGrp->getGroup("fields");
-      SLIC_ERROR_IF(fieldsGrp == nullptr, "Input blueprint mesh lacks the 'fields' Group/Node.");
-      for(auto& group : fieldsGrp->groups())
+      if(fieldsGrp != nullptr)
       {
-        std::string materialName = fieldNameToMaterialName(group.getName());
-        if(!materialName.empty())
+        for(auto& group : fieldsGrp->groups())
         {
-          materialNames.emplace_back(materialName);
+          std::string materialName = fieldNameToMaterialName(group.getName());
+          if(!materialName.empty())
+          {
+            materialNames.emplace_back(materialName);
+          }
         }
       }
     }