diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 68ea4a719b..6f848d8491 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -46,6 +46,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Adds the `AXOM_TEST_NUM_OMP_THREADS` configuration variable to control the default OpenMP thread count for tests. ### Changed +- Evaluation methods for line integrals in `axom::primal` have been generalized, and + `evaluate_scalar_line_integral` has been renamed to `evaluate_line_integral`. - Treatment of materials on strided-structured Blueprint meshes has changed in `axom::mir`. Materials are now expected to be defined only on the valid subset of zones in the mesh. This more closely matches VisIt behavior. diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 64c3600a51..8355885570 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -59,6 +59,7 @@ template class BezierCurve { public: + using NumericType = T; using PointType = Point; using VectorType = Vector; using SegmentType = Segment; diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 4fea324170..d83ceaec46 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -34,31 +34,7 @@ namespace primal 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 +/// \name Type traits for a compile-time flag for the cached object template struct has_cached_data : std::false_type { }; @@ -80,8 +56,7 @@ std::ostream& operator<<(std::ostream& os, const CurvedPolygon& poly) /*! * \class CurvedPolygon * - * \brief Represents a polygon with curved edges defined by BezierCurves - * \tparam T the coordinate type, e.g., double, float, etc. + * \brief Represents a polygon with generic curves for edges * \tparam NDIMS the number of dimensions * \note The component curves should be ordered in a counter clockwise * orientation with respect to the polygon's normal vector @@ -90,8 +65,7 @@ template class CurvedPolygon { public: - using T = typename internal::get_numeric_type::type; - + using NumericType = typename CurveType::NumericType; using PointType = typename CurveType::PointType; using VectorType = typename CurveType::VectorType; using BoundingBoxType = typename CurveType::BoundingBoxType; @@ -164,7 +138,7 @@ class CurvedPolygon void addEdge(CurveType&& c1) { m_edges.push_back(std::move(c1)); } /// Splits an edge "in place" - void splitEdge(int idx, T t) + void splitEdge(int idx, NumericType t) { SLIC_ASSERT(idx >= 0 && idx < static_cast(m_edges.size())); AXOM_STATIC_ASSERT_MSG(!internal::has_cached_data::value, diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index 1beedb8094..d4b8c4c547 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -61,6 +61,7 @@ template class NURBSCurve { public: + using NumericType = T; using PointType = Point; using VectorType = Vector; using SegmentType = Segment; diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index cf1f8c4820..7624d4618a 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -3111,25 +3111,16 @@ class NURBSPatch { SLIC_ASSERT(NDIMS == 3); - // Split the patch along the unique knot values to improve convergence - auto split_patches = extractTrimmedBezier(); - VectorType ret_vec; - for(int n = 0; n < split_patches.size(); ++n) - { - // Integrand for the surface area integral - auto& nPatch = split_patches[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]; - }; - // Find the area of the resulting projection - ret_vec[N] += - evaluate_area_integral(nPatch.getTrimmingCurves(), avg_surface_normal_integrand, npts); - } + // Split the patch along the unique knot values to improve convergence + for(const auto& nPatch : extractTrimmedBezier()) + { + // Integrate the surface normal over the patches + ret_vec += evaluate_area_integral( + nPatch.getTrimmingCurves(), + [&nPatch](Point2D x) -> Vector { return nPatch.normal(x[0], x[1]); }, + npts); } return ret_vec; @@ -3397,11 +3388,9 @@ class NURBSPatch const double sq_tol = 1e-14; const double EPS = 1e-6; - // Extract the Bezier curves of the NURBS curve - auto beziers = curve.extractBezier(); + // Extract the Bezier curves of the NURBS curve, checking each for intersection axom::Array knot_vals = curve.getKnots().getUniqueKnots(); - - // Check each Bezier segment for intersection + const auto beziers = curve.extractBezier(); for(int i = 0; i < beziers.size(); ++i) { axom::Array temp_curve_p; @@ -3849,11 +3838,9 @@ class NURBSPatch const double sq_tol = 1e-14; const double EPS = 1e-6; - // Extract the Bezier curves of the NURBS curve - auto beziers = curve.extractBezier(); + // Extract the Bezier curves of the NURBS curve, and check each for intersection axom::Array knot_vals = curve.getKnots().getUniqueKnots(); - - // Check each Bezier segment for intersection + const auto beziers = curve.extractBezier(); for(int i = 0; i < beziers.size(); ++i) { axom::Array temp_curve_p; diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 804bdb9ecf..bd59d48dc2 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -3,6 +3,24 @@ // // SPDX-License-Identifier: (BSD-3-Clause) +/*! + * \file evaluate_integral.hpp + * + * \brief Consists of methods that evaluate scalar-field integrals on curves and + * regions defined by 2D curves, and vector-field integrals on curves + * + * All integrals are evaluated numerically with Gauss-Legendre quadrature + * + * Scalar-field line integrals and scalar-field area integrals are of form + * int_D f(x) dr, with f : R^n -> R^m, D is a curve or a 2D region bound by curves + * + * Vector-field line integrals are of form int_C f(x) \cdot d\vec{r}, + * with f : R^n -> R^n, C is a curve + * + * 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar + * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. + */ + #ifndef PRIMAL_EVAL_INTEGRAL_IMPL_HPP_ #define PRIMAL_EVAL_INTEGRAL_IMPL_HPP_ @@ -18,6 +36,8 @@ // C++ includes #include +#include +#include namespace axom { @@ -25,27 +45,63 @@ namespace primal { namespace detail { +namespace internal +{ +///@{ +/// \name Type traits to support integrals of functions with general return types, +/// provided it supports addition and scalar multiplication +template +struct has_addition : std::false_type +{ }; + +template +struct has_addition() + std::declval())>> : std::true_type +{ }; + +template +struct has_scalar_multiplication : std::false_type +{ }; + +template +struct has_scalar_multiplication() * std::declval())>> + : std::true_type +{ }; + +template +using is_integrable = std::conjunction, has_scalar_multiplication>; + +template +constexpr bool is_integrable_v = is_integrable::value; +///@} +} // namespace internal + +///@{ +/// \name Evaluates scalar-field line integrals for functions f : R^n -> R^m /*! - * \brief Evaluate a scalar field line integral on a single Bezier curve. + * \brief Evaluate a line integral on a single Bezier curve. * - * Evaluate the scalar field line integral with Gauss-Legendre quadrature + * Evaluate the line integral with Gauss-Legendre quadrature * - * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] c the Bezier curve object - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] 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::BezierCurve& c, - Lambda&& scalar_integrand, +template ::PointType>> +inline LambdaRetType evaluate_line_integral_component(const BezierCurve& c, + Lambda&& integrand, const int npts) { const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); // Store/compute quadrature result - double full_quadrature = 0.0; + LambdaRetType full_quadrature = LambdaRetType {}; for(int q = 0; q < npts; q++) { // Get intermediate quadrature point @@ -53,91 +109,98 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< auto x_q = c.evaluate(quad.node(q)); auto dx_q = c.dt(quad.node(q)); - full_quadrature += quad.weight(q) * scalar_integrand(x_q) * dx_q.norm(); + full_quadrature += quad.weight(q) * integrand(x_q) * dx_q.norm(); } return full_quadrature; } /*! - * \brief Evaluate a scalar field line integral on a single NURBS curve. + * \brief Evaluate a 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 + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam LambdaRetType The return type of Lambda, which must support addition and scalar multiplication * \param [in] n The NURBS curve object - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] 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, +template ::PointType>> +inline LambdaRetType evaluate_line_integral_component(const NURBSCurve& n, + Lambda&& integrand, const int npts) { - const auto beziers = n.extractBezier(); - double total_integral = 0.0; - for(int i = 0; i < beziers.size(); ++i) + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) { total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], - std::forward(scalar_integrand), - npts); + detail::evaluate_line_integral_component(bez, std::forward(integrand), npts); } return total_integral; } /*! - * \brief Evaluate the scalar integral on a single NURBS curve with cached data for GWN evaluation. + * \brief Evaluate an 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 + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] n The NURBS curve object - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] 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, +template ::PointType>> +inline LambdaRetType evaluate_line_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& integrand, const int npts) { - double total_integral = 0.0; + LambdaRetType total_integral = LambdaRetType {}; 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); + total_integral += detail::evaluate_line_integral_component(this_bezier_data.getCurve(), + std::forward(integrand), + npts); } return total_integral; } +//@} + +///@{ +/// \name Evaluates vector-field line integrals for functions f : R^n -> R^n /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * - * Evaluate the scalar field line integral with Gauss-Legendre quadrature + * Evaluate the vector field line integral with Gauss-Legendre quadrature * - * \tparam Lambda A callable type taking an axom::Point type returning an axom::Vector type + * \tparam Lambda A callable type taking an CurveType's PointType and returning its numeric type * \param [in] c the Bezier 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::BezierCurve& c, - Lambda&& vector_integrand, - const int npts) +inline T evaluate_vector_line_integral_component(const primal::BezierCurve& c, + Lambda&& vector_integrand, + const int npts) { const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); // Store/compute quadrature result - double full_quadrature = 0.0; + T full_quadrature = T {}; for(int q = 0; q < npts; q++) { // Get intermediate quadrature point @@ -148,6 +211,7 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< full_quadrature += quad.weight(q) * Vector::dot_product(func_val, dx_q); } + return full_quadrature; } @@ -156,23 +220,22 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< * * 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 + * \tparam Lambda A callable type taking an CurveType's PointType and returning its numeric 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) +inline T 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) + T total_integral = T {}; + for(const auto& bez : n.extractBezier()) { total_integral += - detail::evaluate_vector_line_integral_component(beziers[i], + detail::evaluate_vector_line_integral_component(bez, std::forward(vector_integrand), npts); } @@ -185,18 +248,18 @@ inline double evaluate_vector_line_integral_component(const primal::NURBSCurve -inline double evaluate_vector_line_integral_component(const primal::detail::NURBSCurveGWNCache& nc, - Lambda&& vector_integrand, - const int npts) +template +inline T evaluate_vector_line_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& vector_integrand, + const int npts) { - double total_integral = 0.0; + T total_integral = T {}; for(int i = 0; i < nc.getNumKnotSpans(); ++i) { // Assuming the cache is properly initialized, this operation will never add to the cache @@ -209,6 +272,10 @@ inline double evaluate_vector_line_integral_component(const primal::detail::NURB return total_integral; } +//@} + +///@{ +/// \name Evaluates scalar-field 2D area integrals for functions f : R^2 -> R^m /*! * \brief Evaluate the area integral across one component of the curved polygon. @@ -219,29 +286,29 @@ inline double evaluate_vector_line_integral_component(const primal::detail::NURB * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \tparam Lambda A callable type taking a 2D axom::Point type returning a numeric type + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] c The component Bezier curve - * \param [in] scalar_integrand The lambda function representing the scalar integrand. + * \param [in] 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::BezierCurve& c, - Lambda&& scalar_integrand, - double int_lb, - const int npts_Q, - const int npts_P) +template ::PointType>> +LambdaRetType evaluate_area_integral_component(const primal::BezierCurve& c, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) { const axom::numerics::QuadratureRule& quad_Q = axom::numerics::get_gauss_legendre(npts_Q); const axom::numerics::QuadratureRule& quad_P = axom::numerics::get_gauss_legendre(npts_P); - // Store some intermediate values - double antiderivative = 0.0; - // Store/compute quadrature result - double full_quadrature = 0.0; + LambdaRetType full_quadrature = LambdaRetType {}; for(int q = 0; q < npts_Q; q++) { // Get intermediate quadrature point @@ -254,7 +321,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) * scalar_integrand(x_qxi); + auto antiderivative = quad_P.weight(xi) * (x_q[1] - int_lb) * integrand(x_qxi); full_quadrature += quad_Q.weight(q) * c.dt(quad_Q.node(q))[0] * -antiderivative; } @@ -268,27 +335,29 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, * * 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 + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \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] integrand The lambda function representing the scalar integrand. + * \param [in] int_lb 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) +template ::PointType>> +LambdaRetType 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) + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) { - total_integral += detail::evaluate_area_integral_component(beziers[i], - std::forward(scalar_integrand), + total_integral += detail::evaluate_area_integral_component(bez, + std::forward(integrand), int_lb, npts_Q, npts_P); @@ -301,28 +370,31 @@ double evaluate_area_integral_component(const primal::NURBSCurve& n, * * 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 + * \tparam Lambda A callable type taking an CurveType's PointType and returning an integrable type + * \tparam RetType The return type of Lambda, which must support addition and scalar multiplication * \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] integrand The lambda function representing the scalar integrand. + * \param [in] int_lb 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) +template ::PointType>> +inline RetType evaluate_area_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) { - double total_integral = 0.0; + RetType total_integral = RetType {}; 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), + std::forward(integrand), int_lb, npts_Q, npts_P); @@ -330,6 +402,7 @@ inline double evaluate_area_integral_component(const primal::detail::NURBSCurveG return total_integral; } +//@} } // end namespace detail } // end namespace primal 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 3bc2d67cac..5cfb479b30 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp @@ -100,6 +100,7 @@ template class NURBSCurveGWNCache { public: + using NumericType = T; using PointType = typename NURBSCurve::PointType; using VectorType = typename NURBSCurve::VectorType; using BoundingBoxType = typename NURBSCurve::BoundingBoxType; diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index c000f27ca2..2f151c759b 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -6,11 +6,17 @@ /*! * \file evaluate_integral.hpp * - * \brief Consists of methods that evaluate integrals over regions defined - * 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 + * \brief Consists of methods that evaluate scalar-field integrals on curves and + * regions defined by 2D curves, and vector-field integrals on curves * + * All integrals are evaluated numerically with Gauss-Legendre quadrature + * + * Scalar-field line integrals and scalar-field area integrals are of form + * int_D f(x) dr, with f : R^n -> R^m, D is a curve or a 2D region bound by curves + * + * Vector-field line integrals are of form int_C f(x) \cdot d\vec{r}, + * with f : R^n -> R^n, C is a curve + * * 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. */ @@ -31,66 +37,84 @@ namespace axom { namespace primal { +///@{ +/// \name Evaluates scalar-field line integrals for functions f : R^n -> R^m + /*! * \brief Evaluate a line integral along the boundary of a CurvedPolygon object - * on a scalar field. + * for a function with an arbitrary return type. * * 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. + * together to represent the total integral. The curved polygon need not be connected. * - * Evaluate the scalar field line integral with Gauss-Legendre quadrature + * Evaluate the 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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] cpoly the CurvedPolygon object - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] integrand the lambda function representing the integrand. * \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, - Lambda&& scalar_integrand, +template > +LambdaRetType evaluate_line_integral(const primal::CurvedPolygon cpoly, + Lambda&& integrand, int npts) { - double total_integral = 0.0; + static_assert( + detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda function " + "return type"); + + LambdaRetType total_integral = LambdaRetType {}; for(int i = 0; i < cpoly.numEdges(); i++) { // Compute the line integral along each component. total_integral += - detail::evaluate_scalar_line_integral_component(cpoly[i], - std::forward(scalar_integrand), - npts); + detail::evaluate_line_integral_component(cpoly[i], std::forward(integrand), npts); } return total_integral; } /*! - * \brief Evaluate a line integral on a single generic curve on a scalar field + * \brief Evaluate a line integral along the boundary of a generic curve + * for a function with an arbitrary return type. * - * \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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] c the generic curve object - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] integrand the lambda function representing the integrand. * \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 > +LambdaRetType evaluate_line_integral(const CurveType& c, Lambda&& integrand, int npts) { - return detail::evaluate_scalar_line_integral_component(c, - std::forward(scalar_integrand), - npts); + static_assert( + detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda function " + "return type"); + + return detail::evaluate_line_integral_component(c, std::forward(integrand), npts); } /*! * \brief Evaluate a line integral on an array of NURBS curves on a scalar field + * for a function with an arbitrary return type. * - * \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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning an integrable type + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \param [in] carray The array of generic curve objects - * \param [in] scalar_integrand the lambda function representing the integrand. + * \param [in] integrand the lambda function representing the integrand. * \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 @@ -98,46 +122,54 @@ 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, - Lambda&& scalar_integrand, - int npts) +template > +LambdaRetType evaluate_line_integral(const axom::Array& carray, Lambda&& integrand, int npts) { - double total_integral = 0.0; + static_assert( + detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda function " + "return type"); + + LambdaRetType total_integral = LambdaRetType {}; for(int i = 0; i < carray.size(); i++) { total_integral += - detail::evaluate_scalar_line_integral_component(carray[i], - std::forward(scalar_integrand), - npts); + detail::evaluate_line_integral_component(carray[i], std::forward(integrand), npts); } return total_integral; } +//@} + +///@{ +/// \name Evaluates vector-field line integrals for functions f : R^n -> R^n /*! - * \brief Evaluate a line integral along the boundary of a CurvedPolygon object - * on a vector field. + * \brief Evaluate a vector-field line integral along the boundary of a CurvedPolygon object * * 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 + * Evaluate the vector field line integral with Gauss-Legendre quadrature * - * \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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning its numeric type + * \tparam FuncRetType The CurveType's numeric type * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. * \param [in] npts the number of quadrature points to evaluate the line integral * on each edge of the CurvedPolygon + * \pre Lambda must return the CurveTypes's vector type * \return the value of the integral */ -template -double evaluate_vector_line_integral(const primal::CurvedPolygon cpoly, - Lambda&& vector_integrand, - int npts) +template +FuncRetType evaluate_vector_line_integral(const CurvedPolygon cpoly, + Lambda&& vector_integrand, + int npts) { - double total_integral = 0.0; + FuncRetType total_integral = FuncRetType {}; for(int i = 0; i < cpoly.numEdges(); i++) { // Compute the line integral along each component. @@ -151,17 +183,20 @@ 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 vector-field line integral on a single generic curve * - * \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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning its numeric type + * \tparam FuncRetType The CurveType's numeric type * \param [in] c the generic curve object * \param [in] vector_integrand the lambda function representing the integrand. * \param [in] npts the number of quadrature nodes + * + * \pre Lambda must return the CurveTypes's vector type * \return the value of the integral */ -template -double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integrand, int npts) +template +FuncRetType evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integrand, int npts) { return detail::evaluate_vector_line_integral_component(c, std::forward(vector_integrand), @@ -171,8 +206,9 @@ double evaluate_vector_line_integral(const CurveType& c, Lambda&& vector_integra /*! * \brief Evaluate a line integral on an array of generic curves on a vector field * - * \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 + * \tparam Lambda A callable type taking a CurveType's PointType and returning its numeric type + * \tparam FuncRetType The CurveType's numeric type * \param [in] carray The array of generic curve objects * \param [in] vector_integrand the lambda function representing the integrand. * \param [in] npts the number of quadrature nodes per curve per knot span @@ -182,12 +218,12 @@ 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, - Lambda&& vector_integrand, - int npts) +template +FuncRetType evaluate_vector_line_integral(const axom::Array& carray, + Lambda&& vector_integrand, + int npts) { - double total_integral = 0.0; + FuncRetType total_integral = FuncRetType {}; for(int i = 0; i < carray.size(); i++) { total_integral += @@ -198,6 +234,10 @@ double evaluate_vector_line_integral(const axom::Array& carray, return total_integral; } +//@} + +///@{ +/// \name Evaluates scalar-field 2D area integrals for functions f : R^2 -> R^m /*! * \brief Evaluate an integral on the interior of a CurvedPolygon object. @@ -209,42 +249,52 @@ double evaluate_vector_line_integral(const axom::Array& carray, * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. * - * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \tparam Lambda A callable type taking a CurveType's PointType and returning an integrable type * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the geometry - * \param [in] cs the array of Bezier curve objects that bound the region - * \param [in] scalar_integrand the lambda function representing the integrand. + * \tparam LambdaRetType A type which supports addition and scalar multiplication + * \param [in] cpoly the CurvedPolygon object + * \param [in] 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 * \return the value of the integral */ -template -double evaluate_area_integral(const primal::CurvedPolygon& cpoly, - Lambda&& scalar_integrand, - int npts_Q, - int npts_P = 0) +template > +LambdaRetType evaluate_area_integral(const primal::CurvedPolygon& cpoly, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) { + using T = typename CurveType::NumericType; + + static_assert( + detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda function " + "return type"); + if(npts_P <= 0) { npts_P = npts_Q; } // Use minimum y-coord of control nodes as lower bound for integration - double int_lb = cpoly[0][0][1]; + T lower_bound_y = cpoly[0][0][1]; for(int i = 0; i < cpoly.numEdges(); i++) { for(int j = 1; j < cpoly[i].getOrder() + 1; j++) { - int_lb = std::min(int_lb, cpoly[i][j][1]); + lower_bound_y = std::min(lower_bound_y, cpoly[i][j][1]); } } // Evaluate the antiderivative line integral along each component - double total_integral = 0.0; + LambdaRetType total_integral = LambdaRetType {}; for(int i = 0; i < cpoly.numEdges(); i++) { total_integral += detail::evaluate_area_integral_component(cpoly[i], - std::forward(scalar_integrand), - int_lb, + std::forward(integrand), + lower_bound_y, npts_Q, npts_P); } @@ -257,10 +307,11 @@ double evaluate_area_integral(const primal::CurvedPolygon& cpoly, * * See above definition for details. * - * \tparam Lambda A callable type taking an axom::Point type returning a numeric type + * \tparam Lambda A callable type taking a CurveType's PointType and returning an integrable type * \tparam CurveType The BezierCurve, NURBSCurve, or NURBSCurveGWNCache which represents the geometry + * \tparam LambdaRetType A type which supports addition and scalar multiplication * \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] 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 * @@ -268,12 +319,21 @@ double evaluate_area_integral(const primal::CurvedPolygon& cpoly, * * \return the value of the integral */ -template -double evaluate_area_integral(const axom::Array& carray, - Lambda&& scalar_integrand, - int npts_Q, - int npts_P = 0) +template > +LambdaRetType evaluate_area_integral(const axom::Array& carray, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) { + using T = typename CurveType::NumericType; + + static_assert( + detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda function " + "return type"); + if(npts_P <= 0) { npts_P = npts_Q; @@ -281,38 +341,36 @@ double evaluate_area_integral(const axom::Array& carray, if(carray.empty()) { - return 0.0; + return LambdaRetType {}; } // Use minimum y-coord of control nodes as lower bound for integration - double int_lb = carray[0][0][1]; + T lower_bound_y = carray[0][0][1]; for(int i = 0; i < carray.size(); i++) { for(int j = 1; j < carray[i].getNumControlPoints(); j++) { - int_lb = std::min(int_lb, carray[i][j][1]); + lower_bound_y = std::min(lower_bound_y, carray[i][j][1]); } } // Evaluate the antiderivative line integral along each component - double total_integral = 0.0; + LambdaRetType total_integral = LambdaRetType {}; for(int i = 0; i < carray.size(); i++) { - auto beziers = carray[i].extractBezier(); - - for(int j = 0; j < beziers.size(); j++) + for(const auto& bez : carray[i].extractBezier()) { - total_integral += - detail::evaluate_area_integral_component(beziers[j], - std::forward(scalar_integrand), - int_lb, - npts_Q, - npts_P); + total_integral += detail::evaluate_area_integral_component(bez, + std::forward(integrand), + lower_bound_y, + npts_Q, + npts_P); } } return total_integral; } +//@} } // namespace primal } // end namespace axom diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index dc814e2f13..f129a4c049 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -91,6 +91,74 @@ TEST(primal_integral, evaluate_area_integral) EXPECT_NEAR(evaluate_area_integral(square, const_integrand, npts), 1.0, abs_tol); } +TEST(primal_integral, evaluate_area_integral_aggregate) +{ + using Point2D = primal::Point; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + using ReturnType = primal::Vector; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Define anonymous function for testing component-wise integration + // of an function with a vector return type. + // This is equivalent to evaluating three separate integrands + // without unnecessarily repeating geometric processing. + auto aggregate_integrand = [](Point2D x) -> ReturnType { + return ReturnType {1.0, x[0] * x[1] * x[1], 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); + + 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); + + BCurve triangle_edges[] = {tri1, tri2, tri3}; + CPolygon triangle(triangle_edges, 3); + + // Compare against hand computed/high-precision calculated values + auto observed = evaluate_area_integral(triangle, aggregate_integrand, npts); + auto expected = ReturnType {0.5, 1.0 / 60.0, 0.0415181074232}; + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], 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); + + 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); + + // Compare against hand computed/high-precision calculated values + observed = evaluate_area_integral(parabola_shape, aggregate_integrand, npts); + expected = ReturnType {8.0 / 3.0, 64.0 / 105.0, 0.0}; + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], abs_tol); + } + + // Ensure compatibility with curved polygons + BCurve pedges[2] = {para1, para2}; + CPolygon parabola_polygon(pedges, 2); + observed = evaluate_area_integral(parabola_polygon, aggregate_integrand, npts); + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], abs_tol); + } +} + TEST(primal_integral, evaluate_line_integral_scalar) { using Point2D = primal::Point; @@ -113,14 +181,10 @@ TEST(primal_integral, evaluate_line_integral_scalar) // 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_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), + EXPECT_NEAR(evaluate_line_integral(parabola_segment, poly_integrand, npts), 37.8010703669, abs_tol); + EXPECT_NEAR(evaluate_line_integral(parabola_segment, transc_integrand, npts), 0.495907795678, abs_tol); @@ -140,13 +204,9 @@ TEST(primal_integral, evaluate_line_integral_scalar) 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), + EXPECT_NEAR(evaluate_line_integral(connected_curve, const_integrand, npts), 8.28968500196, abs_tol); + EXPECT_NEAR(evaluate_line_integral(connected_curve, poly_integrand, npts), -5.97565740064, abs_tol); + EXPECT_NEAR(evaluate_line_integral(connected_curve, transc_integrand, npts), -0.574992518405, abs_tol); @@ -154,17 +214,83 @@ TEST(primal_integral, evaluate_line_integral_scalar) 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), + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, const_integrand, npts), 6.05361702446, abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, poly_integrand, npts), + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, poly_integrand, npts), -6.34833539689, abs_tol); - EXPECT_NEAR(evaluate_scalar_line_integral(disconnected_curve, transc_integrand, npts), + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, transc_integrand, npts), -0.914161242161, abs_tol); } +TEST(primal_integral, evaluate_line_integral_aggregate) +{ + using Point2D = primal::Point; + using BCurve = primal::BezierCurve; + using CPolygon = primal::CurvedPolygon; + using ReturnType = primal::Vector; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + // Define anonymous function for testing component-wise integration + // of an function with a vector return type. + // This is equivalent to evaluating three separate integrands + // without unnecessarily repeating geometric processing. + auto aggregate_integrand = [](Point2D x) -> ReturnType { + return ReturnType {1.0, x[0] * x[1] * x[1], 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. + auto observed = evaluate_line_integral(parabola_segment, aggregate_integrand, npts); + auto expected = ReturnType {6.12572661998, 37.8010703669, 0.495907795678}; + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], 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); + + observed = evaluate_line_integral(connected_curve, aggregate_integrand, npts); + expected = ReturnType {8.28968500196, -5.97565740064, -0.574992518405}; + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], abs_tol); + } + + // Test algorithm on disconnected curves + BCurve disconnected_curve_edges[] = {cubic_segment, quadratic_segment}; + CPolygon disconnected_curve(disconnected_curve_edges, 2); + + observed = evaluate_line_integral(disconnected_curve, aggregate_integrand, npts); + expected = ReturnType {6.05361702446, -6.34833539689, -0.914161242161}; + for(int i = 0; i < 3; ++i) + { + EXPECT_NEAR(observed[i], expected[i], abs_tol); + } +} + TEST(primal_integral, evaluate_line_integral_vector) { using Point2D = primal::Point; @@ -256,12 +382,8 @@ TEST(primal_integral, evaluate_integral_3D) }; // 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); + EXPECT_NEAR(evaluate_line_integral(spatial_arc, const_integrand, npts), 4.09193268998, abs_tol); + EXPECT_NEAR(evaluate_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); @@ -309,12 +431,8 @@ TEST(primal_integral, evaluate_integral_rational) 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); + EXPECT_NEAR(evaluate_line_integral(ellipse_arc, const_integrand, npts), 2.42211205514, abs_tol); + EXPECT_NEAR(evaluate_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), @@ -370,12 +488,8 @@ TEST(primal_integral, evaluate_integral_nurbs) 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_line_integral(ellipse_arc, const_integrand, npts), 2.42211205514, abs_tol); + EXPECT_NEAR(evaluate_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, @@ -498,12 +612,8 @@ TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) 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_line_integral(ellipse_arc, const_integrand, npts), 2.42211205514, abs_tol); + EXPECT_NEAR(evaluate_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, diff --git a/src/axom/quest/LinearizeCurves.cpp b/src/axom/quest/LinearizeCurves.cpp index 8f6b0fcc0a..a357ed4246 100644 --- a/src/axom/quest/LinearizeCurves.cpp +++ b/src/axom/quest/LinearizeCurves.cpp @@ -150,7 +150,7 @@ double computeArcLength(const LinearizeCurves::NURBSCurve& nurbs, int npts) for(const auto& bezier : nurbs.extractBezier()) { arcLength += - primal::evaluate_scalar_line_integral(bezier, [](PointType /*x*/) -> double { return 1.; }, npts); + primal::evaluate_line_integral(bezier, [](PointType /*x*/) -> double { return 1.; }, npts); } return arcLength; } diff --git a/src/axom/quest/io/STEPReader.cpp b/src/axom/quest/io/STEPReader.cpp index 9aaafd4590..81fb049a7f 100644 --- a/src/axom/quest/io/STEPReader.cpp +++ b/src/axom/quest/io/STEPReader.cpp @@ -934,13 +934,14 @@ class StepFileProcessor // Ensure that the trimming curves form ccw loops if(patch.isTrimmed()) { - auto area_field = [](PointType x) -> VectorType { - return primal::Vector({-0.5 * x[1], 0.5 * x[0]}); - }; - + // Evaluate the signed area over the patch trimming curves constexpr int n_quad_pts = 20; - auto area = - primal::evaluate_vector_line_integral(patch.getTrimmingCurves(), area_field, n_quad_pts); + auto area = primal::evaluate_vector_line_integral( + patch.getTrimmingCurves(), + [](PointType x) -> VectorType { + return primal::Vector({-0.5 * x[1], 0.5 * x[0]}); + }, + n_quad_pts); // Signed areas should be positive if(area < 0)