diff --git a/include/integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp b/include/integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp index 02e729a..87aa160 100644 --- a/include/integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp +++ b/include/integratorxx/composite_quadratures/pruned_spherical_quadrature.hpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -48,6 +49,9 @@ class RadialGridPartition { template void finalize( const RadialQuad& rq ) { + if( partition_idx_.empty() ) + throw std::runtime_error( + "RadialGridPartition: cannot finalize without any angular quadrature"); if( partition_idx_.back() != rq.npts() ) partition_idx_.emplace_back( rq.npts() ); } @@ -134,6 +138,11 @@ class RadialGridPartition { return iterator( partition_idx_.begin(), quads_.begin() ); } iterator end() { + // partition_idx_ carries one entry past the last partition (the npts + // sentinel added by finalize), so the end index iterator is one before its + // end. An unfinalized partition holds no entries at all, and end()-1 would + // step before begin(); the range is simply empty in that case. + if( partition_idx_.empty() ) return begin(); return iterator( partition_idx_.end()-1, quads_.end() ); } @@ -141,6 +150,7 @@ class RadialGridPartition { return const_iterator( partition_idx_.cbegin(), quads_.cbegin() ); } const_iterator cend() const { + if( partition_idx_.empty() ) return cbegin(); return const_iterator( partition_idx_.cend()-1, quads_.cend() ); } @@ -148,6 +158,7 @@ class RadialGridPartition { return const_iterator( partition_idx_.cbegin(), quads_.cbegin() ); } const_iterator end() const { + if( partition_idx_.empty() ) return begin(); return const_iterator( partition_idx_.cend()-1, quads_.cend() ); } }; @@ -183,13 +194,13 @@ class PrunedSphericalQuadrature : const point_container& sph_points_adaptor() const override { return quad_base_type::points(); } - point_container& sph_points_adaptor() override {; + point_container& sph_points_adaptor() override { return quad_base_type::points(); } - const weight_container& sph_weights_adaptor() const override {; + const weight_container& sph_weights_adaptor() const override { return quad_base_type::weights(); } - weight_container& sph_weights_adaptor() override {; + weight_container& sph_weights_adaptor() override { return quad_base_type::weights(); } diff --git a/include/integratorxx/composite_quadratures/spherical_quadrature.hpp b/include/integratorxx/composite_quadratures/spherical_quadrature.hpp index 09ef168..e7c0c0d 100644 --- a/include/integratorxx/composite_quadratures/spherical_quadrature.hpp +++ b/include/integratorxx/composite_quadratures/spherical_quadrature.hpp @@ -133,13 +133,13 @@ class SphericalQuadrature : const point_container& sph_points_adaptor() const override { return quad_base_type::points(); } - point_container& sph_points_adaptor() override {; + point_container& sph_points_adaptor() override { return quad_base_type::points(); } - const weight_container& sph_weights_adaptor() const override {; + const weight_container& sph_weights_adaptor() const override { return quad_base_type::weights(); } - weight_container& sph_weights_adaptor() override {; + weight_container& sph_weights_adaptor() override { return quad_base_type::weights(); } diff --git a/include/integratorxx/composite_quadratures/sub_quadrature.hpp b/include/integratorxx/composite_quadratures/sub_quadrature.hpp index 0f35375..9ff10f4 100644 --- a/include/integratorxx/composite_quadratures/sub_quadrature.hpp +++ b/include/integratorxx/composite_quadratures/sub_quadrature.hpp @@ -49,16 +49,13 @@ struct quadrature_traits< inline static std::tuple generate( range_idx_type idx_range, const Quadrature& quad ) { const auto [begin_idx, end_idx] = idx_range; - assert( begin_idx >= 0 ); - assert( begin_idx <= quad.npts() ); - assert( end_idx >= 0 ); - assert( end_idx <= quad.npts() ); - assert( begin_idx < end_idx ); - - point_container points( quad.points().begin() + begin_idx, - quad.points().begin() + end_idx ); - point_container weights( quad.weights().begin() + begin_idx, - quad.weights().begin() + end_idx ); + assert( begin_idx < end_idx ); + assert( end_idx <= quad.npts() ); + + point_container points ( quad.points().begin() + begin_idx, + quad.points().begin() + end_idx ); + weight_container weights( quad.weights().begin() + begin_idx, + quad.weights().begin() + end_idx ); return std::make_tuple( points, weights ); } diff --git a/test/composite_quadratures.cxx b/test/composite_quadratures.cxx index 02fdf79..acba24e 100644 --- a/test/composite_quadratures.cxx +++ b/test/composite_quadratures.cxx @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -239,6 +240,29 @@ TEST_CASE( "Pruned Spherical Quadratures", "[sph-quad]" ) { } + SECTION("Unfinalized Partition") { + + size_t nrad = 4; + IntegratorXX::MuraKnowles r(nrad); + IntegratorXX::LebedevLaikov q(6); + IntegratorXX::RadialGridPartition rgp; + + // Nothing has been added, so there is no last index to compare against + // rq.npts(): finalize must say so rather than read back() off an empty + // vector. + REQUIRE_THROWS( rgp.finalize(r) ); + + // The iterators must present an empty range over it. partition_idx_ is + // empty, so the end()-1 that skips the npts sentinel would step before + // begin(). + CHECK( rgp.begin() == rgp.end() ); + CHECK( rgp.cbegin() == rgp.cend() ); + const auto& crgp = rgp; + CHECK( crgp.begin() == crgp.end() ); + + } + + SECTION("Validity") { size_t nrad = 12; @@ -300,9 +324,30 @@ TEST_CASE( "Pruned Spherical Quadratures", "[sph-quad]" ) { } +TEST_CASE( "Sub Quadratures", "[sph-quad]" ) { + // The slice's weights must be built in the parent's weight_container, not in + // its point_container. The two coincide for the usual case, + // which is why nothing caught it; under mixed precision the weights would be + // narrowed to the point type. + IntegratorXX::Becke r(20); + const std::pair range{ 5ul, 12ul }; + IntegratorXX::SubQuadrature sq( range, r ); + static_assert( std::is_same_v< decltype(sq)::weight_container, + decltype(r)::weight_container >, + "sub-quadrature weights must keep the parent's weight type" ); + static_assert( std::is_same_v< decltype(sq)::point_container, + decltype(r)::point_container >, + "sub-quadrature points must keep the parent's point type" ); + REQUIRE( sq.npts() == range.second - range.first ); + // The slice is a copy, so it must agree with the parent exactly + for( size_t i = 0; i < sq.npts(); ++i ) { + CHECK( sq.points()[i] == r.points()[i + range.first] ); + CHECK( sq.weights()[i] == r.weights()[i + range.first] ); + } +}