diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index 9b538aa..f07427e 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -55,6 +55,8 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { point_container points(npts); weight_container weights(npts); diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp index 83c8c17..338c6c3 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp @@ -59,6 +59,8 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { const weight_type oonpp = 1.0 / (npts + 1); diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index 39a8034..bb242d8 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -55,6 +55,8 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_2n_p_1 = M_PI / (2 * npts + 1); diff --git a/include/integratorxx/quadratures/primitive/gausslegendre.hpp b/include/integratorxx/quadratures/primitive/gausslegendre.hpp index 755c305..183fc9b 100644 --- a/include/integratorxx/quadratures/primitive/gausslegendre.hpp +++ b/include/integratorxx/quadratures/primitive/gausslegendre.hpp @@ -44,6 +44,8 @@ struct quadrature_traits< using point_container = std::vector< point_type >; using weight_container = std::vector< weight_type >; + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate( size_t npts ) { point_container points( npts ); diff --git a/include/integratorxx/quadratures/primitive/gausslobatto.hpp b/include/integratorxx/quadratures/primitive/gausslobatto.hpp index 830e83d..f662248 100644 --- a/include/integratorxx/quadratures/primitive/gausslobatto.hpp +++ b/include/integratorxx/quadratures/primitive/gausslobatto.hpp @@ -34,6 +34,8 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + inline static constexpr bool bound_inclusive = true; + inline static std::tuple generate( size_t npts) { point_container points(npts); diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index d106801..a7cc128 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -50,6 +50,23 @@ auto chebyshev_T(int n, T x) { + +// Compose one primitive rule with every radial transform. The pairing is only +// legal if the primitive's traits declare bound_inclusive, so instantiating +// these is itself the regression guard for #116. +template +void compose_with_all_transforms() { + constexpr size_t npts = 16; + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature>(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); +} + TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { // Reference integral for polynomial evaluated over [-1,1] @@ -191,6 +208,50 @@ TEST_CASE( "Becke Quadratures", "[1d-quad]" ) { test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); } + +TEST_CASE( "Radial transforms over the primitive rules", "[1d-quad]" ) { + + // Becke and Treutler-Ahlrichs are fixed aliases over GaussChebyshev2, and + // Mura-Knowles and MHL over UniformTrapezoid, so nothing above pairs a + // transform with any other primitive by hand. + + SECTION("Gauss-Lobatto x Becke") { + // Lobatto is the only bound_inclusive Gauss rule, so this is the + // endpoint-dropping path for those: the base rule is asked for npts+2 + // nodes and the two that map to r = 0 and r = infinity are skipped. + // UniformTrapezoid is bound_inclusive as well, so the Mura-Knowles and + // MHL cases above already cover that path. + IntegratorXX::RadialTransformQuadrature< + IntegratorXX::GaussLobatto, + IntegratorXX::BeckeRadialTraits> quad(350); + REQUIRE( quad.npts() == 350 ); + const auto msg = "Gauss-Lobatto x Becke N = " + std::to_string(quad.npts()); + test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); + } + + SECTION("Gauss-Legendre x Becke") { + IntegratorXX::RadialTransformQuadrature< + IntegratorXX::GaussLegendre, + IntegratorXX::BeckeRadialTraits> quad(350); + REQUIRE( quad.npts() == 350 ); + const auto msg = "Gauss-Legendre x Becke N = " + std::to_string(quad.npts()); + test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); + } + + SECTION("Every primitive composes with every transform") { + // All 28 pairings. Accuracy is not the point here: several of these are + // not sensible rules, but all of them must build. + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + } + +} + TEST_CASE( "Lebedev-Laikov", "[1d-quad]" ) {