Adds Gregory surfaces to primal

RELEASE-NOTES.md

@@ -57,6 +57,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Klee: Adds support for lua-based input decks for shaping
 - Slam: Adds convenience aliases in `axom/slam/Aliases.hpp` for the most common set and relation configurations, 
   including `ArraySet`, `ArrayViewSet`, `VariableRelation`, `ConstantRelation` and their `View` forms.
+- Primal: Adds Gregory surface classes, with a `primal::GregoryPatch` class for cubic quadrilateral surfaces and a `primal::GregoryTriangle` class for hybrid cubic-quartic triangle surfaces.
 
 ### Removed
 - Bump: Removed `axom::bump::views::MultiBufferMaterialView`, which was a view type for an obsolete flavor of Blueprint matset.

src/axom/primal/CMakeLists.txt

@@ -28,6 +28,8 @@ set( primal_headers
     geometry/CoordinateTransformer.hpp
     geometry/Cone.hpp
     geometry/CurvedPolygon.hpp
+    geometry/GregoryPatch.hpp
+    geometry/GregoryTriangle.hpp
     geometry/Hexahedron.hpp
     geometry/KnotVector.hpp
     geometry/Line.hpp

src/axom/primal/geometry/BezierCurve.hpp

@@ -395,6 +395,91 @@ class BezierCurve
     }
   }
 
+  /*!
+   * \brief Degree-elevates this Bezier curve by \a degrees
+   *
+   * Degree elevation increases the polynomial order while preserving the curve geometry.
+   * For a polynomial curve of order n, one elevation step produces an order n+1 curve
+   * with control points:
+   *   P'_0 = P_0
+   *   P'_{n+1} = P_n
+   *   P'_i = (i/(n+1)) P_{i-1} + (1 - i/(n+1)) P_i, for i=1..n
+   *
+   * For a rational curve, the degree elevation is performed in projective space by
+   * elevating both the projective control points (w*P) and the weights (w), then
+   * converting back to Euclidean control points.
+   *
+   * \param [in] degrees Number of elevation steps to apply (must be nonnegative)
+   *
+   * \note This is a no-op for empty curves (order < 0) or when \a degrees == 0
+   */
+  void degreeElevate(int degrees = 1)
+  {
+    SLIC_ASSERT(degrees >= 0);
+    if(degrees == 0)
+    {
+      return;
+    }
+
+    int ord = getOrder();
+    if(ord < 0)
+    {
+      return;
+    }
+
+    for(int step = 0; step < degrees; ++step)
+    {
+      const int n = ord;
+      const int np1 = n + 1;
+
+      axom::Array<PointType> newPts(np1 + 1);
+      newPts[0] = m_controlPoints[0];
+      newPts[np1] = m_controlPoints[n];
+
+      if(!isRational())
+      {
+        for(int i = 1; i <= n; ++i)
+        {
+          const T alpha = static_cast<T>(i) / static_cast<T>(np1);
+          newPts[i] = PointType(alpha * m_controlPoints[i - 1].array() +
+                                (T(1) - alpha) * m_controlPoints[i].array());
+        }
+
+        m_controlPoints = newPts;
+      }
+      else
+      {
+        axom::Array<T> newWts(np1 + 1);
+        newWts[0] = m_weights[0];
+        newWts[np1] = m_weights[n];
+
+        // Projective control points H_i = w_i * P_i
+        axom::Array<PointType> H(n + 1);
+        for(int i = 0; i <= n; ++i)
+        {
+          H[i] = PointType(m_weights[i] * m_controlPoints[i].array());
+        }
+
+        axom::Array<PointType> newH(np1 + 1);
+        newH[0] = H[0];
+        newH[np1] = H[n];
+
+        for(int i = 1; i <= n; ++i)
+        {
+          const T alpha = static_cast<T>(i) / static_cast<T>(np1);
+          newWts[i] = alpha * m_weights[i - 1] + (T(1) - alpha) * m_weights[i];
+          newH[i] = PointType(alpha * H[i - 1].array() + (T(1) - alpha) * H[i].array());
+          newPts[i] = PointType(newH[i].array() / newWts[i]);
+        }
+
+        m_controlPoints = newPts;
+        m_weights = newWts;
+      }
+
+      ord = np1;
+    }
+  }
+
   ///@}
 
   ///@{

