From af6637f3457736e93f9f8080efc0b4d3dcae2a88 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 16 Sep 2025 15:43:23 -0700 Subject: [PATCH 01/21] Refactor CurvedPolygon to take arbitrary type --- src/axom/primal/geometry/BezierCurve.hpp | 3 + src/axom/primal/geometry/CurvedPolygon.hpp | 85 +++++--- src/axom/primal/geometry/NURBSCurve.hpp | 8 + src/axom/primal/operators/compute_moments.hpp | 4 +- .../primal/operators/in_curved_polygon.hpp | 2 +- src/axom/primal/operators/winding_number.hpp | 44 +++- .../primal/tests/primal_curved_polygon.cpp | 196 ++++++++++++------ src/axom/quest/SamplingShaper.hpp | 2 +- .../detail/shaping/WindingNumberSampler.hpp | 10 +- src/axom/quest/io/MFEMReader.cpp | 4 +- 10 files changed, 255 insertions(+), 103 deletions(-) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 2db53b2e2c..1d9a24bb1f 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -256,6 +256,9 @@ class BezierCurve /// Retrieves the control point at index \a idx const PointType& operator[](int idx) const { return m_controlPoints[idx]; } + PointType getInitPoint() const { return m_controlPoints[0]; } + PointType getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + /// Returns a copy of the Bezier curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index a4afe50567..0e95fbb3d9 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -18,6 +18,7 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" #include "axom/primal/geometry/BoundingBox.hpp" #include @@ -27,13 +28,43 @@ namespace axom { namespace primal { -// Forward declare the templated classes and operator functions +// Boilerplate to make the containers work as expected +template +struct get_numeric_type; + +template +struct get_numeric_type> +{ + using type = T; +}; + template +struct get_numeric_type> +{ + using type = T; +}; + +template +struct get_numeric_type> +{ + using type = T; +}; + +template +struct has_cached_data : std::false_type +{ }; + +template +struct has_cached_data> : std::true_type +{ }; + +// Forward declare the templated classes and operator functions +template class CurvedPolygon; /*! \brief Overloaded output operator for polygons */ -template -std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly); +template +std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly); /*! * \class CurvedPolygon @@ -44,15 +75,15 @@ std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly); * \note The component curves should be ordered in a counter clockwise * orientation with respect to the polygon's normal vector */ -template +template class CurvedPolygon { public: - using PointType = Point; - using VectorType = Vector; - using NumArrayType = axom::NumericArray; - using BezierCurveType = BezierCurve; - using BoundingBoxType = typename BezierCurveType::BoundingBoxType; + using T = typename get_numeric_type::type; + + using PointType = typename CurveType::PointType; + using VectorType = typename CurveType::VectorType; + using BoundingBoxType = typename CurveType::BoundingBoxType; public: /// Default constructor for an empty polygon @@ -74,7 +105,7 @@ class CurvedPolygon } /// Constructor from an array of \a nEdges curves - CurvedPolygon(BezierCurveType* curves, int nEdges) + CurvedPolygon(CurveType* curves, int nEdges) { SLIC_ASSERT(curves != nullptr); SLIC_ASSERT(nEdges >= 1); @@ -106,35 +137,39 @@ class CurvedPolygon } /// Appends a BezierCurve to the list of edges - void addEdge(const BezierCurveType& c1) { m_edges.push_back(c1); } + void addEdge(const CurveType& c1) { m_edges.push_back(c1); } /// Splits an edge "in place" void splitEdge(int idx, T t) { SLIC_ASSERT(idx < static_cast(m_edges.size())); + AXOM_STATIC_ASSERT_MSG(!has_cached_data::value, + "splitEdge cannot be called on objects with associated cache data"); m_edges.insert(m_edges.begin() + idx + 1, 1, m_edges[idx]); auto& csplit = m_edges[idx]; csplit.split(t, m_edges[idx], m_edges[idx + 1]); } - axom::Array> getEdges() const { return m_edges; } + axom::Array getEdges() const { return m_edges; } /// @} - /*! Retrieves the Bezier Curve at index idx */ - BezierCurveType& operator[](int idx) { return m_edges[idx]; } + /*! Retrieves the curve at index idx */ + CurveType& operator[](int idx) { return m_edges[idx]; } /*! Retrieves the vertex at index idx */ - const BezierCurveType& operator[](int idx) const { return m_edges[idx]; } + const CurveType& operator[](int idx) const { return m_edges[idx]; } /// Tests equality of two CurvedPolygons - friend inline bool operator==(const CurvedPolygon& lhs, const CurvedPolygon& rhs) + friend inline bool operator==(const CurvedPolygon& lhs, + const CurvedPolygon& rhs) { return lhs.m_edges == rhs.m_edges; } /// Tests inequality of two CurvedPolygons - friend inline bool operator!=(const CurvedPolygon& lhs, const CurvedPolygon& rhs) + friend inline bool operator!=(const CurvedPolygon& lhs, + const CurvedPolygon& rhs) { return !(lhs == rhs); } @@ -149,7 +184,7 @@ class CurvedPolygon { const int sz = numEdges(); - os << "{" << sz << "-sided Bezier polygon:"; + os << "{" << sz << "-sided curved polygon:"; for(int i = 0; i < sz - 1; ++i) { os << m_edges[i] << ","; @@ -185,9 +220,8 @@ class CurvedPolygon // foreach edge: check last vertex of previous edge against first vertex of current edge for(int cur = 0, prev = nEdges - 1; cur < nEdges; prev = cur++) { - const auto ord = m_edges[prev].getOrder(); - const auto& lastPrev = m_edges[prev][ord]; - const auto& firstCur = m_edges[cur][0]; + const auto& lastPrev = m_edges[prev].getEndPoint(); + const auto& firstCur = m_edges[cur].getInitPoint(); if(!isNearlyEqual(squared_distance(lastPrev, firstCur), 0., sq_tol)) { return false; @@ -200,6 +234,9 @@ class CurvedPolygon /// \brief Reverses orientation of a CurvedPolygon void reverseOrientation() { + AXOM_STATIC_ASSERT_MSG(!has_cached_data::value, + "reverseOrientation cannot be called on objects with associated cache data"); + const int ngon = numEdges(); const int mid = ngon >> 1; const bool isOdd = (ngon & 1) != 0; @@ -233,14 +270,14 @@ class CurvedPolygon } private: - axom::Array> m_edges; + axom::Array m_edges; }; //------------------------------------------------------------------------------ /// Free functions implementing Polygon's operators //------------------------------------------------------------------------------ -template -std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly) +template +std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly) { poly.print(os); return os; diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index 628fb4b830..2c24df6057 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -78,6 +78,11 @@ struct SubdivisionData template class NURBSCurveGWNCache { +public: + using PointType = typename NURBSCurve::PointType; + using VectorType = typename NURBSCurve::VectorType; + using BoundingBoxType = typename NURBSCurve::BoundingBoxType; + public: NURBSCurveGWNCache() = default; @@ -759,6 +764,9 @@ class NURBSCurve */ const PointType& operator[](int idx) const { return m_controlPoints[idx]; } + PointType getInitPoint() const { return m_controlPoints[0]; } + PointType getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + /// \brief Returns a copy of the NURBS curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } diff --git a/src/axom/primal/operators/compute_moments.hpp b/src/axom/primal/operators/compute_moments.hpp index ccc3492623..8acdfd6041 100644 --- a/src/axom/primal/operators/compute_moments.hpp +++ b/src/axom/primal/operators/compute_moments.hpp @@ -95,7 +95,7 @@ primal::Point sector_centroid(const primal::BezierCurve& curve) /// \brief Returns the area enclosed by the CurvedPolygon template -T area(const primal::CurvedPolygon& poly, double tol = 1e-8) +T area(const primal::CurvedPolygon>& poly, double tol = 1e-8) { const int ngon = poly.numEdges(); T A = 0.0; @@ -119,7 +119,7 @@ T area(const primal::CurvedPolygon& poly, double tol = 1e-8) /// \brief Returns the centroid of the CurvedPolygon template -primal::Point centroid(const primal::CurvedPolygon& poly, double tol = 1e-8) +primal::Point centroid(const primal::CurvedPolygon>& poly, double tol = 1e-8) { using PointType = primal::Point; diff --git a/src/axom/primal/operators/in_curved_polygon.hpp b/src/axom/primal/operators/in_curved_polygon.hpp index ea8356de0d..9fe1ff7885 100644 --- a/src/axom/primal/operators/in_curved_polygon.hpp +++ b/src/axom/primal/operators/in_curved_polygon.hpp @@ -48,7 +48,7 @@ namespace primal */ template bool in_curved_polygon(const Point& query, - const CurvedPolygon& cpoly, + const CurvedPolygon>& cpoly, bool useNonzeroRule = true, double edge_tol = 1e-8, double EPS = 1e-8) diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index 13d547f039..442e17b2df 100644 --- a/src/axom/primal/operators/winding_number.hpp +++ b/src/axom/primal/operators/winding_number.hpp @@ -181,19 +181,16 @@ double winding_number(const Point& q, * * \return The GWN. */ -template +template double winding_number(const Point& q, - const CurvedPolygon& cpoly, + const CurvedPolygon& cpoly, double edge_tol = 1e-8, double EPS = 1e-8) { - bool dummy_isOnCurve = false, isConvexControlPolygon = false; - double ret_val = 0.0; for(int i = 0; i < cpoly.numEdges(); i++) { - ret_val += - detail::bezier_winding_number(q, cpoly[i], isConvexControlPolygon, dummy_isOnCurve, edge_tol, EPS); + ret_val += winding_number(q, cpoly[i], edge_tol, EPS); } return ret_val; @@ -425,9 +422,9 @@ axom::Array winding_number(const axom::Array>& q_arr, * * \return The GWN. */ -template +template axom::Array winding_number(const axom::Array>& q_arr, - const CurvedPolygon& cpoly, + const CurvedPolygon& cpoly, double edge_tol = 1e-8, double EPS = 1e-8) { @@ -446,6 +443,37 @@ axom::Array winding_number(const axom::Array>& q_arr, return ret_val; } + +/*! + * \brief Computes the GWN for an array of 2D points wrt to a 2D curved polygon + * + * \param [in] q_arr The query point to test + * \param [in] cpoly The CurvedPolygon object + * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable + * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances + * + * Computes the GWN for the curved polygon by summing the GWN for each curved edge + * + * \warning Because the cache isKdiscarded immediately after computation, + * this method is not accelerated by memoization + * + * \return The GWN. + */ +template +axom::Array winding_number(const axom::Array>& q_arr, + const CurvedPolygon>& cpoly, + double edge_tol = 1e-8, + double EPS = 1e-8) +{ + axom::Array ret_val(q_arr.size()); + for(int n = 0; n < q_arr.size(); ++n) + { + ret_val[n] = winding_number(q_arr[n], cache_arr, edge_tol, EPS); + } + + return ret_val; +} + ///@{ //! @name Winding number operations between 3D points and primitives diff --git a/src/axom/primal/tests/primal_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp index 80f1cac48e..ae17502096 100644 --- a/src/axom/primal/tests/primal_curved_polygon.cpp +++ b/src/axom/primal/tests/primal_curved_polygon.cpp @@ -31,7 +31,7 @@ namespace primal = axom::primal; * stored in \a expArea and \a expCentroid. Areas and Moments are computed within tolerance \a eps and checks use \a test_eps. */ template -void checkMoments(const primal::CurvedPolygon& bPolygon, +void checkMoments(const primal::CurvedPolygon>& bPolygon, const CoordType expArea, const primal::Point& expMoment, double eps, @@ -47,19 +47,19 @@ void checkMoments(const primal::CurvedPolygon& bPolygon, } /*! - * Helper function to create a CurvedPolygon from a list of control points and a list of orders of component curves. + * Helper function to create a Bezier CurvedPolygon from a list of control points and a list of orders of component curves. * Control points should be given as a list of Points in order of orientation with no duplicates except that * the first control point should also be the last control point (if the polygon is closed). * Orders should be given as a list of ints in order of orientation, representing the orders of the component curves. */ template -primal::CurvedPolygon createPolygon( +primal::CurvedPolygon> createBezierPolygon( const axom::Array> ControlPoints, const axom::Array orders) { using PointType = primal::Point; - using CurvedPolygonType = primal::CurvedPolygon; using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; const int num_edges = orders.size(); const int num_unique_control_points = ControlPoints.size(); @@ -93,8 +93,8 @@ TEST(primal_curvedpolygon, constructor) { const int DIM = 3; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; { SLIC_INFO("Testing default CurvedPolygon constructor "); @@ -121,9 +121,9 @@ TEST(primal_curvedpolygon, add_edges) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; SLIC_INFO("Test adding edges to empty CurvedPolygon"); @@ -152,12 +152,88 @@ TEST(primal_curvedpolygon, add_edges) } } +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, add_edges_nurbs) +{ + const int DIM = 2; + using CoordType = double; + using CurveType = primal::NURBSCurve; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test adding edges to empty CurvedPolygon"); + + CurvedPolygonType nPolygon; + EXPECT_EQ(0, nPolygon.numEdges()); + + PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; + + CurveType nCurve(controlPoints, 2, 1); + + nPolygon.addEdge(nCurve); + nPolygon.addEdge(nCurve); + + EXPECT_EQ(2, nPolygon.numEdges()); + for(int p = 0; p < nPolygon.numEdges(); ++p) + { + CurveType& nc = nPolygon[p]; + for(int sz = 0; sz <= nc.getDegree(); ++sz) + { + auto& pt = nc[sz]; + for(int i = 0; i < DIM; ++i) + { + EXPECT_DOUBLE_EQ(controlPoints[sz][i], pt[i]); + } + } + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, add_edges_nurbs_cache) +{ + const int DIM = 2; + using CoordType = double; + using CurveType = primal::NURBSCurveGWNCache; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test adding edges to empty CurvedPolygon"); + + CurvedPolygonType nPolygon; + EXPECT_EQ(0, nPolygon.numEdges()); + + PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; + + CurveType nCurve(primal::NURBSCurve(controlPoints, 2, 1)); + + nPolygon.addEdge(nCurve); + nPolygon.addEdge(nCurve); + + EXPECT_EQ(2, nPolygon.numEdges()); + for(int p = 0; p < nPolygon.numEdges(); ++p) + { + CurveType& nc = nPolygon[p]; + auto& init_pt = nc.getInitPoint(); + for(int i = 0; i < DIM; ++i) + { + EXPECT_DOUBLE_EQ(controlPoints[0][i], init_pt[i]); + } + + auto& end_pt = nc.getEndPoint(); + for(int i = 0; i < DIM; ++i) + { + EXPECT_DOUBLE_EQ(controlPoints[1][i], end_pt[i]); + } + } +} + //---------------------------------------------------------------------------------- TEST(primal_curvedpolygon, isClosed) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; + using CurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; SLIC_INFO("Test checking if CurvedPolygon is closed."); @@ -177,12 +253,12 @@ TEST(primal_curvedpolygon, isClosed) { axom::Array subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; axom::Array suborders = {1}; - CurvedPolygonType subPolygon = createPolygon(subCP, suborders); + CurvedPolygonType subPolygon = createBezierPolygon(subCP, suborders); EXPECT_FALSE(subPolygon.isClosed()); } { - CurvedPolygonType bPolygon = createPolygon(CP, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); EXPECT_EQ(3, bPolygon.numEdges()); EXPECT_TRUE(bPolygon.isClosed()); @@ -191,7 +267,7 @@ TEST(primal_curvedpolygon, isClosed) } { - CurvedPolygonType bPolygon = createPolygon(CP, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); bPolygon[1][0][0] = 5; EXPECT_FALSE(bPolygon.isClosed(1e-15)); @@ -199,43 +275,43 @@ TEST(primal_curvedpolygon, isClosed) } //---------------------------------------------------------------------------------- -TEST(primal_curvedpolygon, isClosed_BiGon) -{ - const int DIM = 2; - using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; - - SLIC_INFO("Test checking if CurvedPolygon is closed for a Bi-Gon"); - - CurvedPolygonType bPolygon; - EXPECT_EQ(0, bPolygon.numEdges()); - EXPECT_FALSE(bPolygon.isClosed()); - - // Bi-gon defined by a quadratic edge and a straight line - axom::Array CP = {PointType {0.8, .25}, - PointType {2.0, .50}, - PointType {0.8, .75}, - PointType {0.8, .25}}; - axom::Array orders = {2, 1}; - - CurvedPolygonType poly = createPolygon(CP, orders); - EXPECT_TRUE(poly.isClosed()); - - // modify a vertex of the quadratic and check again - CurvedPolygonType poly2 = poly; - poly2[0][2] = PointType {0.8, 1.0}; - EXPECT_FALSE(poly2.isClosed()); -} +//TEST(primal_curvedpolygon, isClosed_BiGon) +//{ +// const int DIM = 2; +// using CoordType = double; +// using CurvedPolygonType = primal::CurvedPolygon; +// using PointType = primal::Point; +// +// SLIC_INFO("Test checking if CurvedPolygon is closed for a Bi-Gon"); +// +// CurvedPolygonType bPolygon; +// EXPECT_EQ(0, bPolygon.numEdges()); +// EXPECT_FALSE(bPolygon.isClosed()); +// +// // Bi-gon defined by a quadratic edge and a straight line +// axom::Array CP = {PointType {0.8, .25}, +// PointType {2.0, .50}, +// PointType {0.8, .75}, +// PointType {0.8, .25}}; +// axom::Array orders = {2, 1}; +// +// CurvedPolygonType poly = createBezierPolygon(CP, orders); +// EXPECT_TRUE(poly.isClosed()); +// +// // modify a vertex of the quadratic and check again +// CurvedPolygonType poly2 = poly; +// poly2[0][2] = PointType {0.8, 1.0}; +// EXPECT_FALSE(poly2.isClosed()); +//} //---------------------------------------------------------------------------------- TEST(primal_curvedpolygon, split_edge) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; + using CurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; - using BezierCurveType = primal::BezierCurve; SLIC_INFO("Test checking CurvedPolygon edge split."); @@ -245,18 +321,18 @@ TEST(primal_curvedpolygon, split_edge) PointType {0.6, 1.2}}; axom::Array orders32 = {1, 1, 1}; - CurvedPolygonType bPolygon32 = createPolygon(CP, orders32); + CurvedPolygonType bPolygon32 = createBezierPolygon(CP, orders32); // Needs to be std::vector to use .assign axom::Array subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; - BezierCurveType bCurve(subCP, 1); + CurveType bCurve(subCP, 1); //std::cout << "Got here!! " << std::endl; //std::cout << bPolygon32 << std::endl; bPolygon32.splitEdge(0, .5); - BezierCurveType bCurve2; - BezierCurveType bCurve3; + CurveType bCurve2; + CurveType bCurve3; bCurve.split(.5, bCurve2, bCurve3); EXPECT_EQ(bPolygon32.numEdges(), 4); @@ -275,9 +351,9 @@ TEST(primal_curvedpolygon, moments_triangle_degenerate) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; SLIC_INFO("Testing area computation of degenerate triangles"); @@ -311,8 +387,9 @@ TEST(primal_curvedpolygon, moments_triangle_linear) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; SLIC_INFO("Test moment computation of a linear triangle"); axom::Array CP = {PointType {0.6, 1.2}, @@ -321,7 +398,7 @@ TEST(primal_curvedpolygon, moments_triangle_linear) PointType {0.6, 1.2}}; axom::Array orders = {1, 1, 1}; - CurvedPolygonType bPolygon = createPolygon(CP, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); CoordType trueA = -.18; PointType trueC = PointType::make_point(0.3, 1.6); @@ -334,8 +411,9 @@ TEST(primal_curvedpolygon, moments_triangle_quadratic) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; SLIC_INFO("Test moment computation of quadratic triangle"); @@ -348,7 +426,7 @@ TEST(primal_curvedpolygon, moments_triangle_quadratic) PointType {0.6, 1.2}}; axom::Array orders = {2, 2, 2}; - CurvedPolygonType bPolygon = createPolygon(CP, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); CoordType trueA = -0.097333333333333; PointType trueC {.294479452054794, 1.548219178082190}; @@ -361,8 +439,9 @@ TEST(primal_curvedpolygon, moments_triangle_mixed_order) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; SLIC_INFO("Test moment computation for curved triangle with mixed order edges"); @@ -374,7 +453,7 @@ TEST(primal_curvedpolygon, moments_triangle_mixed_order) PointType {0.6, 1.2}}; axom::Array orders = {2, 2, 1}; - CurvedPolygonType bPolygon = createPolygon(CP, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); CoordType trueA = -.0906666666666666666666; PointType trueC {.2970147058823527, 1.55764705882353}; @@ -387,8 +466,9 @@ TEST(primal_curvedpolygon, moments_quad_all_orders) { const int DIM = 2; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; SLIC_INFO("Test moment computation for quads of different orders"); axom::Array CPorig = {PointType {0.0, 0.0}, @@ -398,7 +478,7 @@ TEST(primal_curvedpolygon, moments_quad_all_orders) PointType {0.0, 0.0}}; axom::Array orders = {1, 1, 1, 1}; - CurvedPolygonType bPolygon = createPolygon(CPorig, orders); + CurvedPolygonType bPolygon = createBezierPolygon(CPorig, orders); CoordType trueA = 1.0; PointType trueC = PointType::make_point(0.5, 0.5); @@ -436,7 +516,7 @@ TEST(primal_curvedpolygon, moments_quad_all_orders) // std::cout << CP[i] << std::endl; // } - bPolygon = createPolygon(CP, orders); + bPolygon = createBezierPolygon(CP, orders); checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); } } @@ -447,10 +527,10 @@ TEST(primal_curvedpolygon, reverseOrientation) const int DIM = 2; const int order = 1; using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; using PointType = primal::Point; - using SegmentType = primal::Segment; using BezierCurveType = primal::BezierCurve; + using CurvedPolygonType = primal::CurvedPolygon; + using SegmentType = primal::Segment; // Test several n-gons discretizing the unit circle const int MAX_SEG = 10; @@ -520,7 +600,7 @@ TEST(primal_curved_polygon, containment) // Test that containment procedures are consistent using Point2D = primal::Point; using Bezier = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; // 8th order, closed curve with internal loop Point2D loop_nodes[] = {Point2D {0.0, 0.0}, diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index eb268c1c76..6c2f52ece1 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -961,7 +961,7 @@ class SamplingShaper : public Shaper std::unique_ptr> m_primitiveSampler3D_hip; std::unique_ptr> m_inoutSamplerWN; - axom::Array> m_contours; + axom::Array>> m_contours; std::set m_knownMaterials; diff --git a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp index a3758430b3..a6c63dab32 100644 --- a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp +++ b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp @@ -45,13 +45,8 @@ namespace detail template bool AXOM_HOST_DEVICE checkInside(const ShapeType& shape, const PointType& pt) { - double wn {}; - for(int c = 0; c < shape.numEdges(); c++) - { - wn += axom::primal::winding_number(pt, shape[c]); - } // A point inside the polygon should have non-zero winding number. - return (round(wn) != 0); + return (round(axom::primal::winding_number(pt, shape)) != 0); } } // end namespace detail @@ -68,7 +63,8 @@ class WindingNumberSampler // For now. using ExecSpace = axom::SEQ_EXEC; - using GeometryView = typename axom::ArrayView>; + using CurvedPolygonType = primal::CurvedPolygon>; + using GeometryView = typename axom::ArrayView; using PointType = primal::Point; using GeometricBoundingBox = axom::primal::BoundingBox; using BVH = typename axom::spin::BVH; diff --git a/src/axom/quest/io/MFEMReader.cpp b/src/axom/quest/io/MFEMReader.cpp index 40decb987d..2a03ca6d53 100644 --- a/src/axom/quest/io/MFEMReader.cpp +++ b/src/axom/quest/io/MFEMReader.cpp @@ -127,8 +127,8 @@ int MFEMReader::read(CurvedPolygonArray &curvedPolygons) auto &poly = curvedPolygons[contourId]; for(int zoneId : zoneIds) { - auto curve = axom::quest::internal::segment_to_curve(mesh, zoneId); - poly.addEdge(std::move(curve)); + auto curve = axom::quest::internal::segment_to_nurbs(mesh, zoneId); + poly.addEdge(std::move(axom::primal::NURBSCurveGWNCache(curve))); } } }); From 0b2df6c7a572e9d70144d96dec2f7616ea7b0cde Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 24 Sep 2025 11:08:44 -0700 Subject: [PATCH 02/21] Add extra methods to make compatible with template overloads --- src/axom/primal/geometry/BezierCurve.hpp | 3 + src/axom/primal/geometry/CurvedPolygon.hpp | 23 ++++--- .../detail/winding_number_2d_memoization.hpp | 61 +++++++++++++++++++ 3 files changed, 80 insertions(+), 7 deletions(-) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 1d9a24bb1f..ed5716c227 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -221,6 +221,9 @@ class BezierCurve /// Returns the order of the Bezier Curve int getOrder() const { return static_cast(m_controlPoints.size()) - 1; } + /// Returns the number of control points of the Bezier Curve + int getNumControlPoints() const { return static_cast(m_controlPoints.size()); } + /// Clears the list of control points, make nonrational void clear() { diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 0e95fbb3d9..42517b7417 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -21,6 +21,9 @@ #include "axom/primal/geometry/NURBSCurve.hpp" #include "axom/primal/geometry/BoundingBox.hpp" +// For NURBSCurveGWNCache objects +#include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" + #include #include @@ -28,7 +31,9 @@ namespace axom { namespace primal { -// Boilerplate to make the containers work as expected + +///@{ +/// \name Boilerplate to deduce the numeric type from the curve object template struct get_numeric_type; @@ -45,18 +50,22 @@ struct get_numeric_type> }; template -struct get_numeric_type> +struct get_numeric_type> { using type = T; }; +///@} +///@{ +/// \name Boilerplate for a compile-time flag for the cached object template struct has_cached_data : std::false_type { }; template -struct has_cached_data> : std::true_type +struct has_cached_data> : std::true_type { }; +///@} // Forward declare the templated classes and operator functions template @@ -152,7 +161,6 @@ class CurvedPolygon } axom::Array getEdges() const { return m_edges; } - /// @} /*! Retrieves the curve at index idx */ @@ -212,7 +220,7 @@ class CurvedPolygon const int nEdges = numEdges(); // initial basic check: no edges, or one edge or linear or quadratic order cannot be closed - if(nEdges < 1 || (nEdges == 1 && m_edges[0].getOrder() <= 2)) + if(nEdges < 1 || (nEdges == 1 && m_edges[0].getNumControlPoints() <= 2)) { return false; } @@ -234,8 +242,9 @@ class CurvedPolygon /// \brief Reverses orientation of a CurvedPolygon void reverseOrientation() { - AXOM_STATIC_ASSERT_MSG(!has_cached_data::value, - "reverseOrientation cannot be called on objects with associated cache data"); + AXOM_STATIC_ASSERT_MSG( + !has_cached_data::value, + "reverseOrientation cannot be called on objects with associated cache data"); const int ngon = numEdges(); const int mid = ngon >> 1; diff --git a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp index 580ac6f6e6..1b5ab962c0 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp @@ -59,12 +59,32 @@ struct BezierCurveData auto isConvexControlPolygon() const { return m_isConvexControlPolygon; } auto getBoundingBox() const { return m_boundingBox; } + friend bool operator==(const BezierCurveData& lhs, const BezierCurveData& rhs) + { + // isConvexControlPolygon will be equal if the curves are + return (lhs.m_curve == rhs.m_curve) && (lhs.m_boundingBox == rhs.m_boundingBox) && + (lhs.m_isConvexControlPolygon == rhs.m_isConvexControlPolygon); + } + + friend bool operator!=(const BezierCurveData& lhs, const BezierCurveData& rhs) + { + return !(lhs == rhs); + } + private: BezierCurve m_curve; bool m_isConvexControlPolygon; BoundingBox m_boundingBox; }; +// Forward declare the templated classes and operator functions +template +class NURBSCurveGWNCache; + +/// \brief Overloaded output operator for Cached Curves +template +std::ostream& operator<<(std::ostream& os, const NURBSCurveGWNCache& nCurveCache); + /*! * \class NURBSCurveGWNCache * @@ -78,6 +98,11 @@ struct BezierCurveData template class NURBSCurveGWNCache { +public: + using PointType = typename NURBSCurve::PointType; + using VectorType = typename NURBSCurve::VectorType; + using BoundingBoxType = typename NURBSCurve::BoundingBoxType; + public: NURBSCurveGWNCache() = default; @@ -167,6 +192,35 @@ class NURBSCurveGWNCache auto getInitPoint() const { return m_initPoint; } auto getEndPoint() const { return m_endPoint; } + friend bool operator==(const NURBSCurveGWNCache& lhs, const NURBSCurveGWNCache& rhs) + { + // numControlPoints, degree, and numSpans will be equal if the subdivision maps are + return (lhs.m_bezierSubdivisionMaps == rhs.m_bezierSubdivisionMaps) && + (lhs.m_boundingBox == rhs.m_boundingBox); + } + + friend bool operator!=(const NURBSCurveGWNCache& lhs, const NURBSCurveGWNCache& rhs) + { + return !(lhs == rhs); + } + + std::ostream& print(std::ostream& os) const + { + os << "{ NURBSCurveGWNCache object with " << m_numSpans << " extracted bezier curves: "; + + if(m_numSpans >= 1) + { + os << m_bezierSubdivisionMaps[0][std::make_pair(0, 0)].getCurve(); + } + for(int i = 1; i < m_numSpans; ++i) + { + os << ", " << m_bezierSubdivisionMaps[i][std::make_pair(0, 0)].getCurve(); + } + os << "}"; + + return os; + } + private: BoundingBox m_boundingBox; int m_numControlPoints; @@ -178,6 +232,13 @@ class NURBSCurveGWNCache mutable axom::Array, BezierCurveData>> m_bezierSubdivisionMaps; }; +template +std::ostream& operator<<(std::ostream& os, const NURBSCurveGWNCache& nCurveCache) +{ + nCurveCache.print(os); + return os; +} + } // namespace detail } // namespace primal } // namespace axom From 539248e105e33cf39261c730eb44ded0b53095ee Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 24 Sep 2025 11:09:15 -0700 Subject: [PATCH 03/21] Make use with new template for CurvedPolygon --- src/axom/primal/operators/evaluate_integral.hpp | 12 ++++++------ src/axom/primal/operators/winding_number.hpp | 2 +- src/axom/primal/tests/primal_bezier_curve.cpp | 5 +++++ src/axom/primal/tests/primal_winding_number.cpp | 4 ++-- 4 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index fb116c109e..d9fd6c49dd 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -54,8 +54,8 @@ namespace primal * on each edge of the CurvedPolygon * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, Lambda&& scalar_integrand, int npts) { @@ -176,8 +176,8 @@ double evaluate_scalar_line_integral(const axom::Array -double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, Lambda&& vector_integrand, int npts) { @@ -294,8 +294,8 @@ double evaluate_vector_line_integral(const axom::Array -double evaluate_area_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_area_integral(const primal::CurvedPolygon cpoly, Lambda&& integrand, int npts_Q, int npts_P = 0) diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index 6dfa9fad36..bc5fadc71f 100644 --- a/src/axom/primal/operators/winding_number.hpp +++ b/src/axom/primal/operators/winding_number.hpp @@ -471,7 +471,7 @@ axom::Array winding_number(const axom::Array>& q_arr, */ template axom::Array winding_number(const axom::Array>& q_arr, - const CurvedPolygon>& cpoly, + const CurvedPolygon>& cpoly, double edge_tol = 1e-8, double EPS = 1e-8) { diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index e6df8e9e6b..49e249f2a7 100644 --- a/src/axom/primal/tests/primal_bezier_curve.cpp +++ b/src/axom/primal/tests/primal_bezier_curve.cpp @@ -32,6 +32,7 @@ TEST(primal_beziercurve, sizing_constructors) constexpr int expOrder = -1; EXPECT_EQ(expOrder, bCurve.getOrder()); EXPECT_EQ(expOrder + 1, bCurve.getControlPoints().size()); + EXPECT_EQ(expOrder + 1, bCurve.getNumControlPoints()); EXPECT_EQ(CoordsVec(), bCurve.getControlPoints()); EXPECT_FALSE(bCurve.isRational()); } @@ -41,6 +42,7 @@ TEST(primal_beziercurve, sizing_constructors) { BezierCurveType bCurve(ord); EXPECT_EQ(ord, bCurve.getOrder()); + EXPECT_EQ(ord + 1, bCurve.getNumControlPoints()); EXPECT_EQ(ord + 1, static_cast(bCurve.getControlPoints().size())); EXPECT_FALSE(bCurve.isRational()); } @@ -161,12 +163,14 @@ TEST(primal_beziercurve, set_order) BezierCurveType bCurve; EXPECT_EQ(-1, bCurve.getOrder()); + EXPECT_EQ(0, bCurve.getNumControlPoints()); const int order = 1; PointType controlPoints[2] = {PointType {0.6, 1.2, 1.0}, PointType {0.0, 1.6, 1.8}}; bCurve.setOrder(order); EXPECT_EQ(order, bCurve.getOrder()); + EXPECT_EQ(order - 1, bCurve.getNumControlPoints()); bCurve[0] = controlPoints[0]; bCurve[1] = controlPoints[1]; @@ -178,6 +182,7 @@ TEST(primal_beziercurve, set_order) bCurve.clear(); EXPECT_EQ(-1, bCurve.getOrder()); + EXPECT_EQ(0, bCurve.getNumControlPoints()); EXPECT_FALSE(bCurve.isRational()); } diff --git a/src/axom/primal/tests/primal_winding_number.cpp b/src/axom/primal/tests/primal_winding_number.cpp index c73e3dbd66..08bdc2b174 100644 --- a/src/axom/primal/tests/primal_winding_number.cpp +++ b/src/axom/primal/tests/primal_winding_number.cpp @@ -25,7 +25,7 @@ TEST(primal_winding_number, simple_cases) using Point2D = primal::Point; using Triangle = primal::Triangle; using Bezier = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-8; double edge_tol = 1e-8; @@ -392,7 +392,7 @@ TEST(primal_winding_number, rational_bezier_winding_number) { using Point2D = primal::Point; using Bezier = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-8; double edge_tol = 0; From fd219080c48f1eadcc1d7f62ae22081cd062bf63 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 24 Sep 2025 11:09:22 -0700 Subject: [PATCH 04/21] Improve main test --- .../primal/tests/primal_curved_polygon.cpp | 365 +++++++++--------- 1 file changed, 184 insertions(+), 181 deletions(-) diff --git a/src/axom/primal/tests/primal_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp index ae17502096..0091d9b6f6 100644 --- a/src/axom/primal/tests/primal_curved_polygon.cpp +++ b/src/axom/primal/tests/primal_curved_polygon.cpp @@ -12,14 +12,7 @@ #include "axom/config.hpp" #include "axom/slic.hpp" -#include "axom/primal/geometry/Point.hpp" -#include "axom/primal/geometry/Segment.hpp" -#include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/primal/geometry/OrientationResult.hpp" -#include "axom/primal/operators/intersect.hpp" -#include "axom/primal/operators/compute_moments.hpp" -#include "axom/primal/operators/in_curved_polygon.hpp" -#include "axom/primal/operators/orientation.hpp" +#include "axom/primal.hpp" #include #include @@ -27,7 +20,7 @@ namespace primal = axom::primal; /*! - * Helper function to compute the area and centroid of a curved polygon and to check that they match expectations, + * Helper function to compute the area and centroid of a Bezier CurvedPolygon and to check that they match expectations, * stored in \a expArea and \a expCentroid. Areas and Moments are computed within tolerance \a eps and checks use \a test_eps. */ template @@ -88,160 +81,140 @@ primal::CurvedPolygon> createBezierPolygon( return bPolygon; } -//------------------------------------------------------------------------------ -TEST(primal_curvedpolygon, constructor) +///! Helper function to convert a Bezier curved polygon into a NURBS curved polygon +template +primal::CurvedPolygon> createNURBSPolygon( + const primal::CurvedPolygon>& bPoly) { - const int DIM = 3; - using CoordType = double; - using BezierCurveType = primal::BezierCurve; - using CurvedPolygonType = primal::CurvedPolygon; + using NURBSCurveType = primal::NURBSCurve; + using NURBSPolygonType = primal::CurvedPolygon; - { - SLIC_INFO("Testing default CurvedPolygon constructor "); - CurvedPolygonType bPolygon; - - const int expNumEdges = 0; - EXPECT_EQ(expNumEdges, bPolygon.numEdges()); - EXPECT_EQ(expNumEdges, bPolygon.getEdges().size()); - EXPECT_EQ(axom::Array(), bPolygon.getEdges()); - } + NURBSPolygonType nPoly; + for(int i = 0; i < bPoly.numEdges(); ++i) { - SLIC_INFO("Testing CurvedPolygon numEdges constructor "); - - CurvedPolygonType bPolygon(1); - const int expNumEdges = 1; - EXPECT_EQ(expNumEdges, bPolygon.numEdges()); - EXPECT_EQ(expNumEdges, static_cast(bPolygon.getEdges().size())); + nPoly.addEdge(NURBSCurveType(bPoly[i])); } + + return nPoly; } -//---------------------------------------------------------------------------------- -TEST(primal_curvedpolygon, add_edges) +///! Helper function to convert a Bezier curved polygon into a curved polygon of GWN caches +template +primal::CurvedPolygon> createNURBSCachePolygon( + const primal::CurvedPolygon>& bPoly) { - const int DIM = 2; - using CoordType = double; - using BezierCurveType = primal::BezierCurve; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; + using NURBSCacheType = primal::detail::NURBSCurveGWNCache; + using NURBSPolygonType = primal::CurvedPolygon; - SLIC_INFO("Test adding edges to empty CurvedPolygon"); + NURBSPolygonType ncPoly; - CurvedPolygonType bPolygon; - EXPECT_EQ(0, bPolygon.numEdges()); + for(int i = 0; i < bPoly.numEdges(); ++i) + { + ncPoly.addEdge(NURBSCacheType(bPoly[i])); + } - PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; + return ncPoly; +} - BezierCurveType bCurve(controlPoints, 1); +template +void test_constructor() +{ + using CurvedPolygon = primal::CurvedPolygon; - bPolygon.addEdge(bCurve); - bPolygon.addEdge(bCurve); + { + CurvedPolygon cPolygon; + const int expNumEdges = 0; + EXPECT_EQ(expNumEdges, cPolygon.numEdges()); + EXPECT_EQ(expNumEdges, cPolygon.getEdges().size()); + EXPECT_EQ(axom::Array(), cPolygon.getEdges()); + } - EXPECT_EQ(2, bPolygon.numEdges()); - for(int p = 0; p < bPolygon.numEdges(); ++p) { - BezierCurveType& bc = bPolygon[p]; - for(int sz = 0; sz <= bc.getOrder(); ++sz) - { - auto& pt = bc[sz]; - for(int i = 0; i < DIM; ++i) - { - EXPECT_DOUBLE_EQ(controlPoints[sz][i], pt[i]); - } - } + CurvedPolygon cPolygon(1); + const int expNumEdges = 1; + EXPECT_EQ(expNumEdges, cPolygon.numEdges()); + EXPECT_EQ(expNumEdges, static_cast(cPolygon.getEdges().size())); } } -//---------------------------------------------------------------------------------- -TEST(primal_curvedpolygon, add_edges_nurbs) +//------------------------------------------------------------------------------ +TEST(primal_curvedpolygon, constructor) { - const int DIM = 2; - using CoordType = double; - using CurveType = primal::NURBSCurve; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; - - SLIC_INFO("Test adding edges to empty CurvedPolygon"); - - CurvedPolygonType nPolygon; - EXPECT_EQ(0, nPolygon.numEdges()); + using BezierCurveType = primal::BezierCurve; + using NURBSCurveType = primal::NURBSCurve; + using NURBSCacheType = primal::detail::NURBSCurveGWNCache; + + SLIC_INFO("Testing default and numEdges CurvedPolygon constructors"); + test_constructor(); + test_constructor(); + test_constructor(); +} - PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; +template +void test_add_edges(primal::Point* controlPoints, const CurveType& added_curve) +{ + using CurvedPolygon = primal::CurvedPolygon; - CurveType nCurve(controlPoints, 2, 1); + CurvedPolygon cPolygon; + EXPECT_EQ(0, cPolygon.numEdges()); - nPolygon.addEdge(nCurve); - nPolygon.addEdge(nCurve); + cPolygon.addEdge(added_curve); + cPolygon.addEdge(added_curve); - EXPECT_EQ(2, nPolygon.numEdges()); - for(int p = 0; p < nPolygon.numEdges(); ++p) + EXPECT_EQ(2, cPolygon.numEdges()); + for(int p = 0; p < cPolygon.numEdges(); ++p) { - CurveType& nc = nPolygon[p]; - for(int sz = 0; sz <= nc.getDegree(); ++sz) - { - auto& pt = nc[sz]; - for(int i = 0; i < DIM; ++i) - { - EXPECT_DOUBLE_EQ(controlPoints[sz][i], pt[i]); - } - } + CurveType& curv = cPolygon[p]; + auto& init_pt = curv.getInitPoint(); + EXPECT_DOUBLE_EQ(controlPoints[0][0], init_pt[0]); + EXPECT_DOUBLE_EQ(controlPoints[0][1], init_pt[1]); + + auto& end_pt = curv.getEndPoint(); + EXPECT_DOUBLE_EQ(controlPoints[1][0], end_pt[0]); + EXPECT_DOUBLE_EQ(controlPoints[1][1], end_pt[1]); } } //---------------------------------------------------------------------------------- -TEST(primal_curvedpolygon, add_edges_nurbs_cache) +TEST(primal_curvedpolygon, add_edges) { - const int DIM = 2; - using CoordType = double; - using CurveType = primal::NURBSCurveGWNCache; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; - - SLIC_INFO("Test adding edges to empty CurvedPolygon"); + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using NURBSCurveType = primal::NURBSCurve; + using NURBSCacheType = primal::detail::NURBSCurveGWNCache; - CurvedPolygonType nPolygon; - EXPECT_EQ(0, nPolygon.numEdges()); + SLIC_INFO("Test adding empty edges to empty CurvedPolygons"); PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; - CurveType nCurve(primal::NURBSCurve(controlPoints, 2, 1)); - - nPolygon.addEdge(nCurve); - nPolygon.addEdge(nCurve); - - EXPECT_EQ(2, nPolygon.numEdges()); - for(int p = 0; p < nPolygon.numEdges(); ++p) - { - CurveType& nc = nPolygon[p]; - auto& init_pt = nc.getInitPoint(); - for(int i = 0; i < DIM; ++i) - { - EXPECT_DOUBLE_EQ(controlPoints[0][i], init_pt[i]); - } - - auto& end_pt = nc.getEndPoint(); - for(int i = 0; i < DIM; ++i) - { - EXPECT_DOUBLE_EQ(controlPoints[1][i], end_pt[i]); - } - } + test_add_edges(controlPoints, BezierCurveType(controlPoints, 1)); + test_add_edges(controlPoints, NURBSCurveType(controlPoints, 2, 1)); + test_add_edges(controlPoints, NURBSCacheType(BezierCurveType(controlPoints, 1))); } //---------------------------------------------------------------------------------- TEST(primal_curvedpolygon, isClosed) { - const int DIM = 2; - using CoordType = double; - using CurveType = primal::BezierCurve; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; + using PointType = primal::Point; + using BezierPolygon = primal::CurvedPolygon>; + using NURBSPolygon = primal::CurvedPolygon>; + using NURBSCachePolygon = primal::CurvedPolygon>; SLIC_INFO("Test checking if CurvedPolygon is closed."); { - CurvedPolygonType bPolygon; + BezierPolygon bPolygon; EXPECT_EQ(0, bPolygon.numEdges()); EXPECT_FALSE(bPolygon.isClosed()); + + NURBSPolygon nPolygon; + EXPECT_EQ(0, nPolygon.numEdges()); + EXPECT_FALSE(nPolygon.isClosed()); + + NURBSCachePolygon ncPolygon; + EXPECT_EQ(0, ncPolygon.numEdges()); + EXPECT_FALSE(ncPolygon.isClosed()); } axom::Array CP = {PointType {0.6, 1.2}, @@ -253,65 +226,110 @@ TEST(primal_curvedpolygon, isClosed) { axom::Array subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; axom::Array suborders = {1}; - CurvedPolygonType subPolygon = createBezierPolygon(subCP, suborders); - EXPECT_FALSE(subPolygon.isClosed()); + + BezierPolygon subBezierPolygon = createBezierPolygon(subCP, suborders); + EXPECT_FALSE(subBezierPolygon.isClosed()); + EXPECT_FALSE(createNURBSPolygon(subBezierPolygon).isClosed()); + EXPECT_FALSE(createNURBSCachePolygon(subBezierPolygon).isClosed()); } { - CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); + BezierPolygon bPolygon = createBezierPolygon(CP, orders); EXPECT_EQ(3, bPolygon.numEdges()); EXPECT_TRUE(bPolygon.isClosed()); + NURBSPolygon nPolygon = createNURBSPolygon(bPolygon); + EXPECT_EQ(3, nPolygon.numEdges()); + EXPECT_TRUE(nPolygon.isClosed()); + + NURBSCachePolygon ncPolygon = createNURBSCachePolygon(bPolygon); + EXPECT_EQ(3, ncPolygon.numEdges()); + EXPECT_TRUE(ncPolygon.isClosed()); + + // Manipulate one vertex bPolygon[2][1][0] -= 2e-15; EXPECT_FALSE(bPolygon.isClosed(1e-15)); - } - { - CurvedPolygonType bPolygon = createBezierPolygon(CP, orders); + nPolygon[2][1][0] -= 2e-15; + EXPECT_FALSE(nPolygon.isClosed(1e-15)); - bPolygon[1][0][0] = 5; - EXPECT_FALSE(bPolygon.isClosed(1e-15)); + // Note: Can't access the vertices of NURBSCurveGWNCache objects directly + //nPolygon[2][1][0] -= 2e-15; // Not allowed + ncPolygon = createNURBSCachePolygon(bPolygon); + EXPECT_FALSE(ncPolygon.isClosed(1e-15)); } } //---------------------------------------------------------------------------------- -//TEST(primal_curvedpolygon, isClosed_BiGon) -//{ -// const int DIM = 2; -// using CoordType = double; -// using CurvedPolygonType = primal::CurvedPolygon; -// using PointType = primal::Point; -// -// SLIC_INFO("Test checking if CurvedPolygon is closed for a Bi-Gon"); -// -// CurvedPolygonType bPolygon; -// EXPECT_EQ(0, bPolygon.numEdges()); -// EXPECT_FALSE(bPolygon.isClosed()); -// -// // Bi-gon defined by a quadratic edge and a straight line -// axom::Array CP = {PointType {0.8, .25}, -// PointType {2.0, .50}, -// PointType {0.8, .75}, -// PointType {0.8, .25}}; -// axom::Array orders = {2, 1}; -// -// CurvedPolygonType poly = createBezierPolygon(CP, orders); -// EXPECT_TRUE(poly.isClosed()); -// -// // modify a vertex of the quadratic and check again -// CurvedPolygonType poly2 = poly; -// poly2[0][2] = PointType {0.8, 1.0}; -// EXPECT_FALSE(poly2.isClosed()); -//} +TEST(primal_curvedpolygon, isClosed_BiGon) +{ + using PointType = primal::Point; + using BezierPolygon = primal::CurvedPolygon>; + using NURBSPolygon = primal::CurvedPolygon>; + using NURBSCachePolygon = primal::CurvedPolygon>; + + SLIC_INFO("Test checking if CurvedPolygon is closed for a Bi-Gon"); + + BezierPolygon bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_FALSE(bPolygon.isClosed()); + + // Bi-gon defined by a quadratic edge and a straight line + axom::Array CP = {PointType {0.8, .25}, + PointType {2.0, .50}, + PointType {0.8, .75}, + PointType {0.8, .25}}; + axom::Array orders = {2, 1}; + + BezierPolygon bPoly = createBezierPolygon(CP, orders); + EXPECT_TRUE(bPoly.isClosed()); + EXPECT_TRUE(createNURBSPolygon(bPoly).isClosed()); + EXPECT_TRUE(createNURBSCachePolygon(bPoly).isClosed()); + + // modify a vertex of the quadratic and check again + BezierPolygon bPoly2 = bPoly; + bPoly2[0][2] = PointType {0.8, 1.0}; + + EXPECT_FALSE(bPoly2.isClosed()); + EXPECT_FALSE(createNURBSPolygon(bPoly2).isClosed()); + EXPECT_FALSE(createNURBSCachePolygon(bPoly2).isClosed()); +} + +template +void test_split_edge(primal::CurvedPolygon& curved_polygon) +{ + using PointType = primal::Point; + + // Needs to be std::vector to use .assign + axom::Array subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; + + CurveType nCurve(subCP, 1); + curved_polygon.splitEdge(0, .5); + + CurveType nCurve2; + CurveType nCurve3; + nCurve.split(.5, nCurve2, nCurve3); + + EXPECT_EQ(curved_polygon.numEdges(), 4); + for(int i = 0; i < curved_polygon[0].getNumControlPoints(); ++i) + { + for(int dimi = 0; dimi < 2; ++dimi) + { + EXPECT_EQ(curved_polygon[0][i][dimi], nCurve2[i][dimi]); + EXPECT_EQ(curved_polygon[1][i][dimi], nCurve3[i][dimi]); + } + } +} //---------------------------------------------------------------------------------- TEST(primal_curvedpolygon, split_edge) { - const int DIM = 2; - using CoordType = double; - using CurveType = primal::BezierCurve; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; + using PointType = primal::Point; + using BezierType = primal::BezierCurve; + using BezierPolygon = primal::CurvedPolygon; + + using NURBSType = primal::NURBSCurve; + using NURBSPolygon = primal::CurvedPolygon; SLIC_INFO("Test checking CurvedPolygon edge split."); @@ -321,29 +339,14 @@ TEST(primal_curvedpolygon, split_edge) PointType {0.6, 1.2}}; axom::Array orders32 = {1, 1, 1}; - CurvedPolygonType bPolygon32 = createBezierPolygon(CP, orders32); - - // Needs to be std::vector to use .assign - axom::Array subCP = {PointType {0.6, 1.2}, PointType {0.3, 2.0}}; + BezierPolygon bPolygon32 = createBezierPolygon(CP, orders32); + NURBSPolygon nPolygon32 = createNURBSPolygon(bPolygon32); - CurveType bCurve(subCP, 1); - //std::cout << "Got here!! " << std::endl; - //std::cout << bPolygon32 << std::endl; - bPolygon32.splitEdge(0, .5); + test_split_edge(bPolygon32); + test_split_edge(nPolygon32); - CurveType bCurve2; - CurveType bCurve3; - bCurve.split(.5, bCurve2, bCurve3); - - EXPECT_EQ(bPolygon32.numEdges(), 4); - for(int i = 0; i < bPolygon32[0].getOrder(); ++i) - { - for(int dimi = 0; dimi < DIM; ++dimi) - { - EXPECT_EQ(bPolygon32[0][i][dimi], bCurve2[i][dimi]); - EXPECT_EQ(bPolygon32[1][i][dimi], bCurve3[i][dimi]); - } - } + // Not allowed, since cached objects cannot be changed + //test_split_edge(createNURBSCachePolygon(bPolygon32)); } //---------------------------------------------------------------------------------- From f4190e367f6f472ff1938b06f3e390a638e4e3b2 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 18:39:15 -0700 Subject: [PATCH 05/21] Reorganize some internal methods --- .../detail/evaluate_integral_impl.hpp | 183 +++++++++++++++++- .../primal/operators/evaluate_integral.hpp | 143 ++++---------- src/axom/primal/operators/winding_number.hpp | 6 +- src/axom/primal/tests/primal_integral.cpp | 79 +++++++- 4 files changed, 297 insertions(+), 114 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 312aa70d2f..f4fc929db6 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -11,6 +11,8 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" +#include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" // Use 1D Gauss-Legendre quadrature generated by MFEM if available, // or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) @@ -53,6 +55,62 @@ class QuadratureRuleView #endif }; +/*! + * \brief Evaluate a scalar field line integral on a single NURBS curve. + * + * Decompose the NURBS curve into Bezier segments, then sum the integral on each + * + * \param [in] n The NURBS curve object + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_scalar_line_integral_component(const primal::NURBSCurve& n, + Lambda&& scalar_integrand, + const int npts) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += + detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); + } + return total_integral; +} + +/*! + * \brief Evaluate the scalar integral on a single NURBS curve with cached data for GWN evaluation. + * + * The cache object has already decomposed the NURBS curve into Bezier segments, + * which are used to evaluate the integral over each + * + * \param [in] n The NURBS curve object + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_scalar_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& scalar_integrand, + const int npts) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_scalar_line_integral_component(this_bezier_data.getCurve(), + scalar_integrand, + npts); + } + + return total_integral; +} + /*! * \brief Evaluate a scalar field line integral on a single Bezier curve. * @@ -93,6 +151,60 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< return full_quadrature; } +/*! + * \brief Evaluate a vector field line integral on a single NURBS curve. + * + * Decompose the NURBS curve into Bezier segments, then sum the integral on each + * + * \param [in] n The NURBS curve object + * \param [in] vec_field the lambda function representing the integrand. + * Must accept a 2D point as input and return a 2D axom vector object + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_vector_line_integral_component(const primal::NURBSCurve& n, + Lambda&& vec_field, + const int npts) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += detail::evaluate_vector_line_integral_component(beziers[i], vec_field, npts); + } + return total_integral; +} + +/*! + * \brief Evaluate the vector integral on a single NURBS curve with cached data for GWN evaluation. + * + * The cache object has already decomposed the NURBS curve into Bezier segments, + * which are used to evaluate the integral over each + * + * \param [in] n The NURBS curve object + * \param [in] vec_field the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_vector_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& vec_field, + const int npts) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += + detail::evaluate_vector_line_integral_component(this_bezier_data.getCurve(), vec_field, npts); + } + + return total_integral; +} + /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * @@ -100,7 +212,7 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \param [in] c the Bezier curve object - * \param [in] integrand the lambda function representing the integrand. + * \param [in] vec_field the lambda function representing the integrand. * Must accept a 2D point as input and return a 2D axom vector object * \param [in] npts The number of quadrature points in the rule * \return the value of the integral @@ -133,6 +245,71 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< return full_quadrature; } +/*! + * \brief Evaluate the area integral across one component NURBS of the region. + * + * Intended to be called for each NURBSCurve object bounding a closed 2D region. + * + * \param [in] n The component NURBSCurve object + * \param [in] integrand The lambda function representing the scalar integrand. + * Must accept a 2D point as input and return a double + * \param [in] The lower bound of integration for the antiderivatives + * \param [in] npts_Q The number of quadrature points for the line integral + * \param [in] npts_P The number of quadrature points for the antiderivative + * \return the value of the integral, which is mathematically meaningless. + */ +template +double evaluate_area_integral_component(const primal::NURBSCurve& n, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += + detail::evaluate_area_integral_component(beziers[i], integrand, int_lb, npts_Q, npts_P); + } + return total_integral; +} + +/*! + * \brief Evaluate the area integral across one component NURBSCurveGWNCache of the region. + * + * Intended to be called for each NURBSCurveGWNCache object bounding a closed 2D region. + * + * \param [in] nc The component NURBSCurveGWNCache object + * \param [in] integrand The lambda function representing the scalar integrand. + * Must accept a 2D point as input and return a double + * \param [in] The lower bound of integration for the antiderivatives + * \param [in] npts_Q The number of quadrature points for the line integral + * \param [in] npts_P The number of quadrature points for the antiderivative + * \return the value of the integral, which is mathematically meaningless. + */ +template +inline double evaluate_area_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_area_integral_component(this_bezier_data.getCurve(), + integrand, + int_lb, + npts_Q, + npts_P); + } + + return total_integral; +} + /*! * \brief Evaluate the area integral across one component of the curved polygon. * @@ -142,8 +319,8 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \param [in] cs the array of Bezier curve objects that bound the region - * \param [in] integrand the lambda function representing the integrand. + * \param [in] c The component Bezier curve + * \param [in] integrand The lambda function representing the scalar integrand. * Must accept a 2D point as input and return a double * \param [in] The lower bound of integration for the antiderivatives * \param [in] npts_Q The number of quadrature points for the line integral diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index b2f4fd45be..0d4684a988 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -23,9 +23,6 @@ #include "axom/config.hpp" #include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/primal/geometry/BezierCurve.hpp" -#include "axom/primal/geometry/NURBSCurve.hpp" - #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" // C++ includes @@ -45,6 +42,7 @@ namespace primal * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input and return a double @@ -69,56 +67,26 @@ double evaluate_scalar_line_integral(const primal::CurvedPolygon cpol } /*! - * \brief Evaluate a line integral on a single Bezier curve on a scalar field + * \brief Evaluate a line integral on a single generic curve on a scalar field * - * \param [in] c the Bezier curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] c the generic curve object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const primal::BezierCurve& c, - Lambda&& scalar_integrand, - int npts) +template +double evaluate_scalar_line_integral(const CurveType& c, Lambda&& scalar_integrand, int npts) { return detail::evaluate_scalar_line_integral_component(c, scalar_integrand, npts); } -/*! - * \brief Evaluate a line integral on a single NURBS curve on a scalar field - * - * \param [in] n the NURBS curve object - * \param [in] scalar_integrand the lambda function representing the integrand. - * Must accept a Point as input, and return a double. - * \param [in] npts the number of quadrature nodes per knot span - * - * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts on each segment. - * - * \return the value of the integral - */ -template -double evaluate_scalar_line_integral(const primal::NURBSCurve& n, - Lambda&& scalar_integrand, - int npts) -{ - // Compute the line integral along each component. - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); i++) - { - total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); - } - - return total_integral; -} - /*! * \brief Evaluate a line integral on an array of NURBS curves on a scalar field * - * \param [in] narray the array of NURBS curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] carray The array of generic curve objects * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -128,15 +96,15 @@ double evaluate_scalar_line_integral(const primal::NURBSCurve& n, * * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const axom::Array>& narray, +template +double evaluate_scalar_line_integral(const axom::Array& carray, Lambda&& scalar_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_scalar_line_integral(narray[i], scalar_integrand, npts); + total_integral += evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); } return total_integral; @@ -152,6 +120,7 @@ double evaluate_scalar_line_integral(const axom::Array ~100) * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector @@ -176,56 +145,26 @@ double evaluate_vector_line_integral(const primal::CurvedPolygon cpol } /*! - * \brief Evaluate a line integral on a single Bezier curve on a vector field + * \brief Evaluate a line integral on a single generic curve on a vector field * - * \param [in] c the Bezier curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] c the generic curve object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a Vector. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_vector_line_integral(const primal::BezierCurve& c, - Lambda&& vector_integrand, - int npts) +template +double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integrand, int npts) { return detail::evaluate_vector_line_integral_component(c, vector_integrand, npts); } /*! - * \brief Evaluate a line integral on a single NURBS curve on a vector field + * \brief Evaluate a line integral on an array of generic curves on a vector field * - * \param [in] n the NURBS curve object - * \param [in] vector_integrand the lambda function representing the integrand. - * Must accept a Point as input, and return a Vector. - * \param [in] npts the number of quadrature nodes per knot span - * - * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts on each segment - * - * \return the value of the integral - */ -template -double evaluate_vector_line_integral(const primal::NURBSCurve& n, - Lambda&& vector_integrand, - int npts) -{ - // Compute the line integral along each component. - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); i++) - { - total_integral += - detail::evaluate_vector_line_integral_component(beziers[i], vector_integrand, npts); - } - - return total_integral; -} - -/*! - * \brief Evaluate a line integral on an array of NURBS curves on a vector field - * - * \param [in] narray the array of NURBS curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] carray The array of generic curve objects * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -235,15 +174,15 @@ double evaluate_vector_line_integral(const primal::NURBSCurve& n, * * \return the value of the integral */ -template -double evaluate_vector_line_integral(const axom::Array>& narray, +template +double evaluate_vector_line_integral(const axom::Array& carray, Lambda&& vector_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_vector_line_integral(narray[i], vector_integrand, npts); + total_integral += evaluate_vector_line_integral(carray[i], vector_integrand, npts); } return total_integral; @@ -252,7 +191,12 @@ double evaluate_vector_line_integral(const axom::Array cpoly, } /*! - * \brief Evaluate an integral on the interior of an array of NURBS curves. + * \brief Evaluate an integral on the interior of a region bound by 2D curves * * See above definition for details. * - * \param [in] narray the array of NURBS curve objects that bound the region + * \param [in] carray the array of generic curve objects that bound the region * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] npts_Q the number of quadrature points to evaluate the line integral * \param [in] npts_P the number of quadrature points to evaluate the antiderivative - * - * \note Each NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts_Q * npts_P on each segment * * \note The numerical result is only meaningful if the curves enclose a region * * \return the value of the integral */ -template -double evaluate_area_integral(const axom::Array> narray, +template +double evaluate_area_integral(const axom::Array carray, Lambda&& integrand, int npts_Q, int npts_P = 0) @@ -322,26 +263,26 @@ double evaluate_area_integral(const axom::Array> narray npts_P = npts_Q; } - if(narray.empty()) + if(carray.empty()) { return 0.0; } // Use minimum y-coord of control nodes as lower bound for integration - double int_lb = narray[0][0][1]; - for(int i = 0; i < narray.size(); i++) + double int_lb = carray[0][0][1]; + for(int i = 0; i < carray.size(); i++) { - for(int j = 1; j < narray[i].getNumControlPoints(); j++) + for(int j = 1; j < carray[i].getNumControlPoints(); j++) { - int_lb = std::min(int_lb, narray[i][j][1]); + int_lb = std::min(int_lb, carray[i][j][1]); } } // Evaluate the antiderivative line integral along each component double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - auto beziers = narray[i].extractBezier(); + auto beziers = carray[i].extractBezier(); for(int j = 0; j < beziers.size(); j++) { diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index 2ba00953da..ee1842ac93 100644 --- a/src/axom/primal/operators/winding_number.hpp +++ b/src/axom/primal/operators/winding_number.hpp @@ -170,6 +170,7 @@ double winding_number(const Point& q, /*! * \brief Computes the GWN for a 2D point wrt to a 2D curved polygon * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] q The query point to test * \param [in] cpoly The CurvedPolygon object * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable @@ -386,7 +387,7 @@ axom::Array winding_number(const axom::Array>& query_arr, /*! * \brief Computes the GWN for an array of 2D points wrt an array of generic 2D curves * - * \tparam CurveType The BezierCurve or NURBSCurve which represents the curve + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] query_arr The array of query point to test * \param [in] curve_arr The array of curve objects * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable @@ -416,7 +417,8 @@ axom::Array winding_number(const axom::Array>& query_arr, /*! * \brief Computes the GWN for an array of 2D points wrt to a 2D curved polygon * - * \param [in] query_arr The query point to test + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] q_arr The query point to test * \param [in] cpoly The CurvedPolygon object * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 7548d74d35..c612320d33 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -18,7 +18,7 @@ TEST(primal_integral, evaluate_area_integral) { using Point2D = primal::Point; using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests @@ -92,7 +92,7 @@ TEST(primal_integral, evaluate_line_integral_scalar) { using Point2D = primal::Point; using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests @@ -167,7 +167,7 @@ TEST(primal_integral, evaluate_line_integral_vector) using Point2D = primal::Point; using Vector2D = primal::Vector; using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests @@ -269,7 +269,7 @@ TEST(primal_integral, evaluate_integral_rational) using Point2D = primal::Point; using Vector2D = primal::Vector; using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests @@ -325,6 +325,7 @@ TEST(primal_integral, evaluate_integral_nurbs) using Point2D = primal::Point; using Vector2D = primal::Vector; using NCurve = primal::NURBSCurve; + using CPolygon = primal::CurvedPolygon; double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests @@ -348,10 +349,10 @@ TEST(primal_integral, evaluate_integral_nurbs) NCurve leg2(leg2_nodes, 2, 1); leg2.insertKnot(0.6, 1); - axom::Array quarter_ellipse; - quarter_ellipse.push_back(ellipse_arc); - quarter_ellipse.push_back(leg1); - quarter_ellipse.push_back(leg2); + CPolygon quarter_ellipse; + quarter_ellipse.addEdge(ellipse_arc); + quarter_ellipse.addEdge(leg1); + quarter_ellipse.addEdge(leg2); auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; @@ -445,6 +446,68 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) } } +TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) +{ + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using NCurve = primal::NURBSCurve; + using NCache = primal::NURBSCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Test integrals with same integrand and curves as `evaluate_integral_rational`, + // but with curves added to detail::NURBSCurveGWNCache objects to ensure template compatibility, + // even if there isn't compelling reason to use GWN caches for this purpose + + // Elliptical arc shape + Point2D ellipse_nodes[] = {Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0}}; + double weights[] = {2.0, 1.0, 1.0}; + NCurve ellipse_arc(ellipse_nodes, weights, 3, 2); + ellipse_arc.insertKnot(0.3, 2); + ellipse_arc.insertKnot(0.7, 1); + + Point2D leg1_nodes[] = {Point2D {0.0, 1.0}, {0.0, 0.0}}; + NCurve leg1(leg1_nodes, 2, 1); + leg1.insertKnot(0.4, 1); + + Point2D leg2_nodes[] = {Point2D {0.0, 0.0}, {2.0, 0.0}}; + NCurve leg2(leg2_nodes, 2, 1); + leg2.insertKnot(0.6, 1); + + CPolygon quarter_ellipse; + quarter_ellipse.addEdge(NCache(ellipse_arc)); + quarter_ellipse.addEdge(NCache(leg1)); + quarter_ellipse.addEdge(NCache(leg2)); + + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + + auto area_field = [](Point2D x) -> Vector2D { return Vector2D({-0.5 * x[1], 0.5 * x[0]}); }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); + + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), + 2.42211205514, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), + 1.38837959326, + abs_tol); + + EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); +} + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); From 9147af9df79c73e42e688bfa99a2d64e254082e5 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 18:39:27 -0700 Subject: [PATCH 06/21] Add templated copy constructor --- src/axom/primal/geometry/CurvedPolygon.hpp | 22 ++++++- .../primal/tests/primal_curved_polygon.cpp | 58 ++++--------------- 2 files changed, 31 insertions(+), 49 deletions(-) diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 42517b7417..979793e31e 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -113,6 +113,20 @@ class CurvedPolygon m_edges.resize(nEdges); } + /*! + * \brief Copy constructor for another curve type + */ + template + CurvedPolygon(const CurvedPolygon& other_poly) + { + m_edges.resize(other_poly.numEdges()); + for(int e = 0; e < other_poly.numEdges(); ++e) + { + //this->addEdge(CurveType(other_poly[e])); + m_edges[e] = CurveType(other_poly[e]); + } + } + /// Constructor from an array of \a nEdges curves CurvedPolygon(CurveType* curves, int nEdges) { @@ -145,9 +159,12 @@ class CurvedPolygon m_edges.resize(ngon); } - /// Appends a BezierCurve to the list of edges + /// Appends a curve to the list of edges void addEdge(const CurveType& c1) { m_edges.push_back(c1); } + /// Consumes then appends a curve to the list of edges + void addEdge(CurveType&& c1) { m_edges.push_back(std::move(c1)); } + /// Splits an edge "in place" void splitEdge(int idx, T t) { @@ -155,6 +172,7 @@ class CurvedPolygon AXOM_STATIC_ASSERT_MSG(!has_cached_data::value, "splitEdge cannot be called on objects with associated cache data"); + m_edges.reserve(numEdges() + 1); m_edges.insert(m_edges.begin() + idx + 1, 1, m_edges[idx]); auto& csplit = m_edges[idx]; csplit.split(t, m_edges[idx], m_edges[idx + 1]); @@ -165,7 +183,7 @@ class CurvedPolygon /*! Retrieves the curve at index idx */ CurveType& operator[](int idx) { return m_edges[idx]; } - /*! Retrieves the vertex at index idx */ + /*! Retrieves the curve at index idx */ const CurveType& operator[](int idx) const { return m_edges[idx]; } /// Tests equality of two CurvedPolygons diff --git a/src/axom/primal/tests/primal_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp index 0091d9b6f6..984074cb94 100644 --- a/src/axom/primal/tests/primal_curved_polygon.cpp +++ b/src/axom/primal/tests/primal_curved_polygon.cpp @@ -81,42 +81,6 @@ primal::CurvedPolygon> createBezierPolygon( return bPolygon; } -///! Helper function to convert a Bezier curved polygon into a NURBS curved polygon -template -primal::CurvedPolygon> createNURBSPolygon( - const primal::CurvedPolygon>& bPoly) -{ - using NURBSCurveType = primal::NURBSCurve; - using NURBSPolygonType = primal::CurvedPolygon; - - NURBSPolygonType nPoly; - - for(int i = 0; i < bPoly.numEdges(); ++i) - { - nPoly.addEdge(NURBSCurveType(bPoly[i])); - } - - return nPoly; -} - -///! Helper function to convert a Bezier curved polygon into a curved polygon of GWN caches -template -primal::CurvedPolygon> createNURBSCachePolygon( - const primal::CurvedPolygon>& bPoly) -{ - using NURBSCacheType = primal::detail::NURBSCurveGWNCache; - using NURBSPolygonType = primal::CurvedPolygon; - - NURBSPolygonType ncPoly; - - for(int i = 0; i < bPoly.numEdges(); ++i) - { - ncPoly.addEdge(NURBSCacheType(bPoly[i])); - } - - return ncPoly; -} - template void test_constructor() { @@ -229,8 +193,8 @@ TEST(primal_curvedpolygon, isClosed) BezierPolygon subBezierPolygon = createBezierPolygon(subCP, suborders); EXPECT_FALSE(subBezierPolygon.isClosed()); - EXPECT_FALSE(createNURBSPolygon(subBezierPolygon).isClosed()); - EXPECT_FALSE(createNURBSCachePolygon(subBezierPolygon).isClosed()); + EXPECT_FALSE(NURBSPolygon(subBezierPolygon).isClosed()); + EXPECT_FALSE(NURBSCachePolygon(subBezierPolygon).isClosed()); } { @@ -238,11 +202,11 @@ TEST(primal_curvedpolygon, isClosed) EXPECT_EQ(3, bPolygon.numEdges()); EXPECT_TRUE(bPolygon.isClosed()); - NURBSPolygon nPolygon = createNURBSPolygon(bPolygon); + NURBSPolygon nPolygon = NURBSPolygon(bPolygon); EXPECT_EQ(3, nPolygon.numEdges()); EXPECT_TRUE(nPolygon.isClosed()); - NURBSCachePolygon ncPolygon = createNURBSCachePolygon(bPolygon); + NURBSCachePolygon ncPolygon = NURBSCachePolygon(bPolygon); EXPECT_EQ(3, ncPolygon.numEdges()); EXPECT_TRUE(ncPolygon.isClosed()); @@ -255,7 +219,7 @@ TEST(primal_curvedpolygon, isClosed) // Note: Can't access the vertices of NURBSCurveGWNCache objects directly //nPolygon[2][1][0] -= 2e-15; // Not allowed - ncPolygon = createNURBSCachePolygon(bPolygon); + ncPolygon = NURBSCachePolygon(bPolygon); EXPECT_FALSE(ncPolygon.isClosed(1e-15)); } } @@ -283,16 +247,16 @@ TEST(primal_curvedpolygon, isClosed_BiGon) BezierPolygon bPoly = createBezierPolygon(CP, orders); EXPECT_TRUE(bPoly.isClosed()); - EXPECT_TRUE(createNURBSPolygon(bPoly).isClosed()); - EXPECT_TRUE(createNURBSCachePolygon(bPoly).isClosed()); + EXPECT_TRUE(NURBSPolygon(bPoly).isClosed()); + EXPECT_TRUE(NURBSCachePolygon(bPoly).isClosed()); // modify a vertex of the quadratic and check again BezierPolygon bPoly2 = bPoly; bPoly2[0][2] = PointType {0.8, 1.0}; EXPECT_FALSE(bPoly2.isClosed()); - EXPECT_FALSE(createNURBSPolygon(bPoly2).isClosed()); - EXPECT_FALSE(createNURBSCachePolygon(bPoly2).isClosed()); + EXPECT_FALSE(NURBSPolygon(bPoly2).isClosed()); + EXPECT_FALSE(NURBSCachePolygon(bPoly2).isClosed()); } template @@ -340,11 +304,11 @@ TEST(primal_curvedpolygon, split_edge) axom::Array orders32 = {1, 1, 1}; BezierPolygon bPolygon32 = createBezierPolygon(CP, orders32); - NURBSPolygon nPolygon32 = createNURBSPolygon(bPolygon32); + NURBSPolygon nPolygon32(bPolygon32); test_split_edge(bPolygon32); test_split_edge(nPolygon32); - + // Not allowed, since cached objects cannot be changed //test_split_edge(createNURBSCachePolygon(bPolygon32)); } From 80ec9b2c48788c07b5fe571b6195005c3374e311 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 18:39:55 -0700 Subject: [PATCH 07/21] More specific interface for memoized sampling method --- src/axom/quest/SamplingShaper.hpp | 10 ++++++++-- .../quest/detail/shaping/WindingNumberSampler.hpp | 12 ++++++++++-- src/axom/quest/io/MFEMReader.cpp | 2 +- src/axom/quest/io/MFEMReader.hpp | 2 +- 4 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index 6c2f52ece1..772b1ffc73 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -210,8 +210,14 @@ class SamplingShaper : public Shaper AXOM_UNUSED_VAR(shapeDimension); if(useWindingNumberSampler(shape)) { + using CachedContourType = + typename axom::primal::CurvedPolygon>; + axom::Array contours_cached; + + for(auto& contour : m_contours) contours_cached.push_back(CachedContourType(contour)); + m_inoutSamplerWN = - std::make_unique>(shapeName, m_contours.view()); + std::make_unique>(shapeName, contours_cached.view()); m_inoutSamplerWN->computeBounds(); m_inoutSamplerWN->initSpatialIndex(this->m_vertexWeldThreshold); } @@ -961,7 +967,7 @@ class SamplingShaper : public Shaper std::unique_ptr> m_primitiveSampler3D_hip; std::unique_ptr> m_inoutSamplerWN; - axom::Array>> m_contours; + axom::Array>> m_contours; std::set m_knownMaterials; diff --git a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp index a6c63dab32..7da9da925c 100644 --- a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp +++ b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp @@ -63,8 +63,8 @@ class WindingNumberSampler // For now. using ExecSpace = axom::SEQ_EXEC; - using CurvedPolygonType = primal::CurvedPolygon>; - using GeometryView = typename axom::ArrayView; + using ContourType = primal::CurvedPolygon>; + using GeometryView = typename axom::ArrayView; using PointType = primal::Point; using GeometricBoundingBox = axom::primal::BoundingBox; using BVH = typename axom::spin::BVH; @@ -90,6 +90,14 @@ class WindingNumberSampler // no-op - We do it in initSpatialIndex. } + /*! + * \brief Initialize the NURBSCurveGWNCache objects that store intermediate calculations for GWN evaluation + */ + void initGWNCaches() + { + + } + /*! * \brief Initialize the BVH that is used for queries from the geometry view. */ diff --git a/src/axom/quest/io/MFEMReader.cpp b/src/axom/quest/io/MFEMReader.cpp index 2a03ca6d53..eea18d0a19 100644 --- a/src/axom/quest/io/MFEMReader.cpp +++ b/src/axom/quest/io/MFEMReader.cpp @@ -128,7 +128,7 @@ int MFEMReader::read(CurvedPolygonArray &curvedPolygons) for(int zoneId : zoneIds) { auto curve = axom::quest::internal::segment_to_nurbs(mesh, zoneId); - poly.addEdge(std::move(axom::primal::NURBSCurveGWNCache(curve))); + poly.addEdge(std::move(curve)); } } }); diff --git a/src/axom/quest/io/MFEMReader.hpp b/src/axom/quest/io/MFEMReader.hpp index 39b45fa8db..49e0496b72 100644 --- a/src/axom/quest/io/MFEMReader.hpp +++ b/src/axom/quest/io/MFEMReader.hpp @@ -35,7 +35,7 @@ class MFEMReader using NURBSCurve = axom::primal::NURBSCurve; using CurveArray = axom::Array; - using CurvedPolygon = axom::primal::CurvedPolygon; + using CurvedPolygon = axom::primal::CurvedPolygon; using CurvedPolygonArray = axom::Array; public: From 9daf5164471c11d2dafee28e87693d128f19a777 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 12:58:55 -0700 Subject: [PATCH 08/21] Fix minor issues --- src/axom/quest/detail/shaping/WindingNumberSampler.hpp | 2 +- src/axom/quest/tests/quest_mfem_reader.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp index 7da9da925c..8795d7cf7a 100644 --- a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp +++ b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp @@ -63,7 +63,7 @@ class WindingNumberSampler // For now. using ExecSpace = axom::SEQ_EXEC; - using ContourType = primal::CurvedPolygon>; + using ContourType = primal::CurvedPolygon>; using GeometryView = typename axom::ArrayView; using PointType = primal::Point; using GeometricBoundingBox = axom::primal::BoundingBox; diff --git a/src/axom/quest/tests/quest_mfem_reader.cpp b/src/axom/quest/tests/quest_mfem_reader.cpp index cc0290dc8b..95a8a59465 100644 --- a/src/axom/quest/tests/quest_mfem_reader.cpp +++ b/src/axom/quest/tests/quest_mfem_reader.cpp @@ -61,7 +61,7 @@ TEST(quest_mfem_reader, read_curved_polygon) reader.setFileName(fileName); // Read as 1 CurvedPolygon with 9 edges - axom::Array> polys; + axom::Array>> polys; EXPECT_EQ(reader.read(polys), 0); EXPECT_EQ(polys.size(), 1); EXPECT_EQ(polys[0].numEdges(), 9); From fe14367b2067665613b6e93d300106213ad469d3 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 13:14:58 -0700 Subject: [PATCH 09/21] Make endpoint access const references --- src/axom/primal/geometry/BezierCurve.hpp | 4 ++-- src/axom/primal/geometry/NURBSCurve.hpp | 4 ++-- .../primal/operators/detail/winding_number_2d_memoization.hpp | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index ed5716c227..fd31c42f91 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -259,8 +259,8 @@ class BezierCurve /// Retrieves the control point at index \a idx const PointType& operator[](int idx) const { return m_controlPoints[idx]; } - PointType getInitPoint() const { return m_controlPoints[0]; } - PointType getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + const PointType& getInitPoint() const { return m_controlPoints[0]; } + const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } /// Returns a copy of the Bezier curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index 1263364590..865b76e65d 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -621,8 +621,8 @@ class NURBSCurve */ const PointType& operator[](int idx) const { return m_controlPoints[idx]; } - PointType getInitPoint() const { return m_controlPoints[0]; } - PointType getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + const PointType& getInitPoint() const { return m_controlPoints[0]; } + const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } /// \brief Returns a copy of the NURBS curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } diff --git a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp index f76c7a4b4b..3bc2d67cac 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp @@ -194,8 +194,8 @@ class NURBSCurveGWNCache auto getNumControlPoints() const { return m_numControlPoints; } auto getDegree() const { return m_degree; } - auto getInitPoint() const { return m_initPoint; } - auto getEndPoint() const { return m_endPoint; } + const auto& getInitPoint() const { return m_initPoint; } + const auto& getEndPoint() const { return m_endPoint; } //@} friend bool operator==(const NURBSCurveGWNCache& lhs, const NURBSCurveGWNCache& rhs) From 96230483f5c1c348f5b1a5811a36b9aa86a39f91 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 13:15:13 -0700 Subject: [PATCH 10/21] Fix bad call --- src/axom/primal/operators/winding_number.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index ee1842ac93..209d7f7b9b 100644 --- a/src/axom/primal/operators/winding_number.hpp +++ b/src/axom/primal/operators/winding_number.hpp @@ -476,7 +476,7 @@ axom::Array winding_number(const axom::Array>& q_arr, axom::Array ret_val(q_arr.size()); for(int n = 0; n < q_arr.size(); ++n) { - ret_val[n] = winding_number(query_arr[n], cache_arr, edge_tol, EPS); + ret_val[n] = winding_number(q_arr[n], cpoly, edge_tol, EPS); } return ret_val; From 9e312288e17b5c0e8dd3f521910822ab7406c1a6 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 13:17:40 -0700 Subject: [PATCH 11/21] Match other branch --- .../detail/evaluate_integral_impl.hpp | 183 +----------------- 1 file changed, 3 insertions(+), 180 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 460a858334..0bfc1b2eec 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -11,8 +11,6 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" -#include "axom/primal/geometry/NURBSCurve.hpp" -#include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" #include "axom/core/numerics/quadrature.hpp" @@ -26,62 +24,6 @@ namespace primal namespace detail { -/*! - * \brief Evaluate a scalar field line integral on a single NURBS curve. - * - * Decompose the NURBS curve into Bezier segments, then sum the integral on each - * - * \param [in] n The NURBS curve object - * \param [in] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_scalar_line_integral_component(const primal::NURBSCurve& n, - Lambda&& scalar_integrand, - const int npts) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); - } - return total_integral; -} - -/*! - * \brief Evaluate the scalar integral on a single NURBS curve with cached data for GWN evaluation. - * - * The cache object has already decomposed the NURBS curve into Bezier segments, - * which are used to evaluate the integral over each - * - * \param [in] n The NURBS curve object - * \param [in] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_scalar_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& scalar_integrand, - const int npts) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += detail::evaluate_scalar_line_integral_component(this_bezier_data.getCurve(), - scalar_integrand, - npts); - } - - return total_integral; -} - /*! * \brief Evaluate a scalar field line integral on a single Bezier curve. * @@ -115,67 +57,13 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< return full_quadrature; } -/*! - * \brief Evaluate a vector field line integral on a single NURBS curve. - * - * Decompose the NURBS curve into Bezier segments, then sum the integral on each - * - * \param [in] n The NURBS curve object - * \param [in] vec_field the lambda function representing the integrand. - * Must accept a 2D point as input and return a 2D axom vector object - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_vector_line_integral_component(const primal::NURBSCurve& n, - Lambda&& vec_field, - const int npts) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += detail::evaluate_vector_line_integral_component(beziers[i], vec_field, npts); - } - return total_integral; -} - -/*! - * \brief Evaluate the vector integral on a single NURBS curve with cached data for GWN evaluation. - * - * The cache object has already decomposed the NURBS curve into Bezier segments, - * which are used to evaluate the integral over each - * - * \param [in] n The NURBS curve object - * \param [in] vec_field the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_vector_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& vec_field, - const int npts) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += - detail::evaluate_vector_line_integral_component(this_bezier_data.getCurve(), vec_field, npts); - } - - return total_integral; -} - /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * * Evaluate the vector field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object - * \param [in] vec_field the lambda function representing the integrand. + * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a 2D axom vector object * \param [in] npts The number of quadrature points in the rule * \return the value of the integral @@ -202,71 +90,6 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< return full_quadrature; } -/*! - * \brief Evaluate the area integral across one component NURBS of the region. - * - * Intended to be called for each NURBSCurve object bounding a closed 2D region. - * - * \param [in] n The component NURBSCurve object - * \param [in] integrand The lambda function representing the scalar integrand. - * Must accept a 2D point as input and return a double - * \param [in] The lower bound of integration for the antiderivatives - * \param [in] npts_Q The number of quadrature points for the line integral - * \param [in] npts_P The number of quadrature points for the antiderivative - * \return the value of the integral, which is mathematically meaningless. - */ -template -double evaluate_area_integral_component(const primal::NURBSCurve& n, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += - detail::evaluate_area_integral_component(beziers[i], integrand, int_lb, npts_Q, npts_P); - } - return total_integral; -} - -/*! - * \brief Evaluate the area integral across one component NURBSCurveGWNCache of the region. - * - * Intended to be called for each NURBSCurveGWNCache object bounding a closed 2D region. - * - * \param [in] nc The component NURBSCurveGWNCache object - * \param [in] integrand The lambda function representing the scalar integrand. - * Must accept a 2D point as input and return a double - * \param [in] The lower bound of integration for the antiderivatives - * \param [in] npts_Q The number of quadrature points for the line integral - * \param [in] npts_P The number of quadrature points for the antiderivative - * \return the value of the integral, which is mathematically meaningless. - */ -template -inline double evaluate_area_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += detail::evaluate_area_integral_component(this_bezier_data.getCurve(), - integrand, - int_lb, - npts_Q, - npts_P); - } - - return total_integral; -} - /*! * \brief Evaluate the area integral across one component of the curved polygon. * @@ -276,8 +99,8 @@ inline double evaluate_area_integral_component(const primal::detail::NURBSCurveG * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \param [in] c The component Bezier curve - * \param [in] integrand The lambda function representing the scalar integrand. + * \param [in] cs the array of Bezier curve objects that bound the region + * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] The lower bound of integration for the antiderivatives * \param [in] npts_Q The number of quadrature points for the line integral From a5e63ac6868ac41a9eceeea3aa6f5fbc58605070 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 13:23:33 -0700 Subject: [PATCH 12/21] Match other branch, easier than rebasing --- .../detail/evaluate_integral_impl.hpp | 183 +----------------- .../primal/operators/evaluate_integral.hpp | 157 ++++++++++----- 2 files changed, 111 insertions(+), 229 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 460a858334..0bfc1b2eec 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -11,8 +11,6 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" -#include "axom/primal/geometry/NURBSCurve.hpp" -#include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" #include "axom/core/numerics/quadrature.hpp" @@ -26,62 +24,6 @@ namespace primal namespace detail { -/*! - * \brief Evaluate a scalar field line integral on a single NURBS curve. - * - * Decompose the NURBS curve into Bezier segments, then sum the integral on each - * - * \param [in] n The NURBS curve object - * \param [in] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_scalar_line_integral_component(const primal::NURBSCurve& n, - Lambda&& scalar_integrand, - const int npts) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); - } - return total_integral; -} - -/*! - * \brief Evaluate the scalar integral on a single NURBS curve with cached data for GWN evaluation. - * - * The cache object has already decomposed the NURBS curve into Bezier segments, - * which are used to evaluate the integral over each - * - * \param [in] n The NURBS curve object - * \param [in] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_scalar_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& scalar_integrand, - const int npts) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += detail::evaluate_scalar_line_integral_component(this_bezier_data.getCurve(), - scalar_integrand, - npts); - } - - return total_integral; -} - /*! * \brief Evaluate a scalar field line integral on a single Bezier curve. * @@ -115,67 +57,13 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< return full_quadrature; } -/*! - * \brief Evaluate a vector field line integral on a single NURBS curve. - * - * Decompose the NURBS curve into Bezier segments, then sum the integral on each - * - * \param [in] n The NURBS curve object - * \param [in] vec_field the lambda function representing the integrand. - * Must accept a 2D point as input and return a 2D axom vector object - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_vector_line_integral_component(const primal::NURBSCurve& n, - Lambda&& vec_field, - const int npts) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += detail::evaluate_vector_line_integral_component(beziers[i], vec_field, npts); - } - return total_integral; -} - -/*! - * \brief Evaluate the vector integral on a single NURBS curve with cached data for GWN evaluation. - * - * The cache object has already decomposed the NURBS curve into Bezier segments, - * which are used to evaluate the integral over each - * - * \param [in] n The NURBS curve object - * \param [in] vec_field the lambda function representing the integrand. - * Must accept a 2D point as input and return a double - * \param [in] npts The number of quadrature points in the rule - * \return the value of the integral - */ -template -inline double evaluate_vector_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& vec_field, - const int npts) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += - detail::evaluate_vector_line_integral_component(this_bezier_data.getCurve(), vec_field, npts); - } - - return total_integral; -} - /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * * Evaluate the vector field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object - * \param [in] vec_field the lambda function representing the integrand. + * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a 2D axom vector object * \param [in] npts The number of quadrature points in the rule * \return the value of the integral @@ -202,71 +90,6 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< return full_quadrature; } -/*! - * \brief Evaluate the area integral across one component NURBS of the region. - * - * Intended to be called for each NURBSCurve object bounding a closed 2D region. - * - * \param [in] n The component NURBSCurve object - * \param [in] integrand The lambda function representing the scalar integrand. - * Must accept a 2D point as input and return a double - * \param [in] The lower bound of integration for the antiderivatives - * \param [in] npts_Q The number of quadrature points for the line integral - * \param [in] npts_P The number of quadrature points for the antiderivative - * \return the value of the integral, which is mathematically meaningless. - */ -template -double evaluate_area_integral_component(const primal::NURBSCurve& n, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) - { - total_integral += - detail::evaluate_area_integral_component(beziers[i], integrand, int_lb, npts_Q, npts_P); - } - return total_integral; -} - -/*! - * \brief Evaluate the area integral across one component NURBSCurveGWNCache of the region. - * - * Intended to be called for each NURBSCurveGWNCache object bounding a closed 2D region. - * - * \param [in] nc The component NURBSCurveGWNCache object - * \param [in] integrand The lambda function representing the scalar integrand. - * Must accept a 2D point as input and return a double - * \param [in] The lower bound of integration for the antiderivatives - * \param [in] npts_Q The number of quadrature points for the line integral - * \param [in] npts_P The number of quadrature points for the antiderivative - * \return the value of the integral, which is mathematically meaningless. - */ -template -inline double evaluate_area_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - double total_integral = 0.0; - for(int i = 0; i < nc.getNumKnotSpans(); ++i) - { - // Assuming the cache is properly initialized, this operation will never add to the cache - const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); - total_integral += detail::evaluate_area_integral_component(this_bezier_data.getCurve(), - integrand, - int_lb, - npts_Q, - npts_P); - } - - return total_integral; -} - /*! * \brief Evaluate the area integral across one component of the curved polygon. * @@ -276,8 +99,8 @@ inline double evaluate_area_integral_component(const primal::detail::NURBSCurveG * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \param [in] c The component Bezier curve - * \param [in] integrand The lambda function representing the scalar integrand. + * \param [in] cs the array of Bezier curve objects that bound the region + * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] The lower bound of integration for the antiderivatives * \param [in] npts_Q The number of quadrature points for the line integral diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index a718a2caac..507d28ebe4 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -22,6 +22,9 @@ #include "axom/config.hpp" #include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" + #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" // C++ includes @@ -40,7 +43,6 @@ namespace primal * * Evaluate the scalar field line integral with Gauss-Legendre quadrature * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input and return a double @@ -48,8 +50,8 @@ namespace primal * on each edge of the CurvedPolygon * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, Lambda&& scalar_integrand, int npts) { @@ -65,26 +67,56 @@ double evaluate_scalar_line_integral(const primal::CurvedPolygon cpol } /*! - * \brief Evaluate a line integral on a single generic curve on a scalar field + * \brief Evaluate a line integral on a single Bezier curve on a scalar field * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve - * \param [in] c the generic curve object + * \param [in] c the Bezier curve object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const CurveType& c, Lambda&& scalar_integrand, int npts) +template +double evaluate_scalar_line_integral(const primal::BezierCurve& c, + Lambda&& scalar_integrand, + int npts) { return detail::evaluate_scalar_line_integral_component(c, scalar_integrand, npts); } +/*! + * \brief Evaluate a line integral on a single NURBS curve on a scalar field + * + * \param [in] n the NURBS curve object + * \param [in] scalar_integrand the lambda function representing the integrand. + * Must accept a Point as input, and return a double. + * \param [in] npts the number of quadrature nodes per knot span + * + * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature + * is computed using npts on each segment. + * + * \return the value of the integral + */ +template +double evaluate_scalar_line_integral(const primal::NURBSCurve& n, + Lambda&& scalar_integrand, + int npts) +{ + // Compute the line integral along each component. + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); i++) + { + total_integral += + detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); + } + + return total_integral; +} + /*! * \brief Evaluate a line integral on an array of NURBS curves on a scalar field * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve - * \param [in] carray The array of generic curve objects + * \param [in] narray the array of NURBS curve object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -94,15 +126,15 @@ double evaluate_scalar_line_integral(const CurveType& c, Lambda&& scalar_integra * * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const axom::Array& carray, +template +double evaluate_scalar_line_integral(const axom::Array>& narray, Lambda&& scalar_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < carray.size(); i++) + for(int i = 0; i < narray.size(); i++) { - total_integral += evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); + total_integral += evaluate_scalar_line_integral(narray[i], scalar_integrand, npts); } return total_integral; @@ -117,7 +149,6 @@ double evaluate_scalar_line_integral(const axom::Array& carray, * * Evaluate the vector field line integral with Gauss-Legendre quadrature * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector @@ -125,8 +156,8 @@ double evaluate_scalar_line_integral(const axom::Array& carray, * on each edge of the CurvedPolygon * \return the value of the integral */ -template -double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, Lambda&& vector_integrand, int npts) { @@ -142,26 +173,56 @@ double evaluate_vector_line_integral(const primal::CurvedPolygon cpol } /*! - * \brief Evaluate a line integral on a single generic curve on a vector field + * \brief Evaluate a line integral on a single Bezier curve on a vector field * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve - * \param [in] c the generic curve object + * \param [in] c the Bezier curve object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a Vector. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integrand, int npts) +template +double evaluate_vector_line_integral(const primal::BezierCurve& c, + Lambda&& vector_integrand, + int npts) { return detail::evaluate_vector_line_integral_component(c, vector_integrand, npts); } /*! - * \brief Evaluate a line integral on an array of generic curves on a vector field + * \brief Evaluate a line integral on a single NURBS curve on a vector field * - * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve - * \param [in] carray The array of generic curve objects + * \param [in] n the NURBS curve object + * \param [in] vector_integrand the lambda function representing the integrand. + * Must accept a Point as input, and return a Vector. + * \param [in] npts the number of quadrature nodes per knot span + * + * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature + * is computed using npts on each segment + * + * \return the value of the integral + */ +template +double evaluate_vector_line_integral(const primal::NURBSCurve& n, + Lambda&& vector_integrand, + int npts) +{ + // Compute the line integral along each component. + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); i++) + { + total_integral += + detail::evaluate_vector_line_integral_component(beziers[i], vector_integrand, npts); + } + + return total_integral; +} + +/*! + * \brief Evaluate a line integral on an array of NURBS curves on a vector field + * + * \param [in] narray the array of NURBS curve object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -171,15 +232,15 @@ double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integra * * \return the value of the integral */ -template -double evaluate_vector_line_integral(const axom::Array& carray, +template +double evaluate_vector_line_integral(const axom::Array>& narray, Lambda&& vector_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < carray.size(); i++) + for(int i = 0; i < narray.size(); i++) { - total_integral += evaluate_vector_line_integral(carray[i], vector_integrand, npts); + total_integral += evaluate_vector_line_integral(narray[i], vector_integrand, npts); } return total_integral; @@ -188,12 +249,7 @@ double evaluate_vector_line_integral(const axom::Array& carray, /*! * \brief Evaluate an integral on the interior of a CurvedPolygon object. * - * Evaluates the integral using a Spectral Mesh-Free Quadrature derived from - * Green's theorem, evaluating the area integral as a line integral of the - * antiderivative over each component curve. - * - * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar - * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. + * See above definition for details. * * \param [in] cs the array of Bezier curve objects that bound the region * \param [in] integrand the lambda function representing the integrand. @@ -202,8 +258,8 @@ double evaluate_vector_line_integral(const axom::Array& carray, * \param [in] npts_P the number of quadrature points to evaluate the antiderivative * \return the value of the integral */ -template -double evaluate_area_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_area_integral(const primal::CurvedPolygon cpoly, Lambda&& integrand, int npts_Q, int npts_P = 0) @@ -235,22 +291,25 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, } /*! - * \brief Evaluate an integral on the interior of a region bound by 2D curves + * \brief Evaluate an integral on the interior of an array of NURBS curves. * * See above definition for details. * - * \param [in] carray the array of generic curve objects that bound the region + * \param [in] narray the array of NURBS curve objects that bound the region * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] npts_Q the number of quadrature points to evaluate the line integral * \param [in] npts_P the number of quadrature points to evaluate the antiderivative + * + * \note Each NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature + * is computed using npts_Q * npts_P on each segment * * \note The numerical result is only meaningful if the curves enclose a region * * \return the value of the integral */ -template -double evaluate_area_integral(const axom::Array carray, +template +double evaluate_area_integral(const axom::Array> narray, Lambda&& integrand, int npts_Q, int npts_P = 0) @@ -260,26 +319,26 @@ double evaluate_area_integral(const axom::Array carray, npts_P = npts_Q; } - if(carray.empty()) + if(narray.empty()) { return 0.0; } // Use minimum y-coord of control nodes as lower bound for integration - double int_lb = carray[0][0][1]; - for(int i = 0; i < carray.size(); i++) + double int_lb = narray[0][0][1]; + for(int i = 0; i < narray.size(); i++) { - for(int j = 1; j < carray[i].getNumControlPoints(); j++) + for(int j = 1; j < narray[i].getNumControlPoints(); j++) { - int_lb = std::min(int_lb, carray[i][j][1]); + int_lb = std::min(int_lb, narray[i][j][1]); } } // Evaluate the antiderivative line integral along each component double total_integral = 0.0; - for(int i = 0; i < carray.size(); i++) + for(int i = 0; i < narray.size(); i++) { - auto beziers = carray[i].extractBezier(); + auto beziers = narray[i].extractBezier(); for(int j = 0; j < beziers.size(); j++) { @@ -294,4 +353,4 @@ double evaluate_area_integral(const axom::Array carray, } // namespace primal } // end namespace axom -#endif +#endif \ No newline at end of file From 6fd55f60f6420ddbb651692087e568127de0ca0f Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 13:27:03 -0700 Subject: [PATCH 13/21] Fix style --- .../primal/operators/evaluate_integral.hpp | 2 +- src/axom/primal/tests/primal_integral.cpp | 885 +++++++++--------- 2 files changed, 443 insertions(+), 444 deletions(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 507d28ebe4..4e3a1b8ece 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -353,4 +353,4 @@ double evaluate_area_integral(const axom::Array> narray } // namespace primal } // end namespace axom -#endif \ No newline at end of file +#endif diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 0c25c9d348..31409e9725 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -16,384 +16,384 @@ namespace primal = axom::primal; TEST(primal_integral, evaluate_area_integral) { - using Point2D = primal::Point; - using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; + using Point2D = primal::Point; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 20; + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; - // Define anonymous functions for testing - auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; - auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; - auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + // Define anonymous functions for testing + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; - // Test on triangular domain - Point2D trinodes1[] = { Point2D {0.0, 0.0}, Point2D {1.0, 0.0} }; - BCurve tri1(trinodes1, 1); + // Test on triangular domain + Point2D trinodes1[] = {Point2D {0.0, 0.0}, Point2D {1.0, 0.0}}; + BCurve tri1(trinodes1, 1); - Point2D trinodes2[] = { Point2D {1.0, 0.0}, Point2D {0.0, 1.0} }; - BCurve tri2(trinodes2, 1); + Point2D trinodes2[] = {Point2D {1.0, 0.0}, Point2D {0.0, 1.0}}; + BCurve tri2(trinodes2, 1); - Point2D trinodes3[] = { Point2D {0.0, 1.0}, Point2D {0.0, 0.0} }; - BCurve tri3(trinodes3, 1); + Point2D trinodes3[] = {Point2D {0.0, 1.0}, Point2D {0.0, 0.0}}; + BCurve tri3(trinodes3, 1); - BCurve triangle_edges[] = { tri1, tri2, tri3 }; - CPolygon triangle(triangle_edges, 3); + BCurve triangle_edges[] = {tri1, tri2, tri3}; + CPolygon triangle(triangle_edges, 3); - // Compare against hand computed/high-precision calculated values - EXPECT_NEAR(evaluate_area_integral(triangle, const_integrand, npts), 0.5, abs_tol); - EXPECT_NEAR(evaluate_area_integral(triangle, poly_integrand, npts), 1.0 / 60.0, abs_tol); - EXPECT_NEAR(evaluate_area_integral(triangle, transc_integrand, npts), 0.0415181074232, abs_tol); + // Compare against hand computed/high-precision calculated values + EXPECT_NEAR(evaluate_area_integral(triangle, const_integrand, npts), 0.5, abs_tol); + EXPECT_NEAR(evaluate_area_integral(triangle, poly_integrand, npts), 1.0 / 60.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(triangle, transc_integrand, npts), 0.0415181074232, abs_tol); - // Test on parabolic domain (between f(x) = 1-x^2 and g(x) = x^2-1, shifted to the right 1 unit) - Point2D paranodes1[] = { Point2D {2.0, 0.0}, Point2D {1.0, 2.0}, Point2D {0.0, 0.0} }; - BCurve para1(paranodes1, 2); + // Test on parabolic domain (between f(x) = 1-x^2 and g(x) = x^2-1, shifted to the right 1 unit) + Point2D paranodes1[] = {Point2D {2.0, 0.0}, Point2D {1.0, 2.0}, Point2D {0.0, 0.0}}; + BCurve para1(paranodes1, 2); - Point2D paranodes2[] = { Point2D {0.0, 0.0}, Point2D {1.0, -2.0}, Point2D {2.0, 0.0} }; - BCurve para2(paranodes2, 2); + Point2D paranodes2[] = {Point2D {0.0, 0.0}, Point2D {1.0, -2.0}, Point2D {2.0, 0.0}}; + BCurve para2(paranodes2, 2); - BCurve parabola_edges[] = { para1, para2 }; - CPolygon parabola_shape(parabola_edges, 2); + BCurve parabola_edges[] = {para1, para2}; + CPolygon parabola_shape(parabola_edges, 2); - // Compare against hand computed/high-precision calculated values - EXPECT_NEAR(evaluate_area_integral(parabola_shape, const_integrand, npts), 8.0 / 3.0, abs_tol); - EXPECT_NEAR(evaluate_area_integral(parabola_shape, poly_integrand, npts), 64.0 / 105.0, abs_tol); - EXPECT_NEAR(evaluate_area_integral(parabola_shape, transc_integrand, npts), 0.0, abs_tol); + // Compare against hand computed/high-precision calculated values + EXPECT_NEAR(evaluate_area_integral(parabola_shape, const_integrand, npts), 8.0 / 3.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_shape, poly_integrand, npts), 64.0 / 105.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_shape, transc_integrand, npts), 0.0, abs_tol); - // Ensure compatibility with curved polygons - BCurve pedges[2] = { para1, para2 }; - CPolygon parabola_polygon(pedges, 2); - EXPECT_NEAR(evaluate_area_integral(parabola_polygon, const_integrand, npts), 8.0 / 3.0, abs_tol); - EXPECT_NEAR(evaluate_area_integral(parabola_polygon, poly_integrand, npts), 64.0 / 105.0, abs_tol); - EXPECT_NEAR(evaluate_area_integral(parabola_polygon, transc_integrand, npts), 0.0, abs_tol); + // Ensure compatibility with curved polygons + BCurve pedges[2] = {para1, para2}; + CPolygon parabola_polygon(pedges, 2); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, const_integrand, npts), 8.0 / 3.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, poly_integrand, npts), 64.0 / 105.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, transc_integrand, npts), 0.0, abs_tol); - // Test on a unit square - Point2D squarenodes1[] = { Point2D {0.0, 0.0}, Point2D {1.0, 0.0} }; - BCurve square1(squarenodes1, 1); + // Test on a unit square + Point2D squarenodes1[] = {Point2D {0.0, 0.0}, Point2D {1.0, 0.0}}; + BCurve square1(squarenodes1, 1); - Point2D squarenodes2[] = { Point2D {1.0, 0.0}, Point2D {1.0, 1.0} }; - BCurve square2(squarenodes2, 1); + Point2D squarenodes2[] = {Point2D {1.0, 0.0}, Point2D {1.0, 1.0}}; + BCurve square2(squarenodes2, 1); - Point2D squarenodes3[] = { Point2D {1.0, 1.0}, Point2D {0.0, 1.0} }; - BCurve square3(squarenodes3, 1); + Point2D squarenodes3[] = {Point2D {1.0, 1.0}, Point2D {0.0, 1.0}}; + BCurve square3(squarenodes3, 1); - Point2D squarenodes4[] = { Point2D {0.0, 1.0}, Point2D {0.0, 0.0} }; - BCurve square4(squarenodes4, 1); + Point2D squarenodes4[] = {Point2D {0.0, 1.0}, Point2D {0.0, 0.0}}; + BCurve square4(squarenodes4, 1); - BCurve square_edges[] = { square1, square2, square3, square4 }; - CPolygon square(square_edges, 4); + BCurve square_edges[] = {square1, square2, square3, square4}; + CPolygon square(square_edges, 4); - EXPECT_NEAR(evaluate_area_integral(square, const_integrand, npts), 1.0, abs_tol); + EXPECT_NEAR(evaluate_area_integral(square, const_integrand, npts), 1.0, abs_tol); } TEST(primal_integral, evaluate_line_integral_scalar) { - using Point2D = primal::Point; - using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 30; - - // Define anonymous functions for testing - auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; - auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; - auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; - - // Test on single parabolic segment - Point2D paranodes[] = { Point2D {-1.0, 1.0}, Point2D {0.5, -2.0}, Point2D {2.0, 4.0} }; - BCurve parabola_segment(paranodes, 2); - - // Compare against hand computed/high-precision calculated values. - - // Constant integrand line integral is equivalent to arc-length calculation - EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, const_integrand, npts), - 6.12572661998, - abs_tol); - - EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, poly_integrand, npts), - 37.8010703669, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, transc_integrand, npts), - 0.495907795678, - abs_tol); - - // Test on a collection of Bezier curves - Point2D segnodes1[] = { Point2D {-1.0, -1.0}, - Point2D {-1.0 / 3.0, 1.0}, - Point2D {1.0 / 3.0, -1.0}, - Point2D {1.0, 1.0} }; - BCurve cubic_segment(segnodes1, 3); - - Point2D segnodes2[] = { Point2D {1.0, 1.0}, Point2D {-1.0, 0.0} }; - BCurve linear_segment(segnodes2, 1); - - Point2D segnodes3[] = { Point2D {-1.0, 0.0}, Point2D {-3.0, 1.0}, Point2D {-1.0, 2.0} }; - BCurve quadratic_segment(segnodes3, 2); - - BCurve connected_curve_edges[] = { cubic_segment, linear_segment, quadratic_segment }; - CPolygon connected_curve(connected_curve_edges, 3); - - EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, const_integrand, npts), - 8.28968500196, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, poly_integrand, npts), - -5.97565740064, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, transc_integrand, npts), - -0.574992518405, - abs_tol); - - // Test algorithm on disconnected curves - BCurve disconnected_curve_edges[] = { cubic_segment, quadratic_segment }; - CPolygon disconnected_curve(disconnected_curve_edges, 2); - - EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, const_integrand, npts), - 6.05361702446, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, poly_integrand, npts), - -6.34833539689, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, transc_integrand, npts), - -0.914161242161, - abs_tol); + using Point2D = primal::Point; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + // Define anonymous functions for testing + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + + // Test on single parabolic segment + Point2D paranodes[] = {Point2D {-1.0, 1.0}, Point2D {0.5, -2.0}, Point2D {2.0, 4.0}}; + BCurve parabola_segment(paranodes, 2); + + // Compare against hand computed/high-precision calculated values. + + // Constant integrand line integral is equivalent to arc-length calculation + EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, const_integrand, npts), + 6.12572661998, + abs_tol); + + EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, poly_integrand, npts), + 37.8010703669, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(parabola_segment, transc_integrand, npts), + 0.495907795678, + abs_tol); + + // Test on a collection of Bezier curves + Point2D segnodes1[] = {Point2D {-1.0, -1.0}, + Point2D {-1.0 / 3.0, 1.0}, + Point2D {1.0 / 3.0, -1.0}, + Point2D {1.0, 1.0}}; + BCurve cubic_segment(segnodes1, 3); + + Point2D segnodes2[] = {Point2D {1.0, 1.0}, Point2D {-1.0, 0.0}}; + BCurve linear_segment(segnodes2, 1); + + Point2D segnodes3[] = {Point2D {-1.0, 0.0}, Point2D {-3.0, 1.0}, Point2D {-1.0, 2.0}}; + BCurve quadratic_segment(segnodes3, 2); + + BCurve connected_curve_edges[] = {cubic_segment, linear_segment, quadratic_segment}; + CPolygon connected_curve(connected_curve_edges, 3); + + EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, const_integrand, npts), + 8.28968500196, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, poly_integrand, npts), + -5.97565740064, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(connected_curve, transc_integrand, npts), + -0.574992518405, + abs_tol); + + // Test algorithm on disconnected curves + BCurve disconnected_curve_edges[] = {cubic_segment, quadratic_segment}; + CPolygon disconnected_curve(disconnected_curve_edges, 2); + + EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, const_integrand, npts), + 6.05361702446, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, poly_integrand, npts), + -6.34833539689, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, transc_integrand, npts), + -0.914161242161, + abs_tol); } TEST(primal_integral, evaluate_line_integral_vector) { - using Point2D = primal::Point; - using Vector2D = primal::Vector; - using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 30; - - // Test on a single line segment - auto vec_field = [](Point2D x) -> Vector2D { return Vector2D({ x[1] * x[1], 3 * x[0] - 6 * x[1] }); }; - - Point2D segnodes[] = { Point2D {3.0, 7.0}, Point2D {0.0, 12.0} }; - BCurve linear_segment(segnodes, 1); - - // Compare against hand computed values - EXPECT_NEAR(evaluate_vector_line_integral(linear_segment, vec_field, npts), -1079.0 / 2.0, abs_tol); - - // Test on a closed curve - auto area_field = [](Point2D x) -> Vector2D { return Vector2D({ -0.5 * x[1], 0.5 * x[0] }); }; - auto conservative_field = [](Point2D x) -> Vector2D { - return Vector2D({ 2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1] }); - }; - auto winding_field = [](Point2D x) -> Vector2D { - double denom = 2 * M_PI * (x[0] * x[0] + x[1] * x[1]); - return Vector2D({ -x[1] / denom, x[0] / denom }); - }; - - Point2D paranodes1[] = { Point2D {1.0, 0.0}, Point2D {0.0, 2.0}, Point2D {-1.0, 0.0} }; - BCurve para1(paranodes1, 2); - - Point2D paranodes2[] = { Point2D {-1.0, 0.0}, Point2D {0.0, -2.0}, Point2D {1.0, 0.0} }; - BCurve para2(paranodes2, 2); - - BCurve parabola_shape_edges[] = { para1, para2 }; - CPolygon parabola_shape(parabola_shape_edges, 2); - - // This vector field calculates the area of the region - EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, area_field, npts), 8.0 / 3.0, abs_tol); - - // This vector field is conservative, so it should evaluate to zero - EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, conservative_field, npts), 0.0, abs_tol); - - // This vector field is generated by a in/out query, should return 1 (inside) - EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, winding_field, npts), 1.0, abs_tol); - - // Test algorithm on disconnected curves - Point2D paranodes2_shifted[] = { Point2D {-1.0, -1.0}, Point2D {0.0, -3.0}, Point2D {1.0, -1.0} }; - BCurve para2_shift(paranodes2_shifted, 2); - - BCurve disconnected_parabola_edges[] = { para1, para2_shift }; - CPolygon disconnected_parabola_shape(disconnected_parabola_edges, 2); - - EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, area_field, npts), - 11.0 / 3.0, - abs_tol); - EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, conservative_field, npts), - 0.0, - abs_tol); - EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, winding_field, npts), - 0.75, - abs_tol); + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + // Test on a single line segment + auto vec_field = [](Point2D x) -> Vector2D { return Vector2D({x[1] * x[1], 3 * x[0] - 6 * x[1]}); }; + + Point2D segnodes[] = {Point2D {3.0, 7.0}, Point2D {0.0, 12.0}}; + BCurve linear_segment(segnodes, 1); + + // Compare against hand computed values + EXPECT_NEAR(evaluate_vector_line_integral(linear_segment, vec_field, npts), -1079.0 / 2.0, abs_tol); + + // Test on a closed curve + auto area_field = [](Point2D x) -> Vector2D { return Vector2D({-0.5 * x[1], 0.5 * x[0]}); }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + auto winding_field = [](Point2D x) -> Vector2D { + double denom = 2 * M_PI * (x[0] * x[0] + x[1] * x[1]); + return Vector2D({-x[1] / denom, x[0] / denom}); + }; + + Point2D paranodes1[] = {Point2D {1.0, 0.0}, Point2D {0.0, 2.0}, Point2D {-1.0, 0.0}}; + BCurve para1(paranodes1, 2); + + Point2D paranodes2[] = {Point2D {-1.0, 0.0}, Point2D {0.0, -2.0}, Point2D {1.0, 0.0}}; + BCurve para2(paranodes2, 2); + + BCurve parabola_shape_edges[] = {para1, para2}; + CPolygon parabola_shape(parabola_shape_edges, 2); + + // This vector field calculates the area of the region + EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, area_field, npts), 8.0 / 3.0, abs_tol); + + // This vector field is conservative, so it should evaluate to zero + EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, conservative_field, npts), 0.0, abs_tol); + + // This vector field is generated by a in/out query, should return 1 (inside) + EXPECT_NEAR(evaluate_vector_line_integral(parabola_shape, winding_field, npts), 1.0, abs_tol); + + // Test algorithm on disconnected curves + Point2D paranodes2_shifted[] = {Point2D {-1.0, -1.0}, Point2D {0.0, -3.0}, Point2D {1.0, -1.0}}; + BCurve para2_shift(paranodes2_shifted, 2); + + BCurve disconnected_parabola_edges[] = {para1, para2_shift}; + CPolygon disconnected_parabola_shape(disconnected_parabola_edges, 2); + + EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, area_field, npts), + 11.0 / 3.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, conservative_field, npts), + 0.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(disconnected_parabola_shape, winding_field, npts), + 0.75, + abs_tol); } TEST(primal_integral, evaluate_integral_3D) { - using Point3D = primal::Point; - using Vector3D = primal::Vector; - using BCurve = primal::BezierCurve; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 30; - - const int order = 3; - Point3D data[order + 1] = { Point3D {0.6, 1.2, 1.0}, - Point3D {1.3, 1.6, 1.8}, - Point3D {2.9, 2.4, 2.3}, - Point3D {3.2, 3.5, 3.0} }; - BCurve spatial_arc(data, 3); - - auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; - auto transc_integrand = [](Point3D x) -> double { return std::sin(x[0] * x[1] * x[2]); }; - - auto vector_field = [](Point3D x) -> Vector3D { - return Vector3D({ 4 * x[1] * x[1], 8 * x[0] * x[1], 1.0 }); - }; - - // Test line integral on scalar domain againt values computed with external software - EXPECT_NEAR(evaluate_scalar_line_integral(spatial_arc, const_integrand, npts), - 4.09193268998, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(spatial_arc, transc_integrand, npts), - 0.515093324547, - abs_tol); - - // Test line integral on vector domain againt values computed with external software - EXPECT_NEAR(evaluate_vector_line_integral(spatial_arc, vector_field, npts), 155.344, abs_tol); + using Point3D = primal::Point; + using Vector3D = primal::Vector; + using BCurve = primal::BezierCurve; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + const int order = 3; + Point3D data[order + 1] = {Point3D {0.6, 1.2, 1.0}, + Point3D {1.3, 1.6, 1.8}, + Point3D {2.9, 2.4, 2.3}, + Point3D {3.2, 3.5, 3.0}}; + BCurve spatial_arc(data, 3); + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + auto transc_integrand = [](Point3D x) -> double { return std::sin(x[0] * x[1] * x[2]); }; + + auto vector_field = [](Point3D x) -> Vector3D { + return Vector3D({4 * x[1] * x[1], 8 * x[0] * x[1], 1.0}); + }; + + // Test line integral on scalar domain againt values computed with external software + EXPECT_NEAR(evaluate_scalar_line_integral(spatial_arc, const_integrand, npts), + 4.09193268998, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(spatial_arc, transc_integrand, npts), + 0.515093324547, + abs_tol); + + // Test line integral on vector domain againt values computed with external software + EXPECT_NEAR(evaluate_vector_line_integral(spatial_arc, vector_field, npts), 155.344, abs_tol); } TEST(primal_integral, evaluate_integral_rational) { - using Point2D = primal::Point; - using Vector2D = primal::Vector; - using BCurve = primal::BezierCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 20; - - // Elliptical arc shape - Point2D ellipse_nodes[] = { Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0} }; - double weights[] = { 2.0, 1.0, 1.0 }; - BCurve ellipse_arc(ellipse_nodes, weights, 2); - - Point2D leg1_nodes[] = { Point2D {0.0, 1.0}, {0.0, 0.0} }; - BCurve leg1(leg1_nodes, 1); - - Point2D leg2_nodes[] = { Point2D {0.0, 0.0}, {2.0, 0.0} }; - BCurve leg2(leg2_nodes, 1); - - CPolygon quarter_ellipse; - quarter_ellipse.addEdge(ellipse_arc); - quarter_ellipse.addEdge(leg1); - quarter_ellipse.addEdge(leg2); - - auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; - auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; - - auto area_field = [](Point2D x) -> Vector2D { return Vector2D({ -0.5 * x[1], 0.5 * x[0] }); }; - auto conservative_field = [](Point2D x) -> Vector2D { - return Vector2D({ 2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1] }); - }; - - // Test area integrals with scalar integrand againt values computed with external software - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); - - // Test line integral on scalar domain againt values computed with external software - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), - 2.42211205514, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), - 1.38837959326, - abs_tol); - - // Test line integral on vector domain againt values computed with external software - EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Elliptical arc shape + Point2D ellipse_nodes[] = {Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0}}; + double weights[] = {2.0, 1.0, 1.0}; + BCurve ellipse_arc(ellipse_nodes, weights, 2); + + Point2D leg1_nodes[] = {Point2D {0.0, 1.0}, {0.0, 0.0}}; + BCurve leg1(leg1_nodes, 1); + + Point2D leg2_nodes[] = {Point2D {0.0, 0.0}, {2.0, 0.0}}; + BCurve leg2(leg2_nodes, 1); + + CPolygon quarter_ellipse; + quarter_ellipse.addEdge(ellipse_arc); + quarter_ellipse.addEdge(leg1); + quarter_ellipse.addEdge(leg2); + + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + + auto area_field = [](Point2D x) -> Vector2D { return Vector2D({-0.5 * x[1], 0.5 * x[0]}); }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + + // Test area integrals with scalar integrand againt values computed with external software + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); + + // Test line integral on scalar domain againt values computed with external software + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), + 2.42211205514, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), + 1.38837959326, + abs_tol); + + // Test line integral on vector domain againt values computed with external software + EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); } TEST(primal_integral, evaluate_integral_nurbs) { - using Point2D = primal::Point; - using Vector2D = primal::Vector; - using NCurve = primal::NURBSCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 20; - - // Test integrals with same integrand and curves as `evaluate_integral_rational`, - // but insert some knots to make the Bezier extraction more interesting - - // Elliptical arc shape - Point2D ellipse_nodes[] = { Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0} }; - double weights[] = { 2.0, 1.0, 1.0 }; - NCurve ellipse_arc(ellipse_nodes, weights, 3, 2); - ellipse_arc.insertKnot(0.3, 2); - ellipse_arc.insertKnot(0.7, 1); - - Point2D leg1_nodes[] = { Point2D {0.0, 1.0}, {0.0, 0.0} }; - NCurve leg1(leg1_nodes, 2, 1); - leg1.insertKnot(0.4, 1); - - Point2D leg2_nodes[] = { Point2D {0.0, 0.0}, {2.0, 0.0} }; - NCurve leg2(leg2_nodes, 2, 1); - leg2.insertKnot(0.6, 1); - - CPolygon quarter_ellipse; - quarter_ellipse.addEdge(ellipse_arc); - quarter_ellipse.addEdge(leg1); - quarter_ellipse.addEdge(leg2); - - auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; - auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; - - auto area_field = [](Point2D x) -> Vector2D { return Vector2D({ -0.5 * x[1], 0.5 * x[0] }); }; - auto conservative_field = [](Point2D x) -> Vector2D { - return Vector2D({ 2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1] }); - }; - - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); - - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), - 2.42211205514, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), - 1.38837959326, - abs_tol); - - EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using NCurve = primal::NURBSCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Test integrals with same integrand and curves as `evaluate_integral_rational`, + // but insert some knots to make the Bezier extraction more interesting + + // Elliptical arc shape + Point2D ellipse_nodes[] = {Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0}}; + double weights[] = {2.0, 1.0, 1.0}; + NCurve ellipse_arc(ellipse_nodes, weights, 3, 2); + ellipse_arc.insertKnot(0.3, 2); + ellipse_arc.insertKnot(0.7, 1); + + Point2D leg1_nodes[] = {Point2D {0.0, 1.0}, {0.0, 0.0}}; + NCurve leg1(leg1_nodes, 2, 1); + leg1.insertKnot(0.4, 1); + + Point2D leg2_nodes[] = {Point2D {0.0, 0.0}, {2.0, 0.0}}; + NCurve leg2(leg2_nodes, 2, 1); + leg2.insertKnot(0.6, 1); + + CPolygon quarter_ellipse; + quarter_ellipse.addEdge(ellipse_arc); + quarter_ellipse.addEdge(leg1); + quarter_ellipse.addEdge(leg2); + + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + + auto area_field = [](Point2D x) -> Vector2D { return Vector2D({-0.5 * x[1], 0.5 * x[0]}); }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); + + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), + 2.42211205514, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), + 1.38837959326, + abs_tol); + + EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); } TEST(primal_integral, evaluate_nurbs_surface_normal) { - const int DIM = 3; - using Point2D = primal::Point; - using Point3D = primal::Point; - using NURBSPatchType = primal::NURBSPatch; + const int DIM = 3; + using Point2D = primal::Point; + using Point3D = primal::Point; + using NURBSPatchType = primal::NURBSPatch; - const int npts_u = 5; - const int npts_v = 4; + const int npts_u = 5; + const int npts_v = 4; - const int degree_u = 3; - const int degree_v = 2; + const int degree_u = 3; + const int degree_v = 2; - // clang-format off + // clang-format off Point3D controlPoints[5 * 4] = { Point3D {0, 0, 0}, Point3D {0, 4, 0}, Point3D {0, 8, -3}, Point3D {0, 12, 0}, Point3D {2, 0, 6}, Point3D {2, 4, 0}, Point3D {2, 8, 0}, Point3D {2, 12, 0}, @@ -407,139 +407,138 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) 3.0, 4.0, 5.0, 4.0, 4.0, 5.0, 6.0, 5.0, 5.0, 6.0, 7.0, 6.0 }; - // clang-format on + // clang-format on - NURBSPatchType nPatch(controlPoints, weights, npts_u, npts_v, degree_u, degree_v); - nPatch.makeTriviallyTrimmed(); + NURBSPatchType nPatch(controlPoints, weights, npts_u, npts_v, degree_u, degree_v); + nPatch.makeTriviallyTrimmed(); - int npts = 20; - double abs_tol = 1e-10; + int npts = 20; + double abs_tol = 1e-10; - // Calculate the average surface normal using boundary formulation - auto ueda_formula = nPatch.calculateUntrimmedPatchNormal(npts); + // Calculate the average surface normal using boundary formulation + auto ueda_formula = nPatch.calculateUntrimmedPatchNormal(npts); - // Calculate the average surface normal component by component, - // i.e. \int_S \hat{n}(x) dA, where \hat{n} is the unit normal at surface point x + // Calculate the average surface normal component by component, + // i.e. \int_S \hat{n}(x) dA, where \hat{n} is the unit normal at surface point x - // Because the NURBS surface has discontinuous derivatives at u, v = 0.5, - // this integrand does too, and we need to split to get accurate results - // using Gaussian quadrature + // Because the NURBS surface has discontinuous derivatives at u, v = 0.5, + // this integrand does too, and we need to split to get accurate results + // using Gaussian quadrature - NURBSPatchType split_patch[4]; - nPatch.split(0.5, 0.5, split_patch[0], split_patch[1], split_patch[2], split_patch[3]); + NURBSPatchType split_patch[4]; + nPatch.split(0.5, 0.5, split_patch[0], split_patch[1], split_patch[2], split_patch[3]); - for (int N = 0; N < 3; ++N) - { - auto avg_surface_normal_integrand = [&nPatch, &N](Point2D x) -> double { - return nPatch.normal(x[0], x[1])[N]; - }; + for(int N = 0; N < 3; ++N) + { + auto avg_surface_normal_integrand = [&nPatch, &N](Point2D x) -> double { + return nPatch.normal(x[0], x[1])[N]; + }; - // Integrate over each patch - double direct_formula = 0.0; - for (int i = 0; i < 4; ++i) - { - direct_formula += - evaluate_area_integral(split_patch[i].getTrimmingCurves(), avg_surface_normal_integrand, npts); - } + // Integrate over each patch + double direct_formula = 0.0; + for(int i = 0; i < 4; ++i) + { + direct_formula += + evaluate_area_integral(split_patch[i].getTrimmingCurves(), avg_surface_normal_integrand, npts); + } - EXPECT_NEAR(ueda_formula[N], direct_formula, abs_tol); - } + EXPECT_NEAR(ueda_formula[N], direct_formula, abs_tol); + } } - TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) { - using Point2D = primal::Point; - using Vector2D = primal::Vector; - using NCurve = primal::NURBSCurve; - using NCache = primal::NURBSCurve; - using CPolygon = primal::CurvedPolygon; - double abs_tol = 1e-10; - - // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 20; - - // Test integrals with same integrand and curves as `evaluate_integral_rational`, - // but with curves added to detail::NURBSCurveGWNCache objects to ensure template compatibility, - // even if there isn't compelling reason to use GWN caches for this purpose - - // Elliptical arc shape - Point2D ellipse_nodes[] = { Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0} }; - double weights[] = { 2.0, 1.0, 1.0 }; - NCurve ellipse_arc(ellipse_nodes, weights, 3, 2); - ellipse_arc.insertKnot(0.3, 2); - ellipse_arc.insertKnot(0.7, 1); - - Point2D leg1_nodes[] = { Point2D {0.0, 1.0}, {0.0, 0.0} }; - NCurve leg1(leg1_nodes, 2, 1); - leg1.insertKnot(0.4, 1); - - Point2D leg2_nodes[] = { Point2D {0.0, 0.0}, {2.0, 0.0} }; - NCurve leg2(leg2_nodes, 2, 1); - leg2.insertKnot(0.6, 1); - - CPolygon quarter_ellipse; - quarter_ellipse.addEdge(NCache(ellipse_arc)); - quarter_ellipse.addEdge(NCache(leg1)); - quarter_ellipse.addEdge(NCache(leg2)); - - auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; - auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; - - auto area_field = [](Point2D x) -> Vector2D { return Vector2D({ -0.5 * x[1], 0.5 * x[0] }); }; - auto conservative_field = [](Point2D x) -> Vector2D { - return Vector2D({ 2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1] }); - }; - - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); - - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), - 2.42211205514, - abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), - 1.38837959326, - abs_tol); - - EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), - M_PI * 2 * 1 / 4.0, - abs_tol); - EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using NCurve = primal::NURBSCurve; + using NCache = primal::NURBSCurve; + using CPolygon = primal::CurvedPolygon; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Test integrals with same integrand and curves as `evaluate_integral_rational`, + // but with curves added to detail::NURBSCurveGWNCache objects to ensure template compatibility, + // even if there isn't compelling reason to use GWN caches for this purpose + + // Elliptical arc shape + Point2D ellipse_nodes[] = {Point2D {2.0, 0.0}, Point2D {2.0, 1.0}, Point2D {0.0, 1.0}}; + double weights[] = {2.0, 1.0, 1.0}; + NCurve ellipse_arc(ellipse_nodes, weights, 3, 2); + ellipse_arc.insertKnot(0.3, 2); + ellipse_arc.insertKnot(0.7, 1); + + Point2D leg1_nodes[] = {Point2D {0.0, 1.0}, {0.0, 0.0}}; + NCurve leg1(leg1_nodes, 2, 1); + leg1.insertKnot(0.4, 1); + + Point2D leg2_nodes[] = {Point2D {0.0, 0.0}, {2.0, 0.0}}; + NCurve leg2(leg2_nodes, 2, 1); + leg2.insertKnot(0.6, 1); + + CPolygon quarter_ellipse; + quarter_ellipse.addEdge(NCache(ellipse_arc)); + quarter_ellipse.addEdge(NCache(leg1)); + quarter_ellipse.addEdge(NCache(leg2)); + + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto transc_integrand = [](Point2D x) -> double { return std::sin(x[0] * x[1]); }; + + auto area_field = [](Point2D x) -> Vector2D { return Vector2D({-0.5 * x[1], 0.5 * x[0]}); }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, const_integrand, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(quarter_ellipse, transc_integrand, npts), 0.472951736306, abs_tol); + + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, const_integrand, npts), + 2.42211205514, + abs_tol); + EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts), + 1.38837959326, + abs_tol); + + EXPECT_NEAR(evaluate_vector_line_integral(ellipse_arc, area_field, npts), + M_PI * 2 * 1 / 4.0, + abs_tol); + EXPECT_NEAR(evaluate_vector_line_integral(quarter_ellipse, conservative_field, npts), 0, abs_tol); } #ifdef AXOM_USE_MFEM TEST(primal_integral, check_axom_mfem_quadrature_values) { - const int N = 200; - - for (int npts = 1; npts <= N; ++npts) - { - // Generate the Axom quadrature rule - axom::numerics::QuadratureRule axom_rule = axom::numerics::compute_gauss_legendre(npts); - - // Generate the MFEM quadrature rule - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& mfem_rule = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - - // Check that the nodes and weights are the same between the two rules - for (int j = 0; j < npts; ++j) - { - EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, 1e-16); - EXPECT_NEAR(axom_rule.weight(j), mfem_rule.IntPoint(j).weight, 1e-16); - } - } + const int N = 200; + + for(int npts = 1; npts <= N; ++npts) + { + // Generate the Axom quadrature rule + axom::numerics::QuadratureRule axom_rule = axom::numerics::compute_gauss_legendre(npts); + + // Generate the MFEM quadrature rule + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& mfem_rule = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + // Check that the nodes and weights are the same between the two rules + for(int j = 0; j < npts; ++j) + { + EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, 1e-16); + EXPECT_NEAR(axom_rule.weight(j), mfem_rule.IntPoint(j).weight, 1e-16); + } + } } #endif int main(int argc, char* argv[]) { - ::testing::InitGoogleTest(&argc, argv); + ::testing::InitGoogleTest(&argc, argv); - axom::slic::SimpleLogger logger; + axom::slic::SimpleLogger logger; - int result = RUN_ALL_TESTS(); + int result = RUN_ALL_TESTS(); - return result; -} + return result; +} \ No newline at end of file From a688f90b8d73cc250b92c8ad5f755c2f04e51e49 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 17:37:15 -0700 Subject: [PATCH 14/21] Rearrange some methods --- .../detail/evaluate_integral_impl.hpp | 187 +++++++++++++++++- .../primal/operators/evaluate_integral.hpp | 165 +++++----------- src/axom/primal/tests/primal_bezier_curve.cpp | 2 +- 3 files changed, 237 insertions(+), 117 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 0bfc1b2eec..0d7783889b 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -11,6 +11,8 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" +#include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" #include "axom/core/numerics/quadrature.hpp" @@ -30,7 +32,7 @@ namespace detail * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object - * \param [in] integrand the lambda function representing the integrand. + * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] npts The number of quadrature points in the rule * \return the value of the integral @@ -57,13 +59,69 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< return full_quadrature; } +/*! + * \brief Evaluate a scalar field line integral on a single NURBS curve. + * + * Decompose the NURBS curve into Bezier segments, then sum the integral on each + * + * \param [in] n The NURBS curve object + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_scalar_line_integral_component(const primal::NURBSCurve& n, + Lambda&& scalar_integrand, + const int npts) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += + detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); + } + return total_integral; +} + +/*! + * \brief Evaluate the scalar integral on a single NURBS curve with cached data for GWN evaluation. + * + * The cache object has already decomposed the NURBS curve into Bezier segments, + * which are used to evaluate the integral over each + * + * \param [in] n The NURBS curve object + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_scalar_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& scalar_integrand, + const int npts) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_scalar_line_integral_component(this_bezier_data.getCurve(), + scalar_integrand, + npts); + } + + return total_integral; +} + /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * - * Evaluate the vector field line integral with Gauss-Legendre quadrature + * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object - * \param [in] integrand the lambda function representing the integrand. + * \param [in] vec_field the lambda function representing the integrand. * Must accept a 2D point as input and return a 2D axom vector object * \param [in] npts The number of quadrature points in the rule * \return the value of the integral @@ -90,6 +148,60 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< return full_quadrature; } +/*! + * \brief Evaluate a vector field line integral on a single NURBS curve. + * + * Decompose the NURBS curve into Bezier segments, then sum the integral on each + * + * \param [in] n The NURBS curve object + * \param [in] vec_field the lambda function representing the integrand. + * Must accept a 2D point as input and return a 2D axom vector object + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_vector_line_integral_component(const primal::NURBSCurve& n, + Lambda&& vec_field, + const int npts) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += detail::evaluate_vector_line_integral_component(beziers[i], vec_field, npts); + } + return total_integral; +} + +/*! + * \brief Evaluate the vector integral on a single NURBS curve with cached data for GWN evaluation. + * + * The cache object has already decomposed the NURBS curve into Bezier segments, + * which are used to evaluate the integral over each + * + * \param [in] n The NURBS curve object + * \param [in] vec_field the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts The number of quadrature points in the rule + * \return the value of the integral + */ +template +inline double evaluate_vector_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& vec_field, + const int npts) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += + detail::evaluate_vector_line_integral_component(this_bezier_data.getCurve(), vec_field, npts); + } + + return total_integral; +} + /*! * \brief Evaluate the area integral across one component of the curved polygon. * @@ -99,8 +211,8 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \param [in] cs the array of Bezier curve objects that bound the region - * \param [in] integrand the lambda function representing the integrand. + * \param [in] c The component Bezier curve + * \param [in] integrand The lambda function representing the scalar integrand. * Must accept a 2D point as input and return a double * \param [in] The lower bound of integration for the antiderivatives * \param [in] npts_Q The number of quadrature points for the line integral @@ -143,6 +255,71 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, return full_quadrature; } +/*! + * \brief Evaluate the area integral across one component NURBS of the region. + * + * Intended to be called for each NURBSCurve object bounding a closed 2D region. + * + * \param [in] n The component NURBSCurve object + * \param [in] integrand The lambda function representing the scalar integrand. + * Must accept a 2D point as input and return a double + * \param [in] The lower bound of integration for the antiderivatives + * \param [in] npts_Q The number of quadrature points for the line integral + * \param [in] npts_P The number of quadrature points for the antiderivative + * \return the value of the integral, which is mathematically meaningless. + */ +template +double evaluate_area_integral_component(const primal::NURBSCurve& n, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + auto beziers = n.extractBezier(); + double total_integral = 0.0; + for(int i = 0; i < beziers.size(); ++i) + { + total_integral += + detail::evaluate_area_integral_component(beziers[i], integrand, int_lb, npts_Q, npts_P); + } + return total_integral; +} + +/*! + * \brief Evaluate the area integral across one component NURBSCurveGWNCache of the region. + * + * Intended to be called for each NURBSCurveGWNCache object bounding a closed 2D region. + * + * \param [in] nc The component NURBSCurveGWNCache object + * \param [in] integrand The lambda function representing the scalar integrand. + * Must accept a 2D point as input and return a double + * \param [in] The lower bound of integration for the antiderivatives + * \param [in] npts_Q The number of quadrature points for the line integral + * \param [in] npts_P The number of quadrature points for the antiderivative + * \return the value of the integral, which is mathematically meaningless. + */ +template +inline double evaluate_area_integral_component(const primal::detail::NURBSCurveGWNCache& nc, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + double total_integral = 0.0; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + // Assuming the cache is properly initialized, this operation will never add to the cache + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_area_integral_component(this_bezier_data.getCurve(), + integrand, + int_lb, + npts_Q, + npts_P); + } + + return total_integral; +} + } // end namespace detail } // end namespace primal } // end namespace axom diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 4e3a1b8ece..e1f179154a 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -10,7 +10,7 @@ * by Bezier and NURBS curves, such as 2D area integrals and scalar/vector field line integrals * * Line integrals are evaluated numerically with Gauss-Legendre quadrature - * + * * 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. */ @@ -22,9 +22,6 @@ #include "axom/config.hpp" #include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/primal/geometry/BezierCurve.hpp" -#include "axom/primal/geometry/NURBSCurve.hpp" - #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" // C++ includes @@ -41,8 +38,10 @@ namespace primal * The line integral is evaluated on each curve in the CurvedPolygon, and added * together to represent the total integral. The Polygon need not be connected. * - * Evaluate the scalar field line integral with Gauss-Legendre quadrature + * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, + * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input and return a double @@ -50,8 +49,8 @@ namespace primal * on each edge of the CurvedPolygon * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly, Lambda&& scalar_integrand, int npts) { @@ -67,56 +66,26 @@ double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly } /*! - * \brief Evaluate a line integral on a single Bezier curve on a scalar field + * \brief Evaluate a line integral on a single generic curve on a scalar field * - * \param [in] c the Bezier curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] c the generic curve object * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const primal::BezierCurve& c, - Lambda&& scalar_integrand, - int npts) +template +double evaluate_scalar_line_integral(const CurveType& c, Lambda&& scalar_integrand, int npts) { return detail::evaluate_scalar_line_integral_component(c, scalar_integrand, npts); } -/*! - * \brief Evaluate a line integral on a single NURBS curve on a scalar field - * - * \param [in] n the NURBS curve object - * \param [in] scalar_integrand the lambda function representing the integrand. - * Must accept a Point as input, and return a double. - * \param [in] npts the number of quadrature nodes per knot span - * - * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts on each segment. - * - * \return the value of the integral - */ -template -double evaluate_scalar_line_integral(const primal::NURBSCurve& n, - Lambda&& scalar_integrand, - int npts) -{ - // Compute the line integral along each component. - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); i++) - { - total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); - } - - return total_integral; -} - /*! * \brief Evaluate a line integral on an array of NURBS curves on a scalar field * - * \param [in] narray the array of NURBS curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] carray The array of generic curve objects * \param [in] scalar_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a double. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -126,15 +95,15 @@ double evaluate_scalar_line_integral(const primal::NURBSCurve& n, * * \return the value of the integral */ -template -double evaluate_scalar_line_integral(const axom::Array>& narray, +template +double evaluate_scalar_line_integral(const axom::Array& carray, Lambda&& scalar_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_scalar_line_integral(narray[i], scalar_integrand, npts); + total_integral += evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); } return total_integral; @@ -147,8 +116,10 @@ double evaluate_scalar_line_integral(const axom::Array ~100) * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector @@ -156,8 +127,8 @@ double evaluate_scalar_line_integral(const axom::Array -double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, Lambda&& vector_integrand, int npts) { @@ -173,56 +144,26 @@ double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly } /*! - * \brief Evaluate a line integral on a single Bezier curve on a vector field + * \brief Evaluate a line integral on a single generic curve on a vector field * - * \param [in] c the Bezier curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] c the generic curve object * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input, and return a Vector. * \param [in] npts the number of quadrature nodes * \return the value of the integral */ -template -double evaluate_vector_line_integral(const primal::BezierCurve& c, - Lambda&& vector_integrand, - int npts) +template +double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integrand, int npts) { return detail::evaluate_vector_line_integral_component(c, vector_integrand, npts); } /*! - * \brief Evaluate a line integral on a single NURBS curve on a vector field + * \brief Evaluate a line integral on an array of generic curves on a vector field * - * \param [in] n the NURBS curve object - * \param [in] vector_integrand the lambda function representing the integrand. - * Must accept a Point as input, and return a Vector. - * \param [in] npts the number of quadrature nodes per knot span - * - * \note The NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts on each segment - * - * \return the value of the integral - */ -template -double evaluate_vector_line_integral(const primal::NURBSCurve& n, - Lambda&& vector_integrand, - int npts) -{ - // Compute the line integral along each component. - auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); i++) - { - total_integral += - detail::evaluate_vector_line_integral_component(beziers[i], vector_integrand, npts); - } - - return total_integral; -} - -/*! - * \brief Evaluate a line integral on an array of NURBS curves on a vector field - * - * \param [in] narray the array of NURBS curve object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] carray The array of generic curve objects * \param [in] vector_integrand the lambda function representing the integrand. * Must accept a Point as input and return a Vector. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -232,15 +173,15 @@ double evaluate_vector_line_integral(const primal::NURBSCurve& n, * * \return the value of the integral */ -template -double evaluate_vector_line_integral(const axom::Array>& narray, +template +double evaluate_vector_line_integral(const axom::Array& carray, Lambda&& vector_integrand, int npts) { double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_vector_line_integral(narray[i], vector_integrand, npts); + total_integral += evaluate_vector_line_integral(carray[i], vector_integrand, npts); } return total_integral; @@ -249,7 +190,12 @@ double evaluate_vector_line_integral(const axom::Array -double evaluate_area_integral(const primal::CurvedPolygon cpoly, +template +double evaluate_area_integral(const primal::CurvedPolygon cpoly, Lambda&& integrand, int npts_Q, int npts_P = 0) @@ -291,25 +237,22 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, } /*! - * \brief Evaluate an integral on the interior of an array of NURBS curves. + * \brief Evaluate an integral on the interior of a region bound by 2D curves * * See above definition for details. * - * \param [in] narray the array of NURBS curve objects that bound the region + * \param [in] carray the array of generic curve objects that bound the region * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] npts_Q the number of quadrature points to evaluate the line integral * \param [in] npts_P the number of quadrature points to evaluate the antiderivative - * - * \note Each NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature - * is computed using npts_Q * npts_P on each segment * * \note The numerical result is only meaningful if the curves enclose a region * * \return the value of the integral */ -template -double evaluate_area_integral(const axom::Array> narray, +template +double evaluate_area_integral(const axom::Array carray, Lambda&& integrand, int npts_Q, int npts_P = 0) @@ -319,26 +262,26 @@ double evaluate_area_integral(const axom::Array> narray npts_P = npts_Q; } - if(narray.empty()) + if(carray.empty()) { return 0.0; } // Use minimum y-coord of control nodes as lower bound for integration - double int_lb = narray[0][0][1]; - for(int i = 0; i < narray.size(); i++) + double int_lb = carray[0][0][1]; + for(int i = 0; i < carray.size(); i++) { - for(int j = 1; j < narray[i].getNumControlPoints(); j++) + for(int j = 1; j < carray[i].getNumControlPoints(); j++) { - int_lb = std::min(int_lb, narray[i][j][1]); + int_lb = std::min(int_lb, carray[i][j][1]); } } // Evaluate the antiderivative line integral along each component double total_integral = 0.0; - for(int i = 0; i < narray.size(); i++) + for(int i = 0; i < carray.size(); i++) { - auto beziers = narray[i].extractBezier(); + auto beziers = carray[i].extractBezier(); for(int j = 0; j < beziers.size(); j++) { @@ -353,4 +296,4 @@ double evaluate_area_integral(const axom::Array> narray } // namespace primal } // end namespace axom -#endif +#endif \ No newline at end of file diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index 49e249f2a7..a6bfdd35ed 100644 --- a/src/axom/primal/tests/primal_bezier_curve.cpp +++ b/src/axom/primal/tests/primal_bezier_curve.cpp @@ -170,7 +170,7 @@ TEST(primal_beziercurve, set_order) bCurve.setOrder(order); EXPECT_EQ(order, bCurve.getOrder()); - EXPECT_EQ(order - 1, bCurve.getNumControlPoints()); + EXPECT_EQ(order + 1, bCurve.getNumControlPoints()); bCurve[0] = controlPoints[0]; bCurve[1] = controlPoints[1]; From 7057b7cbf368b66cc0fb78ba6f9fd8d8cacdfbbb Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 11:18:43 -0700 Subject: [PATCH 15/21] Fix some scope issues --- src/axom/core/numerics/quadrature.hpp | 5 +-- src/axom/quest/SamplingShaper.hpp | 8 +--- .../detail/shaping/WindingNumberSampler.hpp | 45 +++++++++---------- 3 files changed, 25 insertions(+), 33 deletions(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index b07224c52f..11a6ba682e 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -39,15 +39,14 @@ class QuadratureRule double weight(size_t idx) const { return m_weights[idx]; }; //! \brief Accessofr for the size of the quadrature rule - int getNumPoints() const { return m_npts; } + int getNumPoints() const { return m_nodes.size(); } private: //! \brief Private constructor for use in compute_() methods. Only allocates memory - QuadratureRule(int npts) : m_nodes(npts), m_weights(npts), m_npts(npts) { }; + QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; axom::Array m_nodes; axom::Array m_weights; - int m_npts; }; /*! diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index 772b1ffc73..84ce4abfe3 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -210,14 +210,8 @@ class SamplingShaper : public Shaper AXOM_UNUSED_VAR(shapeDimension); if(useWindingNumberSampler(shape)) { - using CachedContourType = - typename axom::primal::CurvedPolygon>; - axom::Array contours_cached; - - for(auto& contour : m_contours) contours_cached.push_back(CachedContourType(contour)); - m_inoutSamplerWN = - std::make_unique>(shapeName, contours_cached.view()); + std::make_unique>(shapeName, m_contours.view()); m_inoutSamplerWN->computeBounds(); m_inoutSamplerWN->initSpatialIndex(this->m_vertexWeldThreshold); } diff --git a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp index 8795d7cf7a..1ea45d4786 100644 --- a/src/axom/quest/detail/shaping/WindingNumberSampler.hpp +++ b/src/axom/quest/detail/shaping/WindingNumberSampler.hpp @@ -63,8 +63,12 @@ class WindingNumberSampler // For now. using ExecSpace = axom::SEQ_EXEC; - using ContourType = primal::CurvedPolygon>; - using GeometryView = typename axom::ArrayView; + using CurvedPolygonType = primal::CurvedPolygon>; + using GeometryView = typename axom::ArrayView; + + using ContourCacheType = primal::CurvedPolygon>; + using ContourCacheArray = typename axom::Array; + using PointType = primal::Point; using GeometricBoundingBox = axom::primal::BoundingBox; using BVH = typename axom::spin::BVH; @@ -77,10 +81,13 @@ class WindingNumberSampler * \param geomView A view that contains the shapes being queried. * */ - WindingNumberSampler(const std::string& shapeName, GeometryView geomView) - : m_shapeName(shapeName) - , m_geometryView(geomView) - { } + WindingNumberSampler(const std::string& shapeName, GeometryView geomView) : m_shapeName(shapeName) + { + for(auto & contour : geomView) + { + m_contourCaches.push_back(ContourCacheType(contour)); + } + } ~WindingNumberSampler() = default; @@ -90,14 +97,6 @@ class WindingNumberSampler // no-op - We do it in initSpatialIndex. } - /*! - * \brief Initialize the NURBSCurveGWNCache objects that store intermediate calculations for GWN evaluation - */ - void initGWNCaches() - { - - } - /*! * \brief Initialize the BVH that is used for queries from the geometry view. */ @@ -106,15 +105,15 @@ class WindingNumberSampler AXOM_ANNOTATE_SCOPE("Initialize spatial index"); // Figure out bounding boxes for each geometric object. - const axom::IndexType geometrySize = m_geometryView.size(); + const axom::IndexType geometrySize = m_contourCaches.size(); axom::Array aabbs(geometrySize, geometrySize, axom::execution_space::allocatorID()); auto aabbsView = aabbs.view(); - const auto geometryView = m_geometryView; + const auto contourCaches = m_contourCaches; axom::for_all( geometrySize, - AXOM_LAMBDA(axom::IndexType i) { aabbsView[i] = geometryView[i].boundingBox(); }); + AXOM_LAMBDA(axom::IndexType i) { aabbsView[i] = contourCaches[i].boundingBox(); }); // Initialize the BVH using the bounding boxes. m_bvh.initialize(aabbs, aabbs.size()); @@ -225,7 +224,7 @@ class WindingNumberSampler axom::Array inOutResult(numQueryPoints, numQueryPoints, allocatorID); auto inOutResultView = inOutResult.view(); const auto candidatesView = candidates.view(); - const auto geometryView = m_geometryView; + const auto contourCaches = m_contourCaches; axom::for_all( numQueryPoints, AXOM_LAMBDA(axom::IndexType qpi) { @@ -236,7 +235,7 @@ class WindingNumberSampler for(axom::IndexType ci = 0; ci < numCandidates && in == false; ci++) { const auto candidateIndex = candidatesView[offsetsView[qpi] + ci]; - in |= detail::checkInside(geometryView[candidateIndex], queryPoint); + in |= detail::checkInside(contourCaches[candidateIndex], queryPoint); } inOutResultView[qpi] = in; }); @@ -295,15 +294,15 @@ class WindingNumberSampler PointProjector projector = {}) { AXOM_ANNOTATE_SCOPE("computeVolumeFractionsBaseline"); - const auto geometryView = m_geometryView; + const auto contourCaches = m_contourCaches; auto checkInside = [=](const PointType& pt) -> bool { // TODO: Use m_bvh to limit which curved polygons might contain point. // Check each candidate bool inside = false; - for(axom::IndexType i = 0; i < geometryView.size() && !inside; i++) + for(axom::IndexType i = 0; i < contourCaches.size() && !inside; i++) { - inside |= detail::checkInside(geometryView[i], pt); + inside |= detail::checkInside(contourCaches[i], pt); } return inside; }; @@ -337,7 +336,7 @@ class WindingNumberSampler std::string m_shapeName; GeometricBoundingBox m_bbox {}; - GeometryView m_geometryView; + ContourCacheArray m_contourCaches; BVH m_bvh {}; }; From 612bf92bad3b1cfd9ab4d0520a73a4197ca739f3 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 14:41:20 -0700 Subject: [PATCH 16/21] Add correct template to in_curved_polygon --- src/axom/primal/operators/detail/winding_number_2d_impl.hpp | 5 ----- src/axom/primal/operators/evaluate_integral.hpp | 4 ++-- src/axom/primal/operators/in_curved_polygon.hpp | 5 +++-- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp index d77f7d8aae..e79d9f84be 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp @@ -21,11 +21,6 @@ // C++ includes #include -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#endif - namespace axom { namespace primal diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index e1f179154a..55673ac082 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -103,7 +103,7 @@ double evaluate_scalar_line_integral(const axom::Array& carray, double total_integral = 0.0; for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); + total_integral += detail::evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); } return total_integral; @@ -181,7 +181,7 @@ double evaluate_vector_line_integral(const axom::Array& carray, double total_integral = 0.0; for(int i = 0; i < carray.size(); i++) { - total_integral += evaluate_vector_line_integral(carray[i], vector_integrand, npts); + total_integral += detail::evaluate_vector_line_integral(carray[i], vector_integrand, npts); } return total_integral; diff --git a/src/axom/primal/operators/in_curved_polygon.hpp b/src/axom/primal/operators/in_curved_polygon.hpp index 9fe1ff7885..f523d320fd 100644 --- a/src/axom/primal/operators/in_curved_polygon.hpp +++ b/src/axom/primal/operators/in_curved_polygon.hpp @@ -31,6 +31,7 @@ namespace primal /*! * \brief Robustly determine if query point is interior to a curved polygon * + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] query The query point to test * \param [in] cpoly The CurvedPolygon object to test for containment * \param [in] edge_tol The physical distance level at which objects are @@ -46,9 +47,9 @@ namespace primal * * \return A boolean value indicating containment. */ -template +template bool in_curved_polygon(const Point& query, - const CurvedPolygon>& cpoly, + const CurvedPolygon& cpoly, bool useNonzeroRule = true, double edge_tol = 1e-8, double EPS = 1e-8) From 0d51b839d9a94abb2d511b5af8c6836c1c708c74 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 15:31:10 -0700 Subject: [PATCH 17/21] Style, namespace guards --- src/axom/core/numerics/quadrature.hpp | 8 +++----- src/axom/core/tests/numerics_quadrature.hpp | 11 ----------- src/axom/primal/geometry/CurvedPolygon.hpp | 10 ++++++---- src/axom/primal/operators/evaluate_integral.hpp | 6 ++++-- src/axom/primal/tests/primal_integral.cpp | 8 ++++++-- 5 files changed, 19 insertions(+), 24 deletions(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index e466ea123c..9d261827c1 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -22,7 +22,7 @@ namespace numerics /*! * \class QuadratureRule * - * \brief Stores a fixed array of 1D quadrature points and weights + * \brief Stores fixed views to arrays of 1D quadrature points and weights */ class QuadratureRule { @@ -38,7 +38,7 @@ class QuadratureRule //! \brief Accessor for quadrature weights double weight(size_t idx) const { return m_weights[idx]; }; - //! \brief Accessofr for the size of the quadrature rule + //! \brief Accessor for the size of the quadrature rule int getNumPoints() const { return m_nodes.size(); } private: @@ -61,9 +61,7 @@ class QuadratureRule * \note This method constructs the points by scratch each time, without caching * \sa get_gauss_legendre(int) */ -void compute_gauss_legendre_data(int npts, - axom::Array& nodes, - axom::Array& weights); +void compute_gauss_legendre_data(int npts, axom::Array& nodes, axom::Array& weights); /*! * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index c1d276806f..9345777662 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -88,14 +88,3 @@ TEST(numerics_quadrature, gauss_legendre_math_check) } } } - -TEST(numerics_quadrature, gauss_legendre_cache_check) -{ - // The first two rules are put in the cache - axom::numerics::QuadratureRule first_rule = axom::numerics::get_gauss_legendre(20); - axom::numerics::QuadratureRule second_rule = axom::numerics::get_gauss_legendre(20); - - // Check that the two rules are located in the same place in memory, and the third isn't - // EXPECT_EQ(&first_rule.node(0), &second_rule.node(0)); - // EXPECT_EQ(&first_rule.weight(0), &second_rule.weight(0)); -} diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 979793e31e..624967663c 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -31,7 +31,8 @@ namespace axom { namespace primal { - +namespace internal +{ ///@{ /// \name Boilerplate to deduce the numeric type from the curve object template @@ -66,6 +67,7 @@ template struct has_cached_data> : std::true_type { }; ///@} +} // namespace internal // Forward declare the templated classes and operator functions template @@ -88,7 +90,7 @@ template class CurvedPolygon { public: - using T = typename get_numeric_type::type; + using T = typename internal::get_numeric_type::type; using PointType = typename CurveType::PointType; using VectorType = typename CurveType::VectorType; @@ -169,7 +171,7 @@ class CurvedPolygon void splitEdge(int idx, T t) { SLIC_ASSERT(idx < static_cast(m_edges.size())); - AXOM_STATIC_ASSERT_MSG(!has_cached_data::value, + AXOM_STATIC_ASSERT_MSG(!internal::has_cached_data::value, "splitEdge cannot be called on objects with associated cache data"); m_edges.reserve(numEdges() + 1); @@ -261,7 +263,7 @@ class CurvedPolygon void reverseOrientation() { AXOM_STATIC_ASSERT_MSG( - !has_cached_data::value, + !internal::has_cached_data::value, "reverseOrientation cannot be called on objects with associated cache data"); const int ngon = numEdges(); diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 55673ac082..5836563738 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -103,7 +103,8 @@ double evaluate_scalar_line_integral(const axom::Array& carray, double total_integral = 0.0; for(int i = 0; i < carray.size(); i++) { - total_integral += detail::evaluate_scalar_line_integral(carray[i], scalar_integrand, npts); + total_integral += + detail::evaluate_scalar_line_integral_component(carray[i], scalar_integrand, npts); } return total_integral; @@ -181,7 +182,8 @@ double evaluate_vector_line_integral(const axom::Array& carray, double total_integral = 0.0; for(int i = 0; i < carray.size(); i++) { - total_integral += detail::evaluate_vector_line_integral(carray[i], vector_integrand, npts); + total_integral += + detail::evaluate_vector_line_integral_component(carray[i], vector_integrand, npts); } return total_integral; diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 55d270b21a..50a789aa84 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -10,6 +10,10 @@ #include "gtest/gtest.h" +#ifdef AXOM_USE_MFEM + #include "mfem.hpp" +#endif + namespace primal = axom::primal; TEST(primal_integral, evaluate_area_integral) @@ -509,9 +513,9 @@ TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) #ifdef AXOM_USE_MFEM TEST(primal_integral, check_axom_mfem_quadrature_values) { - const int N = 3; + const int N = 200; - for(int npts = N; npts <= N; ++npts) + for(int npts = 1; npts <= N; ++npts) { // Generate the Axom quadrature rule axom::numerics::QuadratureRule axom_rule = axom::numerics::get_gauss_legendre(npts); From 1145e8353fb5f87438c152c1d4dc13120f017921 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 18:06:14 -0700 Subject: [PATCH 18/21] Fix typo --- src/axom/core/numerics/quadrature.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 2fe3230bea..a828b5e031 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -23,7 +23,7 @@ namespace numerics /*! * \class QuadratureRule * - * \brief Stores a fixed array of 1D quadrature points and weights + * \brief Stores fixed views to arrays of 1D quadrature points and weights */ class QuadratureRule { From 2fc9dc0bab494eabaaa076e2da484b9966cd4357 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 18:14:48 -0700 Subject: [PATCH 19/21] Handful of comments --- src/axom/primal/operators/evaluate_integral.hpp | 10 ++++------ src/axom/primal/tests/primal_integral.cpp | 2 +- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 5836563738..6736146a0e 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -10,7 +10,7 @@ * by Bezier and NURBS curves, such as 2D area integrals and scalar/vector field line integrals * * Line integrals are evaluated numerically with Gauss-Legendre quadrature - * + * * 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. */ @@ -38,8 +38,7 @@ namespace primal * The line integral is evaluated on each curve in the CurvedPolygon, and added * together to represent the total integral. The Polygon need not be connected. * - * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, - * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object @@ -117,9 +116,8 @@ double evaluate_scalar_line_integral(const axom::Array& carray, * The line integral is evaluated on each curve in the CurvedPolygon, and added * together to represent the total integral. The Polygon need not be connected. * - * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, - * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) - * + * Evaluate the scalar field line integral with Gauss-Legendre quadrature + * * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 7e35b95fe7..051d04b13b 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -546,4 +546,4 @@ int main(int argc, char* argv[]) int result = RUN_ALL_TESTS(); return result; -} \ No newline at end of file +} From 90a7d94b4e33ce688cc727db7ac7fda9fb36f858 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 18:38:25 -0700 Subject: [PATCH 20/21] Fix constructors --- src/axom/primal/geometry/CurvedPolygon.hpp | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 624967663c..1cab4f73f1 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -129,19 +129,16 @@ class CurvedPolygon } } - /// Constructor from an array of \a nEdges curves - CurvedPolygon(CurveType* curves, int nEdges) - { - SLIC_ASSERT(curves != nullptr); - SLIC_ASSERT(nEdges >= 1); + ///! \brief Constructor for CurvedPolygon from an ArrayView of curves + CurvedPolygon(axom::ArrayView curves) : m_edges(curves) { } - m_edges.reserve(nEdges); + ///! \brief Constructor from a c-style array of \a nEdges curves + CurvedPolygon(const CurveType* curves, int nEdges) + : CurvedPolygon(axom::ArrayView(curves, nEdges)) + { } - for(int e = 0; e < nEdges; ++e) - { - this->addEdge(curves[e]); - } - } + ///! \brief Constructor from a Axom array of curves + CurvedPolygon(const axom::Array& curves) : CurvedPolygon(curves.view()) { } /// Clears the list of edges void clear() { m_edges.clear(); } From a452c1bd3aaf267437714ccad0c4927d604703b0 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 23 Oct 2025 16:50:49 -0700 Subject: [PATCH 21/21] Minor changes from review --- src/axom/primal/geometry/CurvedPolygon.hpp | 2 +- src/axom/primal/operators/detail/evaluate_integral_impl.hpp | 4 ++-- src/axom/primal/tests/primal_integral.cpp | 2 +- src/axom/quest/detail/shaping/WindingNumberSampler.hpp | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 1cab4f73f1..7fc9a9e461 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -167,7 +167,7 @@ class CurvedPolygon /// Splits an edge "in place" void splitEdge(int idx, T t) { - SLIC_ASSERT(idx < static_cast(m_edges.size())); + SLIC_ASSERT(idx >= 0 && idx < static_cast(m_edges.size())); AXOM_STATIC_ASSERT_MSG(!internal::has_cached_data::value, "splitEdge cannot be called on objects with associated cache data"); diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 0d7783889b..c126f01fbc 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -75,7 +75,7 @@ inline double evaluate_scalar_line_integral_component(const primal::NURBSCurve>; using GeometryView = typename axom::ArrayView; - + using ContourCacheType = primal::CurvedPolygon>; using ContourCacheArray = typename axom::Array; @@ -83,7 +83,7 @@ class WindingNumberSampler */ WindingNumberSampler(const std::string& shapeName, GeometryView geomView) : m_shapeName(shapeName) { - for(auto & contour : geomView) + for(const auto& contour : geomView) { m_contourCaches.push_back(ContourCacheType(contour)); }