From 6be497c82d60b4f37a3fa2e767462ec6d8cb529a Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 30 Sep 2025 16:06:11 -0700 Subject: [PATCH 01/28] First commit of search algorithm --- src/axom/core/CMakeLists.txt | 2 + src/axom/core/numerics/gauss_legendre.cpp | 79 +++++++++++++++++++++ src/axom/core/numerics/gauss_legendre.hpp | 27 +++++++ src/axom/core/tests/CMakeLists.txt | 1 + src/axom/core/tests/core_serial_main.cpp | 1 + src/axom/core/tests/numerics_quadrature.hpp | 23 ++++++ 6 files changed, 133 insertions(+) create mode 100644 src/axom/core/numerics/gauss_legendre.cpp create mode 100644 src/axom/core/numerics/gauss_legendre.hpp create mode 100644 src/axom/core/tests/numerics_quadrature.hpp diff --git a/src/axom/core/CMakeLists.txt b/src/axom/core/CMakeLists.txt index 89fa639774..1c268b5da7 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/gauss_legendre.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/gauss_legendre.cpp Path.cpp Types.cpp diff --git a/src/axom/core/numerics/gauss_legendre.cpp b/src/axom/core/numerics/gauss_legendre.cpp new file mode 100644 index 0000000000..cd99d02583 --- /dev/null +++ b/src/axom/core/numerics/gauss_legendre.cpp @@ -0,0 +1,79 @@ +// 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/config.hpp" +#include "axom/core/numerics/gauss_legendre.hpp" + +#include + +namespace axom +{ +namespace numerics +{ + +void compute_rule(int np, double* nodes, double* weights) +{ + if(np == 1) + { + nodes[0] = 0.5; + weights[0] = 1.0; + return; + } + if(np == 2) + { + nodes[0] = 0.21132486540518711775; + nodes[1] = 0.78867513459481288225; + + weights[0] = weights[1] = 0.5; + return; + } + + const int n = np; + const int m = (n + 1) / 2; + + // Nodes are mirrored across x = 0.5, so only need to evaluate half + for(int i = 1; i <= m; ++i) + { + double z = std::cos(M_PI * (i - 0.25) / (n + 0.5)); + double pp, p1, dz, xi = 0.0; + + bool done = false; + while(true) + { + double p2 = 1.0; + p1 = z; + for(int j = 2; j <= n; ++j) + { + double p3 = p2; + p2 = p1; + p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j; + } + // p1 is Legendre polynomial + + pp = n * (z * p1 - p2) / (z * z - 1); + + if(done) + { + break; + } + + dz = p1 / pp; + + if(std::fabs(dz) < 1e-16) + { + done = true; + xi = ((1 - z) + dz) / 2; + } + + z -= dz; + } + + nodes[i - 1] = xi; + nodes[n - i] = 1.0 - xi; + weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * pp * pp); + } +} +} /* end namespace numerics */ +} /* end namespace axom */ \ No newline at end of file diff --git a/src/axom/core/numerics/gauss_legendre.hpp b/src/axom/core/numerics/gauss_legendre.hpp new file mode 100644 index 0000000000..103f7fcb2e --- /dev/null +++ b/src/axom/core/numerics/gauss_legendre.hpp @@ -0,0 +1,27 @@ +// 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_GAUSS_LEGENDGRE_HPP_ +#define AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ + +#include + +/*! + * \file gauss_legendre.hpp + * The functions declared in this header file find the nodes and weights of + * arbitrary order Gauss-Legendre quadrature rules + */ + +namespace axom +{ +namespace numerics +{ + +void compute_rule(int np, double* nodes, double* weights); + +} /* end namespace numerics */ +} /* end namespace axom */ + +#endif // AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ diff --git a/src/axom/core/tests/CMakeLists.txt b/src/axom/core/tests/CMakeLists.txt index 60bc0779d8..1ed1d6e771 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..9890903610 --- /dev/null +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -0,0 +1,23 @@ +// 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/core/numerics/gauss_legendre.hpp" + +TEST(numerics_quadrature, generate_rules) +{ + const int N = 3; + + double nodes[N]; + double weights[N]; + + axom::numerics::compute_rule(N, nodes, weights); + + for(int i = 0; i < N; ++i) + { + std::cout << nodes[i] << " " << weights[i] << std::endl; + } +} From 67fdfc68724ba8d6239975ccf3d6d75f8f5dc323 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 1 Oct 2025 09:12:23 -0700 Subject: [PATCH 02/28] Add test --- src/axom/core/numerics/gauss_legendre.cpp | 32 +++++++++----- src/axom/core/tests/numerics_quadrature.hpp | 48 +++++++++++++++++++-- 2 files changed, 66 insertions(+), 14 deletions(-) diff --git a/src/axom/core/numerics/gauss_legendre.cpp b/src/axom/core/numerics/gauss_legendre.cpp index cd99d02583..d9930a950b 100644 --- a/src/axom/core/numerics/gauss_legendre.cpp +++ b/src/axom/core/numerics/gauss_legendre.cpp @@ -36,43 +36,55 @@ void compute_rule(int np, double* nodes, double* weights) // 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 distirbuted 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, p1, dz, xi = 0.0; + double Pp_n, P_n, dz, xi = 0.0; bool done = false; while(true) { - double p2 = 1.0; - p1 = z; + // 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 p3 = p2; - p2 = p1; - p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j; + double P_nm2 = P_nm1; + P_nm1 = P_n; + P_n = ((2 * j - 1) * z * P_nm1 - (j - 1) * P_nm2) / j; } - // p1 is Legendre polynomial - pp = n * (z * p1 - p2) / (z * z - 1); + // 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; } - dz = p1 / pp; + // Compute the Newton method step size + dz = P_n / Pp_n; if(std::fabs(dz) < 1e-16) { 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; - weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * pp * pp); + + // Evaluate the weight as w_i = 2 / (1 - z^2) / P'_n(z)^2, with z \in [-1, 1] + weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); } } } /* end namespace numerics */ diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 9890903610..8974060571 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -6,18 +6,58 @@ #include "gtest/gtest.h" #include "axom/core/numerics/gauss_legendre.hpp" +#include TEST(numerics_quadrature, generate_rules) { - const int N = 3; + const int N = 50; double nodes[N]; double weights[N]; - axom::numerics::compute_rule(N, nodes, weights); + // Define a collection of random coefficients for a polynomial + double coeffs[2 * N - 1]; + for(int j = 0; j < 2 * N - 1; ++j) + { + coeffs[j] = axom::utilities::random_real(-5.0, 5.0); + std::cout << coeffs[j] << std::endl; + } - for(int i = 0; i < N; ++i) + // Test that the rules provide exact integration for polynomials of degree 2n - 1 + for(int npts = 1; npts < N; ++npts) { - std::cout << nodes[i] << " " << weights[i] << std::endl; + int degree = 2 * npts - 1; + + // 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 using the quadrature rule + axom::numerics::compute_rule(npts, nodes, weights); + + // 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 += weights[j] * eval_polynomial(nodes[j]); + } + + std::cout << std::setprecision(20); + std::cout << "=====" << std::endl; + std::cout << analytic_result << std::endl; + std::cout << quadrature_result << std::endl; + std::cout << "=====" << std::endl; } } From a0b488f42e8b53260606130232d74fd9020d71b2 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 1 Oct 2025 16:35:38 -0700 Subject: [PATCH 03/28] Add cache methods, rename to quadrature.hpp --- src/axom/core/CMakeLists.txt | 4 +- src/axom/core/numerics/gauss_legendre.hpp | 27 --------- .../{gauss_legendre.cpp => quadrature.cpp} | 59 ++++++++++++------- src/axom/core/numerics/quadrature.hpp | 48 +++++++++++++++ src/axom/core/tests/numerics_quadrature.hpp | 12 ++-- 5 files changed, 93 insertions(+), 57 deletions(-) delete mode 100644 src/axom/core/numerics/gauss_legendre.hpp rename src/axom/core/numerics/{gauss_legendre.cpp => quadrature.cpp} (58%) create mode 100644 src/axom/core/numerics/quadrature.hpp diff --git a/src/axom/core/CMakeLists.txt b/src/axom/core/CMakeLists.txt index 1c268b5da7..7c78bc2ec6 100644 --- a/src/axom/core/CMakeLists.txt +++ b/src/axom/core/CMakeLists.txt @@ -42,7 +42,7 @@ set(core_headers numerics/eigen_solve.hpp numerics/eigen_sort.hpp numerics/floating_point_limits.hpp - numerics/gauss_legendre.hpp + numerics/quadrature.hpp numerics/jacobi_eigensolve.hpp numerics/linear_solve.hpp numerics/matvecops.hpp @@ -107,7 +107,7 @@ set(core_sources ${PROJECT_BINARY_DIR}/axom/core/utilities/About.cpp numerics/polynomial_solvers.cpp - numerics/gauss_legendre.cpp + numerics/quadrature.cpp Path.cpp Types.cpp diff --git a/src/axom/core/numerics/gauss_legendre.hpp b/src/axom/core/numerics/gauss_legendre.hpp deleted file mode 100644 index 103f7fcb2e..0000000000 --- a/src/axom/core/numerics/gauss_legendre.hpp +++ /dev/null @@ -1,27 +0,0 @@ -// 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_GAUSS_LEGENDGRE_HPP_ -#define AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ - -#include - -/*! - * \file gauss_legendre.hpp - * The functions declared in this header file find the nodes and weights of - * arbitrary order Gauss-Legendre quadrature rules - */ - -namespace axom -{ -namespace numerics -{ - -void compute_rule(int np, double* nodes, double* weights); - -} /* end namespace numerics */ -} /* end namespace axom */ - -#endif // AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ diff --git a/src/axom/core/numerics/gauss_legendre.cpp b/src/axom/core/numerics/quadrature.cpp similarity index 58% rename from src/axom/core/numerics/gauss_legendre.cpp rename to src/axom/core/numerics/quadrature.cpp index d9930a950b..0862febec6 100644 --- a/src/axom/core/numerics/gauss_legendre.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -3,9 +3,11 @@ // // SPDX-License-Identifier: (BSD-3-Clause) -#include "axom/config.hpp" -#include "axom/core/numerics/gauss_legendre.hpp" +#include "axom/core/Array.hpp" +#include "axom/core/numerics/quadrature.hpp" +// For math constants and includes +#include "axom/config.hpp" #include namespace axom @@ -13,30 +15,32 @@ namespace axom namespace numerics { -void compute_rule(int np, double* nodes, double* weights) +QuadratureRule compute_gauss_legendre(int npts) { - if(np == 1) + QuadratureRule rule(npts); + + if(npts == 1) { - nodes[0] = 0.5; - weights[0] = 1.0; - return; + rule.m_nodes[0] = 0.5; + rule.m_weights[0] = 1.0; + return rule; } - if(np == 2) + if(npts == 2) { - nodes[0] = 0.21132486540518711775; - nodes[1] = 0.78867513459481288225; + rule.m_nodes[0] = 0.21132486540518711775; + rule.m_nodes[1] = 0.78867513459481288225; - weights[0] = weights[1] = 0.5; - return; + rule.m_weights[0] = rule.m_weights[1] = 0.5; + return rule; } - const int n = np; - const int m = (n + 1) / 2; + 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, + // Each node is the root of a Legendre polynomial, // which are approximately uniformly distirbuted 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)); @@ -47,8 +51,8 @@ void compute_rule(int np, double* nodes, double* weights) { // 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 + 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; @@ -80,12 +84,27 @@ void compute_rule(int np, double* nodes, double* weights) z -= dz; } - nodes[i - 1] = xi; - nodes[n - i] = 1.0 - xi; + rule.m_nodes[i - 1] = xi; + rule.m_nodes[n - i] = 1.0 - xi; // Evaluate the weight as w_i = 2 / (1 - z^2) / P'_n(z)^2, with z \in [-1, 1] - weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); + rule.m_weights[i - 1] = rule.m_weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); } + + return rule; } + +const QuadratureRule& get_gauss_legendre(int npts) +{ + // Define a static map that stores the GL quadrature rule for a given order + static std::map rule_library; + if(rule_library.find(npts) == rule_library.end()) + { + rule_library[npts] = compute_gauss_legendre(npts); + } + + return rule_library[npts]; +} + } /* end namespace numerics */ } /* end namespace axom */ \ No newline at end of file diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp new file mode 100644 index 0000000000..757af89bb9 --- /dev/null +++ b/src/axom/core/numerics/quadrature.hpp @@ -0,0 +1,48 @@ +// 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_GAUSS_LEGENDGRE_HPP_ +#define AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ + +#include "axom/core/Array.hpp" + +/*! + * \file gauss_legendre.hpp + * The functions declared in this header file find the nodes and weights of + * arbitrary order Gauss-Legendre quadrature rules + */ + +namespace axom +{ +namespace numerics +{ + +//! \brief Class to store quadature nodes and weights +class QuadratureRule +{ + // Define friend functions so rules can only be created via compute_rule() methods + friend QuadratureRule compute_gauss_legendre(int npts); + +public: + QuadratureRule() = default; + + double node(size_t idx) const { return m_nodes[idx]; }; + double weight(size_t idx) const { return m_weights[idx]; }; + +private: + QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; + + axom::Array m_nodes; + axom::Array m_weights; +}; + +QuadratureRule compute_gauss_legendre(int npts); + +const QuadratureRule& get_gauss_legendre(int npts); + +} /* end namespace numerics */ +} /* end namespace axom */ + +#endif // AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 8974060571..5e4e54da94 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -5,22 +5,18 @@ #include "gtest/gtest.h" -#include "axom/core/numerics/gauss_legendre.hpp" +#include "axom/core/numerics/quadrature.hpp" #include TEST(numerics_quadrature, generate_rules) { const int N = 50; - double nodes[N]; - double weights[N]; - // Define a collection of random coefficients for a polynomial double coeffs[2 * N - 1]; for(int j = 0; j < 2 * N - 1; ++j) { coeffs[j] = axom::utilities::random_real(-5.0, 5.0); - std::cout << coeffs[j] << std::endl; } // Test that the rules provide exact integration for polynomials of degree 2n - 1 @@ -36,7 +32,7 @@ TEST(numerics_quadrature, generate_rules) } // Evaluate using the quadrature rule - axom::numerics::compute_rule(npts, nodes, weights); + auto rule = axom::numerics::compute_gauss_legendre(npts); // Evaluate the polynomial using Horner's rule auto eval_polynomial = [°ree, &coeffs](double x) { @@ -50,8 +46,8 @@ TEST(numerics_quadrature, generate_rules) double quadrature_result = 0.0; for(int j = 0; j < npts; ++j) - { - quadrature_result += weights[j] * eval_polynomial(nodes[j]); + { + quadrature_result += rule.weight(j) * eval_polynomial(rule.node(j)); } std::cout << std::setprecision(20); From fe2358086e105f2d4df575c566bdb1817e116b61 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 1 Oct 2025 18:18:23 -0700 Subject: [PATCH 04/28] Update comments, test --- src/axom/core/numerics/quadrature.hpp | 46 +++++++++++-- src/axom/core/tests/numerics_quadrature.hpp | 75 ++++++++++++++++----- 2 files changed, 100 insertions(+), 21 deletions(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 757af89bb9..693c259f1b 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -3,15 +3,15 @@ // // SPDX-License-Identifier: (BSD-3-Clause) -#ifndef AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ -#define AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ +#ifndef AXOM_NUMERICS_QUADRATURE_HPP_ +#define AXOM_NUMERICS_QUADRATURE_HPP_ #include "axom/core/Array.hpp" /*! - * \file gauss_legendre.hpp + * \file quadrature.hpp * The functions declared in this header file find the nodes and weights of - * arbitrary order Gauss-Legendre quadrature rules + * arbitrary order quadrature rules */ namespace axom @@ -19,7 +19,11 @@ namespace axom namespace numerics { -//! \brief Class to store quadature nodes and weights +/*! + * \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 compute_rule() methods @@ -28,21 +32,51 @@ class QuadratureRule public: QuadratureRule() = default; + //! \brief Accessor for quadrature nodes double node(size_t idx) const { return m_nodes[idx]; }; + + //! \brief Accessor for quadrature weights double weight(size_t idx) const { return m_weights[idx]; }; private: + //! \brief Private constructor for use in compute_() methods. Only allocates memory QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; axom::Array m_nodes; axom::Array m_weights; }; +/*! + * \brief Computes a 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 This method constructs the points by scratch each time, without caching + * \sa get_gauss_legendre(int) + * + * \return The `QuadratureRule` object which contains axom::Array's of nodes and weights + */ QuadratureRule compute_gauss_legendre(int npts); +/*! + * \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::Array's of nodes and weights + */ const QuadratureRule& get_gauss_legendre(int npts); } /* end namespace numerics */ } /* end namespace axom */ -#endif // AXOM_NUMERICS_GAUSS_LEGENDGRE_HPP_ +#endif // AXOM_NUMERICS_QUADRATURE_HPP_ diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 5e4e54da94..3d40e1faf7 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -8,22 +8,26 @@ #include "axom/core/numerics/quadrature.hpp" #include -TEST(numerics_quadrature, generate_rules) +TEST(numerics_quadrature, gauss_legendre_math_check) { - const int N = 50; + const int N = 200; - // Define a collection of random coefficients for a polynomial double coeffs[2 * N - 1]; - for(int j = 0; j < 2 * N - 1; ++j) - { - coeffs[j] = axom::utilities::random_real(-5.0, 5.0); - } // 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) @@ -31,9 +35,6 @@ TEST(numerics_quadrature, generate_rules) analytic_result += coeffs[j] / (j + 1); } - // Evaluate using the quadrature rule - auto rule = axom::numerics::compute_gauss_legendre(npts); - // Evaluate the polynomial using Horner's rule auto eval_polynomial = [°ree, &coeffs](double x) { double result = coeffs[degree - 1]; @@ -50,10 +51,54 @@ TEST(numerics_quadrature, generate_rules) quadrature_result += rule.weight(j) * eval_polynomial(rule.node(j)); } - std::cout << std::setprecision(20); - std::cout << "=====" << std::endl; - std::cout << analytic_result << std::endl; - std::cout << quadrature_result << std::endl; - std::cout << "=====" << std::endl; + 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); + } } } + +TEST(numerics_quadrature, gauss_legendre_cache_check) +{ + // The first two rules are put in the cache + const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); + const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); + + // The third is not + const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); + + // Check that the two rules are located in the same place in memory, and the third isn't + EXPECT_EQ(&first_rule, &second_rule); + EXPECT_NE(&first_rule, &third_rule); +} \ No newline at end of file From 651fe3fa0847cea415baa34cc7233205aa216147 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 2 Oct 2025 17:53:16 -0700 Subject: [PATCH 05/28] Update integral methods to use MFEM if available --- src/axom/primal/CMakeLists.txt | 11 +- src/axom/primal/geometry/NURBSPatch.hpp | 10 +- .../detail/evaluate_integral_impl.hpp | 106 +++++++++++++----- .../detail/winding_number_3d_impl.hpp | 6 +- .../detail/winding_number_3d_memoization.hpp | 21 +--- .../primal/operators/evaluate_integral.hpp | 92 ++++----------- src/axom/primal/operators/winding_number.hpp | 13 +-- src/axom/primal/tests/CMakeLists.txt | 9 +- 8 files changed, 117 insertions(+), 151 deletions(-) diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index ba99c1b9d5..4dc8180b3d 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 @@ -83,15 +85,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..f214c397d5 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,8 @@ 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 line integrals with MFEM integration rules if available, + * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \return The calculated mean surface normal */ @@ -2572,7 +2572,6 @@ class NURBSPatch ret_vec[1] = -ret_vec[1]; return ret_vec; } -#endif ///@} ///@{ @@ -3092,7 +3091,6 @@ class NURBSPatch return split_patches; } -#ifdef AXOM_USE_MFEM /*! * \brief Calculate the average normal for the trimmed patch * @@ -3101,7 +3099,8 @@ 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 line integrals with MFEM integration rules if available, + * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \return The calculated mean surface normal */ @@ -3132,7 +3131,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..bf676dc908 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -10,11 +10,12 @@ #include "axom/config.hpp" // for compile-time configuration options #include "axom/primal.hpp" -// MFEM includes +// Use 1D Gauss-Legendre quadrature generated by MFEM if available, +// or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) #ifdef AXOM_USE_MFEM #include "mfem.hpp" #else - #error "Primal's integral evaluation functions require mfem library." + #include "axom/core/numerics/quadrature.hpp" #endif // C++ includes @@ -26,32 +27,65 @@ namespace primal { namespace detail { + +class QuadratureRuleView +{ +#ifdef AXOM_USE_MFEM +public: + QuadratureRuleView(const mfem::IntegrationRule& rule) : m_quad_rule(rule) { } + + double node(int idx) const { return m_quad_rule.IntPoint(idx).x; } + double weight(int idx) const { return m_quad_rule.IntPoint(idx).weight; } + +private: + const mfem::IntegrationRule& m_quad_rule; +#else +public: + QuadratureRuleView(const axom::numerics::QuadratureRule& rule) : m_quad_rule(rule) { } + + double node(int idx) const { return m_quad_rule.node(idx); } + double weight(int idx) const { return m_quad_rule.weight(idx); } + +private: + const axom::numerics::QuadratureRule& m_quad_rule; +#endif +}; + /*! * \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 MFEM integration rule if available, + * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \param [in] c the Bezier curve object * \param [in] integrand the lambda function representing the integrand. * 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) { +#ifdef AXOM_USE_MFEM + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + QuadratureRuleView quad(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); +#else + SLIC_WARNING_IF(npts > 100, "Axom's Gauss-Legendre nodes may be less precise for n > 100"); + QuadratureRuleView quad(axom::numerics::get_gauss_legendre(npts)); +#endif + // 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 +94,39 @@ 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 MFEM integration rule if available, + * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \param [in] c the Bezier curve object * \param [in] integrand the lambda function representing the integrand. * 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) { +#ifdef AXOM_USE_MFEM + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + QuadratureRuleView quad(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); +#else + SLIC_WARNING_IF(npts > 100, "Axom's Gauss-Legendre nodes may be less precise for n > 100"); + QuadratureRuleView quad(axom::numerics::get_gauss_legendre(npts)); +#endif + // 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 +144,48 @@ 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) { +#ifdef AXOM_USE_MFEM + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + QuadratureRuleView quad_Q(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); + QuadratureRuleView quad_P(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1)); +#else + SLIC_WARNING_IF(npts_Q > 100 || npts_P > 100, + "Axom's Gauss-Legendre nodes may be less precise for n > 100"); + QuadratureRuleView quad_Q(axom::numerics::get_gauss_legendre(npts_Q)); + QuadratureRuleView quad_P(axom::numerics::get_gauss_legendre(npts_P)); +#endif + // 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 = Point2D({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_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..915ef70f3f 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -7,13 +7,13 @@ * \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 computed with 1D quadrature rules supplied by MFEM if available, + * or by methods in "core/numerics/quadrature" otherwise, which are less robust for + * high numbers of quadrature points (> ~100) + * 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_ @@ -25,13 +25,6 @@ #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#else - #error "Primal's integral evaluation functions require mfem library." -#endif - // C++ includes #include @@ -45,7 +38,9 @@ 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. + * + * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, + * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) * * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. @@ -59,18 +54,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 +79,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 +100,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 +145,9 @@ double evaluate_scalar_line_integral(const axom::Array ~100) * * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. @@ -181,18 +161,12 @@ double evaluate_vector_line_integral(const primal::CurvedPolygon 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} ) From d90bab02803a4a311d11935fc8737b31cd6bb3ff Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 2 Oct 2025 17:59:04 -0700 Subject: [PATCH 06/28] Improve comments --- src/axom/core/numerics/quadrature.cpp | 28 +++++++++++++++++++++ src/axom/core/numerics/quadrature.hpp | 2 ++ src/axom/core/tests/numerics_quadrature.hpp | 2 +- 3 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 0862febec6..0692a8e402 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -15,6 +15,21 @@ namespace axom namespace numerics { +/*! + * \brief Computes a 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 + * + * 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) + * + * \return The `QuadratureRule` object which contains axom::Array's of nodes and weights + */ QuadratureRule compute_gauss_legendre(int npts) { QuadratureRule rule(npts); @@ -94,6 +109,19 @@ QuadratureRule compute_gauss_legendre(int npts) return rule; } +/*! + * \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::Array's of nodes and weights + */ const QuadratureRule& get_gauss_legendre(int npts) { // Define a static map that stores the GL quadrature rule for a given order diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 693c259f1b..4b4ea942c5 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -54,6 +54,8 @@ class QuadratureRule * 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) * diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 3d40e1faf7..659708d656 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -101,4 +101,4 @@ TEST(numerics_quadrature, gauss_legendre_cache_check) // Check that the two rules are located in the same place in memory, and the third isn't EXPECT_EQ(&first_rule, &second_rule); EXPECT_NE(&first_rule, &third_rule); -} \ No newline at end of file +} From b5b199f67af0e208f4c2b63913eebc26c8cc1180 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 2 Oct 2025 18:20:48 -0700 Subject: [PATCH 07/28] Re-add missing headers from CMakeLists.txt --- src/axom/core/numerics/quadrature.cpp | 2 +- src/axom/primal/CMakeLists.txt | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 0692a8e402..775fd8820b 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -135,4 +135,4 @@ const QuadratureRule& get_gauss_legendre(int npts) } } /* end namespace numerics */ -} /* end namespace axom */ \ No newline at end of file +} /* end namespace axom */ diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 4dc8180b3d..bfea4f835d 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -76,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 From ac385b3d76b0871656bb3d1ebecd9fb722d378d8 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 09:11:56 -0700 Subject: [PATCH 08/28] Temporarily revert test --- src/axom/core/tests/numerics_quadrature.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 659708d656..f76f98b6fd 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -99,6 +99,6 @@ TEST(numerics_quadrature, gauss_legendre_cache_check) const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); // Check that the two rules are located in the same place in memory, and the third isn't - EXPECT_EQ(&first_rule, &second_rule); - EXPECT_NE(&first_rule, &third_rule); + //EXPECT_EQ(&first_rule, &second_rule); + //EXPECT_NE(&first_rule, &third_rule); } From e90dc6c1bddd2bc7b22c79dfe728c979be21eddd Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 10:39:56 -0700 Subject: [PATCH 09/28] Move test to other suite --- src/axom/primal/tests/primal_integral.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 7548d74d35..c5df1cef05 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -14,6 +14,20 @@ namespace primal = axom::primal; +TEST(primal_integral, gauss_legendre_cache_check) +{ + // The first two rules are put in the cache + const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); + const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); + + // The third is not + const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); + + // Check that the two rules are located in the same place in memory, and the third isn't + //EXPECT_EQ(&first_rule, &second_rule); + //EXPECT_NE(&first_rule, &third_rule); +} + TEST(primal_integral, evaluate_area_integral) { using Point2D = primal::Point; From 68002f3a57dc97427fb2e2127ce20103c80fbb3c Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 12:31:28 -0700 Subject: [PATCH 10/28] More testing --- src/axom/core/tests/numerics_quadrature.hpp | 2 +- src/axom/primal/tests/primal_integral.cpp | 14 -------------- 2 files changed, 1 insertion(+), 15 deletions(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index f76f98b6fd..c013e5455c 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -93,7 +93,7 @@ TEST(numerics_quadrature, gauss_legendre_cache_check) { // The first two rules are put in the cache const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); - const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); + //const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); // The third is not const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index c5df1cef05..7548d74d35 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -14,20 +14,6 @@ namespace primal = axom::primal; -TEST(primal_integral, gauss_legendre_cache_check) -{ - // The first two rules are put in the cache - const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); - const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); - - // The third is not - const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); - - // Check that the two rules are located in the same place in memory, and the third isn't - //EXPECT_EQ(&first_rule, &second_rule); - //EXPECT_NE(&first_rule, &third_rule); -} - TEST(primal_integral, evaluate_area_integral) { using Point2D = primal::Point; From cc857c6b7c36bf4f8bf5068aff16a4e7efe44fa6 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 14:45:24 -0700 Subject: [PATCH 11/28] Another test check --- src/axom/core/tests/numerics_quadrature.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index c013e5455c..0c275a0c61 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -17,6 +17,8 @@ TEST(numerics_quadrature, gauss_legendre_math_check) // Test that the rules provide exact integration for polynomials of degree 2n - 1 for(int npts = 1; npts < N; ++npts) { + std::cout << "Testing with " << npts << std::endl; + // Evaluate using the quadrature rule auto rule = axom::numerics::get_gauss_legendre(npts); int degree = 2 * npts - 1; @@ -92,11 +94,11 @@ TEST(numerics_quadrature, gauss_legendre_math_check) TEST(numerics_quadrature, gauss_legendre_cache_check) { // The first two rules are put in the cache - const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); + //const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); //const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); // The third is not - const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); + //const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); // Check that the two rules are located in the same place in memory, and the third isn't //EXPECT_EQ(&first_rule, &second_rule); From e300d22d4862e60d93fdc142cf010b689f551c52 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 16:16:06 -0700 Subject: [PATCH 12/28] Try to restore the test with necessary includes --- src/axom/core/tests/numerics_quadrature.hpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 0c275a0c61..9793f45925 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -6,7 +6,7 @@ #include "gtest/gtest.h" #include "axom/core/numerics/quadrature.hpp" -#include +#include "axom/core/utilities/Utilities.hpp" TEST(numerics_quadrature, gauss_legendre_math_check) { @@ -17,8 +17,6 @@ TEST(numerics_quadrature, gauss_legendre_math_check) // Test that the rules provide exact integration for polynomials of degree 2n - 1 for(int npts = 1; npts < N; ++npts) { - std::cout << "Testing with " << npts << std::endl; - // Evaluate using the quadrature rule auto rule = axom::numerics::get_gauss_legendre(npts); int degree = 2 * npts - 1; @@ -94,13 +92,13 @@ TEST(numerics_quadrature, gauss_legendre_math_check) TEST(numerics_quadrature, gauss_legendre_cache_check) { // The first two rules are put in the cache - //const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); - //const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); + const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); + const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); // The third is not - //const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); + const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); // Check that the two rules are located in the same place in memory, and the third isn't - //EXPECT_EQ(&first_rule, &second_rule); - //EXPECT_NE(&first_rule, &third_rule); + EXPECT_EQ(&first_rule, &second_rule); + EXPECT_NE(&first_rule, &third_rule); } From be3c09a21d76ae7132c9a197a0bbb8abdf526fa9 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 3 Oct 2025 16:22:35 -0700 Subject: [PATCH 13/28] Fix some includes --- src/axom/primal/operators/detail/evaluate_integral_impl.hpp | 6 ++++-- src/axom/primal/operators/evaluate_integral.hpp | 5 ++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index bf676dc908..312aa70d2f 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -8,7 +8,9 @@ // 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" // Use 1D Gauss-Legendre quadrature generated by MFEM if available, // or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) @@ -181,7 +183,7 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, 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.node(xi) + int_lb}); + auto x_qxi = Point({x_q[0], (x_q[1] - int_lb) * quad_P.node(xi) + int_lb}); antiderivative = quad_P.weight(xi) * (x_q[1] - int_lb) * integrand(x_qxi); diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 915ef70f3f..fc75f31376 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -21,7 +21,10 @@ // Axom includes #include "axom/config.hpp" -#include "axom/primal.hpp" + +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/NURBSCurve.hpp" #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" From 49677e98c947db563b8476a83eb1e744cc73f0c5 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 11:32:49 -0700 Subject: [PATCH 14/28] Remove the rest of MFEM from evaluate_integral --- RELEASE-NOTES.md | 1 + src/axom/core/numerics/quadrature.cpp | 11 ++++ src/axom/core/numerics/quadrature.hpp | 6 +- src/axom/core/tests/numerics_quadrature.hpp | 4 +- .../detail/evaluate_integral_impl.hpp | 65 ++----------------- .../primal/operators/evaluate_integral.hpp | 11 ++-- src/axom/primal/tests/primal_integral.cpp | 24 +++++++ 7 files changed, 54 insertions(+), 68 deletions(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 8381f94f5c..d4e38e473b 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -75,6 +75,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/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 775fd8820b..a5740b0212 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -48,6 +48,17 @@ QuadratureRule compute_gauss_legendre(int npts) rule.m_weights[0] = rule.m_weights[1] = 0.5; return rule; } + if(npts == 3) + { + rule.m_nodes[0] = 0.11270166537925831148207345; + rule.m_nodes[1] = 0.5; + rule.m_nodes[2] = 0.88729833462074168851792655; + + rule.m_weights[0] = 0.2777777777777777777777778; + rule.m_weights[1] = 0.4444444444444444444444444; + rule.m_weights[2] = 0.2777777777777777777777778; + return rule; + } const int n = npts; const int m = (npts + 1) / 2; diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 4b4ea942c5..b07224c52f 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -38,12 +38,16 @@ class QuadratureRule //! \brief Accessor for quadrature weights double weight(size_t idx) const { return m_weights[idx]; }; + //! \brief Accessofr for the size of the quadrature rule + int getNumPoints() const { return m_npts; } + private: //! \brief Private constructor for use in compute_() methods. Only allocates memory - QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; + QuadratureRule(int npts) : m_nodes(npts), m_weights(npts), m_npts(npts) { }; axom::Array m_nodes; axom::Array m_weights; + int m_npts; }; /*! diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 9793f45925..d1e097e38d 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -6,7 +6,7 @@ #include "gtest/gtest.h" #include "axom/core/numerics/quadrature.hpp" -#include "axom/core/utilities/Utilities.hpp" +#include "axom/core/utilities/Utilities.hpp" TEST(numerics_quadrature, gauss_legendre_math_check) { @@ -15,7 +15,7 @@ TEST(numerics_quadrature, gauss_legendre_math_check) 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) + for(int npts = 1; npts <= N; ++npts) { // Evaluate using the quadrature rule auto rule = axom::numerics::get_gauss_legendre(npts); diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 312aa70d2f..0bfc1b2eec 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -12,13 +12,7 @@ #include "axom/primal/geometry/Vector.hpp" #include "axom/primal/geometry/BezierCurve.hpp" -// Use 1D Gauss-Legendre quadrature generated by MFEM if available, -// or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#else - #include "axom/core/numerics/quadrature.hpp" -#endif +#include "axom/core/numerics/quadrature.hpp" // C++ includes #include @@ -30,34 +24,10 @@ namespace primal namespace detail { -class QuadratureRuleView -{ -#ifdef AXOM_USE_MFEM -public: - QuadratureRuleView(const mfem::IntegrationRule& rule) : m_quad_rule(rule) { } - - double node(int idx) const { return m_quad_rule.IntPoint(idx).x; } - double weight(int idx) const { return m_quad_rule.IntPoint(idx).weight; } - -private: - const mfem::IntegrationRule& m_quad_rule; -#else -public: - QuadratureRuleView(const axom::numerics::QuadratureRule& rule) : m_quad_rule(rule) { } - - double node(int idx) const { return m_quad_rule.node(idx); } - double weight(int idx) const { return m_quad_rule.weight(idx); } - -private: - const axom::numerics::QuadratureRule& m_quad_rule; -#endif -}; - /*! * \brief Evaluate a scalar field line integral on a single Bezier curve. * - * Evaluate the scalar field line integral with MFEM integration rule if available, - * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * 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. @@ -70,13 +40,7 @@ inline double evaluate_scalar_line_integral_component(const primal::BezierCurve< Lambda&& scalar_integrand, const int npts) { -#ifdef AXOM_USE_MFEM - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - QuadratureRuleView quad(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); -#else - SLIC_WARNING_IF(npts > 100, "Axom's Gauss-Legendre nodes may be less precise for n > 100"); - QuadratureRuleView quad(axom::numerics::get_gauss_legendre(npts)); -#endif + const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); // Store/compute quadrature result double full_quadrature = 0.0; @@ -96,8 +60,7 @@ 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 if available, - * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * 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. @@ -110,13 +73,7 @@ inline double evaluate_vector_line_integral_component(const primal::BezierCurve< Lambda&& vec_field, const int npts) { -#ifdef AXOM_USE_MFEM - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - QuadratureRuleView quad(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); -#else - SLIC_WARNING_IF(npts > 100, "Axom's Gauss-Legendre nodes may be less precise for n > 100"); - QuadratureRuleView quad(axom::numerics::get_gauss_legendre(npts)); -#endif + const axom::numerics::QuadratureRule& quad = axom::numerics::get_gauss_legendre(npts); // Store/compute quadrature result double full_quadrature = 0.0; @@ -157,16 +114,8 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, const int npts_Q, const int npts_P) { -#ifdef AXOM_USE_MFEM - static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); - QuadratureRuleView quad_Q(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); - QuadratureRuleView quad_P(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1)); -#else - SLIC_WARNING_IF(npts_Q > 100 || npts_P > 100, - "Axom's Gauss-Legendre nodes may be less precise for n > 100"); - QuadratureRuleView quad_Q(axom::numerics::get_gauss_legendre(npts_Q)); - QuadratureRuleView quad_P(axom::numerics::get_gauss_legendre(npts_P)); -#endif + 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; diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index fc75f31376..4e3a1b8ece 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -9,9 +9,8 @@ * \brief Consists of methods that evaluate integrals over regions defined * by Bezier and NURBS curves, such as 2D area integrals and scalar/vector field line integrals * - * Line integrals are computed with 1D quadrature rules supplied by MFEM if available, - * or by methods in "core/numerics/quadrature" otherwise, which are less robust for - * high numbers of quadrature points (> ~100) + * 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. */ @@ -42,8 +41,7 @@ namespace primal * The line integral is evaluated on each curve in the CurvedPolygon, and added * together to represent the total integral. The Polygon need not be connected. * - * Uses 1D Gauss-Legendre quadrature generated by MFEM if available, - * or core/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * Evaluate the scalar field line integral with Gauss-Legendre quadrature * * \param [in] cpoly the CurvedPolygon object * \param [in] scalar_integrand the lambda function representing the integrand. @@ -149,8 +147,7 @@ double evaluate_scalar_line_integral(const axom::Array ~100) + * Evaluate the vector field line integral with Gauss-Legendre quadrature * * \param [in] cpoly the CurvedPolygon object * \param [in] vector_integrand the lambda function representing the integrand. diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 7548d74d35..f4fc34d60a 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -445,6 +445,30 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) } } +#ifdef AXOM_USE_MFEM +TEST(primal_integral, check_axom_mfem_quadrature_values) +{ + const int N = 200; + + for(int npts = 1; npts <= N; ++npts) + { + // Generate the Axom quadrature rule + axom::numerics::QuadratureRule axom_rule = axom::numerics::compute_gauss_legendre(npts); + + // Generate the MFEM quadrature rule + static mfem::IntegrationRules my_IntRules(0, mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& mfem_rule = my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + // Check that the nodes and weights are the same between the two rules + for(int j = 0; j < npts; ++j) + { + EXPECT_NEAR(axom_rule.node(j), mfem_rule.IntPoint(j).x, 1e-16); + EXPECT_NEAR(axom_rule.weight(j), mfem_rule.IntPoint(j).weight, 1e-16); + } + } +} +#endif + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); From 6c546542d77fb55bea55e21ca72fd29cc5031f24 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 12:03:46 -0700 Subject: [PATCH 15/28] Swap to axom::FlatMap --- src/axom/core/numerics/quadrature.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index a5740b0212..1da368a67a 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -4,6 +4,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include "axom/core/Array.hpp" +#include "axom/core/FlatMap.hpp" #include "axom/core/numerics/quadrature.hpp" // For math constants and includes @@ -136,7 +137,7 @@ QuadratureRule compute_gauss_legendre(int npts) const QuadratureRule& get_gauss_legendre(int npts) { // Define a static map that stores the GL quadrature rule for a given order - static std::map rule_library; + static axom::FlatMap rule_library; if(rule_library.find(npts) == rule_library.end()) { rule_library[npts] = compute_gauss_legendre(npts); From 03b9ebee073a911f1310733dd76975b9dca90346 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 17:43:34 -0700 Subject: [PATCH 16/28] Use new method in LinearizeCurves --- src/axom/quest/LinearizeCurves.cpp | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/axom/quest/LinearizeCurves.cpp b/src/axom/quest/LinearizeCurves.cpp index a471887314..9b576872b1 100644 --- a/src/axom/quest/LinearizeCurves.cpp +++ b/src/axom/quest/LinearizeCurves.cpp @@ -514,13 +514,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 +540,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 +564,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); From 6474280e8b2501e32618029cb70f6f16d0415eb7 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 6 Oct 2025 17:45:33 -0700 Subject: [PATCH 17/28] Remove unnecessary member data --- src/axom/core/numerics/quadrature.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index b07224c52f..11a6ba682e 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -39,15 +39,14 @@ class QuadratureRule double weight(size_t idx) const { return m_weights[idx]; }; //! \brief Accessofr for the size of the quadrature rule - int getNumPoints() const { return m_npts; } + int getNumPoints() const { return m_nodes.size(); } private: //! \brief Private constructor for use in compute_() methods. Only allocates memory - QuadratureRule(int npts) : m_nodes(npts), m_weights(npts), m_npts(npts) { }; + QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; axom::Array m_nodes; axom::Array m_weights; - int m_npts; }; /*! From 85f995e1c59c8371560da994b22c2c1a866a8f5b Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 11:31:16 -0700 Subject: [PATCH 18/28] Fix comments, asserts, numeric limits --- src/axom/core/numerics/quadrature.cpp | 11 +++++++++-- src/axom/primal/geometry/NURBSPatch.hpp | 6 ++---- src/axom/primal/tests/primal_integral.cpp | 8 ++++---- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 1da368a67a..8691fc4f0b 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -5,6 +5,7 @@ #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 @@ -33,6 +34,8 @@ namespace numerics */ QuadratureRule compute_gauss_legendre(int npts) { + SLIC_ASSERT(npts >= 1, "Quadrature rules must have >= 1 point"); + QuadratureRule rule(npts); if(npts == 1) @@ -68,7 +71,7 @@ QuadratureRule compute_gauss_legendre(int npts) for(int i = 1; i <= m; ++i) { // Each node is the root of a Legendre polynomial, - // which are approximately uniformly distirbuted in arccos(xi). + // 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; @@ -99,7 +102,7 @@ QuadratureRule compute_gauss_legendre(int npts) // Compute the Newton method step size dz = P_n / Pp_n; - if(std::fabs(dz) < 1e-16) + if(std::fabs(dz) < axom::numeric_limits::epsilon()) { done = true; @@ -132,10 +135,14 @@ QuadratureRule compute_gauss_legendre(int npts) * \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::Array's of nodes and weights */ const QuadratureRule& get_gauss_legendre(int npts) { + SLIC_ASSERT(npts >= 1, "Quadrature rules must have >= 1 point"); + // Define a static map that stores the GL quadrature rule for a given order static axom::FlatMap rule_library; if(rule_library.find(npts) == rule_library.end()) diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index f214c397d5..7c8949ece9 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -2488,8 +2488,7 @@ class NURBSPatch * then computes the 2D area of that projection to get the corresponding * component of the average surface normal. * - * Evaluates line integrals with MFEM integration rules if available, - * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * Evaluates the integral with Gauss-Legendre quadrature on each boundary curve * * \return The calculated mean surface normal */ @@ -3099,8 +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 * - * Evaluates line integrals with MFEM integration rules if available, - * or axom/numerics/quadrature otherwise, which are less robust for large npts (> ~100) + * Evaluates the integral with Gauss-Legendre quadrature on each boundary curve * * \return The calculated mean surface normal */ diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index f4fc34d60a..30bf1a49ef 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -6,8 +6,6 @@ #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" @@ -462,8 +460,10 @@ 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, 1e-16); - EXPECT_NEAR(axom_rule.weight(j), mfem_rule.IntPoint(j).weight, 1e-16); + 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()); } } } From 55ad9f47b679f5774769a5591a5884326005cd8d Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 12:46:47 -0700 Subject: [PATCH 19/28] Make an ArrayView version --- src/axom/core/numerics/quadrature.cpp | 68 +++++++++++---------- src/axom/core/numerics/quadrature.hpp | 23 ++++--- src/axom/core/tests/numerics_quadrature.hpp | 11 ++-- src/axom/primal/tests/primal_integral.cpp | 6 +- 4 files changed, 55 insertions(+), 53 deletions(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 8691fc4f0b..e04e44ff6c 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -3,6 +3,7 @@ // // 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" @@ -10,6 +11,7 @@ // For math constants and includes #include "axom/config.hpp" + #include namespace axom @@ -21,6 +23,8 @@ 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 @@ -29,39 +33,38 @@ namespace numerics * * \note This method constructs the points by scratch each time, without caching * \sa get_gauss_legendre(int) - * - * \return The `QuadratureRule` object which contains axom::Array's of nodes and weights */ -QuadratureRule compute_gauss_legendre(int npts) +void compute_gauss_legendre_data(int npts, axom::Array& nodes, axom::Array& weights) { - SLIC_ASSERT(npts >= 1, "Quadrature rules must have >= 1 point"); + assert("Quadrature rules must have >= 1 point" && (npts >= 1)); - QuadratureRule rule(npts); + nodes.resize(npts); + weights.resize(npts); if(npts == 1) { - rule.m_nodes[0] = 0.5; - rule.m_weights[0] = 1.0; - return rule; + nodes[0] = 0.5; + weights[0] = 1.0; + return; } if(npts == 2) { - rule.m_nodes[0] = 0.21132486540518711775; - rule.m_nodes[1] = 0.78867513459481288225; + nodes[0] = 0.21132486540518711775; + nodes[1] = 0.78867513459481288225; - rule.m_weights[0] = rule.m_weights[1] = 0.5; - return rule; + weights[0] = weights[1] = 0.5; + return; } if(npts == 3) { - rule.m_nodes[0] = 0.11270166537925831148207345; - rule.m_nodes[1] = 0.5; - rule.m_nodes[2] = 0.88729833462074168851792655; - - rule.m_weights[0] = 0.2777777777777777777777778; - rule.m_weights[1] = 0.4444444444444444444444444; - rule.m_weights[2] = 0.2777777777777777777777778; - return rule; + 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; @@ -114,14 +117,12 @@ QuadratureRule compute_gauss_legendre(int npts) z -= dz; } - rule.m_nodes[i - 1] = xi; - rule.m_nodes[n - i] = 1.0 - xi; + nodes[i - 1] = xi; + nodes[n - i] = 1.0 - xi; // Evaluate the weight as w_i = 2 / (1 - z^2) / P'_n(z)^2, with z \in [-1, 1] - rule.m_weights[i - 1] = rule.m_weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); + weights[i - 1] = weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n); } - - return rule; } /*! @@ -137,20 +138,25 @@ QuadratureRule compute_gauss_legendre(int npts) * * \warning The use of a static variable to store cached nodes makes this method not threadsafe. * - * \return The `QuadratureRule` object which contains axom::Array's of nodes and weights + * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights */ -const QuadratureRule& get_gauss_legendre(int npts) +QuadratureRule get_gauss_legendre(int npts) { - SLIC_ASSERT(npts >= 1, "Quadrature rules must have >= 1 point"); + assert("Quadrature rules must have >= 1 point" && (npts >= 1)); // Define a static map that stores the GL quadrature rule for a given order - static axom::FlatMap rule_library; + static axom::FlatMap, axom::Array>> rule_library; if(rule_library.find(npts) == rule_library.end()) { - rule_library[npts] = compute_gauss_legendre(npts); + rule_library[npts] = std::make_pair(axom::Array(npts), axom::Array(npts)); + compute_gauss_legendre_data(npts, rule_library[npts].first, rule_library[npts].second); } - return rule_library[npts]; + QuadratureRule rule; + rule.m_nodes = rule_library[npts].first.view(); + rule.m_weights = rule_library[npts].second.view(); + + return rule; } } /* end namespace numerics */ diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 11a6ba682e..e466ea123c 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -26,8 +26,8 @@ namespace numerics */ class QuadratureRule { - // Define friend functions so rules can only be created via compute_rule() methods - friend QuadratureRule compute_gauss_legendre(int npts); + // Define friend functions so rules can only be created via get_rule() methods + friend QuadratureRule get_gauss_legendre(int npts); public: QuadratureRule() = default; @@ -42,17 +42,16 @@ class QuadratureRule int getNumPoints() const { return m_nodes.size(); } private: - //! \brief Private constructor for use in compute_() methods. Only allocates memory - QuadratureRule(int npts) : m_nodes(npts), m_weights(npts) { }; - - axom::Array m_nodes; - axom::Array m_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 @@ -61,10 +60,10 @@ class QuadratureRule * * \note This method constructs the points by scratch each time, without caching * \sa get_gauss_legendre(int) - * - * \return The `QuadratureRule` object which contains axom::Array's of nodes and weights */ -QuadratureRule compute_gauss_legendre(int npts); +void compute_gauss_legendre_data(int npts, + axom::Array& nodes, + axom::Array& weights); /*! * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points @@ -77,9 +76,9 @@ QuadratureRule compute_gauss_legendre(int npts); * \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::Array's of nodes and weights + * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights */ -const QuadratureRule& get_gauss_legendre(int npts); +QuadratureRule get_gauss_legendre(int npts); } /* end namespace numerics */ } /* end namespace axom */ diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index d1e097e38d..c1d276806f 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -92,13 +92,10 @@ TEST(numerics_quadrature, gauss_legendre_math_check) TEST(numerics_quadrature, gauss_legendre_cache_check) { // The first two rules are put in the cache - const axom::numerics::QuadratureRule& first_rule = axom::numerics::get_gauss_legendre(20); - const axom::numerics::QuadratureRule& second_rule = axom::numerics::get_gauss_legendre(20); - - // The third is not - const axom::numerics::QuadratureRule& third_rule = axom::numerics::compute_gauss_legendre(20); + axom::numerics::QuadratureRule first_rule = axom::numerics::get_gauss_legendre(20); + axom::numerics::QuadratureRule second_rule = axom::numerics::get_gauss_legendre(20); // Check that the two rules are located in the same place in memory, and the third isn't - EXPECT_EQ(&first_rule, &second_rule); - EXPECT_NE(&first_rule, &third_rule); + // EXPECT_EQ(&first_rule.node(0), &second_rule.node(0)); + // EXPECT_EQ(&first_rule.weight(0), &second_rule.weight(0)); } diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 30bf1a49ef..f3a8986097 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -446,12 +446,12 @@ TEST(primal_integral, evaluate_nurbs_surface_normal) #ifdef AXOM_USE_MFEM TEST(primal_integral, check_axom_mfem_quadrature_values) { - const int N = 200; + const int N = 3; - for(int npts = 1; npts <= N; ++npts) + for(int npts = N; npts <= N; ++npts) { // Generate the Axom quadrature rule - axom::numerics::QuadratureRule axom_rule = axom::numerics::compute_gauss_legendre(npts); + 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); From 095f65835bfc212f959b6607eaaa16c9657b9d17 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 15:54:36 -0700 Subject: [PATCH 20/28] Remove unnecessary mfem include --- src/axom/primal/operators/detail/winding_number_2d_impl.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp index d77f7d8aae..e79d9f84be 100644 --- a/src/axom/primal/operators/detail/winding_number_2d_impl.hpp +++ b/src/axom/primal/operators/detail/winding_number_2d_impl.hpp @@ -21,11 +21,6 @@ // C++ includes #include -// MFEM includes -#ifdef AXOM_USE_MFEM - #include "mfem.hpp" -#endif - namespace axom { namespace primal From 5396af20d1b1e06fadc45aeb6d84c2f1fbabcaac Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 7 Oct 2025 17:05:12 -0700 Subject: [PATCH 21/28] Finish removing mfem from LinearizeCurves.cpp --- src/axom/quest/LinearizeCurves.cpp | 30 ------------------------------ 1 file changed, 30 deletions(-) diff --git a/src/axom/quest/LinearizeCurves.cpp b/src/axom/quest/LinearizeCurves.cpp index 9b576872b1..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; From cf066707d740083d816e862d3dfa82ef516af352 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 11:38:32 -0700 Subject: [PATCH 22/28] Update to use views and allocator IDs --- src/axom/core/numerics/quadrature.cpp | 30 ++++++++------- src/axom/core/numerics/quadrature.hpp | 12 ++++-- src/axom/core/tests/numerics_quadrature.hpp | 42 +++++++++++++++++---- src/axom/primal/tests/primal_integral.cpp | 5 +++ 4 files changed, 65 insertions(+), 24 deletions(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index e04e44ff6c..49aadda0b6 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -34,12 +34,15 @@ namespace numerics * \note This method constructs the points by scratch each time, without caching * \sa get_gauss_legendre(int) */ -void compute_gauss_legendre_data(int npts, axom::Array& nodes, axom::Array& weights) +void compute_gauss_legendre_data(int npts, + axom::Array& nodes, + axom::Array& weights, + int allocatorID) { assert("Quadrature rules must have >= 1 point" && (npts >= 1)); - nodes.resize(npts); - weights.resize(npts); + nodes = axom::Array(npts, npts, allocatorID); + weights = axom::Array(npts, npts, allocatorID); if(npts == 1) { @@ -140,23 +143,24 @@ void compute_gauss_legendre_data(int npts, axom::Array& nodes, axom::Arr * * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights */ -QuadratureRule get_gauss_legendre(int npts) +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 axom::FlatMap, axom::Array>> rule_library; - if(rule_library.find(npts) == rule_library.end()) + 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()) { - rule_library[npts] = std::make_pair(axom::Array(npts), axom::Array(npts)); - compute_gauss_legendre_data(npts, rule_library[npts].first, rule_library[npts].second); + auto& vals = rule_library[key]; + compute_gauss_legendre_data(npts, vals.first, vals.second); + value_it = rule_library.find(key); } - QuadratureRule rule; - rule.m_nodes = rule_library[npts].first.view(); - rule.m_weights = rule_library[npts].second.view(); - - return rule; + return QuadratureRule {value_it->second.first.view(), value_it->second.second.view()}; } } /* end namespace numerics */ diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index e466ea123c..3d0fb4bb6b 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -7,6 +7,7 @@ #define AXOM_NUMERICS_QUADRATURE_HPP_ #include "axom/core/Array.hpp" +#include "axom/core/memory_management.hpp" /*! * \file quadrature.hpp @@ -27,7 +28,7 @@ namespace numerics class QuadratureRule { // Define friend functions so rules can only be created via get_rule() methods - friend QuadratureRule get_gauss_legendre(int npts); + friend QuadratureRule get_gauss_legendre(int, int); public: QuadratureRule() = default; @@ -42,6 +43,10 @@ class QuadratureRule int getNumPoints() const { return m_nodes.size(); } private: + QuadratureRule(axom::ArrayView nodes, axom::ArrayView weights) + : m_nodes(nodes) + , m_weights(weights) { }; + axom::ArrayView m_nodes; axom::ArrayView m_weights; }; @@ -63,7 +68,8 @@ class QuadratureRule */ void compute_gauss_legendre_data(int npts, axom::Array& nodes, - axom::Array& weights); + axom::Array& weights, + int allocatorID = axom::getDefaultAllocatorID()); /*! * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points @@ -78,7 +84,7 @@ void compute_gauss_legendre_data(int npts, * * \return The `QuadratureRule` object which contains axom::ArrayView's of stored nodes and weights */ -QuadratureRule get_gauss_legendre(int npts); +QuadratureRule get_gauss_legendre(int npts, int allocatorID = axom::getDefaultAllocatorID()); } /* end namespace numerics */ } /* end namespace axom */ diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index c1d276806f..7a391117e6 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -19,6 +19,8 @@ TEST(numerics_quadrature, gauss_legendre_math_check) { // Evaluate using the quadrature rule auto rule = axom::numerics::get_gauss_legendre(npts); + std::cout << rule.getNumPoints() << std::endl; + for(int i = 0; i < npts; ++i) std::cout << rule.node(i) << std::endl; int degree = 2 * npts - 1; // Define a collection of random coefficients for a polynomial @@ -89,13 +91,37 @@ TEST(numerics_quadrature, gauss_legendre_math_check) } } -TEST(numerics_quadrature, gauss_legendre_cache_check) +template +struct test_device_quadrature { - // The first two rules are put in the cache - axom::numerics::QuadratureRule first_rule = axom::numerics::get_gauss_legendre(20); - axom::numerics::QuadratureRule second_rule = axom::numerics::get_gauss_legendre(20); + static void test() + { + const int npts = 15; - // Check that the two rules are located in the same place in memory, and the third isn't - // EXPECT_EQ(&first_rule.node(0), &second_rule.node(0)); - // EXPECT_EQ(&first_rule.weight(0), &second_rule.weight(0)); -} + // 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(M_PI * rule.node(i)); }); + + EXPECT_NEAR(quadrature_sum.get(), 0.636619772368, 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/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index f3a8986097..2bdb18c7ff 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -10,6 +10,11 @@ #include "gtest/gtest.h" +// MFEM includes +#ifdef AXOM_USE_MFEM + #include "mfem.hpp" +#endif + namespace primal = axom::primal; TEST(primal_integral, evaluate_area_integral) From 1caf0987465952037841490334ba260f24f1a77d Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 12:01:33 -0700 Subject: [PATCH 23/28] fix test --- src/axom/core/tests/numerics_quadrature.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 7a391117e6..06c8543d7d 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -5,6 +5,7 @@ #include "gtest/gtest.h" +#include "axom/config.hpp" #include "axom/core/numerics/quadrature.hpp" #include "axom/core/utilities/Utilities.hpp" @@ -106,9 +107,9 @@ struct test_device_quadrature axom::ReduceSum quadrature_sum(0.0); axom::for_all( npts, - AXOM_LAMBDA(axom::IndexType i) { quadrature_sum += rule.weight(i) * sin(M_PI * rule.node(i)); }); + AXOM_LAMBDA(axom::IndexType i) { quadrature_sum += rule.weight(i) * sin(rule.node(i)); }); - EXPECT_NEAR(quadrature_sum.get(), 0.636619772368, 1e-6); + EXPECT_NEAR(quadrature_sum.get(), 0.459697694132, 1e-6); } }; From aa242b168d352c8e01eebd140be2eb02a8e2bfa1 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 13:07:54 -0700 Subject: [PATCH 24/28] Correct prints, add decorators --- src/axom/core/numerics/quadrature.hpp | 8 +++++--- src/axom/core/tests/numerics_quadrature.hpp | 2 -- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/axom/core/numerics/quadrature.hpp b/src/axom/core/numerics/quadrature.hpp index 3d0fb4bb6b..2fe3230bea 100644 --- a/src/axom/core/numerics/quadrature.hpp +++ b/src/axom/core/numerics/quadrature.hpp @@ -31,18 +31,20 @@ class QuadratureRule friend QuadratureRule get_gauss_legendre(int, int); public: - QuadratureRule() = default; - //! \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 Accessofr for the size of the quadrature rule + //! \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) { }; diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 06c8543d7d..f0bc2e6ec1 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -20,8 +20,6 @@ TEST(numerics_quadrature, gauss_legendre_math_check) { // Evaluate using the quadrature rule auto rule = axom::numerics::get_gauss_legendre(npts); - std::cout << rule.getNumPoints() << std::endl; - for(int i = 0; i < npts; ++i) std::cout << rule.node(i) << std::endl; int degree = 2 * npts - 1; // Define a collection of random coefficients for a polynomial From 10eac223a61de452a831f8473c404a0d28c46c69 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 14:29:13 -0700 Subject: [PATCH 25/28] Push debugging change --- src/axom/core/tests/numerics_quadrature.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index f0bc2e6ec1..7f93d17b51 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -98,6 +98,7 @@ struct test_device_quadrature const int npts = 15; // Create the rule with the proper allocator + std::cout << axom::execution_space::allocatorID() << std::endl; const auto rule = axom::numerics::get_gauss_legendre(npts, axom::execution_space::allocatorID()); From cdda86976e0e4d9b618d60c2984c9938e2ccfb51 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 14:33:02 -0700 Subject: [PATCH 26/28] Fix a missing allocatorID pass --- src/axom/core/numerics/quadrature.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index 49aadda0b6..ec7f0ca59c 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -156,7 +156,7 @@ QuadratureRule get_gauss_legendre(int npts, int allocatorID) if(value_it == rule_library.end()) { auto& vals = rule_library[key]; - compute_gauss_legendre_data(npts, vals.first, vals.second); + compute_gauss_legendre_data(npts, vals.first, vals.second, allocatorID); value_it = rule_library.find(key); } From 95dab9dbe9b4757926eea009754e55891aa6ae8b Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 10 Oct 2025 16:03:25 -0700 Subject: [PATCH 27/28] Remove debugging statement --- src/axom/core/tests/numerics_quadrature.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/axom/core/tests/numerics_quadrature.hpp b/src/axom/core/tests/numerics_quadrature.hpp index 7f93d17b51..f0bc2e6ec1 100644 --- a/src/axom/core/tests/numerics_quadrature.hpp +++ b/src/axom/core/tests/numerics_quadrature.hpp @@ -98,7 +98,6 @@ struct test_device_quadrature const int npts = 15; // Create the rule with the proper allocator - std::cout << axom::execution_space::allocatorID() << std::endl; const auto rule = axom::numerics::get_gauss_legendre(npts, axom::execution_space::allocatorID()); From 69620c6904b8eec01d09bf1b8268ecb1af5e6107 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 23 Oct 2025 15:20:32 -0700 Subject: [PATCH 28/28] Improve comment quality --- src/axom/core/numerics/quadrature.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp index ec7f0ca59c..9a54819f84 100644 --- a/src/axom/core/numerics/quadrature.cpp +++ b/src/axom/core/numerics/quadrature.cpp @@ -123,7 +123,8 @@ void compute_gauss_legendre_data(int npts, nodes[i - 1] = xi; nodes[n - i] = 1.0 - xi; - // Evaluate the weight as w_i = 2 / (1 - z^2) / P'_n(z)^2, with z \in [-1, 1] + // 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); } }