Add Gauss-Legendre quadrature detached from MFEM include

RELEASE-NOTES.md

@@ -81,6 +81,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Adds `quest::MFEMReader` for reading 1D MFEM contours in 2D space.
 - Adds an option to `quest::SamplingShaper` to allow in/out tests based on winding numbers for MFEM contours.
 - The `shaping_driver` example program can select `--sampling inout` to do the default In/Out sampling and `--sampling windingnumber` to select winding number in/out tests for MFEM data.
+- Adds Axom-native Gauss-Legendre quadrature rules that can be used without an MFEM dependency
 
 ###  Changed
 - Updates blt submodule to [BLT version 0.7.1][https://github.com/LLNL/blt/releases/tag/v0.7.1]

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,167 @@
+// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "axom/core/utilities/Utilities.hpp"
+#include "axom/core/Array.hpp"
+#include "axom/core/FlatMap.hpp"
+#include "axom/core/NumericLimits.hpp"
+#include "axom/core/numerics/quadrature.hpp"
+
+// For math constants and includes
+#include "axom/config.hpp"
+
+#include <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
+ * \param [out] nodes The array of 1D nodes
+ * \param [out] weights The array of weights
+ * 
+ * A Gauss-Legendre rule with \a npts points can exactly integrate
+ *  polynomials of order 2 * npts - 1
+ *
+ * Algorithm adapted from the MFEM implementation in `mfem/fem/intrules.cpp`
+ * 
+ * \note This method constructs the points by scratch each time, without caching
+ * \sa get_gauss_legendre(int)
+ */
+void compute_gauss_legendre_data(int npts,
+                                 axom::Array<double>& nodes,
+                                 axom::Array<double>& weights,
+                                 int allocatorID)
+{
+  assert("Quadrature rules must have >= 1 point" && (npts >= 1));
+
+  nodes = axom::Array<double>(npts, npts, allocatorID);
+  weights = axom::Array<double>(npts, npts, allocatorID);
+
+  if(npts == 1)
+  {
+    nodes[0] = 0.5;
+    weights[0] = 1.0;
+    return;
+  }
+  if(npts == 2)
+  {
+    nodes[0] = 0.21132486540518711775;
+    nodes[1] = 0.78867513459481288225;
+
+    weights[0] = weights[1] = 0.5;
+    return;
+  }
+  if(npts == 3)
+  {
+    nodes[0] = 0.11270166537925831148207345;
+    nodes[1] = 0.5;
+    nodes[2] = 0.88729833462074168851792655;
+
+    weights[0] = 0.2777777777777777777777778;
+    weights[1] = 0.4444444444444444444444444;
+    weights[2] = 0.2777777777777777777777778;
+    return;
+  }
+
+  const int n = npts;
+  const int m = (npts + 1) / 2;
+
+  // Nodes are mirrored across x = 0.5, so only need to evaluate half
+  for(int i = 1; i <= m; ++i)
+  {
+    // Each node is the root of a Legendre polynomial,
+    //  which are approximately uniformly distributed in arccos(xi).
+    // This makes cos a good initial guess for subsequent Newton iterations
+    double z = std::cos(M_PI * (i - 0.25) / (n + 0.5));
+    double Pp_n, P_n, dz, xi = 0.0;
+
+    bool done = false;
+    while(true)
+    {
+      // Recursively evaluate P_n(z) through the recurrence relation
+      //  n * P_n(z) = (2n - 1) * P_{n-1}(z) - (n - 1) * P_{n - 2}(z)
+      double P_nm1 = 1.0;  // P_0(z) = 1
+      P_n = z;             // P_1(z) = z
+      for(int j = 2; j <= n; ++j)
+      {
+        double P_nm2 = P_nm1;
+        P_nm1 = P_n;
+        P_n = ((2 * j - 1) * z * P_nm1 - (j - 1) * P_nm2) / j;
+      }
+
+      // Evaluate P'_n(z) using another recurrence relation
+      //  (z^2 - 1) * P'_n(z) = n * z * P_n(z) - n * P_{n-1}(z)
+      Pp_n = n * (z * P_n - P_nm1) / (z * z - 1);
+
+      if(done)
+      {
+        break;
+      }
+
+      // Compute the Newton method step size
+      dz = P_n / Pp_n;
+
+      if(std::fabs(dz) < axom::numeric_limits<double>::epsilon())
+      {
+        done = true;
+
+        // Scale back to [0, 1]
+        xi = ((1 - z) + dz) / 2;
+      }
+
+      // Take the Newton step, repeat the process
+      z -= dz;
+    }
+
+    nodes[i - 1] = xi;
+    nodes[n - i] = 1.0 - xi;
+
+    // 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);
+  }
+}
+
+/*!
+ * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points 
+ *
+ * \param [in] npts The number of points in the rule
+ * 
+ * A Gauss-Legendre rule with \a npts points can exactly integrate
+ *  polynomials of order 2 * npts - 1
+ *
+ * \note If this method has already been called for a given order, it will reuse the same quadrature points
+ *  without needing to recompute them
+ *
+ * \warning The use of a static variable to store cached nodes makes this method not threadsafe.
+ * 
+ * \return The `QuadratureRule` object which contains axom::ArrayView<double>'s of stored nodes and weights
+ */
+QuadratureRule get_gauss_legendre(int npts, int allocatorID)
+{
+  assert("Quadrature rules must have >= 1 point" && (npts >= 1));
+
+  // Define a static map that stores the GL quadrature rule for a given order
+  static std::map<std::pair<int, int>, std::pair<axom::Array<double>, axom::Array<double>>> rule_library;
+
+  const std::pair<int, int> key = std::make_pair(npts, allocatorID);
+
+  auto value_it = rule_library.find(key);
+  if(value_it == rule_library.end())
+  {
+    auto& vals = rule_library[key];
+    compute_gauss_legendre_data(npts, vals.first, vals.second, allocatorID);
+    value_it = rule_library.find(key);
+  }
+
+  return QuadratureRule {value_it->second.first.view(), value_it->second.second.view()};
+}
+
+} /* end namespace numerics */
+} /* end namespace axom */

