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 { diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 2db53b2e2c..fd31c42f91 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() { @@ -256,6 +259,9 @@ class BezierCurve /// Retrieves the control point at index \a idx const PointType& operator[](int idx) const { return m_controlPoints[idx]; } + const PointType& getInitPoint() const { return m_controlPoints[0]; } + const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + /// Returns a copy of the Bezier curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index a4afe50567..c8f285c950 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -18,8 +18,12 @@ #include "axom/primal/geometry/Point.hpp" #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" #include "axom/primal/geometry/BoundingBox.hpp" +// For NURBSCurveGWNCache objects +#include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" + #include #include @@ -27,13 +31,51 @@ namespace axom { namespace primal { -// Forward declare the templated classes and operator functions +namespace internal +{ +///@{ +/// \name Boilerplate to deduce the numeric type from the curve object +template +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; +}; +///@} + +///@{ +/// \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 +{ }; +///@} +} // namespace internal + +// 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 +86,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 internal::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 @@ -73,20 +115,30 @@ class CurvedPolygon m_edges.resize(nEdges); } - /// Constructor from an array of \a nEdges curves - CurvedPolygon(BezierCurveType* curves, int nEdges) + /*! + * \brief Copy constructor for another curve type + */ + template + CurvedPolygon(const CurvedPolygon& other_poly) { - SLIC_ASSERT(curves != nullptr); - SLIC_ASSERT(nEdges >= 1); - - m_edges.reserve(nEdges); - - for(int e = 0; e < nEdges; ++e) + m_edges.resize(other_poly.numEdges()); + for(int e = 0; e < other_poly.numEdges(); ++e) { - this->addEdge(curves[e]); + m_edges[e] = CurveType(other_poly[e]); } } + ///! \brief Constructor for CurvedPolygon from an ArrayView of curves + CurvedPolygon(axom::ArrayView curves) : m_edges(curves) { } + + ///! \brief Constructor from a c-style array of \a nEdges curves + CurvedPolygon(const CurveType* curves, int nEdges) + : CurvedPolygon(axom::ArrayView(curves, nEdges)) + { } + + ///! \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(); } @@ -105,36 +157,43 @@ class CurvedPolygon m_edges.resize(ngon); } - /// Appends a BezierCurve to the list of edges - void addEdge(const BezierCurveType& c1) { m_edges.push_back(c1); } + /// Appends a curve to the list of edges + void addEdge(const CurveType& c1) { m_edges.push_back(c1); } + + /// Consumes then appends a curve to the list of edges + void addEdge(CurveType&& c1) { m_edges.push_back(std::move(c1)); } /// Splits an edge "in place" void splitEdge(int idx, T t) { - SLIC_ASSERT(idx < static_cast(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"); + m_edges.reserve(numEdges() + 1); m_edges.insert(m_edges.begin() + idx + 1, 1, m_edges[idx]); auto& csplit = m_edges[idx]; csplit.split(t, m_edges[idx], m_edges[idx + 1]); } - axom::Array> 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 vertex at index idx */ - const BezierCurveType& operator[](int idx) const { return m_edges[idx]; } + /*! Retrieves the curve at index idx */ + CurveType& operator[](int idx) { return m_edges[idx]; } + /*! Retrieves the curve at index idx */ + const CurveType& operator[](int idx) const { return m_edges[idx]; } /// Tests equality of two CurvedPolygons - friend inline bool operator==(const CurvedPolygon& 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 +208,7 @@ class CurvedPolygon { const int sz = numEdges(); - os << "{" << sz << "-sided Bezier polygon:"; + os << "{" << sz << "-sided curved polygon:"; for(int i = 0; i < sz - 1; ++i) { os << m_edges[i] << ","; @@ -177,7 +236,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; } @@ -185,9 +244,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 +258,10 @@ class CurvedPolygon /// \brief Reverses orientation of a CurvedPolygon void reverseOrientation() { + AXOM_STATIC_ASSERT_MSG( + !internal::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 +295,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 a7366fdb1f..865b76e65d 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -621,6 +621,9 @@ class NURBSCurve */ const PointType& operator[](int idx) const { return m_controlPoints[idx]; } + const PointType& getInitPoint() const { return m_controlPoints[0]; } + const PointType& getEndPoint() const { return m_controlPoints[m_controlPoints.size() - 1]; } + /// \brief Returns a copy of the NURBS curve's control points CoordsVec getControlPoints() const { return m_controlPoints; } 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/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 0bfc1b2eec..804bdb9ecf 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" @@ -29,13 +31,13 @@ namespace detail * * Evaluate the scalar field line integral with Gauss-Legendre quadrature * + * \tparam Lambda A callable type taking an axom::Point type returning a numeric type * \param [in] c the Bezier curve object - * \param [in] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double + * \param [in] scalar_integrand the lambda function representing the integrand. * \param [in] npts The number of quadrature points in the rule * \return the value of the integral */ -template +template inline double evaluate_scalar_line_integral_component(const primal::BezierCurve& c, Lambda&& scalar_integrand, const int npts) @@ -57,20 +59,79 @@ 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 + * + * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \param [in] n The NURBS curve object + * \param [in] scalar_integrand the lambda function representing the integrand. + * \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) +{ + const 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], + std::forward(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 + * + * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \param [in] n The NURBS curve object + * \param [in] scalar_integrand the lambda function representing the integrand. + * \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(), + std::forward(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 * + * \tparam Lambda A callable type taking an axom::Point type returning an axom::Vector type * \param [in] c the Bezier curve object - * \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] vector_integrand the lambda function representing the integrand. * \param [in] npts The number of quadrature points in the rule * \return the value of the integral */ -template +template inline double evaluate_vector_line_integral_component(const primal::BezierCurve& c, - Lambda&& vec_field, + Lambda&& vector_integrand, const int npts) { const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); @@ -83,13 +144,72 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< // on which to evaluate dot product auto x_q = c.evaluate(quad.node(q)); auto dx_q = c.dt(quad.node(q)); - auto func_val = vec_field(x_q); + auto func_val = vector_integrand(x_q); full_quadrature += quad.weight(q) * Vector::dot_product(func_val, dx_q); } 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 + * + * \tparam Lambda A callable type taking an axom::Point type returning an axom::Vector type + * \param [in] n The NURBS curve object + * \param [in] vector_integrand the lambda function representing the integrand. + * \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&& vector_integrand, + const int npts) +{ + const 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], + std::forward(vector_integrand), + 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 + * + * \tparam Lambda A callable type taking an axom::Point type returning an axom::Vector type + * \param [in] n The NURBS curve object + * \param [in] vector_integrand the lambda function representing the integrand. + * \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&& vector_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_vector_line_integral_component(this_bezier_data.getCurve(), + std::forward(vector_integrand), + npts); + } + + return total_integral; +} + /*! * \brief Evaluate the area integral across one component of the curved polygon. * @@ -99,17 +219,17 @@ 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. - * Must accept a 2D point as input and return a double + * \tparam Lambda A callable type taking a 2D axom::Point type returning a numeric type + * \param [in] c The component Bezier curve + * \param [in] scalar_integrand The lambda function representing the scalar integrand. * \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 +template double evaluate_area_integral_component(const primal::BezierCurve& c, - Lambda&& integrand, + Lambda&& scalar_integrand, double int_lb, const int npts_Q, const int npts_P) @@ -134,7 +254,7 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, // Define interior quadrature points auto x_qxi = Point({x_q[0], (x_q[1] - int_lb) * quad_P.node(xi) + int_lb}); - antiderivative = quad_P.weight(xi) * (x_q[1] - int_lb) * integrand(x_qxi); + antiderivative = quad_P.weight(xi) * (x_q[1] - int_lb) * scalar_integrand(x_qxi); full_quadrature += quad_Q.weight(q) * c.dt(quad_Q.node(q))[0] * -antiderivative; } @@ -143,6 +263,74 @@ 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. + * + * \tparam Lambda A callable type taking a 2D axom::Point type returning a numeric type + * \param [in] n The component NURBSCurve object + * \param [in] scalar_integrand The lambda function representing the scalar integrand. + * \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&& scalar_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], + std::forward(scalar_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. + * + * \tparam Lambda A callable type taking a 2D axom::Point type returning a numeric type + * \param [in] nc The component NURBSCurveGWNCache object + * \param [in] scalar_integrand The lambda function representing the scalar integrand. + * \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&& scalar_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(), + std::forward(scalar_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/detail/winding_number_2d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp index 72b93c885c..3bc2d67cac 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 * @@ -79,6 +99,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; @@ -169,10 +194,39 @@ 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) + { + // 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; @@ -184,6 +238,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 diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 4e3a1b8ece..c000f27ca2 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -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 @@ -43,15 +40,16 @@ namespace primal * * Evaluate the scalar field line integral with Gauss-Legendre quadrature * + * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \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 * \param [in] npts the number of quadrature points to evaluate the line integral * 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) { @@ -60,65 +58,39 @@ double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly { // Compute the line integral along each component. total_integral += - detail::evaluate_scalar_line_integral_component(cpoly[i], scalar_integrand, npts); + detail::evaluate_scalar_line_integral_component(cpoly[i], + std::forward(scalar_integrand), + npts); } return total_integral; } /*! - * \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 Lambda A callable type taking an axom::Point type returning a numeric type + * \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; + return detail::evaluate_scalar_line_integral_component(c, + std::forward(scalar_integrand), + npts); } /*! * \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 Lambda A callable type taking an axom::Point type returning a numeric type + * \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 * * \note Each NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature @@ -126,15 +98,18 @@ 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 += + detail::evaluate_scalar_line_integral_component(carray[i], + std::forward(scalar_integrand), + npts); } return total_integral; @@ -147,17 +122,18 @@ double evaluate_scalar_line_integral(const axom::Array as input and return a Vector * \param [in] npts the number of quadrature points to evaluate the line integral * 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) { @@ -166,65 +142,39 @@ double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly { // Compute the line integral along each component. total_integral += - detail::evaluate_vector_line_integral_component(cpoly[i], vector_integrand, npts); + detail::evaluate_vector_line_integral_component(cpoly[i], + std::forward(vector_integrand), + npts); } return total_integral; } /*! - * \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 Lambda A callable type taking an axom::Point type returning an axom::Vector type + * \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); + return detail::evaluate_vector_line_integral_component(c, + std::forward(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 + * \tparam Lambda A callable type taking an axom::Point type returning an axom::Vector type + * \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 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 * * \note Each NURBS curve is decomposed into Bezier segments, and the Gaussian quadrature @@ -232,15 +182,18 @@ 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 += + detail::evaluate_vector_line_integral_component(carray[i], + std::forward(vector_integrand), + npts); } return total_integral; @@ -249,18 +202,24 @@ double evaluate_vector_line_integral(const axom::Array -double evaluate_area_integral(const primal::CurvedPolygon cpoly, - Lambda&& integrand, +template +double evaluate_area_integral(const primal::CurvedPolygon& cpoly, + Lambda&& scalar_integrand, int npts_Q, int npts_P = 0) { @@ -283,34 +242,35 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, double total_integral = 0.0; for(int i = 0; i < cpoly.numEdges(); i++) { - total_integral += - detail::evaluate_area_integral_component(cpoly[i], integrand, int_lb, npts_Q, npts_P); + total_integral += detail::evaluate_area_integral_component(cpoly[i], + std::forward(scalar_integrand), + int_lb, + npts_Q, + npts_P); } return total_integral; } /*! - * \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] integrand the lambda function representing the integrand. - * Must accept a 2D point as input and return a double + * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the geometry + * \param [in] carray the array of generic curve objects that bound the region + * \param [in] scalar_integrand the lambda function representing the integrand. * \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, - Lambda&& integrand, +template +double evaluate_area_integral(const axom::Array& carray, + Lambda&& scalar_integrand, int npts_Q, int npts_P = 0) { @@ -319,31 +279,35 @@ 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++) { total_integral += - detail::evaluate_area_integral_component(beziers[j], integrand, int_lb, npts_Q, npts_P); + detail::evaluate_area_integral_component(beziers[j], + std::forward(scalar_integrand), + int_lb, + npts_Q, + npts_P); } } @@ -353,4 +317,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/operators/in_curved_polygon.hpp b/src/axom/primal/operators/in_curved_polygon.hpp index ea8356de0d..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) diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index aca97e0174..2939dcff2a 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 @@ -179,18 +180,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; - double ret_val = 0.0; for(int i = 0; i < cpoly.numEdges(); i++) { - ret_val += detail::bezier_winding_number(q, cpoly[i], dummy_isOnCurve, edge_tol, EPS); + ret_val += winding_number(q, cpoly[i], edge_tol, EPS); } return ret_val; @@ -350,7 +349,7 @@ double winding_number(const Point& query, /*! * \brief Computes the GWN for an array of 2D points wrt an array of cached data for 2D NURBS curves * - * \param [in] query_arr The array of query point to test + * \param [in] query_arr The array of query points to test * \param [in] nurbs_curve_arr The array of memoized curve objects * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances @@ -388,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 @@ -418,8 +417,9 @@ 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 - * \param [in] cpoly The CurvedPolygon object + * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the curve + * \param [in] q_arr The array of query points to test + * \param [in] cpoly The CurvedPolygon object of generic curves * \param [in] edge_tol The physical distance level at which objects are considered indistinguishable * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances * @@ -430,9 +430,9 @@ axom::Array winding_number(const axom::Array>& query_arr, * * \return The GWN. */ -template -axom::Array winding_number(const axom::Array>& query_arr, - const CurvedPolygon& cpoly, +template +axom::Array winding_number(const axom::Array>& q_arr, + const CurvedPolygon& cpoly, double edge_tol = 1e-8, double EPS = 1e-8) { @@ -443,14 +443,45 @@ axom::Array winding_number(const axom::Array>& query_arr, cache_arr.emplace_back(detail::NURBSCurveGWNCache(cpoly[i], edge_tol)); } - axom::Array ret_val(query_arr.size()); - for(int n = 0; n < query_arr.size(); ++n) + 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; +} + +/*! + * \brief Computes the GWN for an array of 2D points wrt to a 2D curved polygon + * + * \param [in] q_arr The array of query points to test + * \param [in] cpoly The CurvedPolygon object of NURBS curves with cached GWN data + * \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(query_arr[n], cache_arr, edge_tol, EPS); + ret_val[n] = winding_number(q_arr[n], cpoly, edge_tol, EPS); } return ret_val; } + ///@{ //! @name Winding number operations between 3D points and primitives diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index e6df8e9e6b..a6bfdd35ed 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_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp index 80f1cac48e..db1dfc43f5 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,11 +20,11 @@ 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 -void checkMoments(const primal::CurvedPolygon& bPolygon, +void checkMoments(const primal::CurvedPolygon>& bPolygon, const CoordType expArea, const primal::Point& expMoment, double eps, @@ -47,19 +40,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(); @@ -88,84 +81,104 @@ primal::CurvedPolygon createPolygon( return bPolygon; } -//------------------------------------------------------------------------------ -TEST(primal_curvedpolygon, constructor) +template +void test_constructor() { - const int DIM = 3; - using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; - using BezierCurveType = primal::BezierCurve; + using CurvedPolygon = primal::CurvedPolygon; { - SLIC_INFO("Testing default CurvedPolygon constructor "); - CurvedPolygonType bPolygon; - + CurvedPolygon cPolygon; const int expNumEdges = 0; - EXPECT_EQ(expNumEdges, bPolygon.numEdges()); - EXPECT_EQ(expNumEdges, bPolygon.getEdges().size()); - EXPECT_EQ(axom::Array(), bPolygon.getEdges()); + EXPECT_EQ(expNumEdges, cPolygon.numEdges()); + EXPECT_EQ(expNumEdges, cPolygon.getEdges().size()); + EXPECT_EQ(axom::Array(), cPolygon.getEdges()); } { - SLIC_INFO("Testing CurvedPolygon numEdges constructor "); - - CurvedPolygonType bPolygon(1); + CurvedPolygon cPolygon(1); const int expNumEdges = 1; - EXPECT_EQ(expNumEdges, bPolygon.numEdges()); - EXPECT_EQ(expNumEdges, static_cast(bPolygon.getEdges().size())); + EXPECT_EQ(expNumEdges, cPolygon.numEdges()); + EXPECT_EQ(expNumEdges, static_cast(cPolygon.getEdges().size())); } } -//---------------------------------------------------------------------------------- -TEST(primal_curvedpolygon, add_edges) +//------------------------------------------------------------------------------ +TEST(primal_curvedpolygon, constructor) { - const int DIM = 2; - using CoordType = double; - using CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; - using BezierCurveType = primal::BezierCurve; - - SLIC_INFO("Test adding edges to empty CurvedPolygon"); - - CurvedPolygonType bPolygon; - EXPECT_EQ(0, bPolygon.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; - BezierCurveType bCurve(controlPoints, 1); + CurvedPolygon cPolygon; + EXPECT_EQ(0, cPolygon.numEdges()); - bPolygon.addEdge(bCurve); - bPolygon.addEdge(bCurve); + cPolygon.addEdge(added_curve); + cPolygon.addEdge(added_curve); - EXPECT_EQ(2, bPolygon.numEdges()); - for(int p = 0; p < bPolygon.numEdges(); ++p) + EXPECT_EQ(2, cPolygon.numEdges()); + for(int p = 0; p < cPolygon.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]); - } - } + 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) +{ + using PointType = primal::Point; + using BezierCurveType = primal::BezierCurve; + using NURBSCurveType = primal::NURBSCurve; + using NURBSCacheType = primal::detail::NURBSCurveGWNCache; + + SLIC_INFO("Test adding empty edges to empty CurvedPolygons"); + + PointType controlPoints[2] = {PointType {0.6, 1.2}, PointType {0.0, 1.6}}; + + 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 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}, @@ -177,38 +190,51 @@ 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); - EXPECT_FALSE(subPolygon.isClosed()); + + BezierPolygon subBezierPolygon = createBezierPolygon(subCP, suborders); + EXPECT_FALSE(subBezierPolygon.isClosed()); + EXPECT_FALSE(NURBSPolygon(subBezierPolygon).isClosed()); + EXPECT_FALSE(NURBSCachePolygon(subBezierPolygon).isClosed()); } { - CurvedPolygonType bPolygon = createPolygon(CP, orders); + BezierPolygon bPolygon = createBezierPolygon(CP, orders); EXPECT_EQ(3, bPolygon.numEdges()); EXPECT_TRUE(bPolygon.isClosed()); + NURBSPolygon nPolygon = NURBSPolygon(bPolygon); + EXPECT_EQ(3, nPolygon.numEdges()); + EXPECT_TRUE(nPolygon.isClosed()); + + NURBSCachePolygon ncPolygon = NURBSCachePolygon(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 = createPolygon(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 = NURBSCachePolygon(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; + 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"); - CurvedPolygonType bPolygon; + BezierPolygon bPolygon; EXPECT_EQ(0, bPolygon.numEdges()); EXPECT_FALSE(bPolygon.isClosed()); @@ -219,23 +245,54 @@ TEST(primal_curvedpolygon, isClosed_BiGon) PointType {0.8, .25}}; axom::Array orders = {2, 1}; - CurvedPolygonType poly = createPolygon(CP, orders); - EXPECT_TRUE(poly.isClosed()); + BezierPolygon bPoly = createBezierPolygon(CP, orders); + EXPECT_TRUE(bPoly.isClosed()); + EXPECT_TRUE(NURBSPolygon(bPoly).isClosed()); + EXPECT_TRUE(NURBSCachePolygon(bPoly).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()); + BezierPolygon bPoly2 = bPoly; + bPoly2[0][2] = PointType {0.8, 1.0}; + + EXPECT_FALSE(bPoly2.isClosed()); + EXPECT_FALSE(NURBSPolygon(bPoly2).isClosed()); + EXPECT_FALSE(NURBSCachePolygon(bPoly2).isClosed()); +} + +template +void test_split_edge(primal::CurvedPolygon& curved_polygon) +{ + using PointType = primal::Point; + + 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 CurvedPolygonType = primal::CurvedPolygon; - using PointType = primal::Point; - using BezierCurveType = primal::BezierCurve; + 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."); @@ -245,29 +302,14 @@ TEST(primal_curvedpolygon, split_edge) PointType {0.6, 1.2}}; axom::Array orders32 = {1, 1, 1}; - CurvedPolygonType bPolygon32 = createPolygon(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); - //std::cout << "Got here!! " << std::endl; - //std::cout << bPolygon32 << std::endl; - bPolygon32.splitEdge(0, .5); + BezierPolygon bPolygon32 = createBezierPolygon(CP, orders32); + NURBSPolygon nPolygon32(bPolygon32); - BezierCurveType bCurve2; - BezierCurveType bCurve3; - bCurve.split(.5, bCurve2, bCurve3); + test_split_edge(bPolygon32); + test_split_edge(nPolygon32); - 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)); } //---------------------------------------------------------------------------------- @@ -275,9 +317,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 +353,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 +364,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 +377,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 +392,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 +405,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 +419,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 +432,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 +444,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 +482,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 +493,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 +566,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/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 2bdb18c7ff..dc814e2f13 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -21,7 +21,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 @@ -95,7 +95,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 @@ -170,7 +170,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 @@ -272,7 +272,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 @@ -328,6 +328,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 @@ -351,10 +352,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]); }; @@ -396,19 +397,19 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) const int degree_v = 2; // 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}, - Point3D {4, 0, 0}, Point3D {4, 4, 0}, Point3D {4, 8, 3}, Point3D {4, 12, 0}, - Point3D {6, 0, 0}, Point3D {6, 4, -3}, Point3D {6, 8, 0}, Point3D {6, 12, 0}, - Point3D {8, 0, 0}, Point3D {8, 4, 0}, Point3D {8, 8, 0}, Point3D {8, 12, 0}}; - - double weights[5 * 4] = { - 1.0, 2.0, 3.0, 2.0, - 2.0, 3.0, 4.0, 3.0, - 3.0, 4.0, 5.0, 4.0, - 4.0, 5.0, 6.0, 5.0, - 5.0, 6.0, 7.0, 6.0}; + 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}, + Point3D {4, 0, 0}, Point3D {4, 4, 0}, Point3D {4, 8, 3}, Point3D {4, 12, 0}, + Point3D {6, 0, 0}, Point3D {6, 4, -3}, Point3D {6, 8, 0}, Point3D {6, 12, 0}, + Point3D {8, 0, 0}, Point3D {8, 4, 0}, Point3D {8, 8, 0}, Point3D {8, 12, 0} }; + + double weights[5 * 4] = { + 1.0, 2.0, 3.0, 2.0, + 2.0, 3.0, 4.0, 3.0, + 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 NURBSPatchType nPatch(controlPoints, weights, npts_u, npts_v, degree_u, degree_v); @@ -448,12 +449,74 @@ 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 a 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 = 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); 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; diff --git a/src/axom/quest/SamplingShaper.hpp b/src/axom/quest/SamplingShaper.hpp index a10fe4b54c..fe5f0c1cf3 100644 --- a/src/axom/quest/SamplingShaper.hpp +++ b/src/axom/quest/SamplingShaper.hpp @@ -981,7 +981,7 @@ class SamplingShaper : public Shaper // Holds an instance of the 2D or 3D sampler; only one can be active at a time SamplerVariant m_sampler; - 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..2ad0cd851f 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,12 @@ class WindingNumberSampler // For now. using ExecSpace = axom::SEQ_EXEC; - 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; @@ -81,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(const auto& contour : geomView) + { + m_contourCaches.push_back(ContourCacheType(contour)); + } + } ~WindingNumberSampler() = default; @@ -102,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()); @@ -221,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) { @@ -232,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; }); @@ -291,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; }; @@ -333,7 +336,7 @@ class WindingNumberSampler std::string m_shapeName; GeometricBoundingBox m_bbox {}; - GeometryView m_geometryView; + ContourCacheArray m_contourCaches; BVH m_bvh {}; }; diff --git a/src/axom/quest/io/MFEMReader.cpp b/src/axom/quest/io/MFEMReader.cpp index 40decb987d..eea18d0a19 100644 --- a/src/axom/quest/io/MFEMReader.cpp +++ b/src/axom/quest/io/MFEMReader.cpp @@ -127,7 +127,7 @@ int MFEMReader::read(CurvedPolygonArray &curvedPolygons) auto &poly = curvedPolygons[contourId]; for(int zoneId : zoneIds) { - auto curve = axom::quest::internal::segment_to_curve(mesh, zoneId); + auto curve = axom::quest::internal::segment_to_nurbs(mesh, zoneId); 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: 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);