Expose gauss jacobi quadrature with alpha parameter

basix/_basixcpp.pyi

@@ -2,6 +2,7 @@ from collections.abc import Sequence
 import enum
 from typing import Annotated, overload
 
+from numpy import float64
 from numpy.typing import ArrayLike
 
 
@@ -493,6 +494,8 @@ def index(arg0: int, arg1: int, arg2: int, /) -> int: ...
 
 def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: PolysetType, arg3: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64')], Annotated[ArrayLike, dict(dtype='float64')]]: ...
 
+def gauss_jacobi_rule(arg0: float64, arg1: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64')], Annotated[ArrayLike, dict(dtype='float64')]]: ...
+
 def polynomials_dim(arg0: PolynomialType, arg1: CellType, arg2: int, /) -> int: ...
 
 def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType, /) -> PolysetType: ...

cpp/basix/quadrature.cpp

@@ -4917,6 +4917,19 @@ quadrature::type quadrature::get_default_rule(cell::type celltype, int m)
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
+std::array<std::vector<T>, 2> quadrature::gauss_jacobi_rule(T a, int m)
+{
+  auto [pts, wts] = compute_gauss_jacobi_rule(a, m);
+  for (std::size_t i = 0; i < wts.size(); ++i)
+  {
+    pts[i] += 1.0;
+    pts[i] /= 2.0;
+    wts[i] /= std::pow(2.0, a + 1);
+  }
+  return {std::move(pts), std::move(wts)};
+}
+//-----------------------------------------------------------------------------
+template <std::floating_point T>
 std::array<std::vector<T>, 2>
 quadrature::make_quadrature(quadrature::type rule, cell::type celltype,
                             polyset::type polytype, int m)
@@ -4966,6 +4979,7 @@ std::vector<T> quadrature::get_gll_points(int m)
   return make_gll_line<T>(m)[0];
 }
 //-----------------------------------------------------------------------------
+/// @cond DOXYGEN_SHOULD_SKIP_THIS
 template std::array<std::vector<float>, 2>
 quadrature::make_quadrature(quadrature::type, cell::type, polyset::type, int);
 template std::array<std::vector<double>, 2>
@@ -4976,4 +4990,10 @@ template std::vector<double> quadrature::get_gl_points(int);
 
 template std::vector<float> quadrature::get_gll_points(int);
 template std::vector<double> quadrature::get_gll_points(int);
+
+template std::array<std::vector<float>, 2> quadrature::gauss_jacobi_rule(float,
+                                                                         int);
+template std::array<std::vector<double>, 2>
+quadrature::gauss_jacobi_rule(double, int);
+/// @endcond
 //-----------------------------------------------------------------------------

cpp/basix/quadrature.h

@@ -26,7 +26,17 @@ enum class type
   strang_fix = 22,
 };
 
+/// @brief Get the Gauss-Jacobi rule for the interval for integrating
+/// f(x) * (1-x)^a on the interval [0, 1].
+/// @tparam The floating point type.
+/// @param[in] a The exponent a.
+/// @param[in] m The number of points.
+/// @return Points and weights of a Gauss-Jacobi rule.
+template <std::floating_point T>
+std::array<std::vector<T>, 2> gauss_jacobi_rule(T a, int m);
+
 /// @brief Make a quadrature rule on a reference cell.
+/// @tparam T The floating point type.
 /// @param[in] rule Type of quadrature rule (or use quadrature::Default).
 /// @param[in] celltype Cell type.
 /// @param[in] polytype Polyset type.
@@ -48,12 +58,14 @@ quadrature::type get_default_rule(cell::type celltype, int m);
 
 /// @brief Get Gauss-Lobatto-Legendre (GLL) points on the interval [0,
 /// 1].
+/// @tparam T The floating point type.
 /// @param[in] m Number of points.
 /// @return Array of GLL points.
 template <std::floating_point T>
 std::vector<T> get_gll_points(int m);
 
 /// @brief Get Gauss-Legendre (GL) points on the interval [0, 1].
+/// @tparam T The floating point type.
 /// @param[in] m Number of points
 /// @return Array of GL points.
 template <std::floating_point T>

doc/cpp/Doxyfile

@@ -2057,6 +2057,7 @@ INCLUDE_FILE_PATTERNS  =
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
 PREDEFINED             = protected=private
