Add Gauss-Legendre quadrature detached from MFEM include

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]

src/axom/core/CMakeLists.txt

@@ -42,6 +42,7 @@ set(core_headers
     numerics/eigen_solve.hpp
     numerics/eigen_sort.hpp
     numerics/floating_point_limits.hpp
+    numerics/quadrature.hpp
     numerics/jacobi_eigensolve.hpp
     numerics/linear_solve.hpp
     numerics/matvecops.hpp
@@ -106,6 +107,7 @@ set(core_sources
     ${PROJECT_BINARY_DIR}/axom/core/utilities/About.cpp
 
     numerics/polynomial_solvers.cpp
+    numerics/quadrature.cpp
 
     Path.cpp
     Types.cpp

src/axom/core/numerics/quadrature.cpp

@@ -0,0 +1,150 @@
+// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "axom/core/Array.hpp"
+#include "axom/core/FlatMap.hpp"
+#include "axom/core/numerics/quadrature.hpp"
+
+// For math constants and includes
+#include "axom/config.hpp"
+#include <cmath>
+
+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<double>'s of nodes and weights
+ */
+QuadratureRule compute_gauss_legendre(int npts)
+{
+  QuadratureRule rule(npts);
+
+  if(npts == 1)
+  {
+    rule.m_nodes[0] = 0.5;
+    rule.m_weights[0] = 1.0;
+    return rule;
+  }
+  if(npts == 2)
+  {
+    rule.m_nodes[0] = 0.21132486540518711775;
+    rule.m_nodes[1] = 0.78867513459481288225;
+
+    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;
+
+  // 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_n, P_n, dz, xi = 0.0;
+
+    bool done = false;
+    while(true)
+    {
+      // Recursively evaluate P_n(z) through the recurrence relation
+      //  n * P_n(z) = (2n - 1) * P_{n-1}(z) - (n - 1) * P_{n - 2}(z)
+      double P_nm1 = 1.0;  // P_0(z) = 1
+      P_n = z;             // P_1(z) = z
+      for(int j = 2; j <= n; ++j)
+      {
+        double P_nm2 = P_nm1;
+        P_nm1 = P_n;
+        P_n = ((2 * j - 1) * z * P_nm1 - (j - 1) * P_nm2) / j;
+      }
+
+      // Evaluate P'_n(z) using another recurrence relation
+      //  (z^2 - 1) * P'_n(z) = n * z * P_n(z) - n * P_{n-1}(z)
+      Pp_n = n * (z * P_n - P_nm1) / (z * z - 1);
+
+      if(done)
+      {
+        break;
+      }
+
+      // Compute the Newton method step size
+      dz = P_n / Pp_n;
+
+      if(std::fabs(dz) < 1e-16)
+      {
+        done = true;
+
+        // Scale back to [0, 1]
+        xi = ((1 - z) + dz) / 2;
+      }
+
+      // Take the Newton step, repeat the process
+      z -= dz;
+    }
+
+    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]
+    rule.m_weights[i - 1] = rule.m_weights[n - i] = 1.0 / (4.0 * xi * (1.0 - xi) * Pp_n * Pp_n);
+  }
+
+  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<double>'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
+  static axom::FlatMap<int, QuadratureRule> 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 */

src/axom/core/numerics/quadrature.hpp

@@ -0,0 +1,88 @@
+// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#ifndef AXOM_NUMERICS_QUADRATURE_HPP_
+#define AXOM_NUMERICS_QUADRATURE_HPP_
+
+#include "axom/core/Array.hpp"
+
+/*!
+ * \file quadrature.hpp
+ * The functions declared in this header file find the nodes and weights of 
+ *   arbitrary order quadrature rules
+ */
+
+namespace axom
+{
+namespace numerics
+{
+
+/*!
+ * \class QuadratureRule
+ *
+ * \brief Stores a fixed array of 1D quadrature points and weights
+ */
+class QuadratureRule
+{
+  // Define friend functions so rules can only be created via compute_rule() methods
+  friend QuadratureRule compute_gauss_legendre(int npts);
+
+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]; };
+
+  //! \brief Accessofr for the size of the quadrature rule
+  int getNumPoints() const { return m_npts; }
+
+private:
+  //! \brief Private constructor for use in compute_<rule>() methods. Only allocates memory
+  QuadratureRule(int npts) : m_nodes(npts), m_weights(npts), m_npts(npts) { };
+
+  axom::Array<double> m_nodes;
+  axom::Array<double> m_weights;
+  int m_npts;
+};
+
+/*!
+ * \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<double>'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<double>'s of nodes and weights
+ */
+const QuadratureRule& get_gauss_legendre(int npts);
+
+} /* end namespace numerics */
+} /* end namespace axom */
+
+#endif  // AXOM_NUMERICS_QUADRATURE_HPP_

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

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"

