From c62a81500b19d74976b14fd8eb6e5b5aa82be954 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Sat, 16 May 2026 18:13:24 -0700 Subject: [PATCH 01/16] Initial effort --- src/axom/primal/CMakeLists.txt | 1 + src/axom/primal/geometry/GregoryPatch.hpp | 246 ++++++++++++++++++++++ 2 files changed, 247 insertions(+) create mode 100644 src/axom/primal/geometry/GregoryPatch.hpp diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 905252da93..60a44cbb82 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -28,6 +28,7 @@ set( primal_headers geometry/CoordinateTransformer.hpp geometry/Cone.hpp geometry/CurvedPolygon.hpp + geometry/GregoryPatch.hpp geometry/Hexahedron.hpp geometry/KnotVector.hpp geometry/Line.hpp diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp new file mode 100644 index 0000000000..454304bd9f --- /dev/null +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -0,0 +1,246 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file GregoryPatch.hpp + * + * \brief A Gregory patch primitive + */ + +#ifndef AXOM_PRIMAL_GREGORY_PATCH_HPP_ +#define AXOM_PRIMAL_GREGORY_PATCH_HPP_ + +#include "axom/core.hpp" +#include "axom/slic.hpp" + +#include "axom/core/NumericArray.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Vector.hpp" +#include "axom/primal/geometry/Segment.hpp" +#include "axom/primal/geometry/BezierPatch.hpp" +#include "axom/primal/geometry/BoundingBox.hpp" +#include "axom/primal/geometry/OrientedBoundingBox.hpp" + +#include +#include + +#include "axom/fmt.hpp" + +namespace axom +{ +namespace primal +{ +// Forward declare the templated classes and operator functions +template +class GregoryPatch; + +/*! \brief Overloaded output operator for Gregory Patches*/ +template +std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch); + +/*! + * \class GregoryPatch + * + * \brief Represents a 3D Gregory patch defined by ... + * \tparam T the coordinate type, e.g., double, float, etc. + * \tparam NDIMS the number of dimensions + */ +template +class GregoryPatch +{ +public: + // The number of control points for a bicubic Gregory patch is fixed: + // - 12 exterior control points + // - 2 interior control points for each of 4 boundary curves + static constexpr int NPTS = 20; + + using PointType = Point; + using VectorType = Vector; + + using CoordsVec = axom::StackArray; + + using BoundingBoxType = BoundingBox; + using OrientedBoundingBoxType = OrientedBoundingBox; + + AXOM_STATIC_ASSERT_MSG((NDIMS == 1) || (NDIMS == 2) || (NDIMS == 3), + "A Gregory Patch object may be defined in 1-, 2-, or 3-D"); + + AXOM_STATIC_ASSERT_MSG(std::is_arithmetic::value, + "A Gregory Patch must be defined using an arithmetic type"); + +public: + GregoryPatch() = default; + + explicit GregoryPatch(ArrayView controlPoints) + { + SLIC_ASSERT(controlPoints.size() == NPTS); + for(int i = 0; i < NPTS; ++i) + { + m_controlPoints[i] = controlPoints[i]; + } + } + + GregoryPatch(ArrayView nodePositions, ArrayView nodeVectors) + { + // Store the position and orthogonal unit vector at each corner + SLIC_ASSERT(nodePositions.size() == 4); + SLIC_ASSERT(nodeVectors.size() == 4); + + axom::Array v(4); + for(int i = 0; i < 4; ++i) + { + getCorner(i) = nodePositions[i]; + v[i] = nodeVectors[i].unitVector(); + } + + axom::Array c0(4), c2(4); + axom::Array a0(4), a3(4); + for(int i = 0; i < 4; ++i) + { + const int ip1 = (i + 1) % 4; + const int im1 = (i + 3) % 4; + + const VectorType dx(getCorner(i), getCorner(ip1)); + + c0[i] = (dx - dx.dot(v[i]) * v[i]) / 3; + a0[i] = VectorType::cross_product(v[i], dx).unitVector(); + c2[i] = (dx - dx.dot(v[ip1]) * v[ip1]) / 3; + a3[i] = VectorType::cross_product(v[ip1], dx).unitVector(); + + // Other stuff for adjusted normals + + // Use Chiyokura algorithm to define the interior control points + const PointType x1(getCorner(i).array() + c0[i].array()); + const PointType x2(getCorner(ip1).array() - c2[i].array()); + + getBoundaryPoint(i, 1) = x1; + getBoundaryPoint(i, 2) = x2; + + const VectorType c1(x1, x2); + const VectorType b0 = -c2[im1]; + const VectorType b3 = c0[ip1]; + + const double k0 = a0[i].dot(b0); + const double k1 = a3[i].dot(b3); + const double h0 = c0[i].dot(b0) / c0[i].dot(c0[i]); + const double h1 = c2[i].dot(b3) / c2[i].dot(c2[i]); + + const VectorType b1 = ((k0 + k1) * a0[i] + k0 * a3[i] + 2.0 * h0 * c1 + h1 * c0[i]) / 3.0; + const VectorType b2 = ((k0 + k1) * a3[i] + k1 * a0[i] + 2.0 * h1 * c1 + h0 * c2[i]) / 3.0; + + getTangent(i, 0) = PointType(x1.array() + b1.array()); + getTangent(i, 1) = PointType(x2.array() + b2.array()); + } + } + + PointType& getCorner(int i) { return m_controlPoints[i]; } + const PointType& getCorner(int i) const { return m_controlPoints[i]; } + + PointType& getTangent(int e, int t) { return m_controlPoints[12 + 2 * e + t]; } + const PointType& getTangent(int e, int t) const { return m_controlPoints[12 + 2 * e + t]; } + + PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[m_index_map[e][k]]; } + const PointType& getBoundaryPoint(int e, int k) const + { + return m_controlPoints[m_index_map[e][k]]; + } + + // Evaluate the patch by constructing the equivalent Bezier patch with interior control nodes + // defined in terms of the tangent vectors and the evaluation parameters + PointType evaluate(T u, T v) const + { + auto gregoryBlend = [&](const PointType& a, const PointType& b, T wa, T wb) -> PointType { + const T denom = wa + wb; + if(axom::utilities::isNearlyEqual(denom, T(0))) + { + return a; + } + return PointType((wa * a.array() + wb * b.array()) / denom); + }; + + const T um = T(1) - u; + const T vm = T(1) - v; + + const PointType Q11 = gregoryBlend(getTangent(0, 0), getTangent(3, 1), u, v); + const PointType Q21 = gregoryBlend(getTangent(0, 1), getTangent(1, 0), um, v); + const PointType Q12 = gregoryBlend(getTangent(2, 1), getTangent(3, 0), u, vm); + const PointType Q22 = gregoryBlend(getTangent(2, 0), getTangent(1, 1), um, vm); + + auto bpatch = get_bezier_boundary(); + + bpatch(1, 1) = Q11; + bpatch(2, 1) = Q21; + bpatch(1, 2) = Q12; + bpatch(2, 2) = Q22; + + return bpatch.evaluate(u, v); + } + + void print(std::ostream& os) const + { + os << "GregoryPatch<" << NDIMS << "D>("; + for(int i = 0; i < NPTS; ++i) + { + os << m_controlPoints[i]; + it if(i + 1 < NPTS) { os << ", "; } + } + os << ")"; + } + +private: + // Copies over the boundary points to a BezierPatch object, + // leaving the 4 interior control points uninitialized + primal::BezierPatch get_bezier_boundary() const + { + BezierPatch bpatch(3, 3); + + bpatch(0, 0) = getCorner(0); + bpatch(1, 0) = getBoundaryPoint(0, 1); + bpatch(2, 0) = getBoundaryPoint(0, 2); + bpatch(3, 0) = getCorner(1); + + bpatch(3, 1) = getBoundaryPoint(1, 1); + bpatch(3, 2) = getBoundaryPoint(1, 2); + bpatch(3, 3) = getCorner(2); + + bpatch(2, 3) = getBoundaryPoint(2, 1); + bpatch(1, 3) = getBoundaryPoint(2, 2); + bpatch(0, 3) = getCorner(3); + + bpatch(0, 2) = getBoundaryPoint(3, 1); + bpatch(0, 1) = getBoundaryPoint(3, 2); + + return bpatch; + } + + CoordsVec m_controlPoints; + + // Map of boundary curve control points into internal storage + static constexpr int m_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, + {/*V1*/ 1, /*E11*/ 6, /*E12*/ 7, /*V2*/ 2}, + {/*V2*/ 2, /*E21*/ 8, /*E22*/ 9, /*V3*/ 3}, + {/*V3*/ 3, /*E31*/ 10, /*E32*/ 11, /*V0*/ 0}}; +}; + +//------------------------------------------------------------------------------ +/// Free functions related to GregoryPatch +//------------------------------------------------------------------------------ +template +std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch) +{ + nPatch.print(os); + return os; +} + +} // namespace primal +} // namespace axom + +/// Overload to format a primal::GregoryPatch using fmt +template +struct axom::fmt::formatter> : ostream_formatter +{ }; + +#endif // AXOM_PRIMAL_GREGORY_PATCH_HPP_ From b3d13857618765efffce453553bf9ddcfc8c5bff Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Sat, 16 May 2026 23:52:45 -0700 Subject: [PATCH 02/16] A little more detail --- src/axom/primal/geometry/GregoryPatch.hpp | 341 +++++++++++++++++++--- 1 file changed, 299 insertions(+), 42 deletions(-) diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp index 454304bd9f..d32147c5ea 100644 --- a/src/axom/primal/geometry/GregoryPatch.hpp +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -34,12 +34,12 @@ namespace axom namespace primal { // Forward declare the templated classes and operator functions -template +template class GregoryPatch; /*! \brief Overloaded output operator for Gregory Patches*/ -template -std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch); +template +std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch); /*! * \class GregoryPatch @@ -48,7 +48,7 @@ std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch) * \tparam T the coordinate type, e.g., double, float, etc. * \tparam NDIMS the number of dimensions */ -template +template class GregoryPatch { public: @@ -57,16 +57,13 @@ class GregoryPatch // - 2 interior control points for each of 4 boundary curves static constexpr int NPTS = 20; - using PointType = Point; - using VectorType = Vector; + using PointType = Point; + using VectorType = Vector; using CoordsVec = axom::StackArray; - using BoundingBoxType = BoundingBox; - using OrientedBoundingBoxType = OrientedBoundingBox; - - AXOM_STATIC_ASSERT_MSG((NDIMS == 1) || (NDIMS == 2) || (NDIMS == 3), - "A Gregory Patch object may be defined in 1-, 2-, or 3-D"); + using BoundingBoxType = BoundingBox; + using OrientedBoundingBoxType = OrientedBoundingBox; AXOM_STATIC_ASSERT_MSG(std::is_arithmetic::value, "A Gregory Patch must be defined using an arithmetic type"); @@ -142,60 +139,303 @@ class GregoryPatch PointType& getTangent(int e, int t) { return m_controlPoints[12 + 2 * e + t]; } const PointType& getTangent(int e, int t) const { return m_controlPoints[12 + 2 * e + t]; } - PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[m_index_map[e][k]]; } + void getTangentsByCorner(int i, PointType& v0, PointType& v1) const + { + v0 = getTangent((i + 3) % 4, 1); + v1 = getTangent(i, 0); + } + + PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_index_map[e][k]]; } const PointType& getBoundaryPoint(int e, int k) const { - return m_controlPoints[m_index_map[e][k]]; + return m_controlPoints[s_index_map[e][k]]; } // Evaluate the patch by constructing the equivalent Bezier patch with interior control nodes // defined in terms of the tangent vectors and the evaluation parameters PointType evaluate(T u, T v) const { - auto gregoryBlend = [&](const PointType& a, const PointType& b, T wa, T wb) -> PointType { - const T denom = wa + wb; - if(axom::utilities::isNearlyEqual(denom, T(0))) - { - return a; - } - return PointType((wa * a.array() + wb * b.array()) / denom); - }; + const auto intermediate = setup_intermediate_bezier(u, v, 0); + return intermediate.bpatch.evaluate(u, v); + } - const T um = T(1) - u; - const T vm = T(1) - v; + /*! + * \brief Evaluates all first derivatives of the Gregory patch at (\a u, \a v) + * + * \param [in] u Parameter value at which to evaluate on the first axis + * \param [in] v Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory patch at (u, v) + * \param [out] Du The vector value of S_u(u, v) + * \param [out] Dv The vector value of S_v(u, v) + */ + void evaluateFirstDerivatives(T u, T v, PointType& eval, VectorType& Du, VectorType& Dv) const + { + const auto intermediate = setup_intermediate_bezier(u, v, 1); + intermediate.bpatch.evaluateFirstDerivatives(u, v, eval, Du, Dv); + + // Chain rule correction due to (u,v)-dependent interior control points + axom::StaticArray Bu, dBu, Bv, dBv; + evaluateCubicBernstein(u, Bu, dBu); + evaluateCubicBernstein(v, Bv, dBv); + + const T w11 = Bu[1] * Bv[1]; + const T w21 = Bu[2] * Bv[1]; + const T w12 = Bu[1] * Bv[2]; + const T w22 = Bu[2] * Bv[2]; + + Du += w11 * intermediate.Q_u[0][0] + w21 * intermediate.Q_u[1][0] + + w12 * intermediate.Q_u[0][1] + w22 * intermediate.Q_u[1][1]; + Dv += w11 * intermediate.Q_v[0][0] + w21 * intermediate.Q_v[1][0] + + w12 * intermediate.Q_v[0][1] + w22 * intermediate.Q_v[1][1]; + } - const PointType Q11 = gregoryBlend(getTangent(0, 0), getTangent(3, 1), u, v); - const PointType Q21 = gregoryBlend(getTangent(0, 1), getTangent(1, 0), um, v); - const PointType Q12 = gregoryBlend(getTangent(2, 1), getTangent(3, 0), u, vm); - const PointType Q22 = gregoryBlend(getTangent(2, 0), getTangent(1, 1), um, vm); + /*! + * \brief Evaluates all second derivatives of the Gregory patch at (\a u, \a v) + * + * \param [in] u Parameter value at which to evaluate on the first axis + * \param [in] v Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory patch at (u, v) + * \param [out] Du The vector value of S_u(u, v) + * \param [out] Dv The vector value of S_v(u, v) + * \param [out] DuDu The vector value of S_uu(u, v) + * \param [out] DvDv The vector value of S_vv(u, v) + * \param [out] DuDv The vector value of S_uv(u, v) == S_vu(u, v) + */ + void evaluateSecondDerivatives(T u, + T v, + PointType& eval, + VectorType& Du, + VectorType& Dv, + VectorType& DuDu, + VectorType& DvDv, + VectorType& DuDv) const + { + const auto intermediate = setup_intermediate_bezier(u, v, 2); + intermediate.bpatch.evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + + // Chain rule correction due to (u,v)-dependent interior control points + axom::StaticArray Bu, dBu, Bv, dBv; + evaluateCubicBernstein(u, Bu, dBu); + evaluateCubicBernstein(v, Bv, dBv); + + const T w11 = Bu[1] * Bv[1]; + const T w21 = Bu[2] * Bv[1]; + const T w12 = Bu[1] * Bv[2]; + const T w22 = Bu[2] * Bv[2]; + + const T wu11 = dBu[1] * Bv[1]; + const T wu21 = dBu[2] * Bv[1]; + const T wu12 = dBu[1] * Bv[2]; + const T wu22 = dBu[2] * Bv[2]; + + const T wv11 = Bu[1] * dBv[1]; + const T wv21 = Bu[2] * dBv[1]; + const T wv12 = Bu[1] * dBv[2]; + const T wv22 = Bu[2] * dBv[2]; + + // First derivative corrections + Du += w11 * intermediate.Q_u[0][0] + w21 * intermediate.Q_u[1][0] + + w12 * intermediate.Q_u[0][1] + w22 * intermediate.Q_u[1][1]; + Dv += w11 * intermediate.Q_v[0][0] + w21 * intermediate.Q_v[1][0] + + w12 * intermediate.Q_v[0][1] + w22 * intermediate.Q_v[1][1]; + + // Second derivative corrections + DuDu += T(2) * + (wu11 * intermediate.Q_u[0][0] + wu21 * intermediate.Q_u[1][0] + + wu12 * intermediate.Q_u[0][1] + wu22 * intermediate.Q_u[1][1]) + + (w11 * intermediate.Q_uu[0][0] + w21 * intermediate.Q_uu[1][0] + + w12 * intermediate.Q_uu[0][1] + w22 * intermediate.Q_uu[1][1]); + + DvDv += T(2) * + (wv11 * intermediate.Q_v[0][0] + wv21 * intermediate.Q_v[1][0] + + wv12 * intermediate.Q_v[0][1] + wv22 * intermediate.Q_v[1][1]) + + (w11 * intermediate.Q_vv[0][0] + w21 * intermediate.Q_vv[1][0] + + w12 * intermediate.Q_vv[0][1] + w22 * intermediate.Q_vv[1][1]); + + DuDv += (wu11 * intermediate.Q_v[0][0] + wu21 * intermediate.Q_v[1][0] + + wu12 * intermediate.Q_v[0][1] + wu22 * intermediate.Q_v[1][1]) + + (wv11 * intermediate.Q_u[0][0] + wv21 * intermediate.Q_u[1][0] + + wv12 * intermediate.Q_u[0][1] + wv22 * intermediate.Q_u[1][1]) + + (w11 * intermediate.Q_uv[0][0] + w21 * intermediate.Q_uv[1][0] + + w12 * intermediate.Q_uv[0][1] + w22 * intermediate.Q_uv[1][1]); + } - auto bpatch = get_bezier_boundary(); + VectorType du(T u, T v) const + { + PointType eval; + VectorType Du, Dv; + evaluateFirstDerivatives(u, v, eval, Du, Dv); + return Du; + } + + VectorType dv(T u, T v) const + { + PointType eval; + VectorType Du, Dv; + evaluateFirstDerivatives(u, v, eval, Du, Dv); + return Dv; + } + + VectorType dudu(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DuDu; + } + + VectorType dvdv(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DvDv; + } + + VectorType dudv(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DuDv; + } - bpatch(1, 1) = Q11; - bpatch(2, 1) = Q21; - bpatch(1, 2) = Q12; - bpatch(2, 2) = Q22; + /// \brief Returns an axis-aligned bounding box containing the patch + BoundingBoxType boundingBox() const + { + return BoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); + } - return bpatch.evaluate(u, v); + /// \brief Returns an oriented bounding box containing the patch + OrientedBoundingBoxType orientedBoundingBox() const + { + return OrientedBoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); } void print(std::ostream& os) const { - os << "GregoryPatch<" << NDIMS << "D>("; + os << "GregoryPatch("; for(int i = 0; i < NPTS; ++i) { os << m_controlPoints[i]; - it if(i + 1 < NPTS) { os << ", "; } + if(i + 1 < NPTS) + { + os << ", "; + } } os << ")"; } private: + struct IntermediateBezierDerivatives + { + BezierPatch bpatch; + PointType Q[2][2]; + VectorType Q_u[2][2]; + VectorType Q_v[2][2]; + VectorType Q_uu[2][2]; + VectorType Q_vv[2][2]; + VectorType Q_uv[2][2]; + }; + + IntermediateBezierDerivatives setup_intermediate_bezier(T u, T v, int derivative_order) const + { + IntermediateBezierDerivatives out; + + out.bpatch = get_bezier_boundary(); + + const T um = T(1) - u; + const T vm = T(1) - v; + + // For each interior point (i,j), select which corner's tangents to use. + // This mapping matches the explicit construction used in evaluate(): + // (0,0)->corner0, (1,0)->corner1, (0,1)->corner3, (1,1)->corner2 + static constexpr int corner_of[2][2] = {{0, 3}, {1, 2}}; + + for(int j = 0; j < 2; ++j) + { + const T wb = (j == 0) ? v : vm; + const T wb_v = (j == 0) ? T(1) : T(-1); + + for(int i = 0; i < 2; ++i) + { + const T wa = (i == 0) ? u : um; + const T wa_u = (i == 0) ? T(1) : T(-1); + + const int corner = corner_of[i][j]; + + PointType tPrev, tNext; + getTangentsByCorner(corner, tPrev, tNext); + + // Corner parity determines which tangent is blended with which weight + // (see explicit mapping in evaluate()). + const bool swap = (corner % 2 == 0); // corners 0 and 2 + const PointType& A = swap ? tNext : tPrev; + const PointType& B = swap ? tPrev : tNext; + + const T denom = wa + wb; + if(axom::utilities::isNearlyEqual(denom, T(0))) + { + out.Q[i][j] = A; + out.Q_u[i][j] = VectorType(T(0)); + out.Q_v[i][j] = VectorType(T(0)); + out.Q_uu[i][j] = VectorType(T(0)); + out.Q_vv[i][j] = VectorType(T(0)); + out.Q_uv[i][j] = VectorType(T(0)); + continue; + } + + out.Q[i][j] = PointType((wa * A.array() + wb * B.array()) / denom); + + if(derivative_order >= 1) + { + // Since wa depends only on u and wb depends only on v: + // dQ/du uses only wa_u; dQ/dv uses only wb_v. + out.Q_u[i][j] = VectorType((wa_u * (A.array() - out.Q[i][j].array())) / denom); + out.Q_v[i][j] = VectorType((wb_v * (B.array() - out.Q[i][j].array())) / denom); + } + else + { + out.Q_u[i][j] = VectorType(T(0)); + out.Q_v[i][j] = VectorType(T(0)); + } + + if(derivative_order >= 2) + { + // Second derivatives for linear weights (wa, wb): + // Q_uu = -(2*wa_u/denom) * Q_u + // Q_vv = -(2*wb_v/denom) * Q_v + // Q_uv = -(wa_u*Q_v + wb_v*Q_u)/denom + out.Q_uu[i][j] = (-T(2) * wa_u / denom) * out.Q_u[i][j]; + out.Q_vv[i][j] = (-T(2) * wb_v / denom) * out.Q_v[i][j]; + out.Q_uv[i][j] = (-(wa_u * out.Q_v[i][j] + wb_v * out.Q_u[i][j])) / denom; + } + else + { + out.Q_uu[i][j] = VectorType(T(0)); + out.Q_vv[i][j] = VectorType(T(0)); + out.Q_uv[i][j] = VectorType(T(0)); + } + } + } + + setBezierInterior(out.bpatch, out.Q); + return out; + } + + static void setBezierInterior(BezierPatch& bpatch, const PointType Q[2][2]) + { + bpatch(1, 1) = Q[0][0]; + bpatch(2, 1) = Q[1][0]; + bpatch(1, 2) = Q[0][1]; + bpatch(2, 2) = Q[1][1]; + } + // Copies over the boundary points to a BezierPatch object, // leaving the 4 interior control points uninitialized - primal::BezierPatch get_bezier_boundary() const + BezierPatch get_bezier_boundary() const { - BezierPatch bpatch(3, 3); + BezierPatch bpatch(3, 3); bpatch(0, 0) = getCorner(0); bpatch(1, 0) = getBoundaryPoint(0, 1); @@ -216,10 +456,27 @@ class GregoryPatch return bpatch; } + static void evaluateCubicBernstein(T t, axom::StaticArray& B, axom::StaticArray& dB) + { + const T tm = T(1) - t; + const T tm2 = tm * tm; + const T t2 = t * t; + + B[0] = tm2 * tm; + B[1] = T(3) * t * tm2; + B[2] = T(3) * t2 * tm; + B[3] = t2 * t; + + dB[0] = -T(3) * tm2; + dB[1] = T(3) * tm2 - T(6) * t * tm; + dB[2] = T(6) * t * tm - T(3) * t2; + dB[3] = T(3) * t2; + } + CoordsVec m_controlPoints; // Map of boundary curve control points into internal storage - static constexpr int m_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, + static constexpr int s_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, {/*V1*/ 1, /*E11*/ 6, /*E12*/ 7, /*V2*/ 2}, {/*V2*/ 2, /*E21*/ 8, /*E22*/ 9, /*V3*/ 3}, {/*V3*/ 3, /*E31*/ 10, /*E32*/ 11, /*V0*/ 0}}; @@ -228,8 +485,8 @@ class GregoryPatch //------------------------------------------------------------------------------ /// Free functions related to GregoryPatch //------------------------------------------------------------------------------ -template -std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch) +template +std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch) { nPatch.print(os); return os; @@ -239,8 +496,8 @@ std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch) } // namespace axom /// Overload to format a primal::GregoryPatch using fmt -template -struct axom::fmt::formatter> : ostream_formatter +template +struct axom::fmt::formatter> : ostream_formatter { }; #endif // AXOM_PRIMAL_GREGORY_PATCH_HPP_ From 171a3b72cf4c5bc7e82d3b0b887628e3bc505273 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 9 Jun 2026 08:23:36 -0700 Subject: [PATCH 03/16] Add degree elevation method --- src/axom/primal/geometry/BezierCurve.hpp | 85 +++++++++++++++++++ src/axom/primal/tests/primal_bezier_curve.cpp | 44 ++++++++++ 2 files changed, 129 insertions(+) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 4e98cd56ce..af173ac908 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/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 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(i) / static_cast(np1); + newPts[i] = + PointType(alpha * m_controlPoints[i - 1].array() + (T(1) - alpha) * m_controlPoints[i].array()); + } + + m_controlPoints = newPts; + } + else + { + axom::Array newWts(np1 + 1); + newWts[0] = m_weights[0]; + newWts[np1] = m_weights[n]; + + // Projective control points H_i = w_i * P_i + axom::Array H(n + 1); + for(int i = 0; i <= n; ++i) + { + H[i] = PointType(m_weights[i] * m_controlPoints[i].array()); + } + + axom::Array newH(np1 + 1); + newH[0] = H[0]; + newH[np1] = H[n]; + + for(int i = 1; i <= n; ++i) + { + const T alpha = static_cast(i) / static_cast(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; + } + } + ///@} ///@{ diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index 20be6e8c46..5775d35cbc 100644 --- a/src/axom/primal/tests/primal_bezier_curve.cpp +++ b/src/axom/primal/tests/primal_bezier_curve.cpp @@ -229,6 +229,50 @@ TEST(primal_beziercurve, evaluate) } } +//------------------------------------------------------------------------------ +TEST(primal_beziercurve, degree_elevate_preserves_geometry) +{ + using CoordType = double; + constexpr int DIM = 3; + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + + constexpr CoordType eps = 1e-12; + + PointType data[4] = {PointType {0.6, 1.2, 1.0}, + PointType {1.3, 1.6, 1.8}, + PointType {2.9, 2.4, 2.3}, + PointType {3.2, 3.5, 3.0}}; + + BezierCurveType poly(data, 3); + CoordType weights[4] = {1.0, 0.5, 2.0, 1.25}; + BezierCurveType rat(data, weights, 3); + + auto check_preserve = [&](BezierCurveType curve, int elevate_by) { + const CoordType ts[] = {0.0, 0.1, 0.3, 0.7, 1.0}; + axom::Array before; + before.reserve(5); + for(CoordType t : ts) + { + before.push_back(curve.evaluate(t)); + } + + curve.degreeElevate(elevate_by); + for(int k = 0; k < 5; ++k) + { + const PointType after = curve.evaluate(ts[k]); + EXPECT_NEAR(after[0], before[k][0], eps); + EXPECT_NEAR(after[1], before[k][1], eps); + EXPECT_NEAR(after[2], before[k][2], eps); + } + }; + + check_preserve(poly, 1); + check_preserve(poly, 2); + check_preserve(rat, 1); + check_preserve(rat, 2); +} + //------------------------------------------------------------------------------ TEST(primal_beziercurve_, first_derivatives) { From 21080537e3a448b7c1d9a63f73a812c0e2e39567 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 9 Jun 2026 08:24:38 -0700 Subject: [PATCH 04/16] Add gregory triangle class --- src/axom/primal/CMakeLists.txt | 1 + src/axom/primal/geometry/GregoryPatch.hpp | 25 +- src/axom/primal/geometry/GregoryTriangle.hpp | 576 ++++++++++++++++++ .../primal/tests/primal_bezier_triangle.cpp | 185 ++++++ 4 files changed, 775 insertions(+), 12 deletions(-) create mode 100644 src/axom/primal/geometry/GregoryTriangle.hpp diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 60a44cbb82..de2ac916d5 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -29,6 +29,7 @@ set( primal_headers geometry/Cone.hpp geometry/CurvedPolygon.hpp geometry/GregoryPatch.hpp + geometry/GregoryTriangle.hpp geometry/Hexahedron.hpp geometry/KnotVector.hpp geometry/Line.hpp diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp index d32147c5ea..79efec0301 100644 --- a/src/axom/primal/geometry/GregoryPatch.hpp +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -7,7 +7,7 @@ /*! * \file GregoryPatch.hpp * - * \brief A Gregory patch primitive + * \brief A bicubic Gregory patch primitive */ #ifndef AXOM_PRIMAL_GREGORY_PATCH_HPP_ @@ -44,9 +44,10 @@ std::ostream& operator<<(std::ostream& os, const GregoryPatch& nPatch); /*! * \class GregoryPatch * - * \brief Represents a 3D Gregory patch defined by ... + * \brief Represents a 3D Gregory patch defined by the control points of 4 cubic Bezier curves + * for each boundary, and an additional two "Gregory points" for each edge which + * determine the internal geometry of the surface. * \tparam T the coordinate type, e.g., double, float, etc. - * \tparam NDIMS the number of dimensions */ template class GregoryPatch @@ -145,10 +146,10 @@ class GregoryPatch v1 = getTangent(i, 0); } - PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_index_map[e][k]]; } + PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_edge_index_map[e][k]]; } const PointType& getBoundaryPoint(int e, int k) const { - return m_controlPoints[s_index_map[e][k]]; + return m_controlPoints[s_edge_index_map[e][k]]; } // Evaluate the patch by constructing the equivalent Bezier patch with interior control nodes @@ -327,7 +328,7 @@ class GregoryPatch } private: - struct IntermediateBezierDerivatives + struct IntermediateBlendingDerivatives { BezierPatch bpatch; PointType Q[2][2]; @@ -338,9 +339,9 @@ class GregoryPatch VectorType Q_uv[2][2]; }; - IntermediateBezierDerivatives setup_intermediate_bezier(T u, T v, int derivative_order) const + IntermediateBlendingDerivatives setup_intermediate_bezier(T u, T v, int derivative_order) const { - IntermediateBezierDerivatives out; + IntermediateBlendingDerivatives out; out.bpatch = get_bezier_boundary(); @@ -476,10 +477,10 @@ class GregoryPatch CoordsVec m_controlPoints; // Map of boundary curve control points into internal storage - static constexpr int s_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, - {/*V1*/ 1, /*E11*/ 6, /*E12*/ 7, /*V2*/ 2}, - {/*V2*/ 2, /*E21*/ 8, /*E22*/ 9, /*V3*/ 3}, - {/*V3*/ 3, /*E31*/ 10, /*E32*/ 11, /*V0*/ 0}}; + static constexpr int s_edge_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, + {/*V1*/ 1, /*E11*/ 6, /*E12*/ 7, /*V2*/ 2}, + {/*V2*/ 2, /*E21*/ 8, /*E22*/ 9, /*V3*/ 3}, + {/*V3*/ 3, /*E31*/ 10, /*E32*/ 11, /*V0*/ 0}}; }; //------------------------------------------------------------------------------ diff --git a/src/axom/primal/geometry/GregoryTriangle.hpp b/src/axom/primal/geometry/GregoryTriangle.hpp new file mode 100644 index 0000000000..a8fc538754 --- /dev/null +++ b/src/axom/primal/geometry/GregoryTriangle.hpp @@ -0,0 +1,576 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file GregoryTriangle.hpp + * + * \brief A bicubic Gregory triangle primitive + */ + +#ifndef AXOM_PRIMAL_GREGORY_TRIANGLE_HPP_ +#define AXOM_PRIMAL_GREGORY_TRIANGLE_HPP_ + +#include "axom/core.hpp" +#include "axom/slic.hpp" + +#include "axom/core/NumericArray.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Vector.hpp" +#include "axom/primal/geometry/Segment.hpp" +#include "axom/primal/geometry/BezierTriangle.hpp" +#include "axom/primal/geometry/BoundingBox.hpp" +#include "axom/primal/geometry/OrientedBoundingBox.hpp" + +#include +#include + +#include "axom/fmt.hpp" + +namespace axom +{ +namespace primal +{ +// Forward declare the templated classes and operator functions +template +class GregoryTriangle; + +/*! \brief Overloaded output operator for Gregory Triangles*/ +template +std::ostream& operator<<(std::ostream& os, const GregoryTriangle& nTri); + +/*! + * \class GregoryTriangle + * + * \brief Represents a 3D Gregory triangle defined by the control points of 3 degree-elevated + * cubic Bezier curves (i.e. quartic curves with identical geometry to cubics), and + * an additional two "Gregory points" for each edge which determine the surface. + * + * Degree elevation of the boundary curves is necessary to provide sufficient degrees of freedom + * for the blending of the internal nodes. + * + * \tparam T the coordinate type, e.g., double, float, etc. + */ +template +class GregoryTriangle +{ +public: + // The number of control points for a hybrid quartic-cubic Gregory triangle is fixed: + // - 12 exterior control points for each of three degree-elevated cubic curves + // - 6 interior control points for each of 3 boundary curves + static constexpr int NPTS = 18; + + using PointType = Point; + using VectorType = Vector; + + using CoordsVec = axom::StackArray; + using BezierType = BezierCurve; + + using BoundingBoxType = BoundingBox; + using OrientedBoundingBoxType = OrientedBoundingBox; + + AXOM_STATIC_ASSERT_MSG(std::is_arithmetic::value, + "A Gregory Triangle must be defined using an arithmetic type"); + +public: + GregoryTriangle() = default; + + explicit GregoryTriangle(ArrayView controlPoints) + { + SLIC_ASSERT(controlPoints.size() == NPTS); + for(int i = 0; i < NPTS; ++i) + { + m_controlPoints[i] = controlPoints[i]; + } + } + + // static constexpr int nDeg = 4; + // static constexpr int nVerts = 3; + // static constexpr int nControlPoints = 18; + // typedef std::array NodePositions; + // typedef std::array NodeVectors; + // typedef std::array NodeUnitVectors; + // typedef std::array ControlPoints; + + GregoryTriangle(ArrayView nodePositions, ArrayView nodeVectors) + { + // Store the position and orthogonal unit vector at each corner + SLIC_ASSERT(nodePositions.size() == 3); + SLIC_ASSERT(nodeVectors.size() == 3); + + axom::Array v(4); + for(int i = 0; i < 3; ++i) + { + getCorner(i) = nodePositions[i]; + v[i] = nodeVectors[i].unitVector(); + } + + // Initialize the boundary points for the quartic Bezier triangle + axom::StackArray, 3> cubic_deriv_cp; + for(int k = 0; k < 3; ++k) // Loop over edges + { + const int kp1 = (k + 1) % 3; + const VectorType dx(nodePositions[k], nodePositions[kp1]); + + const VectorType c0 = (dx - dx.dot(v[k]) * v[k]) / 3.0; + const VectorType c2 = (dx - dx.dot(v[kp1]) * v[kp1]) / 3.0; + + // Define the cubic Bezier which represents the boundary of the curve + BezierType cubic(axom::Array {nodePositions[k], + PointType {nodePositions[k].array() + c0.array()}, + PointType {nodePositions[kp1].array() - c2.array()}, + nodePositions[kp1]}, + 3); + + // Store the control points of the derivative of this cubic for later + cubic_deriv_cp[k][0] = VectorType(cubic[0], cubic[1]); + cubic_deriv_cp[k][1] = VectorType(cubic[1], cubic[2]); + cubic_deriv_cp[k][2] = VectorType(cubic[2], cubic[3]); + + // Do degree elevation on the cubic curve, which defines the quartic boundary control points + cubic.degreeElevate(); + getBoundaryPoint(k, 0) = cubic[0]; + getBoundaryPoint(k, 1) = cubic[1]; + getBoundaryPoint(k, 2) = cubic[2]; + getBoundaryPoint(k, 3) = cubic[3]; + // getBoundaryPoint(i, 4) will be set for the next edge + } + + // Define the interior control nodse at each vertex + for(int k = 0; k < 3; ++k) + { + const int km1 = (k + 2) % 3; + const int kp1 = (k + 1) % 3; + + // Get control points at and around the edge's starting vertex + const auto& p0 = getBoundaryPoint(km1, 3); // Vertex - 1 + const auto& q0 = getBoundaryPoint(k, 0); // Vertex + const auto& q1 = getBoundaryPoint(k, 1); // Vertex + 1 + const auto deriv0 = 0.5 * VectorType(p0, q0) + 0.5 * VectorType(p0, q1); + + // Get control points at and around the edge's ending vertex + const auto& q3 = getBoundaryPoint(k, 3); // Vertex - 1 + const auto& q4 = getBoundaryPoint(k, 4); // Vertex + const auto& p3 = getBoundaryPoint(kp1, 1); // Vertex + 1 + const auto deriv1 = 0.5 * VectorType(p3, q3) + 0.5 * VectorType(p3, q4); + + // Compute a boundary cross derivative that varies across the edge + const VectorType dx(getCorner(k), getCorner(kp1)); + const VectorType a0 = VectorType::cross_product(nodeVectors[k], dx).unitVector(); + const VectorType a3 = VectorType::cross_product(nodeVectors[kp1], dx).unitVector(); + + // Elevate the linear cross derivative a(t) = (1-t)a0 + t*a3 to quadratic + const axom::StackArray aHat = {a0, 0.5 * (a0 + a3), a3}; + + // Compute the CPs for blending functions k(t) = (1-t)*k0 + t*k1 and + // h(t) = (1-t)*h- + t*h1 + auto& c = cubic_deriv_cp[k]; + const double k0 = aHat[0].dot(deriv0); + const double k1 = aHat[2].dot(deriv1); + + const double h0 = c[0].dot(deriv0) / c[0].dot(c[0]); + const double h1 = c[2].dot(deriv1) / c[2].dot(c[2]); + + // Compute cross derivatives at edge interior points and use them for interior CP + axom::StackArray deriv; + for(int j = 1; j < 3; ++j) + { + const double fac = j / 3.0; + deriv[j - 1] = + (1. - fac) * (k0 * aHat[j] + h0 * c[j]) + fac * (k1 * aHat[j - 1] + h1 * c[j - 1]); + } + + const auto& q2 = getBoundaryPoint(k, 2); + getTangent(k, 0) = PointType(PointType::lerp(q1, q2, 0.5).array() - deriv[0].array()); + getTangent(k, 1) = PointType(PointType::lerp(q2, q3, 0.5).array() - deriv[1].array()); + } + } + + PointType& getCorner(int i) { return m_controlPoints[i]; } + const PointType& getCorner(int i) const { return m_controlPoints[i]; } + + PointType& getTangent(int e, int t) { return m_controlPoints[12 + 2 * e + t]; } + const PointType& getTangent(int e, int t) const { return m_controlPoints[12 + 2 * e + t]; } + + void getTangentsByCorner(int i, PointType& v0, PointType& v1) const + { + v0 = getTangent((i + 2) % 3, 1); + v1 = getTangent(i, 0); + } + + PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_edge_index_map[e][k]]; } + const PointType& getBoundaryPoint(int e, int k) const + { + return m_controlPoints[s_edge_index_map[e][k]]; + } + + // Evaluate the triangle by constructing the equivalent Bezier triangle with interior control nodes + // defined in terms of the tangent vectors and the evaluation parameters + PointType evaluate(T u0, T v0) const + { + const auto intermediate = setup_intermediate_bezier(u0, v0, 0); + return intermediate.btri.evaluate(u0, v0); + } + + /*! + * \brief Evaluates all first derivatives of the Gregory patch at (\a u, \a v) + * + * \param [in] u Parameter value at which to evaluate on the first axis + * \param [in] v Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory patch at (u, v) + * \param [out] Du The vector value of S_u(u, v) + * \param [out] Dv The vector value of S_v(u, v) + */ + void evaluateFirstDerivatives(T u0, T v0, PointType& eval, VectorType& Du, VectorType& Dv) const + { + const auto intermediate = setup_intermediate_bezier(u0, v0, 1); + intermediate.btri.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); + + // Chain rule correction due to (u0,v0)-dependent interior control points. + // For a quartic Bezier triangle, the basis weight of control point (i,j) + // is: 4!/(i! j! k!) * u^i v^j w^k with k = 4-i-j and barycentric weights + // (lambda0, lambda1, lambda2) = (1-u0-v0, v0, u0) corresponding to vertices + // (0,0), (0,order), (order,0), respectively. + const T u = u0; + const T v = v0; + const T w = T(1) - u0 - v0; + + const T w11 = T(12) * u * v * w * w; // (i,j)=(1,1), k=2 + const T w21 = T(12) * u * u * v * w; // (i,j)=(2,1), k=1 + const T w12 = T(12) * u * v * v * w; // (i,j)=(1,2), k=1 + + Du += w11 * intermediate.Q_u[0] + w21 * intermediate.Q_u[1] + w12 * intermediate.Q_u[2]; + Dv += w11 * intermediate.Q_v[0] + w21 * intermediate.Q_v[1] + w12 * intermediate.Q_v[2]; + } + + /*! + * \brief Evaluates all second derivatives of the Gregory patch at (\a u, \a v) + * + * \param [in] u Parameter value at which to evaluate on the first axis + * \param [in] v Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory patch at (u, v) + * \param [out] Du The vector value of S_u(u, v) + * \param [out] Dv The vector value of S_v(u, v) + * \param [out] DuDu The vector value of S_uu(u, v) + * \param [out] DvDv The vector value of S_vv(u, v) + * \param [out] DuDv The vector value of S_uv(u, v) == S_vu(u, v) + */ + // void evaluateSecondDerivatives(T u, + // T v, + // PointType& eval, + // VectorType& Du, + // VectorType& Dv, + // VectorType& DuDu, + // VectorType& DvDv, + // VectorType& DuDv) const + // { + // const auto intermediate = setup_intermediate_bezier(u, v, 2); + // intermediate.bpatch.evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + + // // Chain rule correction due to (u,v)-dependent interior control points + // axom::StaticArray Bu, dBu, Bv, dBv; + // evaluateCubicBernstein(u, Bu, dBu); + // evaluateCubicBernstein(v, Bv, dBv); + + // const T w11 = Bu[1] * Bv[1]; + // const T w21 = Bu[2] * Bv[1]; + // const T w12 = Bu[1] * Bv[2]; + // const T w22 = Bu[2] * Bv[2]; + + // const T wu11 = dBu[1] * Bv[1]; + // const T wu21 = dBu[2] * Bv[1]; + // const T wu12 = dBu[1] * Bv[2]; + // const T wu22 = dBu[2] * Bv[2]; + + // const T wv11 = Bu[1] * dBv[1]; + // const T wv21 = Bu[2] * dBv[1]; + // const T wv12 = Bu[1] * dBv[2]; + // const T wv22 = Bu[2] * dBv[2]; + + // // First derivative corrections + // Du += w11 * intermediate.Q_u[0][0] + w21 * intermediate.Q_u[1][0] + + // w12 * intermediate.Q_u[0][1] + w22 * intermediate.Q_u[1][1]; + // Dv += w11 * intermediate.Q_v[0][0] + w21 * intermediate.Q_v[1][0] + + // w12 * intermediate.Q_v[0][1] + w22 * intermediate.Q_v[1][1]; + + // // Second derivative corrections + // DuDu += T(2) * + // (wu11 * intermediate.Q_u[0][0] + wu21 * intermediate.Q_u[1][0] + + // wu12 * intermediate.Q_u[0][1] + wu22 * intermediate.Q_u[1][1]) + + // (w11 * intermediate.Q_uu[0][0] + w21 * intermediate.Q_uu[1][0] + + // w12 * intermediate.Q_uu[0][1] + w22 * intermediate.Q_uu[1][1]); + + // DvDv += T(2) * + // (wv11 * intermediate.Q_v[0][0] + wv21 * intermediate.Q_v[1][0] + + // wv12 * intermediate.Q_v[0][1] + wv22 * intermediate.Q_v[1][1]) + + // (w11 * intermediate.Q_vv[0][0] + w21 * intermediate.Q_vv[1][0] + + // w12 * intermediate.Q_vv[0][1] + w22 * intermediate.Q_vv[1][1]); + + // DuDv += (wu11 * intermediate.Q_v[0][0] + wu21 * intermediate.Q_v[1][0] + + // wu12 * intermediate.Q_v[0][1] + wu22 * intermediate.Q_v[1][1]) + + // (wv11 * intermediate.Q_u[0][0] + wv21 * intermediate.Q_u[1][0] + + // wv12 * intermediate.Q_u[0][1] + wv22 * intermediate.Q_u[1][1]) + + // (w11 * intermediate.Q_uv[0][0] + w21 * intermediate.Q_uv[1][0] + + // w12 * intermediate.Q_uv[0][1] + w22 * intermediate.Q_uv[1][1]); + // } + + // VectorType du(T u, T v) const + // { + // PointType eval; + // VectorType Du, Dv; + // evaluateFirstDerivatives(u, v, eval, Du, Dv); + // return Du; + // } + + // VectorType dv(T u, T v) const + // { + // PointType eval; + // VectorType Du, Dv; + // evaluateFirstDerivatives(u, v, eval, Du, Dv); + // return Dv; + // } + + // VectorType dudu(T u, T v) const + // { + // PointType eval; + // VectorType Du, Dv, DuDu, DvDv, DuDv; + // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + // return DuDu; + // } + + // VectorType dvdv(T u, T v) const + // { + // PointType eval; + // VectorType Du, Dv, DuDu, DvDv, DuDv; + // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + // return DvDv; + // } + + // VectorType dudv(T u, T v) const + // { + // PointType eval; + // VectorType Du, Dv, DuDu, DvDv, DuDv; + // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + // return DuDv; + // } + + /// \brief Returns an axis-aligned bounding box containing the patch + BoundingBoxType boundingBox() const + { + return BoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); + } + + /// \brief Returns an oriented bounding box containing the patch + OrientedBoundingBoxType orientedBoundingBox() const + { + return OrientedBoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); + } + + void print(std::ostream& os) const + { + os << "GregoryTriangle("; + for(int i = 0; i < NPTS; ++i) + { + os << m_controlPoints[i]; + if(i + 1 < NPTS) + { + os << ", "; + } + } + os << ")"; + } + +private: + struct IntermediateBlendingDerivatives + { + BezierTriangle btri; + PointType Q[3]; + VectorType Q_u[3]; + VectorType Q_v[3]; + }; + + IntermediateBlendingDerivatives setup_intermediate_bezier(T u0, T v0, int derivative_order) const + { + IntermediateBlendingDerivatives out; + out.btri = get_bezier_boundary(); + + // Parameter convention matches BezierTriangle::evaluate(): + // barycentric weights (lambda0, lambda1, lambda2) are (1-u0-v0, v0, u0) + const T u = u0; + const T v = v0; + const T w = T(1) - u0 - v0; + + auto blend = [&](const PointType& A, + const PointType& B, + T wa, + T wb, + T wa_u0, + T wb_u0, + T wa_v0, + T wb_v0, + PointType& Q, + VectorType& Q_u0, + VectorType& Q_v0) { + const T denom = wa + wb; + if(axom::utilities::isNearlyEqual(denom, T(0))) + { + Q = A; + Q_u0 = VectorType(T(0)); + Q_v0 = VectorType(T(0)); + return; + } + + Q = PointType((wa * A.array() + wb * B.array()) / denom); + + if(derivative_order >= 1) + { + const auto dQ_u0 = + (wa_u0 * (A.array() - Q.array()) + wb_u0 * (B.array() - Q.array())) / denom; + const auto dQ_v0 = + (wa_v0 * (A.array() - Q.array()) + wb_v0 * (B.array() - Q.array())) / denom; + Q_u0 = VectorType(dQ_u0); + Q_v0 = VectorType(dQ_v0); + } + else + { + Q_u0 = VectorType(T(0)); + Q_v0 = VectorType(T(0)); + } + }; + + // These are the three (u0,v0)-dependent interior control points for the + // equivalent quartic Bezier triangle. Each is a blend of the two Gregory + // points adjacent to the corresponding vertex. + // + // By convention, edge k connects vertex k -> vertex (k+1)%3. + // - getTangent(k,0) is the Gregory point near the starting vertex k + // - getTangent(k,1) is the Gregory point near the ending vertex (k+1)%3 + // + // Control net slots: + // Q[0] -> btri(1,1) + // Q[1] -> btri(2,1) + // Q[2] -> btri(1,2) + blend(getTangent(2, 1), + getTangent(0, 0), + u, + v, + T(1), + T(0), + T(0), + T(1), + out.Q[0], + out.Q_u[0], + out.Q_v[0]); + + blend(getTangent(1, 1), + getTangent(2, 0), + v, + w, + T(0), + T(-1), + T(1), + T(-1), + out.Q[1], + out.Q_u[1], + out.Q_v[1]); + + blend(getTangent(0, 1), + getTangent(1, 0), + w, + u, + T(-1), + T(1), + T(-1), + T(0), + out.Q[2], + out.Q_u[2], + out.Q_v[2]); + + setBezierInterior(out.btri, out.Q); + return out; + } + + static void setBezierInterior(BezierTriangle& btri, const PointType Q[3]) + { + btri(1, 1) = Q[0]; + btri(2, 1) = Q[1]; + btri(1, 2) = Q[2]; + } + + // Copies over the boundary points to a BezierPatch object, + // leaving the 4 interior control points uninitialized + BezierTriangle get_bezier_boundary() const + { + BezierTriangle btri(4); + + // Edge 0 + btri(0, 0) = getBoundaryPoint(0, 0); + btri(0, 1) = getBoundaryPoint(0, 1); + btri(0, 2) = getBoundaryPoint(0, 2); + btri(0, 3) = getBoundaryPoint(0, 3); + btri(0, 4) = getBoundaryPoint(0, 4); + + // Edge 1 + // btri(0, 4) = getBoundaryPOint(1, 0); + btri(1, 3) = getBoundaryPoint(1, 1); + btri(2, 2) = getBoundaryPoint(1, 2); + btri(3, 1) = getBoundaryPoint(1, 3); + btri(4, 0) = getBoundaryPoint(1, 4); + + // Edge 2 + // btri(4, 0) = getBoundaryPoint(2, 0); + btri(3, 0) = getBoundaryPoint(2, 1); + btri(2, 0) = getBoundaryPoint(2, 2); + btri(1, 0) = getBoundaryPoint(2, 3); + // btri(0, 0) = getBoundaryPoint(2, 4); + + return btri; + } + + // static void evaluateCubicBernstein(T t, axom::StaticArray& B, axom::StaticArray& dB) + // { + // const T tm = T(1) - t; + // const T tm2 = tm * tm; + // const T t2 = t * t; + + // B[0] = tm2 * tm; + // B[1] = T(3) * t * tm2; + // B[2] = T(3) * t2 * tm; + // B[3] = t2 * t; + + // dB[0] = -T(3) * tm2; + // dB[1] = T(3) * tm2 - T(6) * t * tm; + // dB[2] = T(6) * t * tm - T(3) * t2; + // dB[3] = T(3) * t2; + // } + + CoordsVec m_controlPoints; + + // Map of boundary curve control points into internal storage + static constexpr int s_edge_index_map[3][5] = { + {/*V0*/ 0, /*E01*/ 3, /*E02*/ 4, /*E03*/ 5, /*V1*/ 1}, + {/*V0*/ 1, /*E01*/ 6, /*E02*/ 7, /*E03*/ 8, /*V1*/ 2}, + {/*V0*/ 2, /*E01*/ 9, /*E02*/ 10, /*E03*/ 11, /*V1*/ 0}}; +}; + +//------------------------------------------------------------------------------ +/// Free functions related to GregoryTriangle +//------------------------------------------------------------------------------ +template +std::ostream& operator<<(std::ostream& os, const GregoryTriangle& nPatch) +{ + nPatch.print(os); + return os; +} + +} // namespace primal +} // namespace axom + +/// Overload to format a primal::GregoryTriangle using fmt +template +struct axom::fmt::formatter> : ostream_formatter +{ }; + +#endif // AXOM_PRIMAL_GREGORY_TRIANGLE_HPP_ diff --git a/src/axom/primal/tests/primal_bezier_triangle.cpp b/src/axom/primal/tests/primal_bezier_triangle.cpp index 838aa1cfea..3878ceb137 100644 --- a/src/axom/primal/tests/primal_bezier_triangle.cpp +++ b/src/axom/primal/tests/primal_bezier_triangle.cpp @@ -14,6 +14,8 @@ #include "axom/slic.hpp" #include "axom/primal/geometry/BezierTriangle.hpp" +#include "axom/primal/geometry/GregoryTriangle.hpp" +#include "axom/primal/geometry/GregoryPatch.hpp" #include @@ -270,6 +272,189 @@ TEST(primal_beziertriangle, evaluate_linear) } } +//------------------------------------------------------------------------------ +TEST(primal_beziertriangle, parameter_convention_matches_barycentric) +{ + using CoordType = double; + using BTri = primal::BezierTriangle; + using PointType = BTri::PointType; + using Barycentric = BTri::Barycentric; + + constexpr CoordType eps = 1e-14; + + BTri tri(1); + const PointType pA {1.0, 2.0, 3.0}; // (i,j) = (0,0) + const PointType pB {4.0, -1.0, 0.5}; // (i,j) = (0,1) + const PointType pC {-2.0, 0.25, 7.0}; // (i,j) = (1,0) + + tri(0, 0) = pA; + tri(0, 1) = pB; + tri(1, 0) = pC; + + // Corner mapping implied by the (u0,v0)->barycentric convention + EXPECT_EQ(tri.evaluate(0.0, 0.0), pA); // bary = (1,0,0) + EXPECT_EQ(tri.evaluate(0.0, 1.0), pB); // bary = (0,1,0) + EXPECT_EQ(tri.evaluate(1.0, 0.0), pC); // bary = (0,0,1) + + // Interior point check: evaluate(u0,v0) == lambda0*A + lambda1*B + lambda2*C, + // with (lambda0,lambda1,lambda2) = (1-u0-v0, v0, u0) + const CoordType u0 = 0.3; + const CoordType v0 = 0.2; + const Barycentric bary {1.0 - u0 - v0, v0, u0}; + + const PointType expected = BTri::triInterpolate(pA, pB, pC, bary); + const PointType eval = tri.evaluate(u0, v0); + for(int d = 0; d < 3; ++d) + { + EXPECT_NEAR(eval[d], expected[d], eps); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_beziertriangle, finite_difference_first_derivatives) +{ + using CoordType = double; + using BTri = primal::BezierTriangle; + using PointType = BTri::PointType; + using VectorType = BTri::VectorType; + + constexpr CoordType h = 1e-7; + constexpr CoordType tol = 5e-6; + + BTri tri(4); + fillControlNet(tri); + + const CoordType u0 = 0.23; + const CoordType v0 = 0.31; + + PointType eval; + VectorType Du, Dv; + tri.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); + + const PointType fp_u = tri.evaluate(u0 + h, v0); + const PointType fm_u = tri.evaluate(u0 - h, v0); + const PointType fp_v = tri.evaluate(u0, v0 + h); + const PointType fm_v = tri.evaluate(u0, v0 - h); + + const VectorType Du_fd(fp_u, fm_u); + const VectorType Dv_fd(fp_v, fm_v); + + for(int d = 0; d < 3; ++d) + { + EXPECT_NEAR(Du[d], Du_fd[d] / (2 * h), tol); + EXPECT_NEAR(Dv[d], Dv_fd[d] / (2 * h), tol); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_beziertriangle, gregorytriangle_finite_difference_first_derivatives) +{ + using CoordType = double; + using PointType = primal::Point; + using VectorType = primal::Vector; + using GTri = primal::GregoryTriangle; + + constexpr CoordType h = 1e-7; + constexpr CoordType tol = 5e-5; + + const std::array corners = { + PointType {0.0, 0.0, 0.0}, + PointType {1.0, 0.0, 0.2}, + PointType {0.0, 1.0, -0.1}, + }; + + const std::array normals = { + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + }; + + GTri gtri(axom::ArrayView(corners.data(), 3), + axom::ArrayView(normals.data(), 3)); + + const CoordType u0 = 0.21; + const CoordType v0 = 0.27; + + PointType eval; + VectorType Du, Dv; + gtri.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); + + const PointType fp_u = gtri.evaluate(u0 + h, v0); + const PointType fm_u = gtri.evaluate(u0 - h, v0); + const PointType fp_v = gtri.evaluate(u0, v0 + h); + const PointType fm_v = gtri.evaluate(u0, v0 - h); + + const VectorType Du_fd(fp_u, fm_u); + const VectorType Dv_fd(fp_v, fm_v); + + for(int d = 0; d < 3; ++d) + { + EXPECT_NEAR(Du[d], Du_fd[d] / (2 * h), tol); + EXPECT_NEAR(Dv[d], Dv_fd[d] / (2 * h), tol); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_beziertriangle, gregorytriangle_parameter_convention_matches_beziertriangle) +{ + using CoordType = double; + using PointType = primal::Point; + using VectorType = primal::Vector; + using GTri = primal::GregoryTriangle; + + constexpr CoordType eps = 1e-12; + + const std::array corners = { + PointType {0.0, 0.0, 0.0}, + PointType {1.0, 0.0, 0.0}, + PointType {0.0, 1.0, 0.0}, + }; + + const std::array normals = { + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + }; + + GTri gtri(axom::ArrayView(corners.data(), 3), + axom::ArrayView(normals.data(), 3)); + + // Corner mapping should match BezierTriangle::evaluate convention + EXPECT_TRUE(gtri.evaluate(0.0, 0.0).isNearlyEqual(corners[0], eps)); + EXPECT_TRUE(gtri.evaluate(0.0, 1.0).isNearlyEqual(corners[1], eps)); + EXPECT_TRUE(gtri.evaluate(1.0, 0.0).isNearlyEqual(corners[2], eps)); + + // Edge mapping: verify straight edges for this planar configuration + auto check_near = [&](const PointType& a, const PointType& b) { + EXPECT_NEAR(a[0], b[0], eps); + EXPECT_NEAR(a[1], b[1], eps); + EXPECT_NEAR(a[2], b[2], eps); + }; + + const CoordType tvals[] = {0.0, 0.2, 0.6, 1.0}; + + // edge u0 = 0: from corner0 -> corner1 as v0 goes 0->1 + for(const CoordType t : tvals) + { + const PointType expected = PointType::lerp(corners[0], corners[1], t); + check_near(gtri.evaluate(0.0, t), expected); + } + + // edge v0 = 0: from corner0 -> corner2 as u0 goes 0->1 + for(const CoordType t : tvals) + { + const PointType expected = PointType::lerp(corners[0], corners[2], t); + check_near(gtri.evaluate(t, 0.0), expected); + } + + // edge u0 + v0 = 1: from corner1 -> corner2 as u0 goes 0->1 (v0 goes 1->0) + for(const CoordType t : tvals) + { + const PointType expected = PointType::lerp(corners[1], corners[2], t); + check_near(gtri.evaluate(t, 1.0 - t), expected); + } +} + //------------------------------------------------------------------------------ TEST(primal_beziertriangle, evaluate_quadratic_second_derivatives_constant) { From 40cfe45a36b0901f20b7c074e0feb2357f69cb99 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Jul 2026 12:55:52 -0700 Subject: [PATCH 05/16] tidy implementation a little --- src/axom/primal/geometry/GregoryPatch.hpp | 6 + src/axom/primal/geometry/GregoryTriangle.hpp | 115 +++---------------- 2 files changed, 22 insertions(+), 99 deletions(-) diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp index 79efec0301..211ef94baf 100644 --- a/src/axom/primal/geometry/GregoryPatch.hpp +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -81,6 +81,8 @@ class GregoryPatch } } + // Constructor from corner nodes and corner vectors. + // Currently not configured to handle any connectivity between patches GregoryPatch(ArrayView nodePositions, ArrayView nodeVectors) { // Store the position and orthogonal unit vector at each corner @@ -154,6 +156,10 @@ class GregoryPatch // Evaluate the patch by constructing the equivalent Bezier patch with interior control nodes // defined in terms of the tangent vectors and the evaluation parameters + // + // A cubic Gregory patch can be evaluated by using the 8 internal, blending nodes to construct + // a set of control points which, when used with the Gregory patch edge nodes, produces + // a polynomial bicubic Bezier patch that is equivalent at the requested parameters PointType evaluate(T u, T v) const { const auto intermediate = setup_intermediate_bezier(u, v, 0); diff --git a/src/axom/primal/geometry/GregoryTriangle.hpp b/src/axom/primal/geometry/GregoryTriangle.hpp index a8fc538754..7a56e95624 100644 --- a/src/axom/primal/geometry/GregoryTriangle.hpp +++ b/src/axom/primal/geometry/GregoryTriangle.hpp @@ -86,14 +86,6 @@ class GregoryTriangle } } - // static constexpr int nDeg = 4; - // static constexpr int nVerts = 3; - // static constexpr int nControlPoints = 18; - // typedef std::array NodePositions; - // typedef std::array NodeVectors; - // typedef std::array NodeUnitVectors; - // typedef std::array ControlPoints; - GregoryTriangle(ArrayView nodePositions, ArrayView nodeVectors) { // Store the position and orthogonal unit vector at each corner @@ -245,18 +237,7 @@ class GregoryTriangle Dv += w11 * intermediate.Q_v[0] + w21 * intermediate.Q_v[1] + w12 * intermediate.Q_v[2]; } - /*! - * \brief Evaluates all second derivatives of the Gregory patch at (\a u, \a v) - * - * \param [in] u Parameter value at which to evaluate on the first axis - * \param [in] v Parameter value at which to evaluate on the second axis - * \param [out] eval The point value of the Gregory patch at (u, v) - * \param [out] Du The vector value of S_u(u, v) - * \param [out] Dv The vector value of S_v(u, v) - * \param [out] DuDu The vector value of S_uu(u, v) - * \param [out] DvDv The vector value of S_vv(u, v) - * \param [out] DuDv The vector value of S_uv(u, v) == S_vu(u, v) - */ + // Not yet completed // void evaluateSecondDerivatives(T u, // T v, // PointType& eval, @@ -266,72 +247,25 @@ class GregoryTriangle // VectorType& DvDv, // VectorType& DuDv) const // { - // const auto intermediate = setup_intermediate_bezier(u, v, 2); - // intermediate.bpatch.evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); - - // // Chain rule correction due to (u,v)-dependent interior control points - // axom::StaticArray Bu, dBu, Bv, dBv; - // evaluateCubicBernstein(u, Bu, dBu); - // evaluateCubicBernstein(v, Bv, dBv); - - // const T w11 = Bu[1] * Bv[1]; - // const T w21 = Bu[2] * Bv[1]; - // const T w12 = Bu[1] * Bv[2]; - // const T w22 = Bu[2] * Bv[2]; - - // const T wu11 = dBu[1] * Bv[1]; - // const T wu21 = dBu[2] * Bv[1]; - // const T wu12 = dBu[1] * Bv[2]; - // const T wu22 = dBu[2] * Bv[2]; - - // const T wv11 = Bu[1] * dBv[1]; - // const T wv21 = Bu[2] * dBv[1]; - // const T wv12 = Bu[1] * dBv[2]; - // const T wv22 = Bu[2] * dBv[2]; - - // // First derivative corrections - // Du += w11 * intermediate.Q_u[0][0] + w21 * intermediate.Q_u[1][0] + - // w12 * intermediate.Q_u[0][1] + w22 * intermediate.Q_u[1][1]; - // Dv += w11 * intermediate.Q_v[0][0] + w21 * intermediate.Q_v[1][0] + - // w12 * intermediate.Q_v[0][1] + w22 * intermediate.Q_v[1][1]; - - // // Second derivative corrections - // DuDu += T(2) * - // (wu11 * intermediate.Q_u[0][0] + wu21 * intermediate.Q_u[1][0] + - // wu12 * intermediate.Q_u[0][1] + wu22 * intermediate.Q_u[1][1]) + - // (w11 * intermediate.Q_uu[0][0] + w21 * intermediate.Q_uu[1][0] + - // w12 * intermediate.Q_uu[0][1] + w22 * intermediate.Q_uu[1][1]); - - // DvDv += T(2) * - // (wv11 * intermediate.Q_v[0][0] + wv21 * intermediate.Q_v[1][0] + - // wv12 * intermediate.Q_v[0][1] + wv22 * intermediate.Q_v[1][1]) + - // (w11 * intermediate.Q_vv[0][0] + w21 * intermediate.Q_vv[1][0] + - // w12 * intermediate.Q_vv[0][1] + w22 * intermediate.Q_vv[1][1]); - - // DuDv += (wu11 * intermediate.Q_v[0][0] + wu21 * intermediate.Q_v[1][0] + - // wu12 * intermediate.Q_v[0][1] + wu22 * intermediate.Q_v[1][1]) + - // (wv11 * intermediate.Q_u[0][0] + wv21 * intermediate.Q_u[1][0] + - // wv12 * intermediate.Q_u[0][1] + wv22 * intermediate.Q_u[1][1]) + - // (w11 * intermediate.Q_uv[0][0] + w21 * intermediate.Q_uv[1][0] + - // w12 * intermediate.Q_uv[0][1] + w22 * intermediate.Q_uv[1][1]); // } - // VectorType du(T u, T v) const - // { - // PointType eval; - // VectorType Du, Dv; - // evaluateFirstDerivatives(u, v, eval, Du, Dv); - // return Du; - // } + VectorType du(T u, T v) const + { + PointType eval; + VectorType Du, Dv; + evaluateFirstDerivatives(u, v, eval, Du, Dv); + return Du; + } - // VectorType dv(T u, T v) const - // { - // PointType eval; - // VectorType Du, Dv; - // evaluateFirstDerivatives(u, v, eval, Du, Dv); - // return Dv; - // } + VectorType dv(T u, T v) const + { + PointType eval; + VectorType Du, Dv; + evaluateFirstDerivatives(u, v, eval, Du, Dv); + return Dv; + } + // Not yet finished // VectorType dudu(T u, T v) const // { // PointType eval; @@ -529,23 +463,6 @@ class GregoryTriangle return btri; } - // static void evaluateCubicBernstein(T t, axom::StaticArray& B, axom::StaticArray& dB) - // { - // const T tm = T(1) - t; - // const T tm2 = tm * tm; - // const T t2 = t * t; - - // B[0] = tm2 * tm; - // B[1] = T(3) * t * tm2; - // B[2] = T(3) * t2 * tm; - // B[3] = t2 * t; - - // dB[0] = -T(3) * tm2; - // dB[1] = T(3) * tm2 - T(6) * t * tm; - // dB[2] = T(6) * t * tm - T(3) * t2; - // dB[3] = T(3) * t2; - // } - CoordsVec m_controlPoints; // Map of boundary curve control points into internal storage From 52279613dbe41c280d6e652a0657ed1b2dc1b70d Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 13 Aug 2026 13:03:07 -0700 Subject: [PATCH 06/16] Rest of gregory methods --- src/axom/primal/geometry/GregoryPatch.hpp | 290 +++++++++++++- src/axom/primal/geometry/GregoryTriangle.hpp | 384 +++++++++++-------- 2 files changed, 498 insertions(+), 176 deletions(-) diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp index 211ef94baf..1d62957709 100644 --- a/src/axom/primal/geometry/GregoryPatch.hpp +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -70,19 +70,136 @@ class GregoryPatch "A Gregory Patch must be defined using an arithmetic type"); public: + ///@{ + /** + * @name Constructors for GregoryPatch + * + * The constructors allow initialization from: + * - the 20 Gregory patch control points, + * - a polynomial bicubic Bezier patch, + * - C-style arrays, Axom StackArrays, or Axom ArrayViews, + * - four corner positions with associated corner normal vectors. + * + * The 20-point control net is stored as: + * - indices 0-3: corners, + * - indices 4-11: two boundary control points for each edge, + * - indices 12-19: two Gregory tangent points for each edge. + * + * Boundary edge \a e is directed from corner \a e to corner `(e+1)%4`. + */ + + /*! + * \brief Default constructor for a Gregory patch + * + * The fixed-size control net is default-initialized. + */ GregoryPatch() = default; + /*! + * \brief Constructor from an ArrayView over the control points + * + * \param [in] controlPoints ArrayView of the 20 Gregory patch control points + * \pre \a controlPoints must contain exactly `NPTS` points + */ explicit GregoryPatch(ArrayView controlPoints) { SLIC_ASSERT(controlPoints.size() == NPTS); + SLIC_ASSERT(controlPoints.data() != nullptr); for(int i = 0; i < NPTS; ++i) { m_controlPoints[i] = controlPoints[i]; } } - // Constructor from corner nodes and corner vectors. - // Currently not configured to handle any connectivity between patches + /*! + * \brief Constructor from a non-const ArrayView over the control points + * + * \param [in] controlPoints ArrayView of the 20 Gregory patch control points + * \pre \a controlPoints must contain exactly `NPTS` points + */ + explicit GregoryPatch(ArrayView controlPoints) + : GregoryPatch(ArrayView(controlPoints.data(), controlPoints.size())) + { } + + /*! + * \brief Constructor from a C-style array of control points + * + * \param [in] pts A C-style array of 20 Gregory patch control points + * \pre \a pts must be non-null and contain at least `NPTS` points + */ + explicit GregoryPatch(const PointType* pts) : GregoryPatch(ArrayView(pts, NPTS)) + { } + + /*! + * \brief Constructor from a C-style array of control points + * + * \param [in] pts A C-style array of 20 Gregory patch control points + * \pre \a pts must be non-null and contain at least `NPTS` points + */ + explicit GregoryPatch(PointType* pts) : GregoryPatch(ArrayView(pts, NPTS)) { } + + /*! + * \brief Constructor from an Axom StackArray of control points + * + * \param [in] pts StackArray containing the 20 Gregory patch control points + */ + explicit GregoryPatch(const CoordsVec& pts) + : GregoryPatch(ArrayView(pts.data(), pts.size())) + { } + + /*! + * \brief Constructor from a polynomial bicubic Bezier patch + * + * \param [in] bPatch A polynomial Bezier patch of order (3, 3) + * + * This creates a Gregory patch that exactly reproduces the input bicubic Bezier patch. + * The Gregory tangent pairs are duplicated from the four Bezier interior control points, + * causing the parameter-dependent Gregory blends to collapse to fixed Bezier points. + * + * \pre \a bPatch must have order (3, 3) + * \pre \a bPatch must be polynomial, not rational + */ + explicit GregoryPatch(const BezierPatch& bPatch) + { + SLIC_ASSERT(bPatch.getOrder_u() == 3); + SLIC_ASSERT(bPatch.getOrder_v() == 3); + SLIC_ASSERT(!bPatch.isRational()); + + getCorner(0) = bPatch(0, 0); + getCorner(1) = bPatch(3, 0); + getCorner(2) = bPatch(3, 3); + getCorner(3) = bPatch(0, 3); + + getBoundaryPoint(0, 1) = bPatch(1, 0); + getBoundaryPoint(0, 2) = bPatch(2, 0); + getBoundaryPoint(1, 1) = bPatch(3, 1); + getBoundaryPoint(1, 2) = bPatch(3, 2); + getBoundaryPoint(2, 1) = bPatch(2, 3); + getBoundaryPoint(2, 2) = bPatch(1, 3); + getBoundaryPoint(3, 1) = bPatch(0, 2); + getBoundaryPoint(3, 2) = bPatch(0, 1); + + getTangent(0, 0) = bPatch(1, 1); + getTangent(0, 1) = bPatch(2, 1); + getTangent(1, 0) = bPatch(2, 1); + getTangent(1, 1) = bPatch(2, 2); + getTangent(2, 0) = bPatch(2, 2); + getTangent(2, 1) = bPatch(1, 2); + getTangent(3, 0) = bPatch(1, 2); + getTangent(3, 1) = bPatch(1, 1); + } + + /*! + * \brief Constructor from corner nodes and corner normal vectors + * + * \param [in] nodePositions ArrayView of the four corner positions + * \param [in] nodeVectors ArrayView of the four corner normal vectors + * + * Deterministically compute cubic boundary control points and Gregory tangent points + * using local corner information. + * + * \pre \a nodePositions and \a nodeVectors must each contain exactly 4 entries + */ GregoryPatch(ArrayView nodePositions, ArrayView nodeVectors) { // Store the position and orthogonal unit vector at each corner @@ -110,8 +227,6 @@ class GregoryPatch c2[i] = (dx - dx.dot(v[ip1]) * v[ip1]) / 3; a3[i] = VectorType::cross_product(v[ip1], dx).unitVector(); - // Other stuff for adjusted normals - // Use Chiyokura algorithm to define the interior control points const PointType x1(getCorner(i).array() + c0[i].array()); const PointType x2(getCorner(ip1).array() - c2[i].array()); @@ -136,30 +251,85 @@ class GregoryPatch } } + ///@} + + /// \brief Returns the \a i-th corner point, oriented ccw PointType& getCorner(int i) { return m_controlPoints[i]; } + + /// \brief Returns the \a i-th corner point, oriented ccw const PointType& getCorner(int i) const { return m_controlPoints[i]; } + /*! + * \brief Returns a Gregory tangent point for an edge + * + * \param [in] e Edge index, oriented ccw + * \param [in] t Tangent point index (either 0 or 1) + */ PointType& getTangent(int e, int t) { return m_controlPoints[12 + 2 * e + t]; } + + /*! + * \brief Returns a Gregory tangent point for an edge + * + * \param [in] e Edge index, oriented ccw + * \param [in] t Tangent point index (either 0 or 1) + */ const PointType& getTangent(int e, int t) const { return m_controlPoints[12 + 2 * e + t]; } + /*! + * \brief Returns the two Gregory tangent points adjacent to a corner + * + * \param [in] i Corner index, oriented ccw + * \param [out] v0 Tangent point from the preceding edge + * \param [out] v1 Tangent point from the following edge + */ void getTangentsByCorner(int i, PointType& v0, PointType& v1) const { v0 = getTangent((i + 3) % 4, 1); v1 = getTangent(i, 0); } + /*! + * \brief Returns a control point on a boundary edge + * + * \param [in] e Edge index in `[0, 3]` + * \param [in] k Boundary point index in `[0, 3]` + * + * Values `k=0` and `k=3` are the edge's corner points; values `k=1` and `k=2` + * are the cubic boundary control points. + */ PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_edge_index_map[e][k]]; } + + /*! + * \brief Returns a control point on a boundary edge + * + * \param [in] e Edge index in `[0, 3]` + * \param [in] k Boundary point index in `[0, 3]` + * + * Values `k=0` and `k=3` are the edge's corner points; values `k=1` and `k=2` + * are the cubic boundary control points. + */ const PointType& getBoundaryPoint(int e, int k) const { return m_controlPoints[s_edge_index_map[e][k]]; } - // Evaluate the patch by constructing the equivalent Bezier patch with interior control nodes - // defined in terms of the tangent vectors and the evaluation parameters - // - // A cubic Gregory patch can be evaluated by using the 8 internal, blending nodes to construct - // a set of control points which, when used with the Gregory patch edge nodes, produces - // a polynomial bicubic Bezier patch that is equivalent at the requested parameters + /*! + * \brief Returns a reference to the patch's control points + */ + CoordsVec& getControlPoints() { return m_controlPoints; } + + /// \brief Returns a reference to the patch's control points + const CoordsVec& getControlPoints() const { return m_controlPoints; } + + /*! + * \brief Evaluates the Gregory patch at the given parameter values + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + * + * A cubic Gregory patch is evaluated by constructing the equivalent bicubic Bezier patch + * whose interior control points are blended from the Gregory tangent points at (\a u, \a v). + */ PointType evaluate(T u, T v) const { const auto intermediate = setup_intermediate_bezier(u, v, 0); @@ -182,8 +352,8 @@ class GregoryPatch // Chain rule correction due to (u,v)-dependent interior control points axom::StaticArray Bu, dBu, Bv, dBv; - evaluateCubicBernstein(u, Bu, dBu); - evaluateCubicBernstein(v, Bv, dBv); + evaluate_cubic_bernstein(u, Bu, dBu); + evaluate_cubic_bernstein(v, Bv, dBv); const T w11 = Bu[1] * Bv[1]; const T w21 = Bu[2] * Bv[1]; @@ -222,8 +392,8 @@ class GregoryPatch // Chain rule correction due to (u,v)-dependent interior control points axom::StaticArray Bu, dBu, Bv, dBv; - evaluateCubicBernstein(u, Bu, dBu); - evaluateCubicBernstein(v, Bv, dBv); + evaluate_cubic_bernstein(u, Bu, dBu); + evaluate_cubic_bernstein(v, Bv, dBv); const T w11 = Bu[1] * Bv[1]; const T w21 = Bu[2] * Bv[1]; @@ -267,6 +437,12 @@ class GregoryPatch w12 * intermediate.Q_uv[0][1] + w22 * intermediate.Q_uv[1][1]); } + /*! + * \brief Evaluates the first derivative in the u direction + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + */ VectorType du(T u, T v) const { PointType eval; @@ -275,6 +451,12 @@ class GregoryPatch return Du; } + /*! + * \brief Evaluates the first derivative in the v direction + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + */ VectorType dv(T u, T v) const { PointType eval; @@ -283,6 +465,12 @@ class GregoryPatch return Dv; } + /*! + * \brief Evaluates the second derivative in the u direction + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + */ VectorType dudu(T u, T v) const { PointType eval; @@ -291,6 +479,12 @@ class GregoryPatch return DuDu; } + /*! + * \brief Evaluates the second derivative in the v direction + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + */ VectorType dvdv(T u, T v) const { PointType eval; @@ -299,6 +493,12 @@ class GregoryPatch return DvDv; } + /*! + * \brief Evaluates the mixed second derivative + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + */ VectorType dudv(T u, T v) const { PointType eval; @@ -319,6 +519,12 @@ class GregoryPatch return OrientedBoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); } + /*! + * \brief Simple formatted print of a Gregory Patch instance + * + * \param os The output stream to write to + * \return A reference to the modified ostream + */ void print(std::ostream& os) const { os << "GregoryPatch("; @@ -334,9 +540,20 @@ class GregoryPatch } private: + /*! + * \brief Stores the temporary Bezier patch and blended interior point derivatives + * + * The Gregory patch evaluation converts the control net to a bicubic Bezier patch at + * a specific parameter value. The four interior Bezier points, `Q`, depend on the + * evaluation parameters, so derivative evaluation also requires their first and second + * derivatives. + */ struct IntermediateBlendingDerivatives { + /// \brief Equivalent bicubic Bezier patch for the requested parameter value BezierPatch bpatch; + + /// \brief Blended interior Bezier control points and derivatives PointType Q[2][2]; VectorType Q_u[2][2]; VectorType Q_v[2][2]; @@ -345,10 +562,19 @@ class GregoryPatch VectorType Q_uv[2][2]; }; + /*! + * \brief Constructs the equivalent Bezier patch and parameter-dependent interior data + * + * \param [in] u Parameter value on the first axis + * \param [in] v Parameter value on the second axis + * \param [in] derivative_order Highest derivative order to compute, in `[0, 2]` + * + * The returned bicubic Bezier patch has the Gregory boundary control points copied + * directly and the four interior control points blended from the Gregory tangent points. + */ IntermediateBlendingDerivatives setup_intermediate_bezier(T u, T v, int derivative_order) const { IntermediateBlendingDerivatives out; - out.bpatch = get_bezier_boundary(); const T um = T(1) - u; @@ -426,11 +652,17 @@ class GregoryPatch } } - setBezierInterior(out.bpatch, out.Q); + set_bezier_interior(out.bpatch, out.Q); return out; } - static void setBezierInterior(BezierPatch& bpatch, const PointType Q[2][2]) + /*! + * \brief Assigns the four interior control points of a bicubic Bezier patch + * + * \param [in,out] bpatch The bicubic Bezier patch to update + * \param [in] Q The four interior control points, indexed by interior u and v position + */ + static void set_bezier_interior(BezierPatch& bpatch, const PointType Q[2][2]) { bpatch(1, 1) = Q[0][0]; bpatch(2, 1) = Q[1][0]; @@ -438,8 +670,12 @@ class GregoryPatch bpatch(2, 2) = Q[1][1]; } - // Copies over the boundary points to a BezierPatch object, - // leaving the 4 interior control points uninitialized + /*! + * \brief Copies the Gregory boundary into a bicubic Bezier patch + * + * The returned patch has its 12 exterior control points initialized from the Gregory + * patch boundary. The four interior control points are intentionally left uninitialized. + */ BezierPatch get_bezier_boundary() const { BezierPatch bpatch(3, 3); @@ -463,7 +699,14 @@ class GregoryPatch return bpatch; } - static void evaluateCubicBernstein(T t, axom::StaticArray& B, axom::StaticArray& dB) + /*! + * \brief Evaluates cubic Bernstein basis functions and their first derivatives + * + * \param [in] t Parameter value + * \param [out] B Cubic Bernstein basis values at \a t + * \param [out] dB First derivative values of the cubic Bernstein basis at \a t + */ + static void evaluate_cubic_bernstein(T t, axom::StaticArray& B, axom::StaticArray& dB) { const T tm = T(1) - t; const T tm2 = tm * tm; @@ -482,7 +725,12 @@ class GregoryPatch CoordsVec m_controlPoints; - // Map of boundary curve control points into internal storage + /*! + * \brief Maps boundary curve control point indices to control net storage indices + * + * The first index selects a directed edge. The second index selects one of the four + * cubic boundary control points on that edge. + */ static constexpr int s_edge_index_map[4][4] = {{/*V0*/ 0, /*E01*/ 4, /*E02*/ 5, /*V1*/ 1}, {/*V1*/ 1, /*E11*/ 6, /*E12*/ 7, /*V2*/ 2}, {/*V2*/ 2, /*E21*/ 8, /*E22*/ 9, /*V3*/ 3}, diff --git a/src/axom/primal/geometry/GregoryTriangle.hpp b/src/axom/primal/geometry/GregoryTriangle.hpp index 7a56e95624..baee6dffff 100644 --- a/src/axom/primal/geometry/GregoryTriangle.hpp +++ b/src/axom/primal/geometry/GregoryTriangle.hpp @@ -103,17 +103,18 @@ class GregoryTriangle axom::StackArray, 3> cubic_deriv_cp; for(int k = 0; k < 3; ++k) // Loop over edges { - const int kp1 = (k + 1) % 3; - const VectorType dx(nodePositions[k], nodePositions[kp1]); + const int start = (k + 1) % 3; + const int end = (k + 2) % 3; + const VectorType dx(nodePositions[start], nodePositions[end]); - const VectorType c0 = (dx - dx.dot(v[k]) * v[k]) / 3.0; - const VectorType c2 = (dx - dx.dot(v[kp1]) * v[kp1]) / 3.0; + const VectorType c0 = (dx - dx.dot(v[start]) * v[start]) / 3.0; + const VectorType c2 = (dx - dx.dot(v[end]) * v[end]) / 3.0; // Define the cubic Bezier which represents the boundary of the curve - BezierType cubic(axom::Array {nodePositions[k], - PointType {nodePositions[k].array() + c0.array()}, - PointType {nodePositions[kp1].array() - c2.array()}, - nodePositions[kp1]}, + BezierType cubic(axom::Array {nodePositions[start], + PointType {nodePositions[start].array() + c0.array()}, + PointType {nodePositions[end].array() - c2.array()}, + nodePositions[end]}, 3); // Store the control points of the derivative of this cubic for later @@ -130,28 +131,30 @@ class GregoryTriangle // getBoundaryPoint(i, 4) will be set for the next edge } - // Define the interior control nodse at each vertex + // Define the interior control nodes at each vertex for(int k = 0; k < 3; ++k) { - const int km1 = (k + 2) % 3; - const int kp1 = (k + 1) % 3; + const int start = (k + 1) % 3; + const int end = (k + 2) % 3; + const int prev_edge = (k + 2) % 3; + const int next_edge = (k + 1) % 3; // Get control points at and around the edge's starting vertex - const auto& p0 = getBoundaryPoint(km1, 3); // Vertex - 1 - const auto& q0 = getBoundaryPoint(k, 0); // Vertex - const auto& q1 = getBoundaryPoint(k, 1); // Vertex + 1 + const auto& p0 = getBoundaryPoint(prev_edge, 3); // Previous edge + const auto& q0 = getBoundaryPoint(k, 0); // Edge start + const auto& q1 = getBoundaryPoint(k, 1); // Current edge const auto deriv0 = 0.5 * VectorType(p0, q0) + 0.5 * VectorType(p0, q1); // Get control points at and around the edge's ending vertex - const auto& q3 = getBoundaryPoint(k, 3); // Vertex - 1 - const auto& q4 = getBoundaryPoint(k, 4); // Vertex - const auto& p3 = getBoundaryPoint(kp1, 1); // Vertex + 1 + const auto& q3 = getBoundaryPoint(k, 3); // Current edge + const auto& q4 = getBoundaryPoint(k, 4); // Edge end + const auto& p3 = getBoundaryPoint(next_edge, 1); // Next edge const auto deriv1 = 0.5 * VectorType(p3, q3) + 0.5 * VectorType(p3, q4); // Compute a boundary cross derivative that varies across the edge - const VectorType dx(getCorner(k), getCorner(kp1)); - const VectorType a0 = VectorType::cross_product(nodeVectors[k], dx).unitVector(); - const VectorType a3 = VectorType::cross_product(nodeVectors[kp1], dx).unitVector(); + const VectorType dx(getCorner(start), getCorner(end)); + const VectorType a0 = VectorType::cross_product(nodeVectors[start], dx).unitVector(); + const VectorType a3 = VectorType::cross_product(nodeVectors[end], dx).unitVector(); // Elevate the linear cross derivative a(t) = (1-t)a0 + t*a3 to quadratic const axom::StackArray aHat = {a0, 0.5 * (a0 + a3), a3}; @@ -188,8 +191,8 @@ class GregoryTriangle void getTangentsByCorner(int i, PointType& v0, PointType& v1) const { - v0 = getTangent((i + 2) % 3, 1); - v1 = getTangent(i, 0); + v0 = getTangent((i + 1) % 3, 1); + v1 = getTangent((i + 2) % 3, 0); } PointType& getBoundaryPoint(int e, int k) { return m_controlPoints[s_edge_index_map[e][k]]; } @@ -207,13 +210,13 @@ class GregoryTriangle } /*! - * \brief Evaluates all first derivatives of the Gregory patch at (\a u, \a v) + * \brief Evaluates all first derivatives of the Gregory triangle at (\a u0, \a v0) * - * \param [in] u Parameter value at which to evaluate on the first axis - * \param [in] v Parameter value at which to evaluate on the second axis - * \param [out] eval The point value of the Gregory patch at (u, v) - * \param [out] Du The vector value of S_u(u, v) - * \param [out] Dv The vector value of S_v(u, v) + * \param [in] u0 Parameter value at which to evaluate on the first axis + * \param [in] v0 Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory triangle at (u0, v0) + * \param [out] Du The vector value of S_u(u0, v0) + * \param [out] Dv The vector value of S_v(u0, v0) */ void evaluateFirstDerivatives(T u0, T v0, PointType& eval, VectorType& Du, VectorType& Dv) const { @@ -221,33 +224,72 @@ class GregoryTriangle intermediate.btri.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); // Chain rule correction due to (u0,v0)-dependent interior control points. - // For a quartic Bezier triangle, the basis weight of control point (i,j) - // is: 4!/(i! j! k!) * u^i v^j w^k with k = 4-i-j and barycentric weights - // (lambda0, lambda1, lambda2) = (1-u0-v0, v0, u0) corresponding to vertices - // (0,0), (0,order), (order,0), respectively. - const T u = u0; + // This follows BezierTriangle's barycentric convention: + // {u, v, w} = {1 - u0 - v0, v0, u0} + const T u = T(1) - u0 - v0; const T v = v0; - const T w = T(1) - u0 - v0; + const T w = u0; - const T w11 = T(12) * u * v * w * w; // (i,j)=(1,1), k=2 - const T w21 = T(12) * u * u * v * w; // (i,j)=(2,1), k=1 - const T w12 = T(12) * u * v * v * w; // (i,j)=(1,2), k=1 + axom::StaticArray B, B_u0, B_v0; + evaluate_quartic_interior_basis(u, v, w, B, B_u0, B_v0); - Du += w11 * intermediate.Q_u[0] + w21 * intermediate.Q_u[1] + w12 * intermediate.Q_u[2]; - Dv += w11 * intermediate.Q_v[0] + w21 * intermediate.Q_v[1] + w12 * intermediate.Q_v[2]; + Du += B[0] * intermediate.Q_u[0] + B[1] * intermediate.Q_u[1] + B[2] * intermediate.Q_u[2]; + Dv += B[0] * intermediate.Q_v[0] + B[1] * intermediate.Q_v[1] + B[2] * intermediate.Q_v[2]; } - // Not yet completed - // void evaluateSecondDerivatives(T u, - // T v, - // PointType& eval, - // VectorType& Du, - // VectorType& Dv, - // VectorType& DuDu, - // VectorType& DvDv, - // VectorType& DuDv) const - // { - // } + /*! + * \brief Evaluates all second derivatives of the Gregory triangle at (\a u0, \a v0) + * + * \param [in] u0 Parameter value at which to evaluate on the first axis + * \param [in] v0 Parameter value at which to evaluate on the second axis + * \param [out] eval The point value of the Gregory triangle at (u0, v0) + * \param [out] Du The vector value of S_u(u0, v0) + * \param [out] Dv The vector value of S_v(u0, v0) + * \param [out] DuDu The vector value of S_uu(u0, v0) + * \param [out] DvDv The vector value of S_vv(u0, v0) + * \param [out] DuDv The vector value of S_uv(u0, v0) == S_vu(u0, v0) + */ + void evaluateSecondDerivatives(T u0, + T v0, + PointType& eval, + VectorType& Du, + VectorType& Dv, + VectorType& DuDu, + VectorType& DvDv, + VectorType& DuDv) const + { + const auto intermediate = setup_intermediate_bezier(u0, v0, 2); + intermediate.btri.evaluateSecondDerivatives(u0, v0, eval, Du, Dv, DuDu, DvDv, DuDv); + + const T u = T(1) - u0 - v0; + const T v = v0; + const T w = u0; + + axom::StaticArray B, B_u0, B_v0; + evaluate_quartic_interior_basis(u, v, w, B, B_u0, B_v0); + + // First derivative corrections + Du += B[0] * intermediate.Q_u[0] + B[1] * intermediate.Q_u[1] + B[2] * intermediate.Q_u[2]; + Dv += B[0] * intermediate.Q_v[0] + B[1] * intermediate.Q_v[1] + B[2] * intermediate.Q_v[2]; + + // Second derivative corrections + DuDu += T(2) * + (B_u0[0] * intermediate.Q_u[0] + B_u0[1] * intermediate.Q_u[1] + + B_u0[2] * intermediate.Q_u[2]) + + (B[0] * intermediate.Q_uu[0] + B[1] * intermediate.Q_uu[1] + B[2] * intermediate.Q_uu[2]); + + DvDv += T(2) * + (B_v0[0] * intermediate.Q_v[0] + B_v0[1] * intermediate.Q_v[1] + + B_v0[2] * intermediate.Q_v[2]) + + (B[0] * intermediate.Q_vv[0] + B[1] * intermediate.Q_vv[1] + B[2] * intermediate.Q_vv[2]); + + DuDv += + (B_u0[0] * intermediate.Q_v[0] + B_u0[1] * intermediate.Q_v[1] + + B_u0[2] * intermediate.Q_v[2]) + + (B_v0[0] * intermediate.Q_u[0] + B_v0[1] * intermediate.Q_u[1] + + B_v0[2] * intermediate.Q_u[2]) + + (B[0] * intermediate.Q_uv[0] + B[1] * intermediate.Q_uv[1] + B[2] * intermediate.Q_uv[2]); + } VectorType du(T u, T v) const { @@ -265,30 +307,29 @@ class GregoryTriangle return Dv; } - // Not yet finished - // VectorType dudu(T u, T v) const - // { - // PointType eval; - // VectorType Du, Dv, DuDu, DvDv, DuDv; - // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); - // return DuDu; - // } - - // VectorType dvdv(T u, T v) const - // { - // PointType eval; - // VectorType Du, Dv, DuDu, DvDv, DuDv; - // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); - // return DvDv; - // } - - // VectorType dudv(T u, T v) const - // { - // PointType eval; - // VectorType Du, Dv, DuDu, DvDv, DuDv; - // evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); - // return DuDv; - // } + VectorType dudu(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DuDu; + } + + VectorType dvdv(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DvDv; + } + + VectorType dudv(T u, T v) const + { + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + evaluateSecondDerivatives(u, v, eval, Du, Dv, DuDu, DvDv, DuDv); + return DuDv; + } /// \brief Returns an axis-aligned bounding box containing the patch BoundingBoxType boundingBox() const @@ -323,6 +364,9 @@ class GregoryTriangle PointType Q[3]; VectorType Q_u[3]; VectorType Q_v[3]; + VectorType Q_uu[3]; + VectorType Q_vv[3]; + VectorType Q_uv[3]; }; IntermediateBlendingDerivatives setup_intermediate_bezier(T u0, T v0, int derivative_order) const @@ -331,28 +375,28 @@ class GregoryTriangle out.btri = get_bezier_boundary(); // Parameter convention matches BezierTriangle::evaluate(): - // barycentric weights (lambda0, lambda1, lambda2) are (1-u0-v0, v0, u0) - const T u = u0; + // barycentric coordinates {u, v, w} are {1-u0-v0, v0, u0} + const T u = T(1) - u0 - v0; const T v = v0; - const T w = T(1) - u0 - v0; - - auto blend = [&](const PointType& A, - const PointType& B, - T wa, - T wb, - T wa_u0, - T wb_u0, - T wa_v0, - T wb_v0, + const T w = u0; + + // clang-format off + auto blend = [&](const PointType& A, const PointType& B, // Internal Gregory points + T wa, T wb, // Barycentric coordinates for eval + T wa_u0, T wb_u0, T wa_v0, T wb_v0, // Derivatives of Barycentric coords PointType& Q, - VectorType& Q_u0, - VectorType& Q_v0) { + VectorType& Q_u0, VectorType& Q_v0, + VectorType& Q_u0u0, VectorType& Q_v0v0, VectorType& Q_u0v0) { + // clang-format on const T denom = wa + wb; if(axom::utilities::isNearlyEqual(denom, T(0))) { Q = A; Q_u0 = VectorType(T(0)); Q_v0 = VectorType(T(0)); + Q_u0u0 = VectorType(T(0)); + Q_v0v0 = VectorType(T(0)); + Q_u0v0 = VectorType(T(0)); return; } @@ -372,104 +416,134 @@ class GregoryTriangle Q_u0 = VectorType(T(0)); Q_v0 = VectorType(T(0)); } + + if(derivative_order >= 2) + { + const T denom_u0 = wa_u0 + wb_u0; + const T denom_v0 = wa_v0 + wb_v0; + Q_u0u0 = (-T(2) * denom_u0 / denom) * Q_u0; + Q_v0v0 = (-T(2) * denom_v0 / denom) * Q_v0; + Q_u0v0 = (-(denom_u0 * Q_v0 + denom_v0 * Q_u0)) / denom; + } + else + { + Q_u0u0 = VectorType(T(0)); + Q_v0v0 = VectorType(T(0)); + Q_u0v0 = VectorType(T(0)); + } }; - // These are the three (u0,v0)-dependent interior control points for the - // equivalent quartic Bezier triangle. Each is a blend of the two Gregory - // points adjacent to the corresponding vertex. - // - // By convention, edge k connects vertex k -> vertex (k+1)%3. - // - getTangent(k,0) is the Gregory point near the starting vertex k - // - getTangent(k,1) is the Gregory point near the ending vertex (k+1)%3 - // - // Control net slots: - // Q[0] -> btri(1,1) - // Q[1] -> btri(2,1) - // Q[2] -> btri(1,2) - blend(getTangent(2, 1), - getTangent(0, 0), - u, - v, - T(1), - T(0), - T(0), - T(1), + // Get the three (u0, v0)-dependent interior control points for the equivalent + // quartic Bezier triangle. Each is blended from the two Gregory points adjacent + // to the corresponding vertex, + + // clang-format off + blend(getTangent(1, 1), getTangent(2, 0), + w, v, + T(1), T(0), T(0), T(1), out.Q[0], - out.Q_u[0], - out.Q_v[0]); - - blend(getTangent(1, 1), - getTangent(2, 0), - v, - w, - T(0), - T(-1), - T(1), - T(-1), + out.Q_u[0], out.Q_v[0], + out.Q_uu[0], out.Q_vv[0], out.Q_uv[0]); + + blend(getTangent(0, 1), getTangent(1, 0), + v, u, + T(0), T(-1), T(1), T(-1), out.Q[1], - out.Q_u[1], - out.Q_v[1]); - - blend(getTangent(0, 1), - getTangent(1, 0), - w, - u, - T(-1), - T(1), - T(-1), - T(0), + out.Q_u[1], out.Q_v[1], + out.Q_uu[1], out.Q_vv[1], out.Q_uv[1]); + + blend(getTangent(2, 1), getTangent(0, 0), + u, w, + T(-1), T(1), T(-1), T(0), out.Q[2], - out.Q_u[2], - out.Q_v[2]); + out.Q_u[2], out.Q_v[2], + out.Q_uu[2], out.Q_vv[2], out.Q_uv[2]); + // clang-format on - setBezierInterior(out.btri, out.Q); + set_bezier_interior(out.btri, out.Q); return out; } - static void setBezierInterior(BezierTriangle& btri, const PointType Q[3]) + /*! + * \brief Assigns the three interior control points of a biquartic Bezier triangle + * + * \param [in,out] btri The biquartic Bezier triangle to update + * \param [in] Q The 3 interior control points + */ + static void set_bezier_interior(BezierTriangle& btri, const PointType Q[3]) { btri(1, 1) = Q[0]; btri(2, 1) = Q[1]; btri(1, 2) = Q[2]; } - // Copies over the boundary points to a BezierPatch object, - // leaving the 4 interior control points uninitialized + /*! + * \brief Evaluates triangular Bernstein basis functions and their first derivatives + * + * \param [in] u First standard barycentric coordinate, equal to `1-u0-v0` + * \param [in] v Second standard barycentric coordinate, equal to `v0` + * \param [in] w Third standard barycentric coordinate, equal to `u0` + * \param [out] B Basis values for the three interior control points + * \param [out] B_u0 First derivatives of \a B with respect to `u0` + * \param [out] B_v0 First derivatives of \a B with respect to `v0` + */ + static void evaluate_quartic_interior_basis(T u, + T v, + T w, + axom::StaticArray& B, + axom::StaticArray& B_u0, + axom::StaticArray& B_v0) + { + B[0] = T(12) * w * v * u * u; // (i,j)=(1,1), k=2 + B[1] = T(12) * w * w * v * u; // (i,j)=(2,1), k=1 + B[2] = T(12) * w * v * v * u; // (i,j)=(1,2), k=1 + + B_u0[0] = T(12) * v * u * (u - T(2) * w); + B_u0[1] = T(12) * w * v * (T(2) * u - w); + B_u0[2] = T(12) * v * v * (u - w); + + B_v0[0] = T(12) * w * u * (u - T(2) * v); + B_v0[1] = T(12) * w * w * (u - v); + B_v0[2] = T(12) * w * v * (T(2) * u - v); + } + + // Copies over the boundary points to a BezierTriangle object, + // leaving the 3 interior control points uninitialized BezierTriangle get_bezier_boundary() const { BezierTriangle btri(4); // Edge 0 - btri(0, 0) = getBoundaryPoint(0, 0); - btri(0, 1) = getBoundaryPoint(0, 1); - btri(0, 2) = getBoundaryPoint(0, 2); - btri(0, 3) = getBoundaryPoint(0, 3); - btri(0, 4) = getBoundaryPoint(0, 4); + btri(0, 4) = getBoundaryPoint(0, 0); + btri(1, 3) = getBoundaryPoint(0, 1); + btri(2, 2) = getBoundaryPoint(0, 2); + btri(3, 1) = getBoundaryPoint(0, 3); + btri(4, 0) = getBoundaryPoint(0, 4); // Edge 1 - // btri(0, 4) = getBoundaryPOint(1, 0); - btri(1, 3) = getBoundaryPoint(1, 1); - btri(2, 2) = getBoundaryPoint(1, 2); - btri(3, 1) = getBoundaryPoint(1, 3); - btri(4, 0) = getBoundaryPoint(1, 4); + // btri(4, 0) = getBoundaryPoint(1, 0); + btri(3, 0) = getBoundaryPoint(1, 1); + btri(2, 0) = getBoundaryPoint(1, 2); + btri(1, 0) = getBoundaryPoint(1, 3); + btri(0, 0) = getBoundaryPoint(1, 4); // Edge 2 - // btri(4, 0) = getBoundaryPoint(2, 0); - btri(3, 0) = getBoundaryPoint(2, 1); - btri(2, 0) = getBoundaryPoint(2, 2); - btri(1, 0) = getBoundaryPoint(2, 3); - // btri(0, 0) = getBoundaryPoint(2, 4); + // btri(0, 0) = getBoundaryPoint(2, 0); + btri(0, 1) = getBoundaryPoint(2, 1); + btri(0, 2) = getBoundaryPoint(2, 2); + btri(0, 3) = getBoundaryPoint(2, 3); + // btri(0, 4) = getBoundaryPoint(2, 4); return btri; } CoordsVec m_controlPoints; - // Map of boundary curve control points into internal storage + // Map of BezierTriangle-style boundary curve control points into internal storage static constexpr int s_edge_index_map[3][5] = { - {/*V0*/ 0, /*E01*/ 3, /*E02*/ 4, /*E03*/ 5, /*V1*/ 1}, - {/*V0*/ 1, /*E01*/ 6, /*E02*/ 7, /*E03*/ 8, /*V1*/ 2}, - {/*V0*/ 2, /*E01*/ 9, /*E02*/ 10, /*E03*/ 11, /*V1*/ 0}}; + {/*V1*/ 1, /*E01*/ 6, /*E02*/ 7, /*E03*/ 8, /*V2*/ 2}, + {/*V2*/ 2, /*E11*/ 9, /*E12*/ 10, /*E13*/ 11, /*V0*/ 0}, + {/*V0*/ 0, /*E21*/ 3, /*E22*/ 4, /*E23*/ 5, /*V1*/ 1}}; }; //------------------------------------------------------------------------------ From 029425fccf1f25d16605956c25db346b076385de Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 13 Aug 2026 13:03:15 -0700 Subject: [PATCH 07/16] Move tests to one file --- src/axom/primal/tests/CMakeLists.txt | 1 + .../primal/tests/primal_bezier_triangle.cpp | 111 ---- .../primal/tests/primal_gregory_surfaces.cpp | 516 ++++++++++++++++++ 3 files changed, 517 insertions(+), 111 deletions(-) create mode 100644 src/axom/primal/tests/primal_gregory_surfaces.cpp diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index 99cb2c2a7a..f045ac3808 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -21,6 +21,7 @@ set( primal_tests primal_cone.cpp primal_coord_transform.cpp primal_curved_polygon.cpp + primal_gregory_surfaces.cpp primal_hexahedron.cpp primal_in_sphere.cpp primal_integral.cpp diff --git a/src/axom/primal/tests/primal_bezier_triangle.cpp b/src/axom/primal/tests/primal_bezier_triangle.cpp index 3878ceb137..f1fc5779f6 100644 --- a/src/axom/primal/tests/primal_bezier_triangle.cpp +++ b/src/axom/primal/tests/primal_bezier_triangle.cpp @@ -14,8 +14,6 @@ #include "axom/slic.hpp" #include "axom/primal/geometry/BezierTriangle.hpp" -#include "axom/primal/geometry/GregoryTriangle.hpp" -#include "axom/primal/geometry/GregoryPatch.hpp" #include @@ -346,115 +344,6 @@ TEST(primal_beziertriangle, finite_difference_first_derivatives) } } -//------------------------------------------------------------------------------ -TEST(primal_beziertriangle, gregorytriangle_finite_difference_first_derivatives) -{ - using CoordType = double; - using PointType = primal::Point; - using VectorType = primal::Vector; - using GTri = primal::GregoryTriangle; - - constexpr CoordType h = 1e-7; - constexpr CoordType tol = 5e-5; - - const std::array corners = { - PointType {0.0, 0.0, 0.0}, - PointType {1.0, 0.0, 0.2}, - PointType {0.0, 1.0, -0.1}, - }; - - const std::array normals = { - VectorType {0.0, 0.0, 1.0}, - VectorType {0.0, 0.0, 1.0}, - VectorType {0.0, 0.0, 1.0}, - }; - - GTri gtri(axom::ArrayView(corners.data(), 3), - axom::ArrayView(normals.data(), 3)); - - const CoordType u0 = 0.21; - const CoordType v0 = 0.27; - - PointType eval; - VectorType Du, Dv; - gtri.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); - - const PointType fp_u = gtri.evaluate(u0 + h, v0); - const PointType fm_u = gtri.evaluate(u0 - h, v0); - const PointType fp_v = gtri.evaluate(u0, v0 + h); - const PointType fm_v = gtri.evaluate(u0, v0 - h); - - const VectorType Du_fd(fp_u, fm_u); - const VectorType Dv_fd(fp_v, fm_v); - - for(int d = 0; d < 3; ++d) - { - EXPECT_NEAR(Du[d], Du_fd[d] / (2 * h), tol); - EXPECT_NEAR(Dv[d], Dv_fd[d] / (2 * h), tol); - } -} - -//------------------------------------------------------------------------------ -TEST(primal_beziertriangle, gregorytriangle_parameter_convention_matches_beziertriangle) -{ - using CoordType = double; - using PointType = primal::Point; - using VectorType = primal::Vector; - using GTri = primal::GregoryTriangle; - - constexpr CoordType eps = 1e-12; - - const std::array corners = { - PointType {0.0, 0.0, 0.0}, - PointType {1.0, 0.0, 0.0}, - PointType {0.0, 1.0, 0.0}, - }; - - const std::array normals = { - VectorType {0.0, 0.0, 1.0}, - VectorType {0.0, 0.0, 1.0}, - VectorType {0.0, 0.0, 1.0}, - }; - - GTri gtri(axom::ArrayView(corners.data(), 3), - axom::ArrayView(normals.data(), 3)); - - // Corner mapping should match BezierTriangle::evaluate convention - EXPECT_TRUE(gtri.evaluate(0.0, 0.0).isNearlyEqual(corners[0], eps)); - EXPECT_TRUE(gtri.evaluate(0.0, 1.0).isNearlyEqual(corners[1], eps)); - EXPECT_TRUE(gtri.evaluate(1.0, 0.0).isNearlyEqual(corners[2], eps)); - - // Edge mapping: verify straight edges for this planar configuration - auto check_near = [&](const PointType& a, const PointType& b) { - EXPECT_NEAR(a[0], b[0], eps); - EXPECT_NEAR(a[1], b[1], eps); - EXPECT_NEAR(a[2], b[2], eps); - }; - - const CoordType tvals[] = {0.0, 0.2, 0.6, 1.0}; - - // edge u0 = 0: from corner0 -> corner1 as v0 goes 0->1 - for(const CoordType t : tvals) - { - const PointType expected = PointType::lerp(corners[0], corners[1], t); - check_near(gtri.evaluate(0.0, t), expected); - } - - // edge v0 = 0: from corner0 -> corner2 as u0 goes 0->1 - for(const CoordType t : tvals) - { - const PointType expected = PointType::lerp(corners[0], corners[2], t); - check_near(gtri.evaluate(t, 0.0), expected); - } - - // edge u0 + v0 = 1: from corner1 -> corner2 as u0 goes 0->1 (v0 goes 1->0) - for(const CoordType t : tvals) - { - const PointType expected = PointType::lerp(corners[1], corners[2], t); - check_near(gtri.evaluate(t, 1.0 - t), expected); - } -} - //------------------------------------------------------------------------------ TEST(primal_beziertriangle, evaluate_quadratic_second_derivatives_constant) { diff --git a/src/axom/primal/tests/primal_gregory_surfaces.cpp b/src/axom/primal/tests/primal_gregory_surfaces.cpp new file mode 100644 index 0000000000..4691a90303 --- /dev/null +++ b/src/axom/primal/tests/primal_gregory_surfaces.cpp @@ -0,0 +1,516 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file primal_gregory_surfaces.cpp + * \brief This file tests primal's Gregory surface functionality + */ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" + +#include "axom/primal/geometry/BezierPatch.hpp" +#include "axom/primal/geometry/BezierTriangle.hpp" +#include "axom/primal/geometry/GregoryPatch.hpp" +#include "axom/primal/geometry/GregoryTriangle.hpp" + +#include +#include + +namespace primal = axom::primal; + +namespace +{ +using CoordType = double; +using GregoryPatchType = primal::GregoryPatch; +using BezierPatchType = primal::BezierPatch; +using GregoryTriangleType = primal::GregoryTriangle; +using BezierTriangleType = primal::BezierTriangle; +using PointType = GregoryPatchType::PointType; +using VectorType = GregoryPatchType::VectorType; +using PatchCoordsVec = GregoryPatchType::CoordsVec; +using TriangleCoordsVec = GregoryTriangleType::CoordsVec; + +BezierPatchType make_sample_bezier_patch() +{ + // clang-format off + PointType bezierControlPoints[16] = { + PointType {0, 0, 0}, PointType {0, 2, 1}, PointType {0, 4, -1}, PointType {0, 6, 0}, + PointType {2, 0, 1}, PointType {2, 2, 3}, PointType {2, 4, 2}, PointType {2, 6, 1}, + PointType {4, 0, 0}, PointType {4, 2, 2}, PointType {4, 4, 3}, PointType {4, 6, 0}, + PointType {6, 0, 0}, PointType {6, 2, -1}, PointType {6, 4, 1}, PointType {6, 6, 0}}; + // clang-format on + + return BezierPatchType(bezierControlPoints, 3, 3); +} + +BezierTriangleType make_sample_bezier_triangle() +{ + BezierTriangleType bTri(4); + + for(int i = 0; i <= bTri.getOrder(); ++i) + { + for(int j = 0; j <= bTri.getOrder() - i; ++j) + { + const auto ii = static_cast(i); + const auto jj = static_cast(j); + bTri(i, j) = PointType {ii + 0.25 * jj, jj - 0.1 * ii, 0.5 * ii * ii - 0.25 * jj + ii * jj}; + } + } + + return bTri; +} + +void check_patch_control_points(const GregoryPatchType& patch, const PatchCoordsVec& expected) +{ + for(int i = 0; i < GregoryPatchType::NPTS; ++i) + { + EXPECT_EQ(patch.getControlPoints()[i], expected[i]); + } +} + +void check_triangle_control_points(const GregoryTriangleType& triangle, const TriangleCoordsVec& expected) +{ + for(int i = 0; i < GregoryTriangleType::NPTS; ++i) + { + EXPECT_EQ(triangle.getControlPoints()[i], expected[i]); + } +} + +} // namespace + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, array_constructors) +{ + SLIC_INFO("Testing Gregory patch point array constructors"); + + auto controlPoints = GregoryPatchType(make_sample_bezier_patch()).getControlPoints(); + const auto constControlPoints = controlPoints; + + { + SCOPED_TRACE("Testing C-array constructor"); + GregoryPatchType patch(controlPoints.data()); + check_patch_control_points(patch, controlPoints); + } + + { + SCOPED_TRACE("Testing const C-array constructor"); + GregoryPatchType patch(constControlPoints.data()); + check_patch_control_points(patch, controlPoints); + } + + { + SCOPED_TRACE("Testing StackArray constructor"); + GregoryPatchType patch(controlPoints); + check_patch_control_points(patch, controlPoints); + } + + { + SCOPED_TRACE("Testing ArrayView constructor"); + axom::ArrayView view(controlPoints.data(), GregoryPatchType::NPTS); + GregoryPatchType patch(view); + check_patch_control_points(patch, controlPoints); + } + + { + SCOPED_TRACE("Testing const ArrayView constructor"); + axom::ArrayView view(constControlPoints.data(), GregoryPatchType::NPTS); + GregoryPatchType patch(view); + check_patch_control_points(patch, controlPoints); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, corner_vector_constructor) +{ + SLIC_INFO("Testing Gregory patch corner and vector constructor"); + + PointType corners[4] = {PointType {0.0, 0.0, 0.0}, + PointType {1.0, 0.0, 0.0}, + PointType {1.0, 1.0, 0.0}, + PointType {0.0, 1.0, 0.0}}; + + VectorType normals[4] = {VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}}; + + GregoryPatchType patch(axom::ArrayView(corners, 4), + axom::ArrayView(normals, 4)); + + for(int i = 0; i < 4; ++i) + { + EXPECT_EQ(patch.getCorner(i), corners[i]); + } + + EXPECT_EQ(patch.evaluate(0.0, 0.0), corners[0]); + EXPECT_EQ(patch.evaluate(1.0, 0.0), corners[1]); + EXPECT_EQ(patch.evaluate(1.0, 1.0), corners[2]); + EXPECT_EQ(patch.evaluate(0.0, 1.0), corners[3]); +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, evaluate) +{ + SLIC_INFO("Testing Gregory patch evaluation"); + + const BezierPatchType bPatch = make_sample_bezier_patch(); + const GregoryPatchType gPatch(bPatch); + + const double uValues[] = {0.0, 0.1, 0.25, 0.3, 0.5, 0.7, 0.75, 0.8, 1.0}; + const double vValues[] = {0.0, 0.1, 0.25, 0.3, 0.5, 0.7, 0.75, 0.8, 1.0}; + + for(double u : uValues) + { + for(double v : vValues) + { + const PointType gregoryPoint = gPatch.evaluate(u, v); + const PointType bezierPoint = bPatch.evaluate(u, v); + for(int dim = 0; dim < 3; ++dim) + { + EXPECT_NEAR(gregoryPoint[dim], bezierPoint[dim], 1e-12); + } + } + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, derivatives) +{ + SLIC_INFO("Testing Gregory patch derivative evaluation"); + + const BezierPatchType bPatch = make_sample_bezier_patch(); + const GregoryPatchType gPatch(bPatch); + + const double u = 0.25; + const double v = 0.75; + + PointType gEval, bEval; + VectorType gDu, gDv, bDu, bDv; + VectorType gDuDu, gDvDv, gDuDv, bDuDu, bDvDv, bDuDv; + + gPatch.evaluateSecondDerivatives(u, v, gEval, gDu, gDv, gDuDu, gDvDv, gDuDv); + bPatch.evaluateSecondDerivatives(u, v, bEval, bDu, bDv, bDuDu, bDvDv, bDuDv); + + for(int dim = 0; dim < 3; ++dim) + { + EXPECT_NEAR(gEval[dim], bEval[dim], 1e-12); + EXPECT_NEAR(gPatch.du(u, v)[dim], bDu[dim], 1e-12); + EXPECT_NEAR(gPatch.dv(u, v)[dim], bDv[dim], 1e-12); + EXPECT_NEAR(gPatch.dudu(u, v)[dim], bDuDu[dim], 1e-12); + EXPECT_NEAR(gPatch.dvdv(u, v)[dim], bDvDv[dim], 1e-12); + EXPECT_NEAR(gPatch.dudv(u, v)[dim], bDuDv[dim], 1e-12); + EXPECT_NEAR(gDu[dim], bDu[dim], 1e-12); + EXPECT_NEAR(gDv[dim], bDv[dim], 1e-12); + EXPECT_NEAR(gDuDu[dim], bDuDu[dim], 1e-12); + EXPECT_NEAR(gDvDv[dim], bDvDv[dim], 1e-12); + EXPECT_NEAR(gDuDv[dim], bDuDv[dim], 1e-12); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, bounding_box) +{ + SLIC_INFO("Testing Gregory patch bounding box"); + + const auto controlPoints = GregoryPatchType(make_sample_bezier_patch()).getControlPoints(); + GregoryPatchType patch(controlPoints); + const auto bbox = patch.boundingBox(); + const auto obb = patch.orientedBoundingBox(); + + for(int i = 0; i < GregoryPatchType::NPTS; ++i) + { + EXPECT_TRUE(bbox.contains(controlPoints[i])); + EXPECT_TRUE(obb.contains(controlPoints[i])); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorypatch, print) +{ + SLIC_INFO("Testing Gregory patch output stream"); + + GregoryPatchType patch(make_sample_bezier_patch()); + std::ostringstream oss; + oss << patch; + + EXPECT_NE(std::string::npos, oss.str().find("GregoryPatch(")); +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, array_constructors) +{ + SLIC_INFO("Testing Gregory triangle point array constructors"); + + auto controlPoints = GregoryTriangleType(make_sample_bezier_triangle()).getControlPoints(); + const auto constControlPoints = controlPoints; + + { + SCOPED_TRACE("Testing C-array constructor"); + GregoryTriangleType triangle(controlPoints.data()); + check_triangle_control_points(triangle, controlPoints); + } + + { + SCOPED_TRACE("Testing const C-array constructor"); + GregoryTriangleType triangle(constControlPoints.data()); + check_triangle_control_points(triangle, controlPoints); + } + + { + SCOPED_TRACE("Testing StackArray constructor"); + GregoryTriangleType triangle(controlPoints); + check_triangle_control_points(triangle, controlPoints); + } + + { + SCOPED_TRACE("Testing ArrayView constructor"); + axom::ArrayView view(controlPoints.data(), GregoryTriangleType::NPTS); + GregoryTriangleType triangle(view); + check_triangle_control_points(triangle, controlPoints); + } + + { + SCOPED_TRACE("Testing const ArrayView constructor"); + axom::ArrayView view(constControlPoints.data(), GregoryTriangleType::NPTS); + GregoryTriangleType triangle(view); + check_triangle_control_points(triangle, controlPoints); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, corner_vector_constructor) +{ + SLIC_INFO("Testing Gregory triangle corner and vector constructor"); + + constexpr CoordType eps = 1e-12; + + const std::array corners = { + PointType {0.0, 0.0, 0.0}, + PointType {1.0, 0.0, 0.0}, + PointType {0.0, 1.0, 0.0}, + }; + + const std::array normals = { + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + }; + + GregoryTriangleType triangle(axom::ArrayView(corners.data(), 3), + axom::ArrayView(normals.data(), 3)); + + EXPECT_TRUE(triangle.evaluate(0.0, 0.0).isNearlyEqual(corners[0], eps)); + EXPECT_TRUE(triangle.evaluate(0.0, 1.0).isNearlyEqual(corners[1], eps)); + EXPECT_TRUE(triangle.evaluate(1.0, 0.0).isNearlyEqual(corners[2], eps)); + + EXPECT_EQ(triangle.getBoundaryPoint(0, 0), corners[1]); + EXPECT_EQ(triangle.getBoundaryPoint(0, 4), corners[2]); + EXPECT_EQ(triangle.getBoundaryPoint(1, 0), corners[2]); + EXPECT_EQ(triangle.getBoundaryPoint(1, 4), corners[0]); + EXPECT_EQ(triangle.getBoundaryPoint(2, 0), corners[0]); + EXPECT_EQ(triangle.getBoundaryPoint(2, 4), corners[1]); + + auto check_near = [&](const PointType& a, const PointType& b) { + EXPECT_NEAR(a[0], b[0], eps); + EXPECT_NEAR(a[1], b[1], eps); + EXPECT_NEAR(a[2], b[2], eps); + }; + + const CoordType tvals[] = {0.0, 0.2, 0.6, 1.0}; + for(const CoordType t : tvals) + { + check_near(triangle.evaluate(0.0, t), PointType::lerp(corners[0], corners[1], t)); + check_near(triangle.evaluate(t, 0.0), PointType::lerp(corners[0], corners[2], t)); + check_near(triangle.evaluate(t, 1.0 - t), PointType::lerp(corners[1], corners[2], t)); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, evaluate) +{ + SLIC_INFO("Testing Gregory triangle evaluation"); + + const BezierTriangleType bTri = make_sample_bezier_triangle(); + const GregoryTriangleType gTri(bTri); + + const CoordType tValues[] = {0.0, 0.1, 0.25, 0.3, 0.5, 0.7, 0.75, 0.8, 1.0}; + for(CoordType u = 0.0; u <= 1.0; u += 0.25) + { + for(const CoordType v : tValues) + { + if(u + v <= 1.0) + { + const PointType gregoryPoint = gTri.evaluate(u, v); + const PointType bezierPoint = bTri.evaluate(u, v); + for(int dim = 0; dim < 3; ++dim) + { + EXPECT_NEAR(gregoryPoint[dim], bezierPoint[dim], 1e-12); + } + } + } + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, derivatives) +{ + SLIC_INFO("Testing Gregory triangle derivative evaluation"); + + const BezierTriangleType bTri = make_sample_bezier_triangle(); + const GregoryTriangleType gTri(bTri); + + const CoordType u = 0.25; + const CoordType v = 0.35; + + PointType gEval, bEval; + VectorType gDu, gDv, bDu, bDv; + VectorType gDuDu, gDvDv, gDuDv, bDuDu, bDvDv, bDuDv; + + gTri.evaluateSecondDerivatives(u, v, gEval, gDu, gDv, gDuDu, gDvDv, gDuDv); + bTri.evaluateSecondDerivatives(u, v, bEval, bDu, bDv, bDuDu, bDvDv, bDuDv); + + for(int dim = 0; dim < 3; ++dim) + { + EXPECT_NEAR(gEval[dim], bEval[dim], 1e-12); + EXPECT_NEAR(gTri.du(u, v)[dim], bDu[dim], 1e-12); + EXPECT_NEAR(gTri.dv(u, v)[dim], bDv[dim], 1e-12); + EXPECT_NEAR(gTri.dudu(u, v)[dim], bDuDu[dim], 1e-12); + EXPECT_NEAR(gTri.dvdv(u, v)[dim], bDvDv[dim], 1e-12); + EXPECT_NEAR(gTri.dudv(u, v)[dim], bDuDv[dim], 1e-12); + EXPECT_NEAR(gDu[dim], bDu[dim], 1e-12); + EXPECT_NEAR(gDv[dim], bDv[dim], 1e-12); + EXPECT_NEAR(gDuDu[dim], bDuDu[dim], 1e-12); + EXPECT_NEAR(gDvDv[dim], bDvDv[dim], 1e-12); + EXPECT_NEAR(gDuDv[dim], bDuDv[dim], 1e-12); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, finite_difference_first_derivatives) +{ + SLIC_INFO("Testing Gregory triangle first derivatives with finite differences"); + + constexpr CoordType h = 1e-7; + constexpr CoordType tol = 5e-5; + + const std::array corners = { + PointType {0.0, 0.0, 0.0}, + PointType {1.0, 0.0, 0.2}, + PointType {0.0, 1.0, -0.1}, + }; + + const std::array normals = { + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + VectorType {0.0, 0.0, 1.0}, + }; + + GregoryTriangleType triangle(axom::ArrayView(corners.data(), 3), + axom::ArrayView(normals.data(), 3)); + + const CoordType u0 = 0.21; + const CoordType v0 = 0.27; + + PointType eval; + VectorType Du, Dv; + triangle.evaluateFirstDerivatives(u0, v0, eval, Du, Dv); + + const PointType fp_u = triangle.evaluate(u0 + h, v0); + const PointType fm_u = triangle.evaluate(u0 - h, v0); + const PointType fp_v = triangle.evaluate(u0, v0 + h); + const PointType fm_v = triangle.evaluate(u0, v0 - h); + + const VectorType Du_fd(fp_u, fm_u); + const VectorType Dv_fd(fp_v, fm_v); + + for(int d = 0; d < 3; ++d) + { + EXPECT_NEAR(Du[d], Du_fd[d] / (2 * h), tol); + EXPECT_NEAR(Dv[d], Dv_fd[d] / (2 * h), tol); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, finite_difference_second_derivatives) +{ + SLIC_INFO("Testing Gregory triangle second derivatives with finite differences"); + + constexpr CoordType h = 1e-5; + constexpr CoordType tol = 1e-4; + + const std::array corners = { + PointType {0.0, 0.0, 0.0}, + PointType {1.2, 0.1, 0.3}, + PointType {-0.1, 1.1, -0.2}, + }; + + const std::array normals = { + VectorType {0.0, 0.0, 1.0}, + VectorType {0.2, 0.1, 1.0}, + VectorType {-0.1, 0.2, 1.0}, + }; + + GregoryTriangleType triangle(axom::ArrayView(corners.data(), 3), + axom::ArrayView(normals.data(), 3)); + + const CoordType u0 = 0.23; + const CoordType v0 = 0.31; + + PointType eval; + VectorType Du, Dv, DuDu, DvDv, DuDv; + triangle.evaluateSecondDerivatives(u0, v0, eval, Du, Dv, DuDu, DvDv, DuDv); + + PointType eval_pu, eval_mu, eval_pv, eval_mv; + VectorType Du_pu, Dv_pu, Du_mu, Dv_mu; + VectorType Du_pv, Dv_pv, Du_mv, Dv_mv; + triangle.evaluateFirstDerivatives(u0 + h, v0, eval_pu, Du_pu, Dv_pu); + triangle.evaluateFirstDerivatives(u0 - h, v0, eval_mu, Du_mu, Dv_mu); + triangle.evaluateFirstDerivatives(u0, v0 + h, eval_pv, Du_pv, Dv_pv); + triangle.evaluateFirstDerivatives(u0, v0 - h, eval_mv, Du_mv, Dv_mv); + + for(int d = 0; d < 3; ++d) + { + EXPECT_NEAR(triangle.dudu(u0, v0)[d], DuDu[d], tol); + EXPECT_NEAR(triangle.dvdv(u0, v0)[d], DvDv[d], tol); + EXPECT_NEAR(triangle.dudv(u0, v0)[d], DuDv[d], tol); + + EXPECT_NEAR(DuDu[d], (Du_pu[d] - Du_mu[d]) / (2 * h), tol); + EXPECT_NEAR(DvDv[d], (Dv_pv[d] - Dv_mv[d]) / (2 * h), tol); + EXPECT_NEAR(DuDv[d], (Du_pv[d] - Du_mv[d]) / (2 * h), tol); + EXPECT_NEAR(DuDv[d], (Dv_pu[d] - Dv_mu[d]) / (2 * h), tol); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, bounding_box) +{ + SLIC_INFO("Testing Gregory triangle bounding box"); + + const auto controlPoints = GregoryTriangleType(make_sample_bezier_triangle()).getControlPoints(); + GregoryTriangleType triangle(controlPoints); + const auto bbox = triangle.boundingBox(); + const auto obb = triangle.orientedBoundingBox(); + + for(int i = 0; i < GregoryTriangleType::NPTS; ++i) + { + EXPECT_TRUE(bbox.contains(controlPoints[i])); + EXPECT_TRUE(obb.contains(controlPoints[i])); + } +} + +//------------------------------------------------------------------------------ +TEST(primal_gregorytriangle, print) +{ + SLIC_INFO("Testing Gregory triangle output stream"); + + GregoryTriangleType triangle(make_sample_bezier_triangle()); + std::ostringstream oss; + oss << triangle; + + EXPECT_NE(std::string::npos, oss.str().find("GregoryTriangle(")); +} From efd7fcb1e45c34f7b5e2cf0432c524b105a99008 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 13 Aug 2026 13:03:59 -0700 Subject: [PATCH 08/16] Added missing constructors --- src/axom/primal/geometry/GregoryTriangle.hpp | 44 ++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/axom/primal/geometry/GregoryTriangle.hpp b/src/axom/primal/geometry/GregoryTriangle.hpp index baee6dffff..c1e77d151f 100644 --- a/src/axom/primal/geometry/GregoryTriangle.hpp +++ b/src/axom/primal/geometry/GregoryTriangle.hpp @@ -80,12 +80,53 @@ class GregoryTriangle explicit GregoryTriangle(ArrayView controlPoints) { SLIC_ASSERT(controlPoints.size() == NPTS); + SLIC_ASSERT(controlPoints.data() != nullptr); for(int i = 0; i < NPTS; ++i) { m_controlPoints[i] = controlPoints[i]; } } + explicit GregoryTriangle(ArrayView controlPoints) + : GregoryTriangle(ArrayView(controlPoints.data(), controlPoints.size())) + { } + + explicit GregoryTriangle(const PointType* pts) : GregoryTriangle(ArrayView(pts, NPTS)) + { } + + explicit GregoryTriangle(PointType* pts) : GregoryTriangle(ArrayView(pts, NPTS)) { } + + explicit GregoryTriangle(const CoordsVec& pts) + : GregoryTriangle(ArrayView(pts.data(), pts.size())) + { } + + explicit GregoryTriangle(const BezierTriangle& bTri) + { + SLIC_ASSERT(bTri.getOrder() == 4); + SLIC_ASSERT(!bTri.isRational()); + + getCorner(0) = bTri(0, 0); + getCorner(1) = bTri(0, 4); + getCorner(2) = bTri(4, 0); + + getBoundaryPoint(0, 1) = bTri(1, 3); + getBoundaryPoint(0, 2) = bTri(2, 2); + getBoundaryPoint(0, 3) = bTri(3, 1); + getBoundaryPoint(1, 1) = bTri(3, 0); + getBoundaryPoint(1, 2) = bTri(2, 0); + getBoundaryPoint(1, 3) = bTri(1, 0); + getBoundaryPoint(2, 1) = bTri(0, 1); + getBoundaryPoint(2, 2) = bTri(0, 2); + getBoundaryPoint(2, 3) = bTri(0, 3); + + getTangent(0, 0) = bTri(1, 2); + getTangent(0, 1) = bTri(2, 1); + getTangent(1, 0) = bTri(2, 1); + getTangent(1, 1) = bTri(1, 1); + getTangent(2, 0) = bTri(1, 1); + getTangent(2, 1) = bTri(1, 2); + } + GregoryTriangle(ArrayView nodePositions, ArrayView nodeVectors) { // Store the position and orthogonal unit vector at each corner @@ -201,6 +242,9 @@ class GregoryTriangle return m_controlPoints[s_edge_index_map[e][k]]; } + CoordsVec& getControlPoints() { return m_controlPoints; } + const CoordsVec& getControlPoints() const { return m_controlPoints; } + // Evaluate the triangle by constructing the equivalent Bezier triangle with interior control nodes // defined in terms of the tangent vectors and the evaluation parameters PointType evaluate(T u0, T v0) const From d041ed441d3514b43166d9de9b659e6ad3188440 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 17 Aug 2026 13:49:18 -0700 Subject: [PATCH 09/16] Fix int double bug --- src/axom/primal/geometry/GregoryPatch.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/geometry/GregoryPatch.hpp b/src/axom/primal/geometry/GregoryPatch.hpp index 1d62957709..ef5a3f4d41 100644 --- a/src/axom/primal/geometry/GregoryPatch.hpp +++ b/src/axom/primal/geometry/GregoryPatch.hpp @@ -222,9 +222,9 @@ class GregoryPatch const VectorType dx(getCorner(i), getCorner(ip1)); - c0[i] = (dx - dx.dot(v[i]) * v[i]) / 3; + c0[i] = (dx - dx.dot(v[i]) * v[i]) / 3.0; a0[i] = VectorType::cross_product(v[i], dx).unitVector(); - c2[i] = (dx - dx.dot(v[ip1]) * v[ip1]) / 3; + c2[i] = (dx - dx.dot(v[ip1]) * v[ip1]) / 3.0; a3[i] = VectorType::cross_product(v[ip1], dx).unitVector(); // Use Chiyokura algorithm to define the interior control points From 5dd88077d53a8837490371ef570f9f5cc3626502 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 12:33:10 -0700 Subject: [PATCH 10/16] fix failing tests --- src/axom/primal/tests/primal_bezier_triangle.cpp | 4 ++-- src/axom/primal/tests/primal_gregory_surfaces.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/axom/primal/tests/primal_bezier_triangle.cpp b/src/axom/primal/tests/primal_bezier_triangle.cpp index 75ff6a8ae5..a5c35aa06f 100644 --- a/src/axom/primal/tests/primal_bezier_triangle.cpp +++ b/src/axom/primal/tests/primal_bezier_triangle.cpp @@ -334,8 +334,8 @@ TEST(primal_beziertriangle, finite_difference_first_derivatives) const PointType fp_v = tri.evaluate(u0, v0 + h); const PointType fm_v = tri.evaluate(u0, v0 - h); - const VectorType Du_fd(fp_u, fm_u); - const VectorType Dv_fd(fp_v, fm_v); + const VectorType Du_fd(fm_u, fp_u); + const VectorType Dv_fd(fm_v, fp_v); for(int d = 0; d < 3; ++d) { diff --git a/src/axom/primal/tests/primal_gregory_surfaces.cpp b/src/axom/primal/tests/primal_gregory_surfaces.cpp index 4691a90303..91de6d7463 100644 --- a/src/axom/primal/tests/primal_gregory_surfaces.cpp +++ b/src/axom/primal/tests/primal_gregory_surfaces.cpp @@ -425,8 +425,8 @@ TEST(primal_gregorytriangle, finite_difference_first_derivatives) const PointType fp_v = triangle.evaluate(u0, v0 + h); const PointType fm_v = triangle.evaluate(u0, v0 - h); - const VectorType Du_fd(fp_u, fm_u); - const VectorType Dv_fd(fp_v, fm_v); + const VectorType Du_fd(fm_u, fp_u); + const VectorType Dv_fd(fm_v, fp_v); for(int d = 0; d < 3; ++d) { From edbec60710f0224a96ac5524bdf9575b42e3a6f1 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 12:33:18 -0700 Subject: [PATCH 11/16] tidy up documentation --- src/axom/primal/geometry/BezierTriangle.hpp | 37 +++++++++++---------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/src/axom/primal/geometry/BezierTriangle.hpp b/src/axom/primal/geometry/BezierTriangle.hpp index f08c590ccc..8ace4ef9b3 100644 --- a/src/axom/primal/geometry/BezierTriangle.hpp +++ b/src/axom/primal/geometry/BezierTriangle.hpp @@ -55,23 +55,26 @@ std::ostream& operator<<(std::ostream& os, const BezierTriangle& 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 @@ -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, From 1db395f567aaa48ce0540d81890e911b43b95b07 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 12:39:54 -0700 Subject: [PATCH 12/16] Fix style! --- src/axom/primal/tests/primal_bezier_triangle.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/axom/primal/tests/primal_bezier_triangle.cpp b/src/axom/primal/tests/primal_bezier_triangle.cpp index a5c35aa06f..30c259e58a 100644 --- a/src/axom/primal/tests/primal_bezier_triangle.cpp +++ b/src/axom/primal/tests/primal_bezier_triangle.cpp @@ -281,9 +281,9 @@ TEST(primal_beziertriangle, parameter_convention_matches_barycentric) constexpr CoordType eps = 1e-14; BTri tri(1); - const PointType pA {1.0, 2.0, 3.0}; // (i,j) = (0,0) - const PointType pB {4.0, -1.0, 0.5}; // (i,j) = (0,1) - const PointType pC {-2.0, 0.25, 7.0}; // (i,j) = (1,0) + const PointType pA {1.0, 2.0, 3.0}; // (i,j) = (0,0) + const PointType pB {4.0, -1.0, 0.5}; // (i,j) = (0,1) + const PointType pC {-2.0, 0.25, 7.0}; // (i,j) = (1,0) tri(0, 0) = pA; tri(0, 1) = pB; From 3fae1dca4929b68c4c95890b2db5c0849b06bc61 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 12:53:51 -0700 Subject: [PATCH 13/16] Fix style again --- src/axom/primal/geometry/BezierCurve.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index aa91909191..deb8ee0a36 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -441,8 +441,8 @@ class BezierCurve for(int i = 1; i <= n; ++i) { const T alpha = static_cast(i) / static_cast(np1); - newPts[i] = - PointType(alpha * m_controlPoints[i - 1].array() + (T(1) - alpha) * m_controlPoints[i].array()); + newPts[i] = PointType(alpha * m_controlPoints[i - 1].array() + + (T(1) - alpha) * m_controlPoints[i].array()); } m_controlPoints = newPts; From 8a53cee3a642e5e2fc7f17359b8e97049689dcf6 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 13:35:03 -0700 Subject: [PATCH 14/16] Fix style again! --- src/axom/primal/tests/primal_gregory_surfaces.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/axom/primal/tests/primal_gregory_surfaces.cpp b/src/axom/primal/tests/primal_gregory_surfaces.cpp index 91de6d7463..0195aa5d15 100644 --- a/src/axom/primal/tests/primal_gregory_surfaces.cpp +++ b/src/axom/primal/tests/primal_gregory_surfaces.cpp @@ -73,7 +73,8 @@ void check_patch_control_points(const GregoryPatchType& patch, const PatchCoords } } -void check_triangle_control_points(const GregoryTriangleType& triangle, const TriangleCoordsVec& expected) +void check_triangle_control_points(const GregoryTriangleType& triangle, + const TriangleCoordsVec& expected) { for(int i = 0; i < GregoryTriangleType::NPTS; ++i) { From 162e98abce33f50da4b6cb29315cd369580974a0 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 13:44:51 -0700 Subject: [PATCH 15/16] final style --- src/axom/primal/geometry/GregoryTriangle.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/axom/primal/geometry/GregoryTriangle.hpp b/src/axom/primal/geometry/GregoryTriangle.hpp index c1e77d151f..10e206bb42 100644 --- a/src/axom/primal/geometry/GregoryTriangle.hpp +++ b/src/axom/primal/geometry/GregoryTriangle.hpp @@ -91,10 +91,12 @@ class GregoryTriangle : GregoryTriangle(ArrayView(controlPoints.data(), controlPoints.size())) { } - explicit GregoryTriangle(const PointType* pts) : GregoryTriangle(ArrayView(pts, NPTS)) + explicit GregoryTriangle(const PointType* pts) + : GregoryTriangle(ArrayView(pts, NPTS)) { } - explicit GregoryTriangle(PointType* pts) : GregoryTriangle(ArrayView(pts, NPTS)) { } + explicit GregoryTriangle(PointType* pts) : GregoryTriangle(ArrayView(pts, NPTS)) + { } explicit GregoryTriangle(const CoordsVec& pts) : GregoryTriangle(ArrayView(pts.data(), pts.size())) @@ -327,11 +329,9 @@ class GregoryTriangle B_v0[2] * intermediate.Q_v[2]) + (B[0] * intermediate.Q_vv[0] + B[1] * intermediate.Q_vv[1] + B[2] * intermediate.Q_vv[2]); - DuDv += - (B_u0[0] * intermediate.Q_v[0] + B_u0[1] * intermediate.Q_v[1] + - B_u0[2] * intermediate.Q_v[2]) + - (B_v0[0] * intermediate.Q_u[0] + B_v0[1] * intermediate.Q_u[1] + - B_v0[2] * intermediate.Q_u[2]) + + DuDv += (B_u0[0] * intermediate.Q_v[0] + B_u0[1] * intermediate.Q_v[1] + + B_u0[2] * intermediate.Q_v[2]) + + (B_v0[0] * intermediate.Q_u[0] + B_v0[1] * intermediate.Q_u[1] + B_v0[2] * intermediate.Q_u[2]) + (B[0] * intermediate.Q_uv[0] + B[1] * intermediate.Q_uv[1] + B[2] * intermediate.Q_uv[2]); } From b83f4bbda2e5da8174e9fd38a444d8ded4851037 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 24 Aug 2026 14:24:01 -0700 Subject: [PATCH 16/16] Update release notes --- RELEASE-NOTES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 879640cd67..25eec60ff9 100644 --- a/RELEASE-NOTES.md +++ b/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.