src/axom/core/numerics/quadrature.hpp

@@ -0,0 +1,94 @@
+// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#ifndef AXOM_NUMERICS_QUADRATURE_HPP_
+#define AXOM_NUMERICS_QUADRATURE_HPP_
+
+#include "axom/core/Array.hpp"
+#include "axom/core/memory_management.hpp"
+
+/*!
+ * \file quadrature.hpp
+ * The functions declared in this header file find the nodes and weights of 
+ *   arbitrary order quadrature rules
+ */
+
+namespace axom
+{
+namespace numerics
+{
+
+/*!
+ * \class QuadratureRule
+ *
+ * \brief Stores a fixed array of 1D quadrature points and weights
+ */
+class QuadratureRule
+{
+  // Define friend functions so rules can only be created via get_rule() methods
+  friend QuadratureRule get_gauss_legendre(int, int);
+
+public:
+  //! \brief Accessor for quadrature nodes
+  AXOM_HOST_DEVICE
+  double node(size_t idx) const { return m_nodes[idx]; };
+
+  //! \brief Accessor for quadrature weights
+  AXOM_HOST_DEVICE
+  double weight(size_t idx) const { return m_weights[idx]; };
+
+  //! \brief Accessor for the size of the quadrature rule
+  AXOM_HOST_DEVICE
+  int getNumPoints() const { return m_nodes.size(); }
+
+private:
+  //! \brief Use a private constructor to avoid creation of an invalid rule
+  QuadratureRule(axom::ArrayView<double> nodes, axom::ArrayView<double> weights)
+    : m_nodes(nodes)
+    , m_weights(weights) { };
+
+  axom::ArrayView<double> m_nodes;
+  axom::ArrayView<double> m_weights;
+};
+
+/*!
+ * \brief Computes a 1D quadrature rule of Gauss-Legendre points 
+ *
+ * \param [in] npts The number of points in the rule
+ * \param [out] nodes The array of 1D nodes
+ * \param [out] weights The array of weights
+ * 
+ * A Gauss-Legendre rule with \a npts points can exactly integrate
+ *  polynomials of order 2 * npts - 1
+ *
+ * Algorithm adapted from the MFEM implementation in `mfem/fem/intrules.cpp`
+ * 
+ * \note This method constructs the points by scratch each time, without caching
+ * \sa get_gauss_legendre(int)
+ */
+void compute_gauss_legendre_data(int npts,
+                                 axom::Array<double>& nodes,
+                                 axom::Array<double>& weights,
+                                 int allocatorID = axom::getDefaultAllocatorID());
+
+/*!
+ * \brief Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points 
+ *
+ * \param [in] npts The number of points in the rule
+ * 
+ * A Gauss-Legendre rule with \a npts points can exactly integrate
+ *  polynomials of order 2 * npts - 1
+ *
+ * \note If this method has already been called for a given order, it will reuse the same quadrature points
+ *  without needing to recompute them
+ *
+ * \return The `QuadratureRule` object which contains axom::ArrayView<double>'s of stored nodes and weights
+ */
+QuadratureRule get_gauss_legendre(int npts, int allocatorID = axom::getDefaultAllocatorID());
+
+} /* end namespace numerics */
+} /* end namespace axom */
+
+#endif  // AXOM_NUMERICS_QUADRATURE_HPP_

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"