src/axom/core/tests/numerics_quadrature.hpp

@@ -0,0 +1,104 @@
+// 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/quadrature.hpp"
+#include "axom/core/utilities/Utilities.hpp"
+
+TEST(numerics_quadrature, gauss_legendre_math_check)
+{
+  const int N = 200;
+
+  double coeffs[2 * N - 1];
+
+  // Test that the rules provide exact integration for polynomials of degree 2n - 1
+  for(int npts = 1; npts <= N; ++npts)
+  {
+    // Evaluate using the quadrature rule
+    auto rule = axom::numerics::get_gauss_legendre(npts);
+    int degree = 2 * npts - 1;
+
+    // Define a collection of random coefficients for a polynomial
+    for(int j = 0; j < degree; ++j)
+    {
+      // Seed the random coefficients for consistency in the test
+      coeffs[j] = axom::utilities::random_real(-1.0, 1.0, npts + j);
+    }
+
+    // Evaluate the area under the curve from 0 to 1 analytically
+    double analytic_result = 0.0;
+    for(int j = 0; j < degree; ++j)
+    {
+      analytic_result += coeffs[j] / (j + 1);
+    }
+
+    // Evaluate the polynomial using Horner's rule
+    auto eval_polynomial = [&degree, &coeffs](double x) {
+      double result = coeffs[degree - 1];
+      for(int i = degree - 2; i >= 0; --i)
+      {
+        result = result * x + coeffs[i];
+      }
+      return result;
+    };
+
+    double quadrature_result = 0.0;
+    for(int j = 0; j < npts; ++j)
+    {
+      quadrature_result += rule.weight(j) * eval_polynomial(rule.node(j));
+    }
+
+    EXPECT_NEAR(quadrature_result, analytic_result, 1e-10);
+
+    // Check that nodes aren't duplicated
+    for(int j = 1; j < npts; ++j)
+    {
+      EXPECT_GT(rule.node(j), rule.node(j - 1));
+    }
+
+    // Check that the sum of the weights is 1, and that all are positive
+    double weight_sum = 0.0;
+    for(int j = 0; j < npts; ++j)
+    {
+      EXPECT_GT(rule.weight(j), 0.0);
+      weight_sum += rule.weight(j);
+    }
+
+    EXPECT_NEAR(weight_sum, 1.0, 1e-10);
+
+    // Check that each node is the root of the next Legendre polynomial
+    for(int j = 0; j < npts; ++j)
+    {
+      // Rescale the node to [-1, 1]
+      double z = 2 * rule.node(j) - 1;
+
+      double P_n = z, P_nm1 = 1.0;
+      for(int k = 2; k <= npts; ++k)
+      {
+        double P_nm2 = P_nm1;
+        P_nm1 = P_n;
+        P_n = ((2 * k - 1) * z * P_nm1 - (k - 1) * P_nm2) / k;
+      }
+
+      // Note that this metric does not directly imply that |z - z*| < tol
+      EXPECT_NEAR(P_n, 0.0, 1e-10);
+    }
+  }
+}
+
+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);
+}