src/axom/primal/geometry/BezierTriangle.hpp

@@ -55,23 +55,26 @@ std::ostream& operator<<(std::ostream& os, const BezierTriangle<T, NDIMS>& bTri)
  * A default-constructed triangle will have order -1, and is "invalid".
  * Arrays of nodes and weights will be empty, and most methods are invalid 
  * 
- * \note This triangle uses permuted barycentric coordinates (u0, v0) for evaluation such that, when
- * `getOrder()==1`, the parameter values correspond to the triangle vertices:
+ * \note The two parameters follow (u0, v0) match the indexing of control points, 
+ * such that when `getOrder()==N`, the parameter values correspond to the triangle vertices:
  * - `evaluate(0,0) == (*this)(0,0)`
- * - `evaluate(0,1) == (*this)(0,1)`
- * - `evaluate(1,0) == (*this)(1,0)`
+ * - `evaluate(1,0) == (*this)(N,0)`
+ * - `evaluate(0,1) == (*this)(0,N)`
  * 
- * These are mapped to standard Barycentric coordinates {u,v,w} through (u0, v0) = {1 - u0 - v0, v0, u0}:
+ * These parameters are mapped to standard Barycentric coordinates {u,v,w} through
+ * (u0, v0) = {1 - u0 - v0, v0, u0}:
  *
  *   Parametric (u0, v0):       Barycentric {u,v,w}:
- *        (1, 0)                      {0,0,1}
- *          /\                          /\
- *         /  \                        /  \
- *   ^    /    \         <--->        /    \
- *   |   /      \                    /      \
- *   |  /        \                  /        \
- *  v0 /__________\                /__________\
- * (0, 0) u0 ---> (0, 1)       {1,0,0}        {0,1,0}
+ *
+ *    v0
+ *     ^       (0,1)                     {0,1,0}
+ *     |       | \                       | \
+ *     |       |  \                      |  \
+ *     |       |   \       <--->         |   \
+ *     |       |    \                    |    \
+ *     |       |_____\                   |_____\
+ *     |     (0,0)  (1,0)            {1,0,0}  {0,0,1}
+ *     +-----------> u0
  * 
  */
 template <typename T, int NDIMS>
@@ -444,15 +447,15 @@ class BezierTriangle
    *
    * See overload returning a `BezierTriangle` for the vertex mapping convention.
    *
-   * \param [in] Qa Barycentric coordinates of the first subtriangle vertex `(u,v,w)`
-   * \param [in] Qb Barycentric coordinates of the second subtriangle vertex `(u,v,w)`
-   * \param [in] Qc Barycentric coordinates of the third subtriangle vertex `(u,v,w)`
+   * \param [in] Qa Barycentric coordinates of the first subtriangle vertex `{u,v,w}`
+   * \param [in] Qb Barycentric coordinates of the second subtriangle vertex `{u,v,w}`
+   * \param [in] Qc Barycentric coordinates of the third subtriangle vertex `{u,v,w}`
    * \param [out] out Output restricted Bezier triangle
    *
    * \pre getOrder() >= 0
    *
    * \note The barycentric inputs \a Qa, \a Qb, \a Qc are standard Barycentric coordiantes
-   *       related to parameter convention through (u0 = Qc, v0 = Qb)
+   *       related to parameter convention through (u0 = Qc, v0 = Qb), and Qa = 1 - u0 - v0
    */
   void restrictToSubtriangle(const Barycentric& Qa,
                              const Barycentric& Qb,