diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index dfb3bde57a..ba31e92f11 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -81,6 +81,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Adds `quest::MFEMReader` for reading 1D MFEM contours in 2D space. - Adds an option to `quest::SamplingShaper` to allow in/out tests based on winding numbers for MFEM contours. - The `shaping_driver` example program can select `--sampling inout` to do the default In/Out sampling and `--sampling windingnumber` to select winding number in/out tests for MFEM data. +- Adds Axom-native Gauss-Legendre quadrature rules that can be used without an MFEM dependency ### Changed - Updates blt submodule to [BLT version 0.7.1][https://github.com/LLNL/blt/releases/tag/v0.7.1] diff --git a/src/axom/core/CMakeLists.txt b/src/axom/core/CMakeLists.txt index 77ec820ef2..39ed6116c7 100644 --- a/src/axom/core/CMakeLists.txt +++ b/src/axom/core/CMakeLists.txt @@ -42,6 +42,7 @@ set(core_headers numerics/eigen_solve.hpp numerics/eigen_sort.hpp numerics/floating_point_limits.hpp + numerics/quadrature.hpp numerics/jacobi_eigensolve.hpp numerics/linear_solve.hpp numerics/matvecops.hpp @@ -106,6 +107,7 @@ set(core_sources ${PROJECT_BINARY_DIR}/axom/core/utilities/About.cpp numerics/polynomial_solvers.cpp + numerics/quadrature.cpp memory_management.cpp Path.cpp diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp new file mode 100644 index 0000000000..9a54819f84 --- /dev/null +++ b/src/axom/core/numerics/quadrature.cpp @@ -0,0 +1,168 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "axom/core/utilities/Utilities.hpp" +#include "axom/core/Array.hpp" +#include "axom/core/FlatMap.hpp" +#include "axom/core/NumericLimits.hpp" +#include "axom/core/numerics/quadrature.hpp" + +// For math constants and includes +#include "axom/config.hpp" + +#include + +namespace axom +{ +namespace numerics +{ + +/*! + * \brief Computes a 1D quadrature rule of Gauss-Legendre points + * + * \param [in] npts The number of points in the rule + * \param [out] nodes The array of 1D nodes + * \param [out] weights The array of weights + * + * A Gauss-Legendre rule with \a npts points can exactly integrate + * polynomials of order 2 * npts - 1 + * + * Algorithm adapted from the MFEM implementation in `mfem/fem/intrules.cpp` + * + * \note This method constructs the points by scratch each time, without caching + * \sa get_gauss_legendre(int) + */ +void compute_gauss_legendre_data(int npts, + axom::Array& nodes, + axom::Array& weights, + int allocatorID) +{ + assert("Quadrature rules must have >= 1 point" && (npts >= 1)); + + nodes = axom::Array(npts, npts, allocatorID); + weights = axom::Array(npts, npts, allocatorID); + + if(npts == 1) + { + nodes[0] = 0.5; + weights[0] = 1.0; + return; + } + if(npts == 2) + { + nodes[0] = 0.21132486540518711775; + nodes[1] = 0.78867513459481288225; + + weights[0] = weights[1] = 0.5; + return; + } + if(npts == 3) + { + nodes[0] = 0.11270166537925831148207345; + nodes[1] = 0.5; + nodes[2] = 0.88729833462074168851792655; + + weights[0] = 0.2777777777777777777777778; + weights[1] = 0.4444444444444444444444444; + weights[2] = 0.2777777777777777777777778; + return; + } + + const int n = npts; + const int m = (npts + 1) / 2; + + // Nodes are mirrored across x = 0.5, so only need to evaluate half + for(int i = 1; i <= m; ++i) + { + // Each node is the root of a Legendre polynomial, + // which are approximately uniformly distributed in arccos(xi). + // This makes cos a good initial guess for subsequent Newton iterations + double z = std::cos(M_PI * (i - 0.25) / (n + 0.5)); + double Pp_n, P_n, dz, xi = 0.0; + + bool done = false; + while(true) + { + // Recursively evaluate P_n(z) through the recurrence relation + // n * P_n(z) = (2n - 1) * P_{n-1}(z) - (n - 1) * P_{n - 2}(z) + double P_nm1 = 1.0; // P_0(z) = 1 + P_n = z; // P_1(z) = z + for(int j = 2; j <= n; ++j) + { + double P_nm2 = P_nm1; + P_nm1 = P_n; + P_n = ((2 * j - 1) * z * P_nm1 - (j - 1) * P_nm2) / j; + } + + // Evaluate P'_n(z) using another recurrence relation + // (z^2 - 1) * P'_n(z) = n * z * P_n(z) - n * P_{n-1}(z) + Pp_n = n * (z * P_n - P_nm1) / (z * z - 1); + + if(done) + { + break; + } + + // Compute the Newton method step size + dz = P_n / Pp_n; + + if(std::fabs(dz) < axom::numeric_limits::epsilon()) + { + done = true; + + // Scale back to [0, 1] + xi = ((1 - z) + dz) / 2; + } + + // Take the Newton step, repeat the process + z -= dz; + } + + nodes[i - 1] = xi; + nodes[n - i] = 1.0 - xi; + + // For z \in [-1, 1], w_i = 2 / (1 - z^2) / P'_n(z)^2. + // For nodes[i] = xi = (1 - z)/2 \in [0, 1], weights[i] = 0.5 * w_i + weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); + } +} + +/*! + * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points + * + * \param [in] npts The number of points in the rule + * + * A Gauss-Legendre rule with \a npts points can exactly integrate + * polynomials of order 2 * npts - 1 + * + * \note If this method has already been called for a given order, it will reuse the same quadrature points + * without needing to recompute them + * + * \warning The use of a static variable to store cached nodes makes this method not threadsafe. + * + * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights + */ +QuadratureRule get_gauss_legendre(int npts, int allocatorID) +{ + assert("Quadrature rules must have >= 1 point" && (npts >= 1)); + + // Define a static map that stores the GL quadrature rule for a given order + static std::map, std::pair, axom::Array>> rule_library; + + const std::pair key = std::make_pair(npts, allocatorID); + + auto value_it = rule_library.find(key); + if(value_it == rule_library.end()) + { + auto& vals = rule_library[key]; + compute_gauss_legendre_data(npts, vals.first, vals.second, allocatorID); + value_it = rule_library.find(key); + } + + return QuadratureRule {value_it->second.first.view(), value_it->second.second.view()}; +} + +} /* end namespace numerics */ +} /* end namespace axom */ diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp new file mode 100644 index 0000000000..2fe3230bea --- /dev/null +++ b/src/axom/core/numerics/quadrature.hpp @@ -0,0 +1,94 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#ifndef AXOM_NUMERICS_QUADRATURE_HPP_ +#define AXOM_NUMERICS_QUADRATURE_HPP_ + +#include "axom/core/Array.hpp" +#include "axom/core/memory_management.hpp" + +/*! + * \file quadrature.hpp + * The functions declared in this header file find the nodes and weights of + * arbitrary order quadrature rules + */ + +namespace axom +{ +namespace numerics +{ + +/*! + * \class QuadratureRule + * + * \brief Stores a fixed array of 1D quadrature points and weights + */ +class QuadratureRule +{ + // Define friend functions so rules can only be created via get_rule() methods + friend QuadratureRule get_gauss_legendre(int, int); + +public: + //! \brief Accessor for quadrature nodes + AXOM_HOST_DEVICE + double node(size_t idx) const { return m_nodes[idx]; }; + + //! \brief Accessor for quadrature weights + AXOM_HOST_DEVICE + double weight(size_t idx) const { return m_weights[idx]; }; + + //! \brief Accessor for the size of the quadrature rule + AXOM_HOST_DEVICE + int getNumPoints() const { return m_nodes.size(); } + +private: + //! \brief Use a private constructor to avoid creation of an invalid rule + QuadratureRule(axom::ArrayView nodes, axom::ArrayView weights) + : m_nodes(nodes) + , m_weights(weights) { }; + + axom::ArrayView m_nodes; + axom::ArrayView m_weights; +}; + +/*! + * \brief Computes a 1D quadrature rule of Gauss-Legendre points + * + * \param [in] npts The number of points in the rule + * \param [out] nodes The array of 1D nodes + * \param [out] weights The array of weights + * + * A Gauss-Legendre rule with \a npts points can exactly integrate + * polynomials of order 2 * npts - 1 + * + * Algorithm adapted from the MFEM implementation in `mfem/fem/intrules.cpp` + * + * \note This method constructs the points by scratch each time, without caching + * \sa get_gauss_legendre(int) + */ +void compute_gauss_legendre_data(int npts, + axom::Array& nodes, + axom::Array& weights, + int allocatorID = axom::getDefaultAllocatorID()); + +/*! + * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points + * + * \param [in] npts The number of points in the rule + * + * A Gauss-Legendre rule with \a npts points can exactly integrate + * polynomials of order 2 * npts - 1 + * + * \note If this method has already been called for a given order, it will reuse the same quadrature points + * without needing to recompute them + * + * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights + */ +QuadratureRule get_gauss_legendre(int npts, int allocatorID = axom::getDefaultAllocatorID()); + +} /* end namespace numerics */ +} /* end namespace axom */ + +#endif // AXOM_NUMERICS_QUADRATURE_HPP_ diff --git a/src/axom/core/tests/CMakeLists.txt b/src/axom/core/tests/CMakeLists.txt index 854ca85edc..33658f5370 100644 --- a/src/axom/core/tests/CMakeLists.txt +++ b/src/axom/core/tests/CMakeLists.txt @@ -39,6 +39,7 @@ set(core_serial_tests numerics_eigen_sort.hpp numerics_floating_point_limits.hpp numerics_jacobi_eigensolve.hpp + numerics_quadrature.hpp numerics_linear_solve.hpp numerics_lu.hpp numerics_matrix.hpp diff --git a/src/axom/core/tests/core_serial_main.cpp b/src/axom/core/tests/core_serial_main.cpp index 70442354d6..d60e03495d 100644 --- a/src/axom/core/tests/core_serial_main.cpp +++ b/src/axom/core/tests/core_serial_main.cpp @@ -39,6 +39,7 @@ #include "numerics_matrix.hpp" #include "numerics_matvecops.hpp" #include "numerics_polynomial_solvers.hpp" +#include "numerics_quadrature.hpp" #include "utils_endianness.hpp" #include "utils_fileUtilities.hpp" diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp new file mode 100644 index 0000000000..f0bc2e6ec1 --- /dev/null +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -0,0 +1,126 @@ +// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "gtest/gtest.h" + +#include "axom/config.hpp" +#include "axom/core/numerics/quadrature.hpp" +#include "axom/core/utilities/Utilities.hpp" + +TEST(numerics_quadrature, gauss_legendre_math_check) +{ + const int N = 200; + + double coeffs[2 * N - 1]; + + // Test that the rules provide exact integration for polynomials of degree 2n - 1 + for(int npts = 1; npts <= N; ++npts) + { + // Evaluate using the quadrature rule + auto rule = axom::numerics::get_gauss_legendre(npts); + int degree = 2 * npts - 1; + + // Define a collection of random coefficients for a polynomial + for(int j = 0; j < degree; ++j) + { + // Seed the random coefficients for consistency in the test + coeffs[j] = axom::utilities::random_real(-1.0, 1.0, npts + j); + } + + // Evaluate the area under the curve from 0 to 1 analytically + double analytic_result = 0.0; + for(int j = 0; j < degree; ++j) + { + analytic_result += coeffs[j] / (j + 1); + } + + // Evaluate the polynomial using Horner's rule + auto eval_polynomial = [°ree, &coeffs](double x) { + double result = coeffs[degree - 1]; + for(int i = degree - 2; i >= 0; --i) + { + result = result * x + coeffs[i]; + } + return result; + }; + + double quadrature_result = 0.0; + for(int j = 0; j < npts; ++j) + { + quadrature_result += rule.weight(j) * eval_polynomial(rule.node(j)); + } + + EXPECT_NEAR(quadrature_result, analytic_result, 1e-10); + + // Check that nodes aren't duplicated + for(int j = 1; j < npts; ++j) + { + EXPECT_GT(rule.node(j), rule.node(j - 1)); + } + + // Check that the sum of the weights is 1, and that all are positive + double weight_sum = 0.0; + for(int j = 0; j < npts; ++j) + { + EXPECT_GT(rule.weight(j), 0.0); + weight_sum += rule.weight(j); + } + + EXPECT_NEAR(weight_sum, 1.0, 1e-10); + + // Check that each node is the root of the next Legendre polynomial + for(int j = 0; j < npts; ++j) + { + // Rescale the node to [-1, 1] + double z = 2 * rule.node(j) - 1; + + double P_n = z, P_nm1 = 1.0; + for(int k = 2; k <= npts; ++k) + { + double P_nm2 = P_nm1; + P_nm1 = P_n; + P_n = ((2 * k - 1) * z * P_nm1 - (k - 1) * P_nm2) / k; + } + + // Note that this metric does not directly imply that |z - z*| < tol + EXPECT_NEAR(P_n, 0.0, 1e-10); + } + } +} + +template +struct test_device_quadrature +{ + static void test() + { + const int npts = 15; + + // Create the rule with the proper allocator + const auto rule = + axom::numerics::get_gauss_legendre(npts, axom::execution_space::allocatorID()); + + // Use the rule in a lambda to integrate the volume under std::sin(pi * x) on [0, 1] + axom::ReduceSum quadrature_sum(0.0); + axom::for_all( + npts, + AXOM_LAMBDA(axom::IndexType i) { quadrature_sum += rule.weight(i) * sin(rule.node(i)); }); + + EXPECT_NEAR(quadrature_sum.get(), 0.459697694132, 1e-6); + } +}; + +TEST(numerics_quadrature, get_nodes_seq) { test_device_quadrature::test(); } + +#if defined(AXOM_RUNTIME_POLICY_USE_OPENMP) +TEST(numerics_quadrature, get_nodes_omp) { test_device_quadrature::test(); } +#endif + +#if defined(AXOM_RUNTIME_POLICY_USE_CUDA) +TEST(numerics_quadrature, get_nodes_cuda) { test_device_quadrature>::test(); } +#endif + +#if defined(AXOM_RUNTIME_POLICY_USE_HIP) +TEST(numerics_quadrature, get_nodes_hip) { test_device_quadrature>::test(); } +#endif \ No newline at end of file diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index ba99c1b9d5..bfea4f835d 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -50,6 +50,7 @@ set( primal_headers ## operators operators/clip.hpp operators/closest_point.hpp + operators/evaluate_integral.hpp operators/intersect.hpp operators/intersection_volume.hpp operators/orientation.hpp @@ -66,6 +67,7 @@ set( primal_headers operators/detail/clip_impl.hpp operators/detail/compute_moments_impl.hpp + operators/detail/evaluate_integral_impl.hpp operators/detail/fuzzy_comparators.hpp operators/detail/intersect_bezier_impl.hpp operators/detail/intersect_bounding_box_impl.hpp @@ -74,6 +76,8 @@ set( primal_headers operators/detail/intersect_ray_impl.hpp operators/detail/winding_number_2d_impl.hpp operators/detail/winding_number_2d_memoization.hpp + operators/detail/winding_number_3d_impl.hpp + operators/detail/winding_number_3d_memoization.hpp ## utils utils/ZipIndexable.hpp @@ -83,15 +87,6 @@ set( primal_headers utils/ZipVector.hpp ) -blt_list_append( - TO primal_headers - ELEMENTS operators/evaluate_integral.hpp - operators/detail/evaluate_integral_impl.hpp - operators/detail/winding_number_3d_impl.hpp - operators/detail/winding_number_3d_memoization.hpp - IF MFEM_FOUND - ) - #------------------------------------------------------------------------------ # Build and install the library #------------------------------------------------------------------------------ diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index 266c6eaa7a..7c8949ece9 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -2476,7 +2476,6 @@ class NURBSPatch return VectorType::cross_product(Du, Dv); } -#ifdef AXOM_USE_MFEM /*! * \brief Calculate the average normal for the (untrimmed) patch * @@ -2489,7 +2488,7 @@ class NURBSPatch * then computes the 2D area of that projection to get the corresponding * component of the average surface normal. * - * \note This requires the MFEM third-party library + * Evaluates the integral with Gauss-Legendre quadrature on each boundary curve * * \return The calculated mean surface normal */ @@ -2572,7 +2571,6 @@ class NURBSPatch ret_vec[1] = -ret_vec[1]; return ret_vec; } -#endif ///@} ///@{ @@ -3092,7 +3090,6 @@ class NURBSPatch return split_patches; } -#ifdef AXOM_USE_MFEM /*! * \brief Calculate the average normal for the trimmed patch * @@ -3101,7 +3098,7 @@ class NURBSPatch * Decomposes the NURBS surface into trimmed Bezier components (to ensure differentiability of the integrand) * and evaluates the integral numerically on each component using trimming curves * - * \note This requires the MFEM third-party library + * Evaluates the integral with Gauss-Legendre quadrature on each boundary curve * * \return The calculated mean surface normal */ @@ -3132,7 +3129,6 @@ class NURBSPatch return ret_vec; } -#endif //@} ///@{ diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 1152cfe76d..0bfc1b2eec 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -8,14 +8,11 @@ // Axom includes #include "axom/config.hpp" // for compile-time configuration options -#include "axom/primal.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Vector.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#else - #error "Primal's integral evaluation functions require mfem library." -#endif +#include "axom/core/numerics/quadrature.hpp" // C++ includes #include @@ -26,32 +23,35 @@ namespace primal { namespace detail { + /*! * \brief Evaluate a scalar field line integral on a single Bezier curve. * - * Evaluate the scalar field line integral with MFEM integration rule + * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double - * \param [in] quad the mfem integration rule containing nodes and weights + * \param [in] npts The number of quadrature points in the rule * \return the value of the integral */ template inline double evaluate_scalar_line_integral_component(const primal::BezierCurve& c, Lambda&& scalar_integrand, - const mfem::IntegrationRule& quad) + const int npts) { + const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); + // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad.GetNPoints(); q++) + for(int q = 0; q < npts; q++) { // Get intermediate quadrature point // at which to evaluate tangent vector - auto x_q = c.evaluate(quad.IntPoint(q).x); - auto dx_q = c.dt(quad.IntPoint(q).x); + auto x_q = c.evaluate(quad.node(q)); + auto dx_q = c.dt(quad.node(q)); - full_quadrature += quad.IntPoint(q).weight * scalar_integrand(x_q) * dx_q.norm(); + full_quadrature += quad.weight(q) * scalar_integrand(x_q) * dx_q.norm(); } return full_quadrature; @@ -60,30 +60,32 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< /*! * \brief Evaluate a vector field line integral on a single Bezier curve. * - * Evaluate the vector field line integral with MFEM integration rule + * Evaluate the vector field line integral with Gauss-Legendre quadrature * * \param [in] c the Bezier curve object * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a 2D axom vector object - * \param [in] quad the mfem integration rule containing nodes and weights + * \param [in] npts The number of quadrature points in the rule * \return the value of the integral */ template inline double evaluate_vector_line_integral_component(const primal::BezierCurve& c, Lambda&& vec_field, - const mfem::IntegrationRule& quad) + const int npts) { + const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); + // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad.GetNPoints(); q++) + for(int q = 0; q < npts; q++) { // Get intermediate quadrature point // on which to evaluate dot product - auto x_q = c.evaluate(quad.IntPoint(q).x); - auto dx_q = c.dt(quad.IntPoint(q).x); + auto x_q = c.evaluate(quad.node(q)); + auto dx_q = c.dt(quad.node(q)); auto func_val = vec_field(x_q); - full_quadrature += quad.IntPoint(q).weight * Vector::dot_product(func_val, dx_q); + full_quadrature += quad.weight(q) * Vector::dot_product(func_val, dx_q); } return full_quadrature; } @@ -101,37 +103,40 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< * \param [in] integrand the lambda function representing the integrand. * Must accept a 2D point as input and return a double * \param [in] The lower bound of integration for the antiderivatives - * \param [in] quad_Q the quadrature rule for the line integral - * \param [in] quad_P the quadrature rule for the antiderivative + * \param [in] npts_Q The number of quadrature points for the line integral + * \param [in] npts_P The number of quadrature points for the antiderivative * \return the value of the integral, which is mathematically meaningless. */ template double evaluate_area_integral_component(const primal::BezierCurve& c, Lambda&& integrand, double int_lb, - const mfem::IntegrationRule& quad_Q, - const mfem::IntegrationRule& quad_P) + const int npts_Q, + const int npts_P) { + const axom::numerics::QuadratureRule& quad_Q = axom::numerics::get_gauss_legendre(npts_Q); + const axom::numerics::QuadratureRule& quad_P = axom::numerics::get_gauss_legendre(npts_P); + // Store some intermediate values double antiderivative = 0.0; // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad_Q.GetNPoints(); q++) + for(int q = 0; q < npts_Q; q++) { // Get intermediate quadrature point // on which to evaluate antiderivative - auto x_q = c.evaluate(quad_Q.IntPoint(q).x); + 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 < quad_P.GetNPoints(); xi++) + for(int xi = 0; xi < npts_P; xi++) { // Define interior quadrature points - auto x_qxi = Point2D({x_q[0], (x_q[1] - int_lb) * quad_P.IntPoint(xi).x + int_lb}); + auto x_qxi = Point({x_q[0], (x_q[1] - int_lb) * quad_P.node(xi) + int_lb}); - antiderivative = quad_P.IntPoint(xi).weight * (x_q[1] - int_lb) * integrand(x_qxi); + antiderivative = quad_P.weight(xi) * (x_q[1] - int_lb) * integrand(x_qxi); - full_quadrature += quad_Q.IntPoint(q).weight * c.dt(quad_Q.IntPoint(q).x)[0] * -antiderivative; + full_quadrature += quad_Q.weight(q) * c.dt(quad_Q.node(q))[0] * -antiderivative; } } diff --git a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp index d77f7d8aae..e79d9f84be 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp @@ -21,11 +21,6 @@ // C++ includes #include -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#endif - namespace axom { namespace primal diff --git a/src/axom/primal/operators/detail/winding_number_3d_impl.hpp b/src/axom/primal/operators/detail/winding_number_3d_impl.hpp index d5e7202109..a73b2519ff 100644 --- a/src/axom/primal/operators/detail/winding_number_3d_impl.hpp +++ b/src/axom/primal/operators/detail/winding_number_3d_impl.hpp @@ -16,17 +16,13 @@ #include "axom/primal/operators/squared_distance.hpp" #include "axom/core/numerics/transforms.hpp" +#include "axom/core/numerics/quadrature.hpp" #include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" // C++ includes #include -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#endif - namespace axom { namespace primal diff --git a/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp index bc2262cc4c..f079ef27c4 100644 --- a/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp +++ b/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp @@ -25,19 +25,14 @@ #include "axom/primal/operators/is_convex.hpp" +#include "axom/core/numerics/quadrature.hpp" + #include #include #include #include "axom/fmt.hpp" -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#else - #error "3D GWN evaluation requires mfem library." -#endif - namespace axom { namespace primal @@ -76,9 +71,7 @@ struct TrimmingCurveQuadratureData : m_quad_npts(quad_npts) { // Generate the (cached) quadrature rules in parameter space - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad_rule = - my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * quad_npts - 1); + const numerics::QuadratureRule& gl_rule = numerics::get_gauss_legendre(quad_npts); auto& the_curve = a_patch.getTrimmingCurve(a_curve_index); @@ -93,7 +86,7 @@ struct TrimmingCurveQuadratureData m_quadrature_tangents.resize(m_quad_npts); for(int q = 0; q < m_quad_npts; ++q) { - T quad_x = quad_rule.IntPoint(q).x * m_span_length + curve_min_knot + span_offset; + T quad_x = gl_rule.node(q) * m_span_length + curve_min_knot + span_offset; Point c_eval; Vector c_Dt; @@ -114,10 +107,8 @@ struct TrimmingCurveQuadratureData { // Because the quadrature weights are identical for each trimming curve (up to a scaling factor), // we query the static rule instead of storing redundant weights - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad_rule = - my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * m_quad_npts - 1); - return quad_rule.IntPoint(idx).weight * m_span_length; + const numerics::QuadratureRule& gl_rule = numerics::get_gauss_legendre(m_quad_npts); + return gl_rule.weight(idx) * m_span_length; } double getNumPoints() const { return m_quad_npts; } diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index fb116c109e..4e3a1b8ece 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -7,13 +7,12 @@ * \file evaluate_integral.hpp * * \brief Consists of methods that evaluate integrals over regions defined - * by Bezier curves, such as 2D area integrals and scalar/vector field line integrals + * by Bezier and NURBS curves, such as 2D area integrals and scalar/vector field line integrals * - * Line integrals are computed with 1D quadrature rules supplied - * by MFEM. 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar + * Line integrals are evaluated numerically with Gauss-Legendre quadrature + * + * 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. - * - * \note This requires the MFEM third-party library */ #ifndef PRIMAL_EVAL_INTEGRAL_HPP_ @@ -21,16 +20,12 @@ // Axom includes #include "axom/config.hpp" -#include "axom/primal.hpp" -#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#else - #error "Primal's integral evaluation functions require mfem library." -#endif +#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" // C++ includes #include @@ -45,7 +40,8 @@ namespace primal * * The line integral is evaluated on each curve in the CurvedPolygon, and added * together to represent the total integral. The Polygon need not be connected. - * Uses 1D Gaussian quadrature generated by MFEM. + * + * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. @@ -59,18 +55,12 @@ double evaluate_scalar_line_integral(const primal::CurvedPolygon cpoly Lambda&& scalar_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Use the same one for every curve in the polygon - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - double total_integral = 0.0; for(int i = 0; i < cpoly.numEdges(); i++) { // Compute the line integral along each component. total_integral += - detail::evaluate_scalar_line_integral_component(cpoly[i], scalar_integrand, quad); + detail::evaluate_scalar_line_integral_component(cpoly[i], scalar_integrand, npts); } return total_integral; @@ -90,13 +80,7 @@ double evaluate_scalar_line_integral(const primal::BezierCurve& c, Lambda&& scalar_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Use the same one for every curve in the polygon - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - - return detail::evaluate_scalar_line_integral_component(c, scalar_integrand, quad); + return detail::evaluate_scalar_line_integral_component(c, scalar_integrand, npts); } /*! @@ -117,18 +101,13 @@ double evaluate_scalar_line_integral(const primal::NURBSCurve& n, Lambda&& scalar_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - // Compute the line integral along each component. auto beziers = n.extractBezier(); double total_integral = 0.0; for(int i = 0; i < beziers.size(); i++) { total_integral += - detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, quad); + detail::evaluate_scalar_line_integral_component(beziers[i], scalar_integrand, npts); } return total_integral; @@ -167,7 +146,8 @@ double evaluate_scalar_line_integral(const axom::Array cpoly Lambda&& vector_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Use the same one for every curve in the polygon - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - double total_integral = 0.0; 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], vector_integrand, quad); + detail::evaluate_vector_line_integral_component(cpoly[i], vector_integrand, npts); } return total_integral; @@ -212,12 +186,7 @@ double evaluate_vector_line_integral(const primal::BezierCurve& c, Lambda&& vector_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - - return detail::evaluate_vector_line_integral_component(c, vector_integrand, quad); + return detail::evaluate_vector_line_integral_component(c, vector_integrand, npts); } /*! @@ -238,18 +207,13 @@ double evaluate_vector_line_integral(const primal::NURBSCurve& n, Lambda&& vector_integrand, int npts) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - const mfem::IntegrationRule& quad = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); - // Compute the line integral along each component. auto beziers = n.extractBezier(); double total_integral = 0.0; for(int i = 0; i < beziers.size(); i++) { total_integral += - detail::evaluate_vector_line_integral_component(beziers[i], vector_integrand, quad); + detail::evaluate_vector_line_integral_component(beziers[i], vector_integrand, npts); } return total_integral; @@ -300,20 +264,11 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, int npts_Q, int npts_P = 0) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - if(npts_P <= 0) { npts_P = npts_Q; } - // Get the quadrature for the line integral. - // Quadrature order is equal to 2*N - 1 - const mfem::IntegrationRule& quad_Q = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1); - const mfem::IntegrationRule& quad_P = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1); - // Use minimum y-coord of control nodes as lower bound for integration double int_lb = cpoly[0][0][1]; for(int i = 0; i < cpoly.numEdges(); i++) @@ -329,7 +284,7 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, for(int i = 0; i < cpoly.numEdges(); i++) { total_integral += - detail::evaluate_area_integral_component(cpoly[i], integrand, int_lb, quad_Q, quad_P); + detail::evaluate_area_integral_component(cpoly[i], integrand, int_lb, npts_Q, npts_P); } return total_integral; @@ -359,10 +314,6 @@ double evaluate_area_integral(const axom::Array> narray int npts_Q, int npts_P = 0) { - // Generate quadrature library, defaulting to GaussLegendre quadrature. - // Quadrature order is equal to 2*npts - 1 - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - if(npts_P <= 0) { npts_P = npts_Q; @@ -373,11 +324,6 @@ double evaluate_area_integral(const axom::Array> narray return 0.0; } - // Get the quadrature for the line integral. - // Quadrature order is equal to 2*N - 1 - const mfem::IntegrationRule& quad_Q = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1); - const mfem::IntegrationRule& quad_P = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1); - // Use minimum y-coord of control nodes as lower bound for integration double int_lb = narray[0][0][1]; for(int i = 0; i < narray.size(); i++) @@ -397,7 +343,7 @@ double evaluate_area_integral(const axom::Array> narray for(int j = 0; j < beziers.size(); j++) { total_integral += - detail::evaluate_area_integral_component(beziers[j], integrand, int_lb, quad_Q, quad_P); + detail::evaluate_area_integral_component(beziers[j], integrand, int_lb, npts_Q, npts_P); } } diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp index 13024c2b54..aca97e0174 100644 --- a/src/axom/primal/operators/winding_number.hpp +++ b/src/axom/primal/operators/winding_number.hpp @@ -33,16 +33,12 @@ #include "axom/primal/operators/detail/winding_number_2d_impl.hpp" #include "axom/primal/operators/detail/winding_number_2d_memoization.hpp" +#include "axom/primal/operators/detail/winding_number_3d_impl.hpp" +#include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" + // C++ includes #include -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" - #include "axom/primal/operators/detail/winding_number_3d_impl.hpp" - #include "axom/primal/operators/detail/winding_number_3d_memoization.hpp" -#endif - namespace axom { namespace primal @@ -672,8 +668,6 @@ int winding_number(const Point& q, return std::lround(wn); } -#ifdef AXOM_USE_MFEM - /*! * \brief Computes the GWN for a 3D point wrt a 3D NURBS patch with precomputed data * @@ -866,7 +860,6 @@ axom::Array winding_number(const axom::Array>& query_arr, return winding_number(query_arr, nurbs_cache_arr, edge_tol, ls_tol, quad_tol, disk_size, EPS); } -#endif ///@} } // namespace primal diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index 845d3cf1b0..462cfe3c6f 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -21,6 +21,7 @@ set( primal_tests primal_curved_polygon.cpp primal_hexahedron.cpp primal_in_sphere.cpp + primal_integral.cpp primal_intersect.cpp primal_intersect_impl.cpp primal_knot_vector.cpp @@ -36,6 +37,7 @@ set( primal_tests primal_ray_intersect.cpp primal_segment.cpp primal_sphere.cpp + primal_solid_angle.cpp primal_split.cpp primal_squared_distance.cpp primal_surface_intersect.cpp @@ -51,13 +53,6 @@ set( primal_tests primal_clip_perf.cpp ) -blt_list_append( - TO primal_tests - ELEMENTS primal_integral.cpp - primal_solid_angle.cpp - IF MFEM_FOUND - ) - set(primal_test_depends primal fmt gtest) foreach ( test ${primal_tests} ) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 7548d74d35..2bdb18c7ff 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -6,12 +6,15 @@ #include "axom/primal.hpp" #include "axom/slic.hpp" #include "axom/fmt.hpp" -#include "axom/primal/operators/evaluate_integral.hpp" -#include "axom/primal/operators/winding_number.hpp" #include #include "gtest/gtest.h" +// MFEM includes +#ifdef AXOM_USE_MFEM + #include "mfem.hpp" +#endif + namespace primal = axom::primal; TEST(primal_integral, evaluate_area_integral) @@ -445,6 +448,32 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) } } +#ifdef AXOM_USE_MFEM +TEST(primal_integral, check_axom_mfem_quadrature_values) +{ + const int N = 3; + + for(int npts = N; npts <= N; ++npts) + { + // Generate the Axom quadrature rule + axom::numerics::QuadratureRule axom_rule = axom::numerics::get_gauss_legendre(npts); + + // Generate the MFEM quadrature rule + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& mfem_rule = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + // Check that the nodes and weights are the same between the two rules + for(int j = 0; j < npts; ++j) + { + EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, axom::numeric_limits::epsilon()); + EXPECT_NEAR(axom_rule.weight(j), + mfem_rule.IntPoint(j).weight, + axom::numeric_limits::epsilon()); + } + } +} +#endif + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); diff --git a/src/axom/quest/LinearizeCurves.cpp b/src/axom/quest/LinearizeCurves.cpp index a471887314..95c2782bc8 100644 --- a/src/axom/quest/LinearizeCurves.cpp +++ b/src/axom/quest/LinearizeCurves.cpp @@ -143,9 +143,7 @@ static void appendPoints(LinearizeCurves::SegmentMesh* mesh, SegmentArray& S, do } } -#ifdef AXOM_USE_MFEM /// Compute the arc length of a curve by numerical quadrature -// This currently uses mfem's quadrature weights double computeArcLength(const LinearizeCurves::NURBSCurve& nurbs, int npts) { double arcLength = 0.; @@ -156,29 +154,6 @@ double computeArcLength(const LinearizeCurves::NURBSCurve& nurbs, int npts) } return arcLength; } -#else -/// Compute the arc length of a curve by discretizing into linear segments and adding their lengths -double computeArcLength(const LinearizeCurves::NURBSCurve& nurbs, int nSamples) -{ - using axom::utilities::lerp; - - // Get the contour start/end parameters. - const auto knots = nurbs.getKnots(); - const double u0 = knots[0]; - const double u1 = knots[knots.getNumKnots() - 1]; - - double arcLength = 0.; - PointType prev = nurbs.evaluate(u0); - for(int i = 1; i <= nSamples; ++i) - { - const double u = lerp(u0, u1, i / static_cast(nSamples)); - PointType cur = nurbs.evaluate(u); - arcLength += sqrt(primal::squared_distance(prev, cur)); - axom::utilities::swap(prev, cur); - } - return arcLength; -} -#endif //--------------------------------------------------------------------------- #ifdef AXOM_DEBUG_WRITE_LINES @@ -408,13 +383,8 @@ void LinearizeCurves::getLinearMeshNonUniform(LinearizeCurves::CurveArrayView cu if(knots.getNumKnotSpans() > 1) { constexpr int MAX_NUMBER_OF_SAMPLES = 2000; - -#ifdef AXOM_USE_MFEM constexpr int quadrature_order = 30; const double hiCurveLen = internal::computeArcLength(nurbs, quadrature_order); -#else - const double hiCurveLen = internal::computeArcLength(nurbs, MAX_NUMBER_OF_SAMPLES); -#endif // The initial curve length. double curveLength = first.length; @@ -514,13 +484,8 @@ double LinearizeCurves::revolvedVolume(const LinearizeCurves::NURBSCurve& nurbs, using PointType = axom::primal::Point; using VectorType = axom::primal::Vector; - // Use 5-point Gauss Quadrature. - const double X[] = {-0.906179845938664, -0.538469310105683, 0, 0.538469310105683, 0.906179845938664}; - const double W[] = {0.23692688505618908, - 0.47862867049936647, - 0.5688888888888889, - 0.47862867049936647, - 0.23692688505618908}; + // Use 5-point Gauss-Legendre Quadrature. + numerics::QuadratureRule quad_rule = numerics::get_gauss_legendre(5); // Make a transform with no translation. We use this to transform // the derivative since we want to permit scaling and rotation but @@ -545,12 +510,12 @@ double LinearizeCurves::revolvedVolume(const LinearizeCurves::NURBSCurve& nurbs, #endif // Approximate the integral "Int pi*x'(u)*y(u)^2du" using quadrature. - constexpr double scale = M_PI * ((bd - ad) / 2.); + constexpr double scale = M_PI * (bd - ad); double sum = 0.; for(size_t i = 0; i < 5; i++) { - // Map quad point x value [-1,1] to [a,b]. - const double u = X[i] * ((bd - ad) / 2.) + ((bd + ad) / 2.); + // Map quad point x value [0,1] to [a,b]. + const double u = quad_rule.node(i) * (bd - ad) + ad; // Compute y(u) to get radius PointType eval; @@ -569,7 +534,7 @@ double LinearizeCurves::revolvedVolume(const LinearizeCurves::NURBSCurve& nurbs, #endif // Accumulate weight times dx*r^2. - sum += W[i] * xp * (r * r); + sum += quad_rule.weight(i) * xp * (r * r); } // Guard against volumes being negative (if the curve went the wrong way) vol += axom::utilities::abs(scale * sum);