Merge pull request #1470 from LLNL/feature/spainhour/nurbs_patches

Adds untrimmed NURBS patches to primal

Author: jcs15c

RELEASE-NOTES.md

@@ -24,6 +24,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Sina C++ library is now a component of Axom
 - Adds optional dependency on [Open CASCADE](https://dev.opencascade.org). The initial intention is 
 to use Open CASCADE's file I/O capabilities in support of Quest applications.
+- Adds `primal::NURBSCurve` and `primal::NURBSPatch` classes, supported by `primal::KnotVector`.
 
 ###  Changed
 - Importing Conduit array data into `sidre::View` now allocates destination

src/axom/core/numerics/matvecops.hpp

@@ -533,10 +533,11 @@ inline bool linspace(const T& x0, const T& x1, T* v, int N)
 
   const T h = (x1 - x0) / static_cast<T>(N - 1);
 
-  for(int i = 0; i < N; ++i)
+  for(int i = 0; i < N - 1; ++i)
   {
     v[i] = x0 + i * h;
   }
+  v[N - 1] = x1;
 
   return true;
 }

src/axom/primal/CMakeLists.txt

@@ -30,6 +30,7 @@ set( primal_headers
     geometry/OrientationResult.hpp
     geometry/NumericArray.hpp
     geometry/NURBSCurve.hpp
+    geometry/NURBSPatch.hpp
     geometry/Plane.hpp
     geometry/Point.hpp
     geometry/Polygon.hpp

src/axom/primal/geometry/KnotVector.hpp

@@ -227,19 +227,28 @@ class KnotVector
    * For knot vector {u_0, ..., u_n}, returns i such that u_i <= t < u_i+1
    *  if t == u_n, returns i such that u_i < t <= u_i+1 (i.e. i = n - degree - 1)
    * 
-   * \pre Assumes that the input t is within the knot vector
-   * 
    * Implementation adapted from Algorithm A2.1 on page 68 of "The NURBS Book"
    * 
+   * \pre Assumes that the input t is in the span [u_0, u_n] (up to some tolerance)
+   * 
+   * \note If t is outside the knot span up to this tolerance, it is clamped to the span
+   * 
    * \return The index of the knot span containing t
    */
   axom::IndexType findSpan(T t) const
   {
-    SLIC_ASSERT(t >= m_knots[0] && t <= m_knots[m_knots.size() - 1]);
+    SLIC_ASSERT(isValidParameter(t));
 
     const axom::IndexType nkts = m_knots.size();
 
-    if(t == m_knots[nkts - 1])
+    // Handle cases where t is outside the knot span within a tolerance
+    //  by implicitly clamping it to the nearest span
+    if(t <= m_knots[0])
+    {
+      return m_deg;
+    }
+
+    if(t >= m_knots[nkts - 1])
     {
       return nkts - m_deg - 2;
     }
@@ -266,19 +275,29 @@ class KnotVector
    * For knot vector {u_0, ..., u_n}, returns i such that u_i <= t < u_i+1
    *  if t == u_n, returns i such that u_i < t <= u_i+1 (i.e. i = n - degree - 1)
    * 
-   * \pre Assumes that the input t is within the knot vector
+   * \pre Assumes that the input t is within the knot vector (up to some tolerance)
+   * 
+   * \note If t is outside the knot span up to this tolerance, the returned multiplicity
+   *  will be equal to the degree + 1 (required for clamped curves)
    * 
    * \return The index of the knot span containing t
    */
   axom::IndexType findSpan(T t, int& multiplicity) const
   {
-    SLIC_ASSERT(t >= m_knots[0] && t <= m_knots[m_knots.size() - 1]);
+    SLIC_ASSERT(isValidParameter(t));
 
     const auto nkts = m_knots.size();
     const auto span = findSpan(t);
 
+    // Early exit for known multiplicities
+    if(t <= m_knots[0] || t >= m_knots[nkts - 1])
+    {
+      multiplicity = m_deg + 1;
+      return span;
+    }
+
     multiplicity = 0;
-    for(auto i = (t == m_knots[nkts - 1]) ? nkts - 1 : span; i >= 0; --i)
+    for(auto i = span; i >= 0; --i)
     {
       if(m_knots[i] == t)
       {
@@ -356,20 +375,22 @@ class KnotVector
    * \param [in] t The value of the knot to insert
    * \param [in] target_mutliplicity The number of times the knot will be present
    *
-   * \pre Assumes that the input t is within the knot vector
+   * \pre Assumes that the input t is within the knot vector (up to some tolerance)
    * 
    * \note If the knot is already present, it will be inserted
    *  up to the given multiplicity, or the maximum permitted by the degree
    */
   void insertKnot(T t, int target_multiplicity)
   {
-    SLIC_ASSERT(t >= m_knots[0] && t <= m_knots[m_knots.size() - 1]);
+    SLIC_ASSERT(isValidParameter(t));
 
     int multiplicity;
     auto span = findSpan(t, multiplicity);
 
-    int r =
-      axom::utilities::clampVal(target_multiplicity - multiplicity, 0, m_deg);
+    // Compute how many knots should be inserted
+    int r = axom::utilities::clampVal(target_multiplicity - multiplicity,
+                                      0,
+                                      m_deg - multiplicity);
 
     insertKnotBySpan(span, t, r);
   }
@@ -424,8 +445,8 @@ class KnotVector
       k2.normalize();
     }
 
-    // SLIC_ASSERT(k1.isValid());
-    // SLIC_ASSERT(k2.isValid());
+    SLIC_ASSERT(k1.isValid());
+    SLIC_ASSERT(k2.isValid());
   }
 
   /*!
@@ -436,11 +457,11 @@ class KnotVector
    * \param [out] k2 The second knot vector
    * \param [in] normalize Whether to normalize the output knot vectors
    * 
-   * \pre Assumes that the input t is within the knot vector
+   * \pre Assumes that the input t is *interior* to the knot vector
    */
   void split(T t, KnotVector& k1, KnotVector& k2, bool normalize = false) const
   {
-    SLIC_ASSERT(t > m_knots[0] && t < m_knots[m_knots.size() - 1]);
+    SLIC_ASSERT(isValidInteriorParameter(t));
 
     int multiplicity;
     axom::IndexType span = findSpan(t, multiplicity);
@@ -457,7 +478,7 @@ class KnotVector
    * \param [in] span The span in which to evaluate the basis functions
    * \param [in] t The parameter value
    * 
-   * \pre Assumes that the input t is within the knot vector and that the span is valid
+   * \pre Assumes that the input t is within the correct span
    * Implementation adapted from Algorithm A2.2 on page 70 of "The NURBS Book".
    * 
    * \return An array of the `m_deg + 1` non-zero basis functions evaluated at t
@@ -496,13 +517,13 @@ class KnotVector
    * 
    * \param [in] t The parameter value
    * 
-   * \pre Assumes that the input t is within the knot vector
+   * \pre Assumes that the input t is within the knot vector (up to a tolerance)
    * 
    * \return An array of the `m_deg + 1` non-zero basis functions evaluated at t
    */
   axom::Array<T> calculateBasisFunctions(T t) const
   {
-    SLIC_ASSERT(t >= m_knots[0] && t <= m_knots[m_knots.size() - 1]);
+    SLIC_ASSERT(isValidParameter(t));
     return calculateBasisFunctionsBySpan(findSpan(t), t);
   }
 
@@ -512,19 +533,23 @@ class KnotVector
    * \param [in] span The span in which to evaluate the basis functions
    * \param [in] t The parameter value
    * \param [in] n The number of derivatives to compute
-   * \param [out] ders An array of the `n + 1` derivatives evaluated at t
    * 
    * Implementation adapted from Algorithm A2.2 on page 70 of "The NURBS Book".
+   *
+   * \pre Assumes that the input t is within the provided knot span
+   * 
+   * \return An array of the `n + 1` derivatives evaluated at t
    */
-  void derivativeBasisFunctionsBySpan(axom::IndexType span,
-                                      T t,
-                                      int n,
-                                      axom::Array<axom::Array<T>>& ders) const
+  axom::Array<axom::Array<T>> derivativeBasisFunctionsBySpan(axom::IndexType span,
+                                                             T t,
+                                                             int n) const
   {
     SLIC_ASSERT(isValidSpan(span, t));
 
     const int m_deg = getDegree();
 
+    axom::Array<axom::Array<T>> ders(n + 1);
+
     axom::Array<axom::Array<T>> ndu(m_deg + 1), a(2);
     axom::Array<T> left(m_deg + 1), right(m_deg + 1);
     for(int j = 0; j <= m_deg; j++)
@@ -596,6 +621,7 @@ class KnotVector
         std::swap(s1, s2);
       }
     }
+
     // Multiply through by the correct factors (Eq. [2.9])
     T r = static_cast<T>(m_deg);
     for(int k = 1; k <= n; k++)
@@ -606,21 +632,24 @@ class KnotVector
       }
       r *= static_cast<T>(m_deg - k);
     }
+
+    return ders;
   }
 
   /*!
    * \brief Evaluates the NURBS basis functions and derivatives for parameter value t
    * 
    * \param [in] t The parameter value
    * \param [in] n The number of derivatives to compute
-   * \param [out] ders An array of the `n + 1` derivatives evaluated at t
    * 
-   * \pre Assumes that the input t is within the knot vector
+   * \pre Assumes that the input t is within the knot vector (up to a tolerance)
+   * 
+   * \return An array of the `n + 1` derivatives evaluated at t
    */
-  void derivativeBasisFunctions(T t, int n, axom::Array<axom::Array<T>>& ders) const
+  axom::Array<axom::Array<T>> derivativeBasisFunctions(T t, int n) const
   {
-    SLIC_ASSERT(t >= m_knots[0] && t <= m_knots[m_knots.size() - 1]);
-    derivativeBasisFunctionsBySpan(findSpan(t), t, n, ders);
+    SLIC_ASSERT(isValidParameter(t));
+    return derivativeBasisFunctionsBySpan(findSpan(t), t, n);
   }
 
   /// \brief Reverse the knot vector
@@ -755,6 +784,18 @@ class KnotVector
     return !(lhs == rhs);
   }
 
+  /// \brief Checks if given parameter is in knot span (to a tolerance)
+  bool isValidParameter(T t, T EPS = 1e-5) const
+  {
+    return t >= m_knots[0] - EPS && t <= m_knots[m_knots.size() - 1] + EPS;
+  }
+
+  /// \brief Checks if given parameter is *interior* to knot span (to a tolerance)
+  bool isValidInteriorParameter(T t) const
+  {
+    return t > m_knots[0] && t < m_knots[m_knots.size() - 1];
+  }
+
   /*!
    * \brief Simple formatted print of a knot vector instance
    *