diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 4b0daf21cc..0c94b673d5 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -19,6 +19,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ ## [Unreleased] - Release date yyyy-mm-dd ### Added +- Primal: Functions to evaluate surface and volume integrals over collections of `BezierPatch` and `NURBSPatch` objects. +- Quest: An example to evaluate surface and volume integrals over a STEP model + and to fit a ellipsoid or oriented bounding box to a model based on its low order moments. ### Removed @@ -54,7 +57,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Core: the `ArrayView` class was modified so it defers initializing an internal allocator id via Umpire, if present. This prevents excessive calls to Umpire, which are not needed in all use cases. -- Quest: A compilation problem with `-DAXOM_NO_INT64_T=1` was fixed. +- Quest: Fixed a compilation problem with `-DAXOM_NO_INT64_T=1` - Quest: `STEPReader` now catches additional edge cases related to orientation of OpenCascade primitives. ## [Version 0.13.0] - Release date 2026-02-05 diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index ca3e591070..6fb85e2205 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -54,6 +54,8 @@ set( primal_headers operators/clip.hpp operators/closest_point.hpp operators/evaluate_integral.hpp + operators/evaluate_integral_curve.hpp + operators/evaluate_integral_surface.hpp operators/intersect.hpp operators/intersection_volume.hpp operators/orientation.hpp @@ -71,6 +73,8 @@ set( primal_headers operators/detail/clip_impl.hpp operators/detail/compute_moments_impl.hpp operators/detail/evaluate_integral_impl.hpp + operators/detail/evaluate_integral_curve_impl.hpp + operators/detail/evaluate_integral_surface_impl.hpp operators/detail/fuzzy_comparators.hpp operators/detail/intersect_bezier_impl.hpp operators/detail/intersect_bounding_box_impl.hpp diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp index 52a61218b8..f543d3106c 100644 --- a/src/axom/primal/geometry/CurvedPolygon.hpp +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -23,7 +23,7 @@ #include "axom/primal/geometry/BoundingBox.hpp" // For NURBSCurveGWNCache objects -#include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" +#include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" #include #include diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index eed10b824d..bf4e8b233d 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -26,6 +26,7 @@ #include "axom/primal/geometry/OrientedBoundingBox.hpp" #include "axom/primal/operators/squared_distance.hpp" +#include "axom/primal/operators/evaluate_integral_curve.hpp" #include "axom/primal/operators/detail/winding_number_2d_impl.hpp" #include "axom/primal/operators/detail/intersect_bezier_impl.hpp" @@ -4192,4 +4193,4 @@ template struct axom::fmt::formatter> : ostream_formatter { }; -#endif // AXOM_PRIMAL_NURBSPATCH_HPP_ \ No newline at end of file +#endif // AXOM_PRIMAL_NURBSPATCH_HPP_ diff --git a/src/axom/primal/operators/detail/evaluate_integral_curve_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_curve_impl.hpp new file mode 100644 index 0000000000..69b7ddf0e9 --- /dev/null +++ b/src/axom/primal/operators/detail/evaluate_integral_curve_impl.hpp @@ -0,0 +1,296 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file evaluate_integral_curve_impl.hpp + * + * \brief Implementation helpers for curve/area integral evaluation. + * + * \note This header intentionally avoids including surface/volume (patch) + * dependencies to prevent circular include chains (e.g. with NURBSPatch). + */ + +#ifndef PRIMAL_EVAL_INTEGRAL_CURVE_IMPL_HPP_ +#define PRIMAL_EVAL_INTEGRAL_CURVE_IMPL_HPP_ + +// Axom includes +#include "axom/core.hpp" +#include "axom/config.hpp" +#include "axom/slic.hpp" + +#include "axom/core/utilities/Utilities.hpp" +#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" + +// C++ includes +#include +#include +#include + +namespace axom +{ +namespace primal +{ +namespace detail +{ +namespace internal +{ +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 Scalar-field line integrals + +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); + + LambdaRetType full_quadrature = LambdaRetType {}; + for(int q = 0; q < npts; q++) + { + auto x_q = c.evaluate(quad.node(q)); + auto dx_q = c.dt(quad.node(q)); + + full_quadrature += quad.weight(q) * integrand(x_q) * dx_q.norm(); + } + + return full_quadrature; +} + +template ::PointType>> +inline LambdaRetType evaluate_line_integral_component(const NURBSCurve& n, + Lambda&& integrand, + const int npts) +{ + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) + { + total_integral += + detail::evaluate_line_integral_component(bez, std::forward(integrand), npts); + } + return total_integral; +} + +template ::PointType>> +inline LambdaRetType evaluate_line_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& integrand, + const int npts) +{ + LambdaRetType total_integral = LambdaRetType {}; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_line_integral_component(this_bezier_data.getCurve(), + std::forward(integrand), + npts); + } + + return total_integral; +} +//@} + +///@{ +/// \name Vector-field line integrals + +template +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); + + T full_quadrature = T {}; + for(int q = 0; q < npts; q++) + { + auto x_q = c.evaluate(quad.node(q)); + auto dx_q = c.dt(quad.node(q)); + auto func_val = vector_integrand(x_q); + + full_quadrature += quad.weight(q) * Vector::dot_product(func_val, dx_q); + } + + return full_quadrature; +} + +template +inline T evaluate_vector_line_integral_component(const primal::NURBSCurve& n, + Lambda&& vector_integrand, + const int npts) +{ + T total_integral = T {}; + for(const auto& bez : n.extractBezier()) + { + total_integral += + detail::evaluate_vector_line_integral_component(bez, + std::forward(vector_integrand), + npts); + } + return total_integral; +} + +template +inline T evaluate_vector_line_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& vector_integrand, + const int npts) +{ + T total_integral = T {}; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + 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; +} +//@} + +///@{ +/// \name Scalar-field 2D area integrals + +template ::PointType>> +inline 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); + + LambdaRetType full_quadrature = LambdaRetType {}; + for(int q = 0; q < npts_Q; q++) + { + auto x_q = c.evaluate(quad_Q.node(q)); + + for(int xi = 0; xi < npts_P; xi++) + { + auto x_qxi = Point({x_q[0], (x_q[1] - int_lb) * quad_P.node(xi) + int_lb}); + + 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; + } + } + + return full_quadrature; +} + +template ::PointType>> +inline LambdaRetType evaluate_area_integral_component(const primal::NURBSCurve& n, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) + { + total_integral += detail::evaluate_area_integral_component(bez, + std::forward(integrand), + int_lb, + npts_Q, + npts_P); + } + return total_integral; +} + +template ::PointType>> +inline RetType evaluate_area_integral_component(const NURBSCurveGWNCache& nc, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P) +{ + RetType total_integral = RetType {}; + for(int i = 0; i < nc.getNumKnotSpans(); ++i) + { + const auto& this_bezier_data = nc.getSubdivisionData(i, 0, 0); + total_integral += detail::evaluate_area_integral_component(this_bezier_data.getCurve(), + std::forward(integrand), + int_lb, + npts_Q, + npts_P); + } + + return total_integral; +} +//@} + +///@{ +/// \name Helper routines for patch-based integrals in 3D + +template +inline typename CurveType::NumericType curve_array_lower_bound_y(const axom::Array& carray) +{ + using T = typename CurveType::NumericType; + + SLIC_ASSERT(!carray.empty() && carray[0].getNumControlPoints() > 0); + + T lower_bound_y = carray[0][0][1]; + for(int i = 0; i < carray.size(); ++i) + { + for(int j = 0; j < carray[i].getNumControlPoints(); ++j) + { + lower_bound_y = axom::utilities::min(lower_bound_y, carray[i][j][1]); + } + } + + return lower_bound_y; +} + +//@} + +} // end namespace detail +} // end namespace primal +} // end namespace axom + +#endif diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 253af353ba..99e69eebbb 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -5,408 +5,18 @@ // SPDX-License-Identifier: (BSD-3-Clause) /*! - * \file evaluate_integral.hpp + * \file evaluate_integral_impl.hpp * - * \brief Consists of methods that evaluate scalar-field integrals on curves and - * regions defined by 2D curves, and vector-field integrals on curves + * \brief Header to include implementation files for evaluating integrals * - * 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. + * Splitting the implementation avoids circular include dependencies when patches + * internally evaluate curve-based integrals (e.g. via trimming curves). */ #ifndef PRIMAL_EVAL_INTEGRAL_IMPL_HPP_ #define PRIMAL_EVAL_INTEGRAL_IMPL_HPP_ -// Axom includes -#include "axom/config.hpp" // for compile-time configuration options -#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/primal/operators/detail/evaluate_integral_curve_impl.hpp" +#include "axom/primal/operators/detail/evaluate_integral_surface_impl.hpp" -#include "axom/core/numerics/quadrature.hpp" - -// C++ includes -#include -#include -#include - -namespace axom -{ -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 line integral on a single Bezier curve. - * - * Evaluate the line integral with Gauss-Legendre quadrature - * - * \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] 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 ::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 - LambdaRetType full_quadrature = LambdaRetType {}; - for(int q = 0; q < npts; q++) - { - // Get intermediate quadrature point - // at which to evaluate tangent vector - auto x_q = c.evaluate(quad.node(q)); - auto dx_q = c.dt(quad.node(q)); - - full_quadrature += quad.weight(q) * integrand(x_q) * dx_q.norm(); - } - - return full_quadrature; -} - -/*! - * \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 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] 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 ::PointType>> -inline LambdaRetType evaluate_line_integral_component(const NURBSCurve& n, - Lambda&& integrand, - const int npts) -{ - LambdaRetType total_integral = LambdaRetType {}; - for(const auto& bez : n.extractBezier()) - { - total_integral += - detail::evaluate_line_integral_component(bez, std::forward(integrand), npts); - } - return total_integral; -} - -/*! - * \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 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] 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 ::PointType>> -inline LambdaRetType evaluate_line_integral_component(const NURBSCurveGWNCache& nc, - Lambda&& integrand, - const int npts) -{ - 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_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 vector field line integral with Gauss-Legendre quadrature - * - * \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 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 - T full_quadrature = T {}; - for(int q = 0; q < npts; q++) - { - // Get intermediate quadrature point - // 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 = 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 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 T evaluate_vector_line_integral_component(const primal::NURBSCurve& n, - Lambda&& vector_integrand, - const int npts) -{ - T total_integral = T {}; - for(const auto& bez : n.extractBezier()) - { - total_integral += - detail::evaluate_vector_line_integral_component(bez, - 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 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 T evaluate_vector_line_integral_component(const NURBSCurveGWNCache& nc, - Lambda&& vector_integrand, - const int npts) -{ - 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 - 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; -} -//@} - -///@{ -/// \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. - * - * Intended to be called for each BezierCurve object in a curved polygon. - * Uses a Spectral Mesh-Free Quadrature derived from Green's theorem, evaluating - * the area integral as a line integral of the antiderivative over the curve. - * 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 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] 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 ::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/compute quadrature result - LambdaRetType full_quadrature = LambdaRetType {}; - for(int q = 0; q < npts_Q; q++) - { - // Get intermediate quadrature point - // on which to evaluate antiderivative - auto x_q = c.evaluate(quad_Q.node(q)); - - // Evaluate the antiderivative at x_q, add it to full quadrature - for(int xi = 0; xi < npts_P; xi++) - { - // Define interior quadrature points - auto x_qxi = Point({x_q[0], (x_q[1] - int_lb) * quad_P.node(xi) + int_lb}); - - 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; - } - } - - 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 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] 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 ::PointType>> -LambdaRetType evaluate_area_integral_component(const primal::NURBSCurve& n, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - LambdaRetType total_integral = LambdaRetType {}; - for(const auto& bez : n.extractBezier()) - { - total_integral += detail::evaluate_area_integral_component(bez, - std::forward(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 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] 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 ::PointType>> -inline RetType evaluate_area_integral_component(const NURBSCurveGWNCache& nc, - Lambda&& integrand, - double int_lb, - const int npts_Q, - const int npts_P) -{ - 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(integrand), - int_lb, - npts_Q, - npts_P); - } - - return total_integral; -} -//@} - -} // end namespace detail -} // end namespace primal -} // end namespace axom - -#endif \ No newline at end of file +#endif diff --git a/src/axom/primal/operators/detail/evaluate_integral_surface_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_surface_impl.hpp new file mode 100644 index 0000000000..78f9268de9 --- /dev/null +++ b/src/axom/primal/operators/detail/evaluate_integral_surface_impl.hpp @@ -0,0 +1,216 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file evaluate_integral_surface_impl.hpp + * + * \brief Implementation helpers for surface/volume integral evaluation. + */ + +#ifndef PRIMAL_EVAL_INTEGRAL_SURFACE_IMPL_HPP_ +#define PRIMAL_EVAL_INTEGRAL_SURFACE_IMPL_HPP_ + +// Axom includes +#include "axom/core.hpp" +#include "axom/config.hpp" + +#include "axom/core/utilities/Utilities.hpp" +#include "axom/primal/geometry/BezierPatch.hpp" +#include "axom/primal/geometry/NURBSPatch.hpp" +#include "axom/primal/operators/detail/evaluate_integral_curve_impl.hpp" + +// C++ includes +#include +#include + +namespace axom +{ +namespace primal +{ +namespace detail +{ +///@{ +/// \name Scalar-field surface/volume integral components + +template ::PointType>> +inline LambdaRetType evaluate_surface_integral_component(const primal::BezierPatch& b, + Lambda&& integrand, + const int npts) +{ + const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); + + LambdaRetType full_quadrature = LambdaRetType {}; + for(int qu = 0; qu < npts; ++qu) + { + for(int qv = 0; qv < npts; ++qv) + { + const auto x_q = b.evaluate(quad.node(qu), quad.node(qv)); + const auto n_q = b.normal(quad.node(qu), quad.node(qv)); + + full_quadrature += quad.weight(qu) * quad.weight(qv) * integrand(x_q) * n_q.norm(); + } + } + + return full_quadrature; +} + +template ::PointType>> +inline LambdaRetType evaluate_surface_integral_component(const primal::NURBSPatch& n, + Lambda&& integrand, + const int npts_Q, + const int npts_P) +{ + if(!n.isTrimmed()) + { + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) + { + total_integral += detail::evaluate_surface_integral_component(bez, integrand, npts_Q); + } + + return total_integral; + } + + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& split_patch : n.extractTrimmedBezier()) + { + const auto& curves = split_patch.getTrimmingCurves(); + if(curves.empty()) + { + continue; + } + + // clamp the lower bound to lie within the parametric space of the patch + const auto lower_bound_y = axom::utilities::clampVal(detail::curve_array_lower_bound_y(curves), + split_patch.getMinKnot_v(), + split_patch.getMaxKnot_v()); + for(int i = 0; i < curves.size(); ++i) + { + total_integral += detail::evaluate_area_integral_component( + curves[i], + [&split_patch, &integrand](Point uv) -> LambdaRetType { + const auto x_q = split_patch.evaluate(uv[0], uv[1]); + const auto n_q = split_patch.normal(uv[0], uv[1]); + return integrand(x_q) * n_q.norm(); + }, + lower_bound_y, + npts_Q, + npts_P); + } + } + + return total_integral; +} + +template ::PointType>> +inline LambdaRetType evaluate_volume_integral_component(const primal::BezierPatch& b, + Lambda&& integrand, + double int_lb, + const int npts_uv, + const int npts_z) +{ + const axom::numerics::QuadratureRule& quad_uv = axom::numerics::get_gauss_legendre(npts_uv); + const axom::numerics::QuadratureRule& quad_z = axom::numerics::get_gauss_legendre(npts_z); + + LambdaRetType full_quadrature = LambdaRetType {}; + for(int qu = 0; qu < npts_uv; ++qu) + { + for(int qv = 0; qv < npts_uv; ++qv) + { + const auto x_q = b.evaluate(quad_uv.node(qu), quad_uv.node(qv)); + const auto n_q = b.normal(quad_uv.node(qu), quad_uv.node(qv)); + + LambdaRetType antiderivative = LambdaRetType {}; + const T z_scale = x_q[2] - int_lb; + for(int qz = 0; qz < npts_z; ++qz) + { + const auto x_qz = Point({x_q[0], x_q[1], z_scale * quad_z.node(qz) + int_lb}); + antiderivative += quad_z.weight(qz) * z_scale * integrand(x_qz); + } + + full_quadrature += quad_uv.weight(qu) * quad_uv.weight(qv) * antiderivative * n_q[2]; + } + } + + return full_quadrature; +} + +template ::PointType>> +inline LambdaRetType evaluate_volume_integral_component(const primal::NURBSPatch& n, + Lambda&& integrand, + double int_lb, + const int npts_Q, + const int npts_P, + const int npts_Z) +{ + if(!n.isTrimmed()) + { + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& bez : n.extractBezier()) + { + total_integral += + detail::evaluate_volume_integral_component(bez, integrand, int_lb, npts_Q, npts_Z); + } + + return total_integral; + } + + const axom::numerics::QuadratureRule& quad_z = axom::numerics::get_gauss_legendre(npts_Z); + + LambdaRetType total_integral = LambdaRetType {}; + for(const auto& split_patch : n.extractTrimmedBezier()) + { + const auto& curves = split_patch.getTrimmingCurves(); + if(curves.empty()) + { + continue; + } + + const auto lower_bound_y = axom::utilities::clampVal(detail::curve_array_lower_bound_y(curves), + split_patch.getMinKnot_v(), + split_patch.getMaxKnot_v()); + for(int i = 0; i < curves.size(); ++i) + { + total_integral += detail::evaluate_area_integral_component( + curves[i], + [&split_patch, &integrand, &int_lb, &quad_z, &npts_Z](Point uv) -> LambdaRetType { + const auto x_q = split_patch.evaluate(uv[0], uv[1]); + const auto n_q = split_patch.normal(uv[0], uv[1]); + + LambdaRetType antiderivative = LambdaRetType {}; + const T z_scale = x_q[2] - int_lb; + for(int qz = 0; qz < npts_Z; ++qz) + { + const auto x_qz = Point({x_q[0], x_q[1], z_scale * quad_z.node(qz) + int_lb}); + antiderivative += quad_z.weight(qz) * z_scale * integrand(x_qz); + } + + return antiderivative * n_q[2]; + }, + lower_bound_y, + npts_Q, + npts_P); + } + } + + return total_integral; +} + +//@} + +} // end namespace detail +} // end namespace primal +} // end namespace axom + +#endif diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index e059d17fe8..1bc47377a4 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -7,374 +7,20 @@ /*! * \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 + * \brief Header for including integral evaluation functions * - * 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. + * The API is split into: + * - evaluate_integral_curve.hpp (line/area) + * - evaluate_integral_surface.hpp (surface/volume) + * + * This split avoids circular include dependencies when patches internally + * evaluate curve-based integrals (e.g. via trimming curves). */ #ifndef PRIMAL_EVAL_INTEGRAL_HPP_ #define PRIMAL_EVAL_INTEGRAL_HPP_ -// Axom includes -#include "axom/config.hpp" - -#include "axom/core/utilities/Utilities.hpp" -#include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" - -// C++ includes -#include - -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 - * 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 curved polygon need not be connected. - * - * Evaluate the line integral with Gauss-Legendre quadrature - * - * \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] 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 > -LambdaRetType evaluate_line_integral(const primal::CurvedPolygon cpoly, - Lambda&& integrand, - int npts) -{ - 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_line_integral_component(cpoly[i], std::forward(integrand), npts); - } - - return total_integral; -} - -/*! - * \brief Evaluate a line integral along the boundary of a generic curve - * for a function with an arbitrary return 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] integrand the lambda function representing the integrand. - * \param [in] npts the number of quadrature nodes - * \return the value of the integral - */ -template > -LambdaRetType evaluate_line_integral(const CurveType& c, Lambda&& integrand, int 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 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] 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 - * is computed using npts on each segment - * - * \return the value of the integral - */ -template > -LambdaRetType evaluate_line_integral(const axom::Array& carray, Lambda&& integrand, int npts) -{ - 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_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 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 vector field line integral with Gauss-Legendre quadrature - * - * \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 -FuncRetType evaluate_vector_line_integral(const CurvedPolygon cpoly, - Lambda&& vector_integrand, - int npts) -{ - FuncRetType total_integral = FuncRetType {}; - for(int i = 0; i < cpoly.numEdges(); i++) - { - // Compute the line integral along each component. - total_integral += - detail::evaluate_vector_line_integral_component(cpoly[i], - std::forward(vector_integrand), - npts); - } - - return total_integral; -} - -/*! - * \brief Evaluate a vector-field line integral on a single generic curve - * - * \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 -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), - npts); -} - -/*! - * \brief Evaluate a line integral on an array of generic curves on a vector field - * - * \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 - * - * \note Each 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 -FuncRetType evaluate_vector_line_integral(const axom::Array& carray, - Lambda&& vector_integrand, - int npts) -{ - FuncRetType total_integral = FuncRetType {}; - for(int i = 0; i < carray.size(); i++) - { - total_integral += - detail::evaluate_vector_line_integral_component(carray[i], - std::forward(vector_integrand), - npts); - } - - 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. - * - * Evaluates the integral using a Spectral Mesh-Free Quadrature derived from - * Green's theorem, evaluating the area integral as a line integral of the - * antiderivative over each component curve. - * - * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar - * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. - * - * \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] 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 > -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 - 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++) - { - lower_bound_y = axom::utilities::min(lower_bound_y, cpoly[i][j][1]); - } - } - - // Evaluate the antiderivative line integral along each component - LambdaRetType total_integral = LambdaRetType {}; - for(int i = 0; i < cpoly.numEdges(); i++) - { - total_integral += detail::evaluate_area_integral_component(cpoly[i], - std::forward(integrand), - lower_bound_y, - npts_Q, - npts_P); - } - - return total_integral; -} - -/*! - * \brief Evaluate an integral on the interior of a region bound by 2D curves - * - * See above definition for details. - * - * \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] 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 The numerical result is only meaningful if the curves enclose a region - * - * \return the value of the integral - */ -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; - } - - if(carray.empty()) - { - return LambdaRetType {}; - } - - // Use minimum y-coord of control nodes as lower bound for integration - 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++) - { - lower_bound_y = axom::utilities::min(lower_bound_y, carray[i][j][1]); - } - } - - // Evaluate the antiderivative line integral along each component - LambdaRetType total_integral = LambdaRetType {}; - for(int i = 0; i < carray.size(); i++) - { - for(const auto& bez : carray[i].extractBezier()) - { - 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 +#include "axom/primal/operators/evaluate_integral_curve.hpp" +#include "axom/primal/operators/evaluate_integral_surface.hpp" #endif diff --git a/src/axom/primal/operators/evaluate_integral_curve.hpp b/src/axom/primal/operators/evaluate_integral_curve.hpp new file mode 100644 index 0000000000..0bbd2be4d0 --- /dev/null +++ b/src/axom/primal/operators/evaluate_integral_curve.hpp @@ -0,0 +1,367 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file evaluate_integral_curve.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 the 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 D. Gunderman et al. + * https://doi.org/10.1016/j.cad.2020.102944 + */ + +#ifndef PRIMAL_EVAL_INTEGRAL_CURVE_HPP_ +#define PRIMAL_EVAL_INTEGRAL_CURVE_HPP_ + +// Axom includes +#include "axom/core.hpp" +#include "axom/config.hpp" + +#include "axom/core/utilities/Utilities.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/operators/detail/evaluate_integral_curve_impl.hpp" + +// C++ includes +#include + +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 + * 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 curved polygon need not be connected. + * + * Evaluate the line integral with Gauss-Legendre quadrature + * + * \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] 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 > +LambdaRetType evaluate_line_integral(const primal::CurvedPolygon& cpoly, + Lambda&& integrand, + int npts) +{ + 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++) + { + total_integral += + detail::evaluate_line_integral_component(cpoly[i], std::forward(integrand), npts); + } + + return total_integral; +} + +/*! + * \brief Evaluate a line integral along the boundary of a generic curve + * for a function with an arbitrary return 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] integrand the lambda function representing the integrand. + * \param [in] npts the number of quadrature nodes + * \return the value of the integral + */ +template > +LambdaRetType evaluate_line_integral(const CurveType& c, Lambda&& integrand, int 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 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] 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 + * is computed using npts on each segment + * + * \return the value of the integral + */ +template > +LambdaRetType evaluate_line_integral(const axom::Array& carray, Lambda&& integrand, int npts) +{ + 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_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 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 vector field line integral with Gauss-Legendre quadrature + * + * \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 +FuncRetType evaluate_vector_line_integral(const CurvedPolygon& cpoly, + Lambda&& vector_integrand, + int npts) +{ + FuncRetType total_integral = FuncRetType {}; + for(int i = 0; i < cpoly.numEdges(); i++) + { + total_integral += + detail::evaluate_vector_line_integral_component(cpoly[i], + std::forward(vector_integrand), + npts); + } + + return total_integral; +} + +/*! + * \brief Evaluate a vector-field line integral on a single generic curve + * + * \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 +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), + npts); +} + +/*! + * \brief Evaluate a line integral on an array of generic curves on a vector field + * + * \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 + * + * \note Each 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 +FuncRetType evaluate_vector_line_integral(const axom::Array& carray, + Lambda&& vector_integrand, + int npts) +{ + FuncRetType total_integral = FuncRetType {}; + for(int i = 0; i < carray.size(); i++) + { + total_integral += + detail::evaluate_vector_line_integral_component(carray[i], + std::forward(vector_integrand), + npts); + } + + 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. + * + * Evaluates the integral using a Spectral Mesh-Free Quadrature derived from + * Green's theorem, evaluating the area integral as a line integral of the + * antiderivative over each component curve. + * + * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar + * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. + * + * \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] 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 > +LambdaRetType evaluate_area_integral(const primal::CurvedPolygon& cpoly, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + 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; + } + + LambdaRetType total_integral = LambdaRetType {}; + if(cpoly.numEdges() == 0) + { + return total_integral; + } + + auto lower_bound_y = cpoly[0][0][1]; + for(int i = 0; i < cpoly.numEdges(); ++i) + { + for(int j = 0; j < cpoly[i].getNumControlPoints(); ++j) + { + lower_bound_y = axom::utilities::min(lower_bound_y, cpoly[i][j][1]); + } + } + + for(int i = 0; i < cpoly.numEdges(); ++i) + { + total_integral += detail::evaluate_area_integral_component(cpoly[i], + std::forward(integrand), + lower_bound_y, + npts_Q, + npts_P); + } + + return total_integral; +} + +/*! + * \brief Evaluate an integral on the interior of a region bound by 2D curves + * + * See above definition for details. + * + * \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] 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 The numerical result is only meaningful if the curves enclose a region + * + * \return the value of the integral + */ +template > +LambdaRetType evaluate_area_integral(const axom::Array& carray, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + 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; + } + + LambdaRetType total_integral = LambdaRetType {}; + if(carray.empty()) + { + return total_integral; + } + + const auto lower_bound_y = detail::curve_array_lower_bound_y(carray); + + for(int i = 0; i < carray.size(); ++i) + { + total_integral += detail::evaluate_area_integral_component(carray[i], + std::forward(integrand), + lower_bound_y, + npts_Q, + npts_P); + } + + return total_integral; +} +//@} + +} // namespace primal +} // end namespace axom + +#endif diff --git a/src/axom/primal/operators/evaluate_integral_surface.hpp b/src/axom/primal/operators/evaluate_integral_surface.hpp new file mode 100644 index 0000000000..8ade0b431e --- /dev/null +++ b/src/axom/primal/operators/evaluate_integral_surface.hpp @@ -0,0 +1,409 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file evaluate_integral_surface.hpp + * + * \brief Consists of methods that evaluate surface and volume integrals on + * surface patches and regions defined by collections of patches. + * + * All integrals are evaluated numerically with Gauss-Legendre quadrature + * + * 3D integrals computed with "High-accuracy mesh-free quadrature + * for trimmed parametric surfaces and volumes" by D. Gunderman et al. + * https://doi.org/10.1016/j.cad.2021.103093 + */ + +#ifndef PRIMAL_EVAL_INTEGRAL_SURFACE_HPP_ +#define PRIMAL_EVAL_INTEGRAL_SURFACE_HPP_ + +// Axom includes +#include "axom/core.hpp" +#include "axom/config.hpp" + +#include "axom/core/utilities/Utilities.hpp" +#include "axom/primal/geometry/BezierPatch.hpp" +#include "axom/primal/geometry/NURBSPatch.hpp" +#include "axom/primal/operators/detail/evaluate_integral_surface_impl.hpp" + +namespace axom +{ +namespace primal +{ +///@{ +/// \name Evaluates scalar-field surface integrals for functions f : R^3 -> R^m + +/*! + * \brief Evaluate a scalar surface integral on a single Bezier patch. + * + * Uses tensor-product Gauss-Legendre quadrature in the patch parameter space. + * + * \param [in] patch the Bezier patch + * \param [in] integrand callable representing the integrand + * \param [in] npts the number of quadrature points in each parametric direction + * + * \pre The patch parameterization must be valid on its full parameter domain. + */ +template ::PointType>> +LambdaRetType evaluate_surface_integral(const primal::BezierPatch& patch, + Lambda&& integrand, + int npts) +{ + static_assert(detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda " + "function return type"); + + return detail::evaluate_surface_integral_component(patch, std::forward(integrand), npts); +} + +/*! + * \brief Evaluate a scalar surface integral on a single NURBS patch. + * + * Untrimmed patches are integrated by Bezier extraction followed by tensor-product + * Gauss-Legendre quadrature. Trimmed patches are integrated by reducing the + * parameter-space area integral to line integrals over the trimming curves. + * + * \param [in] patch the NURBS patch + * \param [in] integrand callable representing the integrand + * \param [in] npts_Q the number of quadrature points on each trimming curve or + * in each parametric direction for untrimmed Bezier pieces + * \param [in] npts_P the number of quadrature points used for numerical + * antidifferentiation in parameter space + * + * \pre The patch parameterization must be valid on its full parameter domain. + * \pre If the patch is trimmed, its trimming curves must bound the intended + * interior region in parameter space. + */ +template ::PointType>> +LambdaRetType evaluate_surface_integral(const primal::NURBSPatch& patch, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + 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; + } + + return detail::evaluate_surface_integral_component(patch, + std::forward(integrand), + npts_Q, + npts_P); +} + +/*! + * \brief Evaluate a scalar surface integral on a collection of Bezier patches. + * + * The result is the sum of the surface integrals over each patch in the array. + * + * \param [in] patches the patch collection + * \param [in] integrand callable representing the integrand + * \param [in] npts the number of quadrature points in each parametric direction + * + * \pre Each patch parameterization must be valid on its full parameter domain. + */ +template ::PointType>> +LambdaRetType evaluate_surface_integral(const axom::Array>& patches, + Lambda&& integrand, + int npts) +{ + 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 < patches.size(); ++i) + { + total_integral += detail::evaluate_surface_integral_component(patches[i], integrand, npts); + } + + return total_integral; +} + +/*! + * \brief Evaluate a scalar surface integral on a collection of NURBS patches. + * + * The result is the sum of the surface integrals over each patch in the array. + * + * \param [in] patches the patch collection + * \param [in] integrand callable representing the integrand + * \param [in] npts_Q the number of quadrature points on each trimming curve or + * in each parametric direction for untrimmed Bezier pieces + * \param [in] npts_P the number of quadrature points used for numerical + * antidifferentiation in parameter space + * + * \pre Each patch parameterization must be valid on its full parameter domain. + * \pre Any trimmed patch in the array must have trimming curves that bound its + * intended interior region in parameter space. + */ +template ::PointType>> +LambdaRetType evaluate_surface_integral(const axom::Array>& patches, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + 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; + } + + LambdaRetType total_integral = LambdaRetType {}; + for(int i = 0; i < patches.size(); ++i) + { + total_integral += + detail::evaluate_surface_integral_component(patches[i], integrand, npts_Q, npts_P); + } + + return total_integral; +} +//@} + +///@{ +/// \name Evaluates scalar-field volume integrals for functions f : R^3 -> R^m + +/*! + * \brief Evaluate a scalar volume-integral contribution from a single Bezier patch. + * + * This applies the Stokes-based reduction used for the full volume algorithm to + * one patch using a z-directed numerical antiderivative. + * + * \param [in] patch the Bezier patch + * \param [in] integrand callable representing the integrand + * \param [in] lower_bound_z the shared lower integration bound used for the + * z-directed antiderivative across the full boundary + * \param [in] npts_uv the number of quadrature points in each patch parameter direction + * \param [in] npts_z the number of quadrature points used for numerical + * antidifferentiation in z + * + * \pre The patch parameterization must be valid on its full parameter domain. + * \pre The returned value is geometrically meaningful as a volume contribution only + * when this patch is interpreted as part of a closed, consistently oriented + * boundary that uses the same lower integration bound. + */ +template ::PointType>> +LambdaRetType evaluate_volume_integral(const primal::BezierPatch& patch, + Lambda&& integrand, + T lower_bound_z, + int npts_uv, + int npts_z = 0) +{ + static_assert(detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda " + "function return type"); + + if(npts_z <= 0) + { + npts_z = npts_uv; + } + + return detail::evaluate_volume_integral_component(patch, + std::forward(integrand), + lower_bound_z, + npts_uv, + npts_z); +} + +/*! + * \brief Evaluate a scalar volume-integral contribution from a single NURBS patch. + * + * Trimmed patches use the same Green/Stokes reduction as the surface-integral + * algorithm, combined with a z-directed numerical antiderivative for the volume + * reduction. + * + * \param [in] patch the NURBS patch + * \param [in] integrand callable representing the integrand + * \param [in] lower_bound_z the shared lower integration bound used for the + * z-directed antiderivative across the full boundary + * \param [in] npts_Q the number of quadrature points on each trimming curve or + * in each parametric direction for untrimmed Bezier pieces + * \param [in] npts_P the number of quadrature points used for numerical + * antidifferentiation in parameter space + * \param [in] npts_Z the number of quadrature points used for numerical + * antidifferentiation in z + * + * \pre The patch parameterization must be valid on its full parameter domain. + * \pre If the patch is trimmed, its trimming curves must bound the intended + * interior region in parameter space. + * \pre The returned value is geometrically meaningful as a volume contribution only + * when this patch is interpreted as part of a closed, consistently oriented + * boundary that uses the same lower integration bound. + */ +template ::PointType>> +LambdaRetType evaluate_volume_integral(const primal::NURBSPatch& patch, + Lambda&& integrand, + T lower_bound_z, + int npts_Q, + int npts_P = 0, + int npts_Z = 0) +{ + 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; + } + if(npts_Z <= 0) + { + npts_Z = npts_Q; + } + + return detail::evaluate_volume_integral_component(patch, + std::forward(integrand), + lower_bound_z, + npts_Q, + npts_P, + npts_Z); +} + +/*! + * \brief Evaluate a scalar volume integral over a collection of Bezier patches. + * + * The result is obtained by summing the Stokes-based contribution from each + * patch in the collection. + * + * \param [in] patches the patch collection + * \param [in] integrand callable representing the integrand + * \param [in] npts_uv the number of quadrature points in each patch parameter direction + * \param [in] npts_z the number of quadrature points used for numerical + * antidifferentiation in z + * + * \pre Each patch parameterization must be valid on its full parameter domain. + * \pre The patch collection must represent a closed, consistently oriented + * boundary of the target volume. + */ +template ::PointType>> +LambdaRetType evaluate_volume_integral(const axom::Array>& patches, + Lambda&& integrand, + int npts_uv, + int npts_z = 0) +{ + static_assert(detail::internal::is_integrable_v, + "evaluate_integral methods require addition and scalar multiplication for lambda " + "function return type"); + + if(npts_z <= 0) + { + npts_z = npts_uv; + } + + if(patches.empty()) + { + return LambdaRetType {}; + } + + T lower_bound_z = patches[0].boundingBox().getMin()[2]; + for(int i = 1; i < patches.size(); ++i) + { + lower_bound_z = axom::utilities::min(lower_bound_z, patches[i].boundingBox().getMin()[2]); + } + + LambdaRetType total_integral = LambdaRetType {}; + for(int i = 0; i < patches.size(); ++i) + { + total_integral += + detail::evaluate_volume_integral_component(patches[i], integrand, lower_bound_z, npts_uv, npts_z); + } + + return total_integral; +} + +/*! + * \brief Evaluate a scalar volume integral over a collection of NURBS patches. + * + * The result is obtained by summing the Stokes-based contribution from each + * patch in the collection. + * + * \param [in] patches the patch collection + * \param [in] integrand callable representing the integrand + * \param [in] npts_Q the number of quadrature points on each trimming curve or + * in each parametric direction for untrimmed Bezier pieces + * \param [in] npts_P the number of quadrature points used for numerical + * antidifferentiation in parameter space + * \param [in] npts_Z the number of quadrature points used for numerical + * antidifferentiation in z + * + * \pre Each patch parameterization must be valid on its full parameter domain. + * \pre Any trimmed patch in the array must have trimming curves that bound its + * intended interior region in parameter space. + * \pre The patch collection must represent a closed, consistently oriented + * boundary of the target volume. + */ +template ::PointType>> +LambdaRetType evaluate_volume_integral(const axom::Array>& patches, + Lambda&& integrand, + int npts_Q, + int npts_P = 0, + int npts_Z = 0) +{ + 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; + } + if(npts_Z <= 0) + { + npts_Z = npts_Q; + } + + if(patches.empty()) + { + return LambdaRetType {}; + } + + T lower_bound_z = patches[0].boundingBox().getMin()[2]; + for(int i = 1; i < patches.size(); ++i) + { + lower_bound_z = axom::utilities::min(lower_bound_z, patches[i].boundingBox().getMin()[2]); + } + + LambdaRetType total_integral = LambdaRetType {}; + for(int i = 0; i < patches.size(); ++i) + { + total_integral += detail::evaluate_volume_integral_component(patches[i], + integrand, + lower_bound_z, + npts_Q, + npts_P, + npts_Z); + } + + return total_integral; +} +//@} + +} // namespace primal +} // end namespace axom + +#endif diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 8724c9c4a1..82e3cdbbd4 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -5,6 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include "axom/primal.hpp" +#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" #include "axom/slic.hpp" #include "axom/fmt.hpp" #include @@ -564,6 +565,140 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) } } +TEST(primal_integral, evaluate_patch_surface_and_volume_integrals) +{ + using Point3D = primal::Point; + using BPatch = primal::BezierPatch; + + auto make_bilinear_patch = + [](const Point3D& p00, const Point3D& p10, const Point3D& p01, const Point3D& p11) { + Point3D control_points[] = {p00, p01, p10, p11}; + return BPatch(control_points, 1, 1); + }; + + axom::Array cube_faces(6); + cube_faces[0] = make_bilinear_patch({0., 0., 0.}, {0., 1., 0.}, {1., 0., 0.}, {1., 1., 0.}); + cube_faces[1] = make_bilinear_patch({0., 0., 1.}, {1., 0., 1.}, {0., 1., 1.}, {1., 1., 1.}); + cube_faces[2] = make_bilinear_patch({0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, {0., 1., 1.}); + cube_faces[3] = make_bilinear_patch({1., 0., 0.}, {1., 1., 0.}, {1., 0., 1.}, {1., 1., 1.}); + cube_faces[4] = make_bilinear_patch({0., 0., 0.}, {1., 0., 0.}, {0., 0., 1.}, {1., 0., 1.}); + cube_faces[5] = make_bilinear_patch({0., 1., 0.}, {0., 1., 1.}, {1., 1., 0.}, {1., 1., 1.}); + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + auto z_integrand = [](Point3D x) -> double { return x[2]; }; + + constexpr int npts = 6; + constexpr double abs_tol = 1e-10; + constexpr double lower_bound_z = 0.0; + + double patchwise_volume = 0.0; + double patchwise_z_moment = 0.0; + for(int i = 0; i < cube_faces.size(); ++i) + { + patchwise_volume += evaluate_volume_integral(cube_faces[i], const_integrand, lower_bound_z, npts); + patchwise_z_moment += evaluate_volume_integral(cube_faces[i], z_integrand, lower_bound_z, npts); + } + + EXPECT_NEAR(evaluate_surface_integral(cube_faces, const_integrand, npts), 6.0, abs_tol); + EXPECT_NEAR(evaluate_volume_integral(cube_faces, const_integrand, npts), 1.0, abs_tol); + EXPECT_NEAR(evaluate_volume_integral(cube_faces, z_integrand, npts), 0.5, abs_tol); + EXPECT_NEAR(patchwise_volume, 1.0, abs_tol); + EXPECT_NEAR(patchwise_z_moment, 0.5, abs_tol); + EXPECT_NEAR(evaluate_volume_integral(cube_faces[1], const_integrand, lower_bound_z, npts), + 1.0, + abs_tol); +} + +TEST(primal_integral, evaluate_trimmed_nurbs_patch_surface_integral) +{ + using Point2D = primal::Point; + using Point3D = primal::Point; + using NPatch = primal::NURBSPatch; + using TrimmingCurve = primal::NURBSCurve; + + Point3D control_points[] = {Point3D {0.0, 0.0, 0.0}, + Point3D {0.0, 1.0, 0.0}, + Point3D {1.0, 0.0, 0.0}, + Point3D {1.0, 1.0, 0.0}}; + + NPatch patch(control_points, 2, 2, 1, 1); + + axom::Array trim_pts {Point2D {0.25, 0.25}, + Point2D {0.75, 0.25}, + Point2D {0.75, 0.75}, + Point2D {0.25, 0.75}, + Point2D {0.25, 0.25}}; + patch.addTrimmingCurve(TrimmingCurve(trim_pts, 1)); + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + auto linear_integrand = [](Point3D x) -> double { return x[0] + x[1]; }; + + constexpr int npts = 8; + constexpr double abs_tol = 1e-10; + + EXPECT_NEAR(evaluate_surface_integral(patch, const_integrand, npts), 0.25, abs_tol); + EXPECT_NEAR(evaluate_surface_integral(patch, linear_integrand, npts), 0.25, abs_tol); +} + +TEST(primal_integral, evaluate_trimmed_nurbs_patch_surface_integral_clamps_bounds) +{ + using Point2D = primal::Point; + using Point3D = primal::Point; + using NPatch = primal::NURBSPatch; + using TrimmingCurve = primal::NURBSCurve; + + Point3D control_points[] = {Point3D {0.0, 0.0, 0.0}, + Point3D {0.0, 1.0, 0.0}, + Point3D {1.0, 0.0, 0.0}, + Point3D {1.0, 1.0, 0.0}}; + + NPatch patch(control_points, 2, 2, 1, 1); + + // Closed trimming loop with a control point below v-min (vmin == 0 for this patch), + // but whose geometric curve remains inside the patch. + axom::Array bottom_pts {Point2D {0.25, 0.25}, Point2D {0.50, -0.05}, Point2D {0.75, 0.25}}; + patch.addTrimmingCurve(TrimmingCurve(bottom_pts, 2)); + + Point2D right_pts[] = {Point2D {0.75, 0.25}, Point2D {0.75, 0.75}}; + patch.addTrimmingCurve(TrimmingCurve(right_pts, 2, 1)); + + Point2D top_pts[] = {Point2D {0.75, 0.75}, Point2D {0.25, 0.75}}; + patch.addTrimmingCurve(TrimmingCurve(top_pts, 2, 1)); + + Point2D left_pts[] = {Point2D {0.25, 0.75}, Point2D {0.25, 0.25}}; + patch.addTrimmingCurve(TrimmingCurve(left_pts, 2, 1)); + + bool saw_clamped_node = false; + const double vmin = patch.getMinKnot_v(); + auto integrand = [&saw_clamped_node, vmin](Point3D x) -> double { + if(x[1] == vmin) + { + saw_clamped_node = true; + } + return x[1]; + }; + + constexpr int npts = 8; + (void)evaluate_surface_integral(patch, integrand, npts); + + EXPECT_FALSE(saw_clamped_node); +} + +TEST(primal_integral, trimming_curve_array_lower_bound_includes_all_first_control_points) +{ + using Point2D = primal::Point; + using TrimmingCurve = primal::NURBSCurve; + + axom::Array curve0_pts {Point2D {0.0, 0.5}, Point2D {1.0, 0.6}}; + axom::Array curve1_pts {Point2D {0.0, -1.0}, Point2D {1.0, 1.0}}; + + axom::Array curves; + curves.push_back(TrimmingCurve(curve0_pts, 1)); + curves.push_back(TrimmingCurve(curve1_pts, 1)); + + EXPECT_DOUBLE_EQ(primal::detail::curve_array_lower_bound_y(curves), -1.0); +} + TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) { using Point2D = primal::Point; @@ -625,8 +760,12 @@ TEST(primal_integral, evaluate_integral_nurbs_gwn_cache) #ifdef AXOM_USE_MFEM TEST(primal_integral, check_axom_mfem_quadrature_values) { - const int N = 200; + // MFEM and Axom both generate Gauss-Legendre rules, but in builds that enable + // `-march=native` we can see ULP-level differences in the computed nodes/weights + // even though the rules are equivalent for integration purposes. + const double fp_tol = 8 * axom::numeric_limits::epsilon(); + constexpr int N = 200; for(int npts = 1; npts <= N; ++npts) { // Generate the Axom quadrature rule @@ -639,12 +778,8 @@ TEST(primal_integral, check_axom_mfem_quadrature_values) // Check that the nodes and weights are the same between the two rules for(int j = 0; j < npts; ++j) { - EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, axom::numeric_limits::epsilon()); - - // Relax tolerance slightly for intel-oneapi - EXPECT_NEAR(axom_rule.weight(j), - mfem_rule.IntPoint(j).weight, - 10 * axom::numeric_limits::epsilon()); + EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, fp_tol); + EXPECT_NEAR(axom_rule.weight(j), mfem_rule.IntPoint(j).weight, fp_tol); } } } diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index fd2215ee43..b9574c2cee 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -728,6 +728,14 @@ endif() if(OPENCASCADE_FOUND) + axom_add_executable( + NAME quest_step_moments_ex + SOURCES quest_step_moments.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON axom cli11 fmt opencascade + FOLDER axom/quest/examples + ) + axom_add_executable( NAME quest_step_file_ex SOURCES quest_step_file.cpp @@ -742,6 +750,15 @@ if(OPENCASCADE_FOUND) set(_num_ranks 2) endif() + axom_add_test( + NAME quest_step_moments_revolved_sphere + COMMAND quest_step_moments_ex -f ${quest_data_dir}/step/revolved_sphere.step + --integral both + --order 1 + --npts 16) + set_tests_properties(quest_step_moments_revolved_sphere PROPERTIES + PASS_REGULAR_EXPRESSION "SUMMARY surface_m0=3\\.14159265358979[0-9]*e\\+02 volume_m0=5\\.23598775598298[0-9]*e\\+02") + axom_add_test( NAME quest_step_file_sliced_cylinder COMMAND quest_step_file_ex -f ${quest_data_dir}/step/sliced_cylinder.step @@ -771,3 +788,4 @@ if(OPENCASCADE_FOUND) endif() endif() + diff --git a/src/axom/quest/examples/quest_step_moments.cpp b/src/axom/quest/examples/quest_step_moments.cpp new file mode 100644 index 0000000000..7d3ad7d0fd --- /dev/null +++ b/src/axom/quest/examples/quest_step_moments.cpp @@ -0,0 +1,1418 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file quest_step_moments.cpp + * \brief Example that reads a STEP file and computes geometric moments of the + * corresponding trimmed NURBS surface and/or enclosed volume. + * + * The example loads a BRep with Quest's STEP reader, then evaluates + * monomial moments \f$\int x^i y^j z^k \, dS\f$ and/or + * \f$\int x^i y^j z^k \, dV\f$ for all nonnegative exponent triples with + * total degree \f$i + j + k \le n\f$, where \a n is user supplied. + * It also derives centroids, inertia tensors, and volume-based principal-axis + * fit proxies from the first- and second-order moments. + * + * \note Volume moments are only geometrically meaningful when the STEP model + * represents a closed, consistently oriented boundary. + */ + +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/primal.hpp" +#include "axom/mint.hpp" +#include "axom/quest.hpp" + +#include "axom/CLI11.hpp" +#include "axom/fmt.hpp" + +#ifdef AXOM_USE_CONDUIT + #include "conduit.hpp" +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mint = axom::mint; +namespace primal = axom::primal; +namespace slic = axom::slic; + +namespace +{ + +using Point3D = primal::Point; +using Vector3D = primal::Vector; +using OBB3D = primal::OrientedBoundingBox; +using TriMesh = mint::UnstructuredMesh; +using PatchArray = axom::quest::STEPReader::PatchArray; +using MomentKey = std::tuple; + +enum class IntegralMode +{ + SURFACE, + VOLUME, + BOTH +}; + +const std::map s_validIntegralModes {{"surface", IntegralMode::SURFACE}, + {"volume", IntegralMode::VOLUME}, + {"both", IntegralMode::BOTH}}; + +constexpr double eps = 1e-14; + +double clamp_near_zero(double value, double tolerance) +{ + return std::abs(value) <= tolerance ? 0.0 : value; +} + +const char* integral_mode_name(IntegralMode mode) +{ + switch(mode) + { + case IntegralMode::SURFACE: + return "surface"; + case IntegralMode::VOLUME: + return "volume"; + case IntegralMode::BOTH: + return "both"; + } + + return "unknown"; +} + +bool should_compute_surface(IntegralMode mode) +{ + return mode == IntegralMode::SURFACE || mode == IntegralMode::BOTH; +} + +bool should_compute_volume(IntegralMode mode) +{ + return mode == IntegralMode::VOLUME || mode == IntegralMode::BOTH; +} + +struct MomentIndex +{ + int degree {}; + int px {}; + int py {}; + int pz {}; +}; + +struct MomentEntry +{ + MomentIndex index {}; + double value {}; +}; + +struct MomentSet +{ + std::vector requested_entries; + std::map values; +}; + +struct InertiaTensor +{ + double xx {}; + double yy {}; + double zz {}; + double xy {}; + double xz {}; + double yz {}; +}; + +struct MassProperties +{ + bool valid {false}; + bool orientation_reversed {false}; + double measure {}; + double signed_measure {}; + Point3D centroid {}; + InertiaTensor inertia_origin {}; + InertiaTensor inertia_centroid {}; +}; + +struct PrincipalFrame +{ + bool valid {false}; + double measure {}; + Point3D centroid {}; + double principal_inertia[3] {}; + double principal_second_moments[3] {}; + Vector3D axes[3]; +}; + +struct EllipsoidFit +{ + bool valid {false}; + Point3D centroid {}; + Vector3D axes[3]; + double radii[3] {}; + double volume {}; +}; + +struct ObbFit +{ + bool valid {false}; + OBB3D box {}; + double volume {}; +}; + +std::vector enumerate_moment_indices(int max_degree) +{ + std::vector moments; + for(int degree = 0; degree <= max_degree; ++degree) + { + for(int px = 0; px <= degree; ++px) + { + for(int py = 0; py <= degree - px; ++py) + { + const int pz = degree - px - py; + moments.push_back(MomentIndex {degree, px, py, pz}); + } + } + } + + return moments; +} + +double ipow(double base, int exponent) +{ + double result = 1.0; + for(int i = 0; i < exponent; ++i) + { + result *= base; + } + return result; +} + +double evaluate_monomial(const Point3D& x, const MomentIndex& idx) +{ + return ipow(x[0], idx.px) * ipow(x[1], idx.py) * ipow(x[2], idx.pz); +} + +template +MomentSet compute_moments(int requested_max_degree, int computed_max_degree, Integrator&& integrator) +{ + MomentSet moments; + const auto indices = enumerate_moment_indices(computed_max_degree); + moments.requested_entries.reserve(indices.size()); + + for(const auto& idx : indices) + { + const double value = integrator(idx); + moments.values[MomentKey {idx.px, idx.py, idx.pz}] = value; + + if(idx.degree <= requested_max_degree) + { + moments.requested_entries.push_back(MomentEntry {idx, value}); + } + } + + return moments; +} + +double get_moment(const MomentSet& moments, int px, int py, int pz) +{ + const auto it = moments.values.find(MomentKey {px, py, pz}); + SLIC_ASSERT(it != moments.values.end()); + return it->second; +} + +InertiaTensor shift_to_centroid(const InertiaTensor& inertia_origin, + double measure, + const Point3D& centroid) +{ + InertiaTensor shifted = inertia_origin; + + shifted.xx -= measure * (centroid[1] * centroid[1] + centroid[2] * centroid[2]); + shifted.yy -= measure * (centroid[0] * centroid[0] + centroid[2] * centroid[2]); + shifted.zz -= measure * (centroid[0] * centroid[0] + centroid[1] * centroid[1]); + shifted.xy += measure * centroid[0] * centroid[1]; + shifted.xz += measure * centroid[0] * centroid[2]; + shifted.yz += measure * centroid[1] * centroid[2]; + + return shifted; +} + +MassProperties compute_mass_properties(const MomentSet& moments, bool normalize_orientation = false) +{ + MassProperties props; + props.signed_measure = get_moment(moments, 0, 0, 0); + if(std::abs(props.signed_measure) <= eps) + { + return props; + } + + // Closed STEP shells can be oriented inward, which flips the sign of every + // volume moment without changing the underlying geometry. When requested, + // normalize that sign here so downstream centroid/inertia fits use a + // positive enclosed-volume convention. + const double orientation_sign = normalize_orientation && props.signed_measure < 0.0 ? -1.0 : 1.0; + + props.valid = true; + props.orientation_reversed = orientation_sign < 0.0; + props.measure = orientation_sign * props.signed_measure; + + const double mx = orientation_sign * get_moment(moments, 1, 0, 0); + const double my = orientation_sign * get_moment(moments, 0, 1, 0); + const double mz = orientation_sign * get_moment(moments, 0, 0, 1); + + props.centroid[0] = mx / props.measure; + props.centroid[1] = my / props.measure; + props.centroid[2] = mz / props.measure; + + const double xx = orientation_sign * get_moment(moments, 2, 0, 0); + const double yy = orientation_sign * get_moment(moments, 0, 2, 0); + const double zz = orientation_sign * get_moment(moments, 0, 0, 2); + const double xy = orientation_sign * get_moment(moments, 1, 1, 0); + const double xz = orientation_sign * get_moment(moments, 1, 0, 1); + const double yz = orientation_sign * get_moment(moments, 0, 1, 1); + + props.inertia_origin.xx = yy + zz; + props.inertia_origin.yy = xx + zz; + props.inertia_origin.zz = xx + yy; + props.inertia_origin.xy = -xy; + props.inertia_origin.xz = -xz; + props.inertia_origin.yz = -yz; + + props.inertia_centroid = shift_to_centroid(props.inertia_origin, props.measure, props.centroid); + + return props; +} + +PrincipalFrame compute_principal_frame(const MassProperties& props) +{ + constexpr int ndims = 3; + + PrincipalFrame frame; + if(!props.valid || std::abs(props.measure) <= eps) + { + return frame; + } + + axom::numerics::Matrix inertia(ndims, ndims); + inertia(0, 0) = props.inertia_centroid.xx; + inertia(0, 1) = props.inertia_centroid.xy; + inertia(0, 2) = props.inertia_centroid.xz; + inertia(1, 0) = props.inertia_centroid.xy; + inertia(1, 1) = props.inertia_centroid.yy; + inertia(1, 2) = props.inertia_centroid.yz; + inertia(2, 0) = props.inertia_centroid.xz; + inertia(2, 1) = props.inertia_centroid.yz; + inertia(2, 2) = props.inertia_centroid.zz; + + axom::numerics::Matrix eigenvectors(ndims, ndims); + double eigenvalues[ndims] {}; + const int rc = axom::numerics::jacobi_eigensolve(inertia, eigenvectors, eigenvalues); + if(rc != axom::numerics::JACOBI_EIGENSOLVE_SUCCESS) + { + return frame; + } + + frame.valid = true; + frame.measure = props.measure; + frame.centroid = props.centroid; + + double inertia_trace = 0.0; + for(int i = 0; i < ndims; ++i) + { + frame.principal_inertia[i] = eigenvalues[i]; + inertia_trace += eigenvalues[i]; + + frame.axes[i][0] = eigenvectors(0, i); + frame.axes[i][1] = eigenvectors(1, i); + frame.axes[i][2] = eigenvectors(2, i); + } + + // In the principal basis, the centroidal inertia tensor diagonal entries are + // pairwise sums of the centroidal second moments. Solve that 3x3 system to + // recover the second moments used by the ellipsoid and OBB proxies. + for(int i = 0; i < ndims; ++i) + { + frame.principal_second_moments[i] = 0.5 * (inertia_trace - 2.0 * frame.principal_inertia[i]); + } + + return frame; +} + +bool compute_shape_dimensions(const PrincipalFrame& frame, double factor, double (&dims)[3]) +{ + constexpr int ndims = 3; + constexpr double min_scale = 1.0; + + if(!frame.valid || std::abs(frame.measure) <= eps) + { + return false; + } + + const double scale = std::max({min_scale, + std::abs(frame.measure), + std::abs(frame.principal_second_moments[0]), + std::abs(frame.principal_second_moments[1]), + std::abs(frame.principal_second_moments[2])}); + const double tol = 1e-12 * scale; + + for(int i = 0; i < ndims; ++i) + { + // The factor selects the proxy family: 5 for a uniform ellipsoid semiaxis + // and 3 for a box half-extent that matches the same centroidal moments. + const double dim_sq = factor * frame.principal_second_moments[i] / frame.measure; + if(dim_sq < -tol) + { + return false; + } + + dims[i] = std::sqrt(std::max(dim_sq, 0.0)); + } + + return true; +} + +double compute_ellipsoid_volume(const double (&radii)[3]) +{ + return (4.0 / 3.0) * M_PI * radii[0] * radii[1] * radii[2]; +} + +double compute_obb_volume(const Vector3D& extents) +{ + return 8.0 * extents[0] * extents[1] * extents[2]; +} + +double compute_volume_ratio(double fit_volume, double reference_volume) +{ + if(std::abs(reference_volume) <= eps) + { + return 0.0; + } + + return fit_volume / reference_volume; +} + +double compute_effective_density(double mass, double geometric_volume) +{ + if(std::abs(geometric_volume) <= eps) + { + return 0.0; + } + + return mass / geometric_volume; +} + +EllipsoidFit compute_inertia_matched_ellipsoid(const PrincipalFrame& frame) +{ + EllipsoidFit fit; + double radii[3] {}; + if(!compute_shape_dimensions(frame, 5.0, radii)) + { + return fit; + } + + fit.valid = true; + fit.centroid = frame.centroid; + for(int i = 0; i < 3; ++i) + { + fit.axes[i] = frame.axes[i]; + fit.radii[i] = radii[i]; + } + fit.volume = compute_ellipsoid_volume(fit.radii); + + return fit; +} + +EllipsoidFit scale_ellipsoid_to_volume(const EllipsoidFit& fit, double target_volume) +{ + EllipsoidFit scaled; + if(!fit.valid || std::abs(fit.volume) <= eps || target_volume <= 0.0) + { + return scaled; + } + + const double scale = std::cbrt(target_volume / fit.volume); + scaled.valid = true; + scaled.centroid = fit.centroid; + for(int i = 0; i < 3; ++i) + { + scaled.axes[i] = fit.axes[i]; + scaled.radii[i] = scale * fit.radii[i]; + } + scaled.volume = compute_ellipsoid_volume(scaled.radii); + + return scaled; +} + +ObbFit compute_inertia_matched_obb(const PrincipalFrame& frame) +{ + ObbFit fit; + double extents_data[3] {}; + if(!compute_shape_dimensions(frame, 3.0, extents_data)) + { + return fit; + } + + Vector3D extents; + for(int i = 0; i < 3; ++i) + { + extents[i] = extents_data[i]; + } + + fit.box = OBB3D(frame.centroid, frame.axes, extents); + fit.valid = fit.box.isValid(); + fit.volume = compute_obb_volume(extents); + + return fit; +} + +ObbFit scale_obb_to_volume(const ObbFit& fit, double target_volume) +{ + ObbFit scaled; + if(!fit.valid || std::abs(fit.volume) <= eps || target_volume <= 0.0) + { + return scaled; + } + + const double scale = std::cbrt(target_volume / fit.volume); + const auto& centroid = fit.box.getCentroid(); + const auto* axes = fit.box.getAxes(); + const auto& extents = fit.box.getExtents(); + + Vector3D scaled_axes[3]; + Vector3D scaled_extents; + for(int i = 0; i < 3; ++i) + { + scaled_axes[i] = axes[i]; + scaled_extents[i] = scale * extents[i]; + } + + scaled.box = OBB3D(centroid, scaled_axes, scaled_extents); + scaled.valid = scaled.box.isValid(); + scaled.volume = compute_obb_volume(scaled_extents); + + return scaled; +} + +std::string make_ellipsoid_vtk_filename(const std::string& prefix, const std::string& variant) +{ + return axom::fmt::format("{}_{}.vtk", axom::utilities::string::removeSuffix(prefix, ".vtk"), variant); +} + +bool write_ellipsoid_vtk(const EllipsoidFit& fit, + const std::string& file_path, + int theta_resolution = 48, + int phi_resolution = 24) +{ + if(!fit.valid || file_path.empty()) + { + return false; + } + + SLIC_ASSERT(theta_resolution >= 3); + SLIC_ASSERT(phi_resolution >= 4); + + const int num_rings = phi_resolution - 2; + const int total_nodes = 2 + theta_resolution * num_rings; + const int total_cells = 2 * theta_resolution * (phi_resolution - 2); + + TriMesh mesh(3, mint::TRIANGLE); + mesh.reserve(total_nodes, total_cells); + + axom::IndexType next_node_id = 0; + + auto append_node = [&mesh, &fit, &next_node_id](double x, double y, double z) { + // Map a point on the unit sphere into the ellipsoid's principal-axis frame. + Point3D point = fit.centroid; + const double local_coords[3] {x, y, z}; + for(int axis = 0; axis < 3; ++axis) + { + for(int dim = 0; dim < 3; ++dim) + { + point[dim] += fit.radii[axis] * local_coords[axis] * fit.axes[axis][dim]; + } + } + + mesh.appendNode(point[0], point[1], point[2]); + return next_node_id++; + }; + + auto append_ring = [&append_node, theta_resolution](std::vector& ring, double phi) { + ring.clear(); + + const double sin_phi = std::sin(phi); + const double cos_phi = std::cos(phi); + if(std::abs(sin_phi) <= eps) + { + ring.push_back(append_node(0.0, 0.0, cos_phi)); + return; + } + + ring.reserve(theta_resolution); + for(int i = 0; i < theta_resolution; ++i) + { + const double theta = 2.0 * M_PI * static_cast(i) / theta_resolution; + ring.push_back(append_node(std::cos(theta) * sin_phi, std::sin(theta) * sin_phi, cos_phi)); + } + }; + + auto append_ring_triangles = [&mesh](const std::vector& lower_ring, + const std::vector& upper_ring) { + SLIC_ASSERT(!lower_ring.empty()); + SLIC_ASSERT(!upper_ring.empty()); + + axom::IndexType cell[3]; + if(lower_ring.size() == 1) + { + const axom::IndexType pole = lower_ring.front(); + const int ring_size = static_cast(upper_ring.size()); + for(int i = 0; i < ring_size; ++i) + { + const int next_i = (i + 1) % ring_size; + cell[0] = pole; + cell[1] = upper_ring[next_i]; + cell[2] = upper_ring[i]; + mesh.appendCell(cell); + } + return; + } + + if(upper_ring.size() == 1) + { + const axom::IndexType pole = upper_ring.front(); + const int ring_size = static_cast(lower_ring.size()); + for(int i = 0; i < ring_size; ++i) + { + const int next_i = (i + 1) % ring_size; + cell[0] = pole; + cell[1] = lower_ring[i]; + cell[2] = lower_ring[next_i]; + mesh.appendCell(cell); + } + return; + } + + SLIC_ASSERT(lower_ring.size() == upper_ring.size()); + const int ring_size = static_cast(lower_ring.size()); + for(int i = 0; i < ring_size; ++i) + { + const int next_i = (i + 1) % ring_size; + cell[0] = lower_ring[i]; + cell[1] = upper_ring[i]; + cell[2] = lower_ring[next_i]; + mesh.appendCell(cell); + + cell[0] = lower_ring[next_i]; + cell[1] = upper_ring[i]; + cell[2] = upper_ring[next_i]; + mesh.appendCell(cell); + } + }; + + std::vector lower_ring; + std::vector upper_ring; + lower_ring.reserve(theta_resolution); + upper_ring.reserve(theta_resolution); + + append_ring(lower_ring, 0.0); + for(int j = 1; j <= num_rings; ++j) + { + const double phi = M_PI * static_cast(j) / (phi_resolution - 1); + append_ring(upper_ring, phi); + append_ring_triangles(lower_ring, upper_ring); + lower_ring.swap(upper_ring); + } + + append_ring(upper_ring, M_PI); + append_ring_triangles(lower_ring, upper_ring); + + return mint::write_vtk(&mesh, file_path) == 0; +} + +std::string join_path(const std::string& prefix, const std::string& key) +{ + return prefix.empty() ? key : axom::fmt::format("{}/{}", prefix, key); +} + +class ResultsStore +{ +public: + explicit ResultsStore(double zero_tolerance = 0.0) : m_zero_tolerance(zero_tolerance) { } + + void set(const std::string& path, const std::array& value) { set_vec3(path, value); } + + void add(const std::string& prefix, const std::string& key, const std::string& value) + { + set_string(join_path(prefix, key), value); + } + + void add(const std::string& prefix, const std::string& key, const char* value) + { + set_string(join_path(prefix, key), value); + } + + void add(const std::string& prefix, const std::string& key, double value) + { + set_real(join_path(prefix, key), value); + } + + void add(const std::string& prefix, const std::string& key, const std::array& value) + { + set_vec3(join_path(prefix, key), value); + } + + void add_unclamped_real(const std::string& prefix, const std::string& key, double value) + { + set_real_unclamped(join_path(prefix, key), value); + } + + void add(const std::string& prefix, const std::string& key, int value) + { + set_integer(join_path(prefix, key), value); + } + + void add(const std::string& prefix, const std::string& key, bool value) + { + set_boolean(join_path(prefix, key), value); + } + + std::string to_yaml() const + { +#ifdef AXOM_USE_CONDUIT + return m_root.to_yaml(); +#else + std::string output; + append_yaml(output, m_root, 0); + return output; +#endif + } + +private: + double sanitize_real(double value) const { return clamp_near_zero(value, m_zero_tolerance); } + + const double m_zero_tolerance; + +#ifdef AXOM_USE_CONDUIT + conduit::Node m_root; + + void set_string(const std::string& path, const std::string& value) { m_root[path] = value; } + + void set_string(const std::string& path, const char* value) { m_root[path].set_string(value); } + + void set_real(const std::string& path, double value) { m_root[path] = sanitize_real(value); } + + void set_vec3(const std::string& path, const std::array& value) + { + const std::array sanitized {sanitize_real(value[0]), + sanitize_real(value[1]), + sanitize_real(value[2])}; + m_root[path].set(sanitized.data(), 3); + } + + void set_real_unclamped(const std::string& path, double value) { m_root[path] = value; } + + void set_integer(const std::string& path, int value) { m_root[path] = value; } + + void set_boolean(const std::string& path, bool value) + { + m_root[path].set_string(value ? "true" : "false"); + } +#else + using ScalarValue = std::variant>; + + struct TreeNode + { + bool has_scalar {false}; + ScalarValue scalar_value {std::string {}}; + std::unordered_map children; + std::vector child_order; + + TreeNode& fetch_child(const std::string& key) + { + auto it = children.find(key); + if(it == children.end()) + { + child_order.push_back(key); + it = children.emplace(key, TreeNode {}).first; + } + + return it->second; + } + }; + + TreeNode m_root; + + TreeNode& fetch_path(const std::string& path) + { + TreeNode* node = &m_root; + + for(const auto& key : axom::utilities::string::split(path, '/')) + { + if(!key.empty()) + { + node = &node->fetch_child(key); + } + } + + return *node; + } + + void set_string(const std::string& path, const std::string& value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = value; + } + + void set_string(const std::string& path, const char* value) + { + set_string(path, std::string(value)); + } + + void set_real(const std::string& path, double value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = sanitize_real(value); + } + + void set_vec3(const std::string& path, const std::array& value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = std::array {sanitize_real(value[0]), + sanitize_real(value[1]), + sanitize_real(value[2])}; + } + + void set_real_unclamped(const std::string& path, double value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = value; + } + + void set_integer(const std::string& path, int value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = value; + } + + void set_boolean(const std::string& path, bool value) + { + TreeNode& node = fetch_path(path); + node.has_scalar = true; + node.scalar_value = value; + } + + static std::string format_scalar(const ScalarValue& value) + { + if(const auto* string_value = std::get_if(&value)) + { + std::string escaped; + escaped.reserve(string_value->size()); + for(char ch : *string_value) + { + escaped += ch; + if(ch == '\'') + { + escaped += '\''; + } + } + + return axom::fmt::format("'{}'", escaped); + } + if(const auto* integer_value = std::get_if(&value)) + { + return axom::fmt::format("{}", *integer_value); + } + if(const auto* real_value = std::get_if(&value)) + { + return axom::fmt::format("{:.16e}", *real_value); + } + if(const auto* vec3_value = std::get_if>(&value)) + { + return axom::fmt::format("[{:.16e}, {:.16e}, {:.16e}]", + (*vec3_value)[0], + (*vec3_value)[1], + (*vec3_value)[2]); + } + + return std::get(value) ? "true" : "false"; + } + + static void append_yaml(std::string& output, const TreeNode& node, int indent) + { + const auto append_line = [&output](int line_indent, + const std::string& key, + const std::string* value = nullptr) { + axom::fmt::format_to(std::back_inserter(output), "{}{}:", std::string(line_indent * 2, ' '), key); + if(value != nullptr) + { + axom::fmt::format_to(std::back_inserter(output), " {}", *value); + } + output += '\n'; + }; + + for(const auto& key : node.child_order) + { + const auto it = node.children.find(key); + SLIC_ASSERT(it != node.children.end()); + const TreeNode& child = it->second; + + if(child.children.empty()) + { + if(child.has_scalar) + { + const std::string value = format_scalar(child.scalar_value); + append_line(indent, key, &value); + } + else + { + append_line(indent, key); + } + continue; + } + + append_line(indent, key); + if(child.has_scalar) + { + const std::string value = format_scalar(child.scalar_value); + append_line(indent + 1, "value", &value); + } + append_yaml(output, child, indent + 1); + } + } +#endif +}; + +class ResultsNode +{ +public: + ResultsNode(ResultsStore& results, std::string prefix) + : m_results(results) + , m_prefix(std::move(prefix)) + { } + + ResultsNode child(const std::string& key) const + { + return ResultsNode {m_results, join_path(m_prefix, key)}; + } + + void add(const std::string& key, const std::string& value) const + { + m_results.add(m_prefix, key, value); + } + + void add(const std::string& key, const char* value) const { m_results.add(m_prefix, key, value); } + + void add(const std::string& key, double value) const { m_results.add(m_prefix, key, value); } + + void add(const std::string& key, const std::array& value) const + { + m_results.add(m_prefix, key, value); + } + + void set(const std::array& value) const { m_results.set(m_prefix, value); } + + void add(const std::string& key, int value) const { m_results.add(m_prefix, key, value); } + + void add(const std::string& key, bool value) const { m_results.add(m_prefix, key, value); } + +private: + ResultsStore& m_results; + std::string m_prefix; +}; + +void write_xyz_triplet(const ResultsNode& node, double x, double y, double z) +{ + node.set(std::array {x, y, z}); +} + +void write(const ResultsNode& node, const Point3D& point) +{ + write_xyz_triplet(node, point[0], point[1], point[2]); +} + +void write(const ResultsNode& node, const Vector3D& vector) +{ + write_xyz_triplet(node, vector[0], vector[1], vector[2]); +} + +void write_axis_scalars(const ResultsNode& node, double a, double b, double c) +{ + node.add("axis_0", a); + node.add("axis_1", b); + node.add("axis_2", c); +} + +void write_axes(const ResultsNode& node, const Vector3D* axes) +{ + for(int i = 0; i < 3; ++i) + { + write(node.child(axom::fmt::format("axis_{}", i)), axes[i]); + } +} + +void write(const ResultsNode& node, const InertiaTensor& tensor) +{ + node.add("xx", tensor.xx); + node.add("yy", tensor.yy); + node.add("zz", tensor.zz); + node.add("xy", tensor.xy); + node.add("xz", tensor.xz); + node.add("yz", tensor.yz); +} + +void write(const ResultsNode& node, const MomentSet& moments) +{ + node.add("measure", get_moment(moments, 0, 0, 0)); + + const ResultsNode entries = node.child("entries"); + entries.add("count", static_cast(moments.requested_entries.size())); + for(const auto& entry : moments.requested_entries) + { + entries.add(axom::fmt::format("m_{}_{}_{}", entry.index.px, entry.index.py, entry.index.pz), + entry.value); + } +} + +void write(const ResultsNode& node, const MassProperties& props) +{ + if(!props.valid) + { + node.add("available", props.valid); + node.add("reason", "measure_is_zero"); + return; + } + + if(props.orientation_reversed) + { + node.add("signed_measure", props.signed_measure); + node.add("orientation_reversed", props.orientation_reversed); + } + + write(node.child("centroid"), props.centroid); + write(node.child("inertia_origin"), props.inertia_origin); + write(node.child("inertia_centroid"), props.inertia_centroid); +} + +void write(const ResultsNode& node, const PrincipalFrame& frame) +{ + if(!frame.valid) + { + node.add("available", frame.valid); + node.add("reason", "eigensolve_failed_or_measure_is_zero"); + return; + } + + write_axis_scalars(node.child("principal_inertia"), + frame.principal_inertia[0], + frame.principal_inertia[1], + frame.principal_inertia[2]); + write_axis_scalars(node.child("principal_second_moments"), + frame.principal_second_moments[0], + frame.principal_second_moments[1], + frame.principal_second_moments[2]); + write_axes(node.child("axes"), &frame.axes[0]); +} + +void write_assumptions(const ResultsNode& node, std::initializer_list assumptions) +{ + int idx = 0; + for(const char* assumption : assumptions) + { + node.add(axom::fmt::format("assumption_{}", idx++), assumption); + } +} + +void write_ellipsoid_fit(const ResultsNode& node, + const EllipsoidFit& fit, + double reference_volume, + const char* assumption, + const std::string& vtk_file) +{ + if(!fit.valid) + { + node.add("available", fit.valid); + node.add("reason", "invalid_second_moments"); + return; + } + + if(assumption != nullptr) + { + node.add("assumption", assumption); + } + if(!vtk_file.empty()) + { + node.add("vtk_file", vtk_file); + } + + write(node.child("center"), fit.centroid); + write_axis_scalars(node.child("semiaxes"), fit.radii[0], fit.radii[1], fit.radii[2]); + write_axes(node.child("axes"), &fit.axes[0]); + node.add("geometric_volume", fit.volume); + node.add("volume_ratio", compute_volume_ratio(fit.volume, reference_volume)); + node.add("effective_density", compute_effective_density(reference_volume, fit.volume)); +} + +void write_obb_fit(const ResultsNode& node, + const ObbFit& fit, + double reference_volume, + const char* assumption) +{ + if(!fit.valid) + { + node.add("available", fit.valid); + node.add("reason", "invalid_second_moments"); + return; + } + + if(assumption != nullptr) + { + node.add("assumption", assumption); + } + + const auto& centroid = fit.box.getCentroid(); + const auto& extents = fit.box.getExtents(); + const auto* axes = fit.box.getAxes(); + write(node.child("center"), centroid); + write_xyz_triplet(node.child("extents"), extents[0], extents[1], extents[2]); + write_axes(node.child("axes"), axes); + node.add("geometric_volume", fit.volume); + node.add("volume_ratio", compute_volume_ratio(fit.volume, reference_volume)); + node.add("effective_density", compute_effective_density(reference_volume, fit.volume)); +} + +void write_results(ResultsStore& results, + const std::string& input_file, + int patch_count, + const std::string& units, + int requested_order, + int computed_order, + int quadrature_order, + double zero_tolerance, + IntegralMode integral_mode, + bool has_surface, + const MomentSet& surface_moments, + const MassProperties& surface_props, + bool has_volume, + const MomentSet& volume_moments, + const MassProperties& volume_props, + const PrincipalFrame& volume_frame, + const EllipsoidFit& inertia_matched_ellipsoid, + const EllipsoidFit& same_volume_ellipsoid, + const ObbFit& inertia_matched_obb, + const ObbFit& same_volume_obb, + const std::string& inertia_matched_ellipsoid_vtk_file, + const std::string& same_volume_ellipsoid_vtk_file) +{ + const ResultsNode root {results, "results"}; + + const ResultsNode model = root.child("model"); + model.add("file", input_file); + model.add("patches", patch_count); + model.add("units", units); + + const ResultsNode config = root.child("config"); + config.add("requested_order", requested_order); + config.add("computed_order", computed_order); + config.add("quadrature_order", quadrature_order); + results.add_unclamped_real("results/config", "zero_tolerance", zero_tolerance); + config.add("integral", integral_mode_name(integral_mode)); + + const ResultsNode summary = root.child("summary"); + if(has_surface) + { + summary.add("surface_measure", surface_props.measure); + } + if(has_volume) + { + summary.add("volume_measure", volume_props.measure); + } + + if(has_surface) + { + const ResultsNode surface = root.child("surface"); + write(surface.child("raw_moments"), surface_moments); + write(surface.child("derived"), surface_props); + } + + if(has_volume) + { + const ResultsNode volume = root.child("volume"); + write_assumptions( + volume.child("assumptions"), + {"volume properties assume the STEP model is a closed, consistently oriented boundary", + "inertia_matched fits reproduce the 0th, 1st, and 2nd moments but may not match the " + "geometric volume at unit density", + "same_volume_scaled fits uniformly scale the inertia_matched shape to match the geometric " + "volume while preserving principal directions and aspect ratios", + "oriented_boxes reported here are moment-based proxies and are not guaranteed to bound the " + "model"}); + if(volume_props.orientation_reversed) + { + volume.add("orientation_note", + "raw volume moments indicated inward orientation; derived properties and fit " + "proxies were reoriented to use a positive enclosed volume"); + } + write(volume.child("raw_moments"), volume_moments); + write(volume.child("derived"), volume_props); + write(volume.child("principal_frame"), volume_frame); + + const ResultsNode fit_proxies = volume.child("fit_proxies"); + const ResultsNode ellipsoids = fit_proxies.child("ellipsoids"); + write_ellipsoid_fit(ellipsoids.child("inertia_matched"), + inertia_matched_ellipsoid, + volume_props.measure, + nullptr, + inertia_matched_ellipsoid_vtk_file); + write_ellipsoid_fit(ellipsoids.child("same_volume_scaled"), + same_volume_ellipsoid, + volume_props.measure, + "uniformly scaled from the inertia_matched ellipsoid", + same_volume_ellipsoid_vtk_file); + + const ResultsNode oriented_boxes = fit_proxies.child("oriented_boxes"); + write_obb_fit(oriented_boxes.child("inertia_matched"), + inertia_matched_obb, + volume_props.measure, + nullptr); + write_obb_fit(oriented_boxes.child("same_volume_scaled"), + same_volume_obb, + volume_props.measure, + "uniformly scaled from the inertia_matched oriented box"); + } +} +} // namespace + +int main(int argc, char** argv) +{ + std::string input_file; + int max_degree {1}; + int quadrature_order {0}; + bool verbose {false}; + bool validate_model {false}; + double zero_tolerance {0.0}; + std::string ellipsoid_vtk_prefix; + std::string annotationMode {"none"}; + IntegralMode integral_mode {IntegralMode::BOTH}; + + axom::CLI::App app { + "Load a STEP model and compute geometric moments, centroids, inertia tensors, and volume-based " + "fit proxies."}; + app.add_option("-f,--file", input_file) + ->description("Input STEP file") + ->required() + ->check(axom::CLI::ExistingFile); + app.add_option("-n,--order", max_degree) + ->description("Maximum total degree n for explicit monomial output x^i y^j z^k with i+j+k<=n") + ->capture_default_str() + ->check(axom::CLI::NonNegativeNumber); + app.add_option("--npts", quadrature_order) + ->description("Quadrature order for each numerical integration stage; defaults to max(8, n+4)") + ->capture_default_str() + ->check(axom::CLI::PositiveNumber); + app.add_option("--integral", integral_mode) + ->description("Which measures to compute: 'surface', 'volume', or 'both'") + ->capture_default_str() + ->transform(axom::CLI::CheckedTransformer(s_validIntegralModes)); + app.add_option("--zero-tolerance", zero_tolerance) + ->description("Clamp reported floating-point values with abs(value) <= tol to exactly zero") + ->capture_default_str() + ->check(axom::CLI::NonNegativeNumber); + app.add_flag("-v,--verbose", verbose, "Enable verbose output")->capture_default_str(); + app.add_flag("--validate", validate_model, "Run STEP model validation checks")->capture_default_str(); + app.add_option("--ellipsoid-vtk-prefix", ellipsoid_vtk_prefix) + ->description( + "Write the two volume ellipsoid fits to '_inertia_matched_ellipsoid.vtk' and " + "'_same_volume_scaled_ellipsoid.vtk'"); +#ifdef AXOM_USE_CALIPER + app.add_option("--caliper", annotationMode) + ->description( + "caliper annotation mode. Valid options include 'none' and 'report'. " + "See Axom's Caliper support for additional modes.") + ->capture_default_str() + ->check(axom::utilities::ValidCaliperMode); +#endif + app.get_formatter()->column_width(44); + + try + { + app.parse(argc, argv); + } + catch(const axom::CLI::ParseError& e) + { + return app.exit(e); + } + + axom::slic::SimpleLogger logger(axom::slic::message::Info); +#ifdef AXOM_USE_CALIPER + axom::utilities::raii::AnnotationsWrapper annotation_raii_wrapper(annotationMode); +#endif + AXOM_ANNOTATE_SCOPE("quest step moments example"); + + if(quadrature_order <= 0) + { + quadrature_order = axom::utilities::max(8, max_degree + 4); + } + + const int computed_max_degree = axom::utilities::max(max_degree, 2); + + SLIC_INFO(axom::fmt::format("Reading STEP file '{}'", input_file)); + + axom::quest::STEPReader reader; + reader.setFileName(input_file); + reader.setVerbosity(verbose); + + { + AXOM_ANNOTATE_SCOPE("read step"); + const int read_status = reader.read(validate_model); + if(read_status != 0) + { + SLIC_ERROR("Failed to read STEP file."); + return 1; + } + } + + const PatchArray& patches = reader.getPatchArray(); + if(patches.empty()) + { + SLIC_ERROR("STEP file did not contain any patches."); + return 1; + } + + if(verbose) + { + SLIC_INFO(axom::fmt::format("STEP file units: '{}'", reader.getFileUnits())); + SLIC_INFO(reader.getBRepStats()); + } + + MomentSet surface_moments; + MomentSet volume_moments; + MassProperties surface_props; + MassProperties volume_props; + PrincipalFrame volume_frame; + EllipsoidFit inertia_matched_ellipsoid; + EllipsoidFit same_volume_ellipsoid; + ObbFit inertia_matched_obb; + ObbFit same_volume_obb; + std::string inertia_matched_ellipsoid_vtk_file; + std::string same_volume_ellipsoid_vtk_file; + + if(should_compute_surface(integral_mode)) + { + AXOM_ANNOTATE_SCOPE("compute surface properties"); + surface_moments = compute_moments(max_degree, computed_max_degree, [&](const MomentIndex& idx) { + auto integrand = [idx](const Point3D& x) -> double { return evaluate_monomial(x, idx); }; + return primal::evaluate_surface_integral(patches, integrand, quadrature_order); + }); + surface_props = compute_mass_properties(surface_moments); + } + + if(should_compute_volume(integral_mode)) + { + AXOM_ANNOTATE_SCOPE("compute volume properties"); + volume_moments = compute_moments(max_degree, computed_max_degree, [&](const MomentIndex& idx) { + auto integrand = [idx](const Point3D& x) -> double { return evaluate_monomial(x, idx); }; + return primal::evaluate_volume_integral(patches, integrand, quadrature_order); + }); + // Normalize inward-oriented closed shells before deriving principal axes and + // fit proxies so the reported volume-based quantities use the same sign + // convention as outward-oriented solids. + volume_props = compute_mass_properties(volume_moments, true); + + if(volume_props.valid) + { + AXOM_ANNOTATE_SCOPE("compute volume fit proxies"); + volume_frame = compute_principal_frame(volume_props); + inertia_matched_ellipsoid = compute_inertia_matched_ellipsoid(volume_frame); + same_volume_ellipsoid = + scale_ellipsoid_to_volume(inertia_matched_ellipsoid, volume_props.measure); + inertia_matched_obb = compute_inertia_matched_obb(volume_frame); + same_volume_obb = scale_obb_to_volume(inertia_matched_obb, volume_props.measure); + } + } + + if(!ellipsoid_vtk_prefix.empty()) + { + AXOM_ANNOTATE_SCOPE("write ellipsoid vtk"); + if(!should_compute_volume(integral_mode)) + { + SLIC_WARNING("Ellipsoid VTK export requested, but volume integrals were not computed."); + } + else if(!inertia_matched_ellipsoid.valid || !same_volume_ellipsoid.valid) + { + SLIC_WARNING( + "Ellipsoid VTK export requested, but the volume ellipsoid fits are unavailable."); + } + else + { + inertia_matched_ellipsoid_vtk_file = + make_ellipsoid_vtk_filename(ellipsoid_vtk_prefix, "inertia_matched_ellipsoid"); + same_volume_ellipsoid_vtk_file = + make_ellipsoid_vtk_filename(ellipsoid_vtk_prefix, "same_volume_scaled_ellipsoid"); + + const bool wrote_inertia = + write_ellipsoid_vtk(inertia_matched_ellipsoid, inertia_matched_ellipsoid_vtk_file); + const bool wrote_same_volume = + write_ellipsoid_vtk(same_volume_ellipsoid, same_volume_ellipsoid_vtk_file); + + if(!wrote_inertia || !wrote_same_volume) + { + SLIC_WARNING("Failed to write one or more ellipsoid VTK files."); + } + } + } + + { + AXOM_ANNOTATE_SCOPE("log results"); + std::string summary = "SUMMARY"; + if(should_compute_surface(integral_mode)) + { + summary += axom::fmt::format(" surface_m0={:.16e}", + clamp_near_zero(surface_props.measure, zero_tolerance)); + } + if(should_compute_volume(integral_mode)) + { + summary += axom::fmt::format(" volume_m0={:.16e}", + clamp_near_zero(volume_props.measure, zero_tolerance)); + } + SLIC_INFO(summary); + + ResultsStore results(zero_tolerance); + write_results(results, + input_file, + static_cast(patches.size()), + reader.getFileUnits(), + max_degree, + computed_max_degree, + quadrature_order, + zero_tolerance, + integral_mode, + should_compute_surface(integral_mode), + surface_moments, + surface_props, + should_compute_volume(integral_mode), + volume_moments, + volume_props, + volume_frame, + inertia_matched_ellipsoid, + same_volume_ellipsoid, + inertia_matched_obb, + same_volume_obb, + inertia_matched_ellipsoid_vtk_file, + same_volume_ellipsoid_vtk_file); + SLIC_INFO(results.to_yaml()); + } + + return 0; +} diff --git a/src/axom/quest/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index 95ad186df4..d30e7f40ab 100644 --- a/src/axom/quest/tests/CMakeLists.txt +++ b/src/axom/quest/tests/CMakeLists.txt @@ -30,6 +30,7 @@ blt_list_append(TO quest_tests ELEMENTS quest_meshtester.cpp) if(OPENCASCADE_FOUND AND AXOM_DATA_DIR) + list(APPEND quest_tests quest_step_quadrature.cpp) list(APPEND quest_tests quest_step_reader.cpp) endif() diff --git a/src/axom/quest/tests/quest_step_quadrature.cpp b/src/axom/quest/tests/quest_step_quadrature.cpp new file mode 100644 index 0000000000..8a32ad1ed9 --- /dev/null +++ b/src/axom/quest/tests/quest_step_quadrature.cpp @@ -0,0 +1,257 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other +// Axom Project Contributors. See top-level LICENSE and COPYRIGHT +// files for dates and other details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "axom/quest/io/STEPReader.hpp" +#include "axom/primal/operators/evaluate_integral.hpp" +#include "axom/core/numerics/transforms.hpp" + +#include "gtest/gtest.h" + +#include +#include + +namespace primal = axom::primal; +namespace quest = axom::quest; + +namespace +{ +using Point3D = primal::Point; +using Matrix = axom::numerics::Matrix; +using PatchArray = quest::STEPReader::PatchArray; + +PatchArray read_step_patches(const std::string& file) +{ + quest::STEPReader reader; + reader.setFileName(file); + EXPECT_EQ(reader.read(false), 0); + return reader.getPatchArray(); +} + +Matrix embed_linear_transform(const Matrix& linear) +{ + EXPECT_EQ(linear.getNumRows(), 3); + EXPECT_EQ(linear.getNumColumns(), 3); + + Matrix affine = Matrix::identity(4); + for(int i = 0; i < 3; ++i) + { + for(int j = 0; j < 3; ++j) + { + affine(i, j) = linear(i, j); + } + } + + return affine; +} + +void transform_patches(PatchArray& patches, const Matrix& transform) +{ + for(int p = 0; p < patches.size(); ++p) + { + auto& control_points = patches[p].getControlPoints(); + for(axom::IndexType ui = 0; ui < control_points.shape()[0]; ++ui) + { + for(axom::IndexType vi = 0; vi < control_points.shape()[1]; ++vi) + { + control_points(ui, vi) = primal::transform_point(control_points(ui, vi), transform); + } + } + } +} + +template +void expect_matching_surface_values(const PatchArray& patches_a, + const PatchArray& patches_b, + const Lambda& integrand, + int npts, + double expected, + double abs_tol) +{ + const double observed_a = primal::evaluate_surface_integral(patches_a, integrand, npts); + const double observed_b = primal::evaluate_surface_integral(patches_b, integrand, npts); + + EXPECT_NEAR(observed_a, expected, abs_tol); + EXPECT_NEAR(observed_b, expected, abs_tol); + EXPECT_NEAR(observed_a, observed_b, abs_tol); +} + +template +void expect_matching_volume_values(const PatchArray& patches_a, + const PatchArray& patches_b, + const Lambda& integrand, + int npts, + double expected, + double abs_tol) +{ + const double observed_a = primal::evaluate_volume_integral(patches_a, integrand, npts); + const double observed_b = primal::evaluate_volume_integral(patches_b, integrand, npts); + + EXPECT_NEAR(observed_a, expected, abs_tol); + EXPECT_NEAR(observed_b, expected, abs_tol); + EXPECT_NEAR(observed_a, observed_b, abs_tol); +} + +} // namespace + +TEST(quest_step_quadrature, tet_surface_and_volume) +{ + auto patches = read_step_patches(std::string(AXOM_DATA_DIR) + "/quest/step/tet.step"); + EXPECT_EQ(patches.size(), 4); + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + + constexpr int npts = 6; + constexpr double abs_tol = 1e-10; + + EXPECT_NEAR(primal::evaluate_surface_integral(patches, const_integrand, npts), + 8.0 * std::sqrt(3.0), + abs_tol); + EXPECT_NEAR(primal::evaluate_volume_integral(patches, const_integrand, npts), 8.0 / 3.0, abs_tol); +} + +TEST(quest_step_quadrature, sphere_models_match_unit_sphere_moments) +{ + const std::string step_dir = std::string(AXOM_DATA_DIR) + "/quest/step/"; + + auto revolved = read_step_patches(step_dir + "revolved_sphere.step"); + auto biquartic = read_step_patches(step_dir + "biquartic_sphere.step"); + + EXPECT_EQ(revolved.size(), 1); + EXPECT_FALSE(biquartic.empty()); + EXPECT_EQ(biquartic.size(), 6); + + // revolved_sphere.step is a radius-5 sphere, while the biquartic sphere patch is unit-radius. + transform_patches(revolved, axom::numerics::transforms::scale(1.0 / 5.0, 4)); + + constexpr int npts = 16; + constexpr double abs_tol = 5e-5; + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + auto x_integrand = [](Point3D x) -> double { return x[0]; }; + auto y_integrand = [](Point3D x) -> double { return x[1]; }; + auto z_integrand = [](Point3D x) -> double { return x[2]; }; + auto xx_integrand = [](Point3D x) -> double { return x[0] * x[0]; }; + auto yy_integrand = [](Point3D x) -> double { return x[1] * x[1]; }; + auto zz_integrand = [](Point3D x) -> double { return x[2] * x[2]; }; + + const double sphere_area = 4.0 * M_PI; + const double sphere_volume = 4.0 * M_PI / 3.0; + const double surface_second_moment = 4.0 * M_PI / 3.0; + const double volume_second_moment = 4.0 * M_PI / 15.0; + + expect_matching_surface_values(revolved, biquartic, const_integrand, npts, sphere_area, abs_tol); + expect_matching_surface_values(revolved, biquartic, x_integrand, npts, 0.0, abs_tol); + expect_matching_surface_values(revolved, biquartic, y_integrand, npts, 0.0, abs_tol); + expect_matching_surface_values(revolved, biquartic, z_integrand, npts, 0.0, abs_tol); + expect_matching_surface_values(revolved, biquartic, xx_integrand, npts, surface_second_moment, abs_tol); + expect_matching_surface_values(revolved, biquartic, yy_integrand, npts, surface_second_moment, abs_tol); + expect_matching_surface_values(revolved, biquartic, zz_integrand, npts, surface_second_moment, abs_tol); + + expect_matching_volume_values(revolved, biquartic, const_integrand, npts, sphere_volume, abs_tol); + expect_matching_volume_values(revolved, biquartic, x_integrand, npts, 0.0, abs_tol); + expect_matching_volume_values(revolved, biquartic, y_integrand, npts, 0.0, abs_tol); + expect_matching_volume_values(revolved, biquartic, z_integrand, npts, 0.0, abs_tol); + expect_matching_volume_values(revolved, biquartic, xx_integrand, npts, volume_second_moment, abs_tol); + expect_matching_volume_values(revolved, biquartic, yy_integrand, npts, volume_second_moment, abs_tol); + expect_matching_volume_values(revolved, biquartic, zz_integrand, npts, volume_second_moment, abs_tol); +} + +TEST(quest_step_quadrature, transformed_sphere_models_match_expected_moments) +{ + const std::string step_dir = std::string(AXOM_DATA_DIR) + "/quest/step/"; + + auto revolved = read_step_patches(step_dir + "revolved_sphere.step"); + auto biquartic = read_step_patches(step_dir + "biquartic_sphere.step"); + + transform_patches(revolved, axom::numerics::transforms::scale(1.0 / 5.0, 4)); + + constexpr double sphere_scale = 1.7; + const Point3D center {1.25, -0.75, 2.0}; + + const Matrix scale = axom::numerics::transforms::scale(sphere_scale, 4); + const Matrix rotation = + embed_linear_transform(axom::numerics::transforms::axisRotation(M_PI / 4.0, 1.0, 2.0, 3.0)); + const Matrix translation = axom::numerics::transforms::translate(center[0], center[1], center[2]); + + transform_patches(revolved, scale); + transform_patches(revolved, rotation); + transform_patches(revolved, translation); + + transform_patches(biquartic, scale); + transform_patches(biquartic, rotation); + transform_patches(biquartic, translation); + + constexpr int npts = 16; + constexpr double abs_tol = 1e-4; + + auto const_integrand = [](Point3D /*x*/) -> double { return 1.0; }; + auto x_integrand = [](Point3D x) -> double { return x[0]; }; + auto y_integrand = [](Point3D x) -> double { return x[1]; }; + auto z_integrand = [](Point3D x) -> double { return x[2]; }; + auto centered_xx_integrand = [center](Point3D x) -> double { + const double dx = x[0] - center[0]; + return dx * dx; + }; + auto centered_yy_integrand = [center](Point3D x) -> double { + const double dy = x[1] - center[1]; + return dy * dy; + }; + auto centered_zz_integrand = [center](Point3D x) -> double { + const double dz = x[2] - center[2]; + return dz * dz; + }; + + const double sphere_area = 4.0 * M_PI * sphere_scale * sphere_scale; + const double sphere_volume = (4.0 * M_PI / 3.0) * sphere_scale * sphere_scale * sphere_scale; + const double surface_second_moment = (4.0 * M_PI / 3.0) * std::pow(sphere_scale, 4); + const double volume_second_moment = (4.0 * M_PI / 15.0) * std::pow(sphere_scale, 5); + + expect_matching_surface_values(revolved, biquartic, const_integrand, npts, sphere_area, abs_tol); + expect_matching_surface_values(revolved, biquartic, x_integrand, npts, sphere_area * center[0], abs_tol); + expect_matching_surface_values(revolved, biquartic, y_integrand, npts, sphere_area * center[1], abs_tol); + expect_matching_surface_values(revolved, biquartic, z_integrand, npts, sphere_area * center[2], abs_tol); + expect_matching_surface_values(revolved, + biquartic, + centered_xx_integrand, + npts, + surface_second_moment, + abs_tol); + expect_matching_surface_values(revolved, + biquartic, + centered_yy_integrand, + npts, + surface_second_moment, + abs_tol); + expect_matching_surface_values(revolved, + biquartic, + centered_zz_integrand, + npts, + surface_second_moment, + abs_tol); + + expect_matching_volume_values(revolved, biquartic, const_integrand, npts, sphere_volume, abs_tol); + expect_matching_volume_values(revolved, biquartic, x_integrand, npts, sphere_volume * center[0], abs_tol); + expect_matching_volume_values(revolved, biquartic, y_integrand, npts, sphere_volume * center[1], abs_tol); + expect_matching_volume_values(revolved, biquartic, z_integrand, npts, sphere_volume * center[2], abs_tol); + expect_matching_volume_values(revolved, + biquartic, + centered_xx_integrand, + npts, + volume_second_moment, + abs_tol); + expect_matching_volume_values(revolved, + biquartic, + centered_yy_integrand, + npts, + volume_second_moment, + abs_tol); + expect_matching_volume_values(revolved, + biquartic, + centered_zz_integrand, + npts, + volume_second_moment, + abs_tol); +}