+PREDEFINED             = DOXYGEN_SHOULD_SKIP_THIS
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
 # tag can be used to specify a list of macro names that should be expanded. The

python/basix/_basixcpp.pyi

@@ -3,6 +3,7 @@ import enum
 from typing import Annotated, overload
 
 from numpy.typing import ArrayLike
+from numpy import float64
 
 
 class CellType(enum.IntEnum):
@@ -493,6 +494,8 @@ def index(arg0: int, arg1: int, arg2: int, /) -> int: ...
 
 def make_quadrature(arg0: QuadratureType, arg1: CellType, arg2: PolysetType, arg3: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64')], Annotated[ArrayLike, dict(dtype='float64')]]: ...
 
+def gauss_jacobi_rule(arg0: float64, arg1: int, /) -> tuple[Annotated[ArrayLike, dict(dtype='float64')], Annotated[ArrayLike, dict(dtype='float64')]]: ...
+
 def polynomials_dim(arg0: PolynomialType, arg1: CellType, arg2: int, /) -> int: ...
 
 def restriction(arg0: PolysetType, arg1: CellType, arg2: CellType, /) -> PolysetType: ...

python/basix/quadrature.py

@@ -5,10 +5,12 @@
 # SPDX-License-Identifier:    MIT
 """Functions to manipulate quadrature types."""
 
+import numpy as _np
 import numpy.typing as _npt
 
 from basix._basixcpp import QuadratureType
 from basix._basixcpp import make_quadrature as _mq
+from basix._basixcpp import gauss_jacobi_rule as _gjr
 from basix.cell import CellType
 from basix.polynomials import PolysetType
 
@@ -56,3 +58,20 @@ def make_quadrature(
         Quadrature points and weights.
     """
     return _mq(rule, cell, polyset_type, degree)
+
+
+def gauss_jacobi_rule(
+    alpha: _np.floating,
+    npoints: int,
+) -> tuple[_npt.ArrayLike, _npt.ArrayLike]:
+    """Create a Gauss-Jacobi quadrature rule for integrating f(x)*(1-x)**alpha
+    on the interval [0, 1].
+
+    Args:
+        alpha: The exponent alpha
+        npoints: Number of points
+
+    Returns:
+        Quadrature points and weights.
+    """
+    return _gjr(_np.float64(alpha), npoints)

python/wrapper.cpp

@@ -695,6 +695,16 @@ NB_MODULE(_basixcpp, m)
                          as_nbarray(std::move(w)));
       });
 
+  m.def(
+      "gauss_jacobi_rule",
+      [](double a, int m)
+      {
+        auto [pts, w]
+            = quadrature::gauss_jacobi_rule<double>(a, m);
+        return std::pair(as_nbarray(std::move(pts)),
+                         as_nbarray(std::move(w)));
+      });
+
   m.def("index", nb::overload_cast<int>(&basix::indexing::idx));
   m.def("index", nb::overload_cast<int, int>(&basix::indexing::idx));
   m.def("index", nb::overload_cast<int, int, int>(&basix::indexing::idx));

test/test_quadrature.py

@@ -137,8 +137,6 @@ def test_qorder_pyramid(m, scheme):
     polyset += [x**a * y**b / (1 - z) ** min(a, b) for a in range(m) for b in range(m + 1 - a, m)]
     for f in polyset:
         q = sympy.integrate(f, (x, 0, 1 - z), (y, 0, 1 - z), (z, 0, 1))
-        print(Qpts)
-        print(Qwts)
         s = 0.0
         for pt, wt in zip(Qpts, Qwts):
             s += wt * f.subs([(x, pt[0]), (y, pt[1]), (z, pt[2])])
@@ -195,3 +193,13 @@ def test_gll():
     assert np.allclose(wts, ref_wts3)
     assert np.isclose((pts * wts.reshape(-1, 1)).sum(), 0)
     assert np.isclose(sum(wts), 8)
+
+
+@pytest.mark.parametrize("alpha", [0, 0.0, 1, 1.0, 1.5, 2.0, 3.0, 4.2])
+@pytest.mark.parametrize("degree", range(6))
+def test_gauss_jacobi_rule(alpha, degree):
+    pts, wts = basix.quadrature.gauss_jacobi_rule(alpha, degree + 1)
+    integral = sum(w * (1 - p) ** degree for p, w in zip(pts, wts))
+
+    expected = 1 / (alpha + degree + 1)
+    assert np.isclose(integral, expected)