Merge pull request #1672 from LLNL/feature/spainhour/sampling_shaper_memoization

Enhance CurvedPolygon so WN memoization works with sampling shaper

Author: jcs15c

src/axom/core/numerics/quadrature.hpp

@@ -23,7 +23,7 @@ namespace numerics
 /*!
  * \class QuadratureRule
  *
- * \brief Stores a fixed array of 1D quadrature points and weights
+ * \brief Stores fixed views to arrays of 1D quadrature points and weights
  */
 class QuadratureRule
 {

src/axom/primal/geometry/BezierCurve.hpp

@@ -221,6 +221,9 @@ class BezierCurve
   /// Returns the order of the Bezier Curve
   int getOrder() const { return static_cast<int>(m_controlPoints.size()) - 1; }
 
+  /// Returns the number of control points of the Bezier Curve
+  int getNumControlPoints() const { return static_cast<int>(m_controlPoints.size()); }
+
   /// Clears the list of control points, make nonrational
   void clear()
   {
@@ -256,6 +259,9 @@ class BezierCurve
   /// Retrieves the control point at index \a idx
   const PointType& operator[](int idx) const { return m_controlPoints[idx]; }
 
+  const PointType& getInitPoint() const { return m_controlPoints[0]; }
+  const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; }
+
   /// Returns a copy of the Bezier curve's control points
   CoordsVec getControlPoints() const { return m_controlPoints; }
 

src/axom/primal/geometry/CurvedPolygon.hpp

@@ -18,22 +18,64 @@
 #include "axom/primal/geometry/Point.hpp"
 #include "axom/primal/geometry/Vector.hpp"
 #include "axom/primal/geometry/BezierCurve.hpp"
+#include "axom/primal/geometry/NURBSCurve.hpp"
 #include "axom/primal/geometry/BoundingBox.hpp"
 
+// For NURBSCurveGWNCache objects
+#include "axom/primal/operators/detail/winding_number_3d_memoization.hpp"
+
 #include <vector>
 #include <ostream>
 
 namespace axom
 {
 namespace primal
 {
-// Forward declare the templated classes and operator functions
+namespace internal
+{
+///@{
+/// \name Boilerplate to deduce the numeric type from the curve object
+template <typename U>
+struct get_numeric_type;
+
 template <typename T, int NDIMS>
+struct get_numeric_type<BezierCurve<T, NDIMS>>
+{
+  using type = T;
+};
+
+template <typename T, int NDIMS>
+struct get_numeric_type<NURBSCurve<T, NDIMS>>
+{
+  using type = T;
+};
+
+template <typename T>
+struct get_numeric_type<detail::NURBSCurveGWNCache<T>>
+{
+  using type = T;
+};
+///@}
+
+///@{
+/// \name Boilerplate for a compile-time flag for the cached object
+template <typename U>
+struct has_cached_data : std::false_type
+{ };
+
+template <typename T>
+struct has_cached_data<detail::NURBSCurveGWNCache<T>> : std::true_type
+{ };
+///@}
+}  // namespace internal
+
+// Forward declare the templated classes and operator functions
+template <typename CurveType>
 class CurvedPolygon;
 
 /*! \brief Overloaded output operator for polygons */
-template <typename T, int NDIMS>
-std::ostream& operator<<(std::ostream& os, const CurvedPolygon<T, NDIMS>& poly);
+template <typename CurveType>
+std::ostream& operator<<(std::ostream& os, const CurvedPolygon<CurveType>& poly);
 
 /*!
  * \class CurvedPolygon
@@ -44,15 +86,15 @@ std::ostream& operator<<(std::ostream& os, const CurvedPolygon<T, NDIMS>& poly);
  * \note The component curves should be ordered in a counter clockwise
  *       orientation with respect to the polygon's normal vector
  */
-template <typename T, int NDIMS>
+template <typename CurveType>
 class CurvedPolygon
 {
 public:
-  using PointType = Point<T, NDIMS>;
-  using VectorType = Vector<T, NDIMS>;
-  using NumArrayType = axom::NumericArray<T, NDIMS>;
-  using BezierCurveType = BezierCurve<T, NDIMS>;
-  using BoundingBoxType = typename BezierCurveType::BoundingBoxType;
+  using T = typename internal::get_numeric_type<CurveType>::type;
+
+  using PointType = typename CurveType::PointType;
+  using VectorType = typename CurveType::VectorType;
+  using BoundingBoxType = typename CurveType::BoundingBoxType;
 
 public:
   /// Default constructor for an empty polygon
@@ -73,20 +115,30 @@ class CurvedPolygon
     m_edges.resize(nEdges);
   }
 
-  /// Constructor from an array of \a nEdges curves
-  CurvedPolygon(BezierCurveType* curves, int nEdges)
+  /*!
+   * \brief Copy constructor for another curve type
+   */
+  template <typename OtherCurveType>
+  CurvedPolygon(const CurvedPolygon<OtherCurveType>& other_poly)
   {
-    SLIC_ASSERT(curves != nullptr);
-    SLIC_ASSERT(nEdges >= 1);
-
-    m_edges.reserve(nEdges);
-
-    for(int e = 0; e < nEdges; ++e)
+    m_edges.resize(other_poly.numEdges());
+    for(int e = 0; e < other_poly.numEdges(); ++e)
     {
-      this->addEdge(curves[e]);
+      m_edges[e] = CurveType(other_poly[e]);
     }
   }
 
+  ///! \brief Constructor for CurvedPolygon from an ArrayView of curves
+  CurvedPolygon(axom::ArrayView<const CurveType> curves) : m_edges(curves) { }
+
+  ///! \brief Constructor from a c-style array of \a nEdges curves
+  CurvedPolygon(const CurveType* curves, int nEdges)
+    : CurvedPolygon(axom::ArrayView<const CurveType>(curves, nEdges))
+  { }
+
+  ///! \brief Constructor from a Axom array of curves
+  CurvedPolygon(const axom::Array<CurveType>& curves) : CurvedPolygon(curves.view()) { }
+
   /// Clears the list of edges
   void clear() { m_edges.clear(); }
 
@@ -105,36 +157,43 @@ class CurvedPolygon
     m_edges.resize(ngon);
   }
 
-  /// Appends a BezierCurve to the list of edges
-  void addEdge(const BezierCurveType& c1) { m_edges.push_back(c1); }
+  /// Appends a curve to the list of edges
+  void addEdge(const CurveType& c1) { m_edges.push_back(c1); }
+
+  /// Consumes then appends a curve to the list of edges
+  void addEdge(CurveType&& c1) { m_edges.push_back(std::move(c1)); }
 
   /// Splits an edge "in place"
   void splitEdge(int idx, T t)
   {
-    SLIC_ASSERT(idx < static_cast<int>(m_edges.size()));
+    SLIC_ASSERT(idx >= 0 && idx < static_cast<int>(m_edges.size()));
+    AXOM_STATIC_ASSERT_MSG(!internal::has_cached_data<CurveType>::value,
+                           "splitEdge cannot be called on objects with associated cache data");
 
+    m_edges.reserve(numEdges() + 1);
     m_edges.insert(m_edges.begin() + idx + 1, 1, m_edges[idx]);
     auto& csplit = m_edges[idx];
     csplit.split(t, m_edges[idx], m_edges[idx + 1]);
   }
 
-  axom::Array<BezierCurve<T, NDIMS>> getEdges() const { return m_edges; }
-
+  axom::Array<CurveType> getEdges() const { return m_edges; }
   /// @}
 
-  /*! Retrieves the Bezier Curve at index idx */
-  BezierCurveType& operator[](int idx) { return m_edges[idx]; }
-  /*! Retrieves the vertex at index idx */
-  const BezierCurveType& operator[](int idx) const { return m_edges[idx]; }
+  /*! Retrieves the curve at index idx */
+  CurveType& operator[](int idx) { return m_edges[idx]; }
+  /*! Retrieves the curve at index idx */
+  const CurveType& operator[](int idx) const { return m_edges[idx]; }
 
   /// Tests equality of two CurvedPolygons
-  friend inline bool operator==(const CurvedPolygon<T, NDIMS>& lhs, const CurvedPolygon<T, NDIMS>& rhs)
+  friend inline bool operator==(const CurvedPolygon<CurveType>& lhs,
+                                const CurvedPolygon<CurveType>& rhs)
   {
     return lhs.m_edges == rhs.m_edges;
   }
 
   /// Tests inequality of two CurvedPolygons
-  friend inline bool operator!=(const CurvedPolygon<T, NDIMS>& lhs, const CurvedPolygon<T, NDIMS>& rhs)
+  friend inline bool operator!=(const CurvedPolygon<CurveType>& lhs,
+                                const CurvedPolygon<CurveType>& rhs)
   {
     return !(lhs == rhs);
   }
@@ -149,7 +208,7 @@ class CurvedPolygon
   {
     const int sz = numEdges();
 
-    os << "{" << sz << "-sided Bezier polygon:";
+    os << "{" << sz << "-sided curved polygon:";
     for(int i = 0; i < sz - 1; ++i)
     {
       os << m_edges[i] << ",";
@@ -177,17 +236,16 @@ class CurvedPolygon
     const int nEdges = numEdges();
 
     // initial basic check: no edges, or one edge or linear or quadratic order cannot be closed
-    if(nEdges < 1 || (nEdges == 1 && m_edges[0].getOrder() <= 2))
+    if(nEdges < 1 || (nEdges == 1 && m_edges[0].getNumControlPoints() <= 2))
     {
       return false;
     }
 
     // foreach edge: check last vertex of previous edge against first vertex of current edge
     for(int cur = 0, prev = nEdges - 1; cur < nEdges; prev = cur++)
     {
-      const auto ord = m_edges[prev].getOrder();
-      const auto& lastPrev = m_edges[prev][ord];
-      const auto& firstCur = m_edges[cur][0];
+      const auto& lastPrev = m_edges[prev].getEndPoint();
+      const auto& firstCur = m_edges[cur].getInitPoint();
       if(!isNearlyEqual(squared_distance(lastPrev, firstCur), 0., sq_tol))
       {
         return false;
@@ -200,6 +258,10 @@ class CurvedPolygon
   /// \brief Reverses orientation of a CurvedPolygon
   void reverseOrientation()
   {
+    AXOM_STATIC_ASSERT_MSG(
+      !internal::has_cached_data<CurveType>::value,
+      "reverseOrientation cannot be called on objects with associated cache data");
+
     const int ngon = numEdges();
     const int mid = ngon >> 1;
     const bool isOdd = (ngon & 1) != 0;
@@ -233,14 +295,14 @@ class CurvedPolygon
   }
 
 private:
-  axom::Array<BezierCurve<T, NDIMS>> m_edges;
+  axom::Array<CurveType> m_edges;
 };
 
 //------------------------------------------------------------------------------
 /// Free functions implementing Polygon's operators
 //------------------------------------------------------------------------------
-template <typename T, int NDIMS>
-std::ostream& operator<<(std::ostream& os, const CurvedPolygon<T, NDIMS>& poly)
+template <typename CurveType>
+std::ostream& operator<<(std::ostream& os, const CurvedPolygon<CurveType>& poly)
 {
   poly.print(os);
   return os;

src/axom/primal/geometry/NURBSCurve.hpp

@@ -621,6 +621,9 @@ class NURBSCurve
    */
   const PointType& operator[](int idx) const { return m_controlPoints[idx]; }
 
+  const PointType& getInitPoint() const { return m_controlPoints[0]; }
+  const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; }
+
   /// \brief Returns a copy of the NURBS curve's control points
   CoordsVec getControlPoints() const { return m_controlPoints; }
 

src/axom/primal/operators/compute_moments.hpp

@@ -95,7 +95,7 @@ primal::Point<T, 2> sector_centroid(const primal::BezierCurve<T, 2>& curve)
 
 /// \brief Returns the area enclosed by the CurvedPolygon
 template <typename T>
-T area(const primal::CurvedPolygon<T, 2>& poly, double tol = 1e-8)
+T area(const primal::CurvedPolygon<BezierCurve<T, 2>>& poly, double tol = 1e-8)
 {
   const int ngon = poly.numEdges();
   T A = 0.0;
@@ -119,7 +119,7 @@ T area(const primal::CurvedPolygon<T, 2>& poly, double tol = 1e-8)
 
 /// \brief Returns the centroid of the CurvedPolygon
 template <typename T>
-primal::Point<T, 2> centroid(const primal::CurvedPolygon<T, 2>& poly, double tol = 1e-8)
+primal::Point<T, 2> centroid(const primal::CurvedPolygon<BezierCurve<T, 2>>& poly, double tol = 1e-8)
 {
   using PointType = primal::Point<T, 2>;