diff --git a/CMakeLists.txt b/CMakeLists.txt index 45337ce..1188a67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,19 +4,33 @@ cmake_minimum_required( VERSION 3.17 ) # Require CMake 3.17+ project( IntegratorXX VERSION 1.2.0 LANGUAGES CXX ) -add_library( integratorxx INTERFACE ) -add_library( IntegratorXX::IntegratorXX ALIAS integratorxx ) #option( INTEGRATORXX_ENABLE_TESTS "Enable Tets" ON ) +option( INTEGRATORXX_ENABLE_C_API "Enable C Bindings" ON ) +if(INTEGRATORXX_ENABLE_C_API) + set( INTEGRATORXX_HEADER_ONLY FALSE ) + enable_language(C) + add_subdirectory(src/c_api) +else() + # Header only if no C-API + set( INTEGRATORXX_HEADER_ONLY TRUE ) + add_library( integratorxx INTERFACE ) +endif() +add_library( IntegratorXX::IntegratorXX ALIAS integratorxx ) # Target features -target_compile_features( integratorxx INTERFACE cxx_std_17 ) +if(INTEGRATORXX_HEADER_ONLY) + set(INTEGRATORXX_PROPERTY_TYPE INTERFACE) +else() + set(INTEGRATORXX_PROPERTY_TYPE PUBLIC) +endif() +target_compile_features( integratorxx ${INTEGRATORXX_PROPERTY_TYPE} cxx_std_17 ) target_include_directories( integratorxx - INTERFACE + ${INTEGRATORXX_PROPERTY_TYPE} $ $ ) @@ -24,7 +38,7 @@ target_include_directories( integratorxx include(CheckCXXCompilerFlag) check_cxx_compiler_flag("-Wno-missing-braces" INTEGRATORXX_HAS_NO_MISSING_BRACES ) if( INTEGRATORXX_HAS_NO_MISSING_BRACES ) - target_compile_options( integratorxx INTERFACE $<$: -Wno-missing-braces> ) + target_compile_options( integratorxx ${INTEGRATORXX_PROPERTY_TYPE} $<$: -Wno-missing-braces> ) endif() diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h new file mode 100644 index 0000000..a6e1609 --- /dev/null +++ b/include/integratorxx/c_api.h @@ -0,0 +1,158 @@ +#pragma once + +/*** Error codes ***/ +#define INTXX_SUCCESS 0 +#define INTXX_INVALID_QUAD -1 +#define INTXX_NULL_QUADPTR -2 +#define INTXX_NULL_INFOPTR -3 +#define INTXX_INVALID_OPT -4 +#define INTXX_INVALID_ARG -5 +#define INTXX_INVALID_OUT -6 + +/*** Quadrature Classes ***/ +#define INTXX_PRM_QUAD 1 +#define INTXX_RAD_QUAD 2 +#define INTXX_ANG_QUAD 3 +#define INDXX_SPH_QUAD 4 + +/*** Primitive Quadratures ***/ +#define INTXX_PRMQ_MASK 0x0000FF +#define INTXX_PRMQ_UNIFORM 0x000001 +#define INTXX_PRMQ_GAUSSLEG 0x000002 +#define INTXX_PRMQ_GAUSSCHEB_1 0x000003 +#define INTXX_PRMQ_GAUSSCHEB_2 0x000004 +#define INTXX_PRMQ_GAUSSCHEB_2MOD 0x000005 +#define INTXX_PRMQ_GAUSSCHEB_3 0x000006 +#define INTXX_PRMQ_GAUSSLOB 0x000007 + +/*** Radial Quadratures ***/ +#define INTXX_RADQ_MASK 0x00FF00 +#define INTXX_RADQ_BECKE 0x000100 // Becke +#define INTXX_RADQ_MHL 0x000200 // Murray-Handy-Laming +#define INTXX_RADQ_TA 0x000300 // Treuter-Ahlrichs +#define INTXX_RADQ_MK 0x000400 // Mura-Knowles + +/*** Angular (S2) Quadratures ***/ +#define INTXX_ANGQ_MASK 0xFF0000 +#define INTXX_ANGQ_LEB 0x010000 // Lebedev-Laikov +#define INTXX_ANGQ_DEL 0x020000 // Delley +#define INTXX_ANGQ_AB 0x030000 // Ahrens-Beylkin +#define INTXX_ANGQ_WOM 0x040000 // Womersley + +#ifdef __cplusplus +extern "C" { +#endif + +// Forward defs of types +struct intxx_quad_type; + +typedef struct { + int n; ///< Number of parameters + + const char** names; ///< Names of the params + const char** descriptions; ///< Long descriptions of params + + int (*set_name)(struct intxx_quad_type* p, const char* name, double v); + ///< Set parameter function + int (*get_name)(struct intxx_quad_type* p, const char* name, double* v); + ///< Get parameter function + int (*generate)(struct intxx_quad_type* p); + int (*destroy) (struct intxx_quad_type* p); +} intxx_quad_params_type; + +typedef struct { + int number; ///< Quadrature identifier + int kind; ///< Type of quadrature (PRM, RAD, ANG, SPH) + int dim; ///< Dimension of the grid (1D or 3D) + + const char* name; ///< Name of the functional, e.g. "Becke" + // TODO References + + intxx_quad_params_type ext_params; ///< External params + + int (*generate)(struct intxx_quad_type* p); + int (*destroy) (struct intxx_quad_type* p); + int (*set_npts)(struct intxx_quad_type* p, int); +} intxx_quad_info_type; + +struct intxx_quad_type { + int npoints; + const intxx_quad_info_type* info; + void* _state_quad; ///< Internal state, NOT FOR EXTERNAL USE + void* _state_parm; ///< Internal state, NOT FOR EXTERNAL USE +}; + +typedef struct intxx_quad_type intxx_quad_type; + + + +/** + * @brief Initialize a specified quadrature + * + * See error code returns for how to interpret failures. + * + * @param[out] p Quadrature instance of the specified type + * @param[in] quad Type of the specified quadrature + * + * @returns INTXX_SUCCESS: no errors were encountered + * INTXX_NULL_QUADPTR: `p == NULL` + * INTXX_INVALID_QUAD: `quad` is invalid + */ +int intxx_quad_init(intxx_quad_type* p, int quad); + +/// Frees a quadrature instance. No throw gurantee +void intxx_quad_end (intxx_quad_type* p); + +/** + * @brief Set the number of quadrature points + * + * Sets the number of nodes for the passed quadrature + * instance. Only sensible for quadratures in which + * the size is a free parameter (most). + * + * See error code returns for how to interpret failures. + * + * @param[out] p The quadrature for which to set npts. + * @param[in] npts The number of points + * + * @returns INTXX_SUCCESS: no errors were encountered + * INTXX_NULL_QUADPTR: `p == NULL` + * INTXX_NULL_INFOPTR: `p->info == NULL` (uninit) + * INTXX_INVALID_ARG: invalid npts (e.g. npts < 0) + */ +int intxx_quad_set_npts(intxx_quad_type* p, int npts); + +/** + * @brief Retrieve the number of points for a quadrature + * + * Retreieves the number of grid points for a specified + * quadrature. In necessicary, this will be initialized + * from the default parameters, otherwise it will be + * read from p->npoints. + * + * See error code returns for how to interpret failures. + * + * @param[in] p The quadrature for which to get npts. + * @param[out] npts The number of points + * + * @returns INTXX_SUCCESS: no errors were encountered + * INTXX_NULL_QUADPTR: `p == NULL` + * INTXX_NULL_INFOPTR: `p->info == NULL` (uninit) + * INTXX_INVALID_ARG: `npts == NULL` + * INTXX_INVALID_OUT: `npts` is not valid (npts < 0) + */ +int intxx_quad_get_npts(intxx_quad_type* p, int* npts); + + +int intxx_generate_quad(intxx_quad_type* p); +int intxx_destroy_quad(intxx_quad_type* p); + +int intxx_get_ext_param_name(intxx_quad_type* p, const char* name, double* v); +int intxx_set_ext_param_name(intxx_quad_type* p, const char* name, double v); + +int intxx_get_rad_scal(intxx_quad_type* p, double *R); +int intxx_set_rad_scal(intxx_quad_type* p, double R); + +#ifdef __cplusplus +} +#endif diff --git a/include/integratorxx/quadratures/radial/becke.hpp b/include/integratorxx/quadratures/radial/becke.hpp index 1c711ae..7c56fa0 100644 --- a/include/integratorxx/quadratures/radial/becke.hpp +++ b/include/integratorxx/quadratures/radial/becke.hpp @@ -47,6 +47,9 @@ class BeckeRadialTraits { inline auto radial_jacobian(PointType x) const noexcept { return R_ * 2.0 / std::pow(1.0 - x, 2); } + + /// Return radial scaling factor + auto R() const { return R_; } }; /** diff --git a/include/integratorxx/quadratures/radial/mhl.hpp b/include/integratorxx/quadratures/radial/mhl.hpp index a45fe44..91565d8 100644 --- a/include/integratorxx/quadratures/radial/mhl.hpp +++ b/include/integratorxx/quadratures/radial/mhl.hpp @@ -47,6 +47,8 @@ class MurrayHandyLamingRadialTraits { return R_ * M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); } + /// Return radial scaling factor + auto R() const { return R_; } }; diff --git a/include/integratorxx/quadratures/radial/muraknowles.hpp b/include/integratorxx/quadratures/radial/muraknowles.hpp index 23a26ff..cf00f0b 100644 --- a/include/integratorxx/quadratures/radial/muraknowles.hpp +++ b/include/integratorxx/quadratures/radial/muraknowles.hpp @@ -155,6 +155,8 @@ class MuraKnowlesRadialTraits { return R_ * 3.0 * x2 / (1.0 - x2 * x); } + /// Return radial scaling factor + auto R() const { return R_; } }; template diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index 93f8096..418d00b 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial/radial_transform.hpp @@ -12,6 +12,9 @@ struct RadialTransformQuadrature : public: + using quad_type = BaseQuad; + using traits_type = RadialTraits; + using point_type = typename base_type::point_type; using weight_type = typename base_type::weight_type; using point_container = typename base_type::point_container; diff --git a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp index 9aa5e96..46c41a6 100644 --- a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp @@ -64,6 +64,8 @@ class TreutlerAhlrichsRadialTraits { return R_ * pow_term / ln_2 * ( alpha_ * log_term / (a+x) + (1./(1.-x)) ); } + /// Return radial scaling factor + auto R() const { return R_; } }; diff --git a/src/c_api/CMakeLists.txt b/src/c_api/CMakeLists.txt new file mode 100644 index 0000000..06c0cbb --- /dev/null +++ b/src/c_api/CMakeLists.txt @@ -0,0 +1 @@ +add_library( integratorxx c_api.c c_primitive.cxx c_radial.cxx c_misc.cxx c_angular.cxx ) diff --git a/src/c_api/c_angular.cxx b/src/c_api/c_angular.cxx new file mode 100644 index 0000000..76d572f --- /dev/null +++ b/src/c_api/c_angular.cxx @@ -0,0 +1,117 @@ +#include + +#include "c_internal.h" +#include "c_api_util.hpp" + + +extern "C" { + +using leb_quad_type = IntegratorXX::LebedevLaikov; +using delley_quad_type = IntegratorXX::Delley; +using wom_quad_type = IntegratorXX::Womersley; +using ab_quad_type = IntegratorXX::AhrensBeylkin; + + +/*******************************************************/ +/****** Forward declare C-API for ANG Quadratures ******/ +/*******************************************************/ + +#define FWD_ANG(cname) \ +int intxx_get_##cname##_info(intxx_quad_info_type* p); \ +int intxx_generate_##cname(intxx_quad_type* p); \ +int intxx_destroy_##cname(intxx_quad_type* p); \ + +FWD_ANG(leb) +FWD_ANG(delley) +FWD_ANG(wom) +FWD_ANG(ab) + +#undef FWD_ANG + +/*******************************************************/ +/****** Runtime Generator for ANG Quadrature Info ******/ +/*******************************************************/ + +int intxx_get_angq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_ANGQ_MASK; + p->number = quad; + p->kind = INTXX_ANG_QUAD; + p->dim = 3; + switch(quad) { + case INTXX_ANGQ_LEB: + p->name = "LEBEDEV_LAIKOV"; + return intxx_get_leb_info(p); + case INTXX_ANGQ_DEL: + p->name = "DELLEY"; + return intxx_get_delley_info(p); + case INTXX_ANGQ_AB: + p->name = "AHRENS_BEYLKIN"; + return intxx_get_ab_info(p); + case INTXX_ANGQ_WOM: + p->name = "WOMERSLEY"; + return intxx_get_wom_info(p); + default: + return INTXX_INVALID_QUAD; + } + return INTXX_SUCCESS; +} + +/**********************************************/ +/****** NPTS Setters for ANG Quadratures ******/ +/**********************************************/ + +#define INTXX_ANG_SET_NPTS(cname, nsp) \ +int intxx_##cname##_set_npts(intxx_quad_type* p, int npts) {\ + using namespace IntegratorXX::detail::nsp;\ + if(algebraic_order_by_npts(npts) > 0) {\ + p->npoints = npts;\ + return INTXX_SUCCESS;\ + } else {\ + p->npoints = -1;\ + return INTXX_INVALID_ARG;\ + }\ +} + +INTXX_ANG_SET_NPTS(leb, lebedev) +INTXX_ANG_SET_NPTS(delley, delley) +INTXX_ANG_SET_NPTS(ab, ahrensbeylkin) +INTXX_ANG_SET_NPTS(wom, womersley) + +/*************************************************/ +/****** Info generation for ANG Quadratures ******/ +/*************************************************/ + +#define INTXX_NOPARAM_GET_INFO_IMPL(cname) \ +int intxx_get_##cname##_info(intxx_quad_info_type* p) { \ + return intxx_noparam_info(p, &intxx_generate_##cname, \ + &intxx_destroy_##cname, intxx_##cname##_set_npts); \ +} + +INTXX_NOPARAM_GET_INFO_IMPL(leb); +INTXX_NOPARAM_GET_INFO_IMPL(delley); +INTXX_NOPARAM_GET_INFO_IMPL(ab); +INTXX_NOPARAM_GET_INFO_IMPL(wom); + +#undef INTXX_NOPARAM_GET_INFO_IMPL + +/*************************************************************/ +/****** Allocate/Deallocate ANG Quadrature States ******/ +/*************************************************************/ + +#define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ +int intxx_generate_##cname(intxx_quad_type* p) { \ + return intxx_generate_noparam_impl(p, ##__VA_ARGS__); \ +} \ +int intxx_destroy_##cname(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ +} + +INTXX_GD_BASIC_IMPL(leb, leb_quad_type); +INTXX_GD_BASIC_IMPL(delley, delley_quad_type); +INTXX_GD_BASIC_IMPL(ab, ab_quad_type); +INTXX_GD_BASIC_IMPL(wom, wom_quad_type); + +#undef INTXX_GD_BASIC_IMPL +} diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c new file mode 100644 index 0000000..abe3e34 --- /dev/null +++ b/src/c_api/c_api.c @@ -0,0 +1,157 @@ +#include "c_internal.h" +#include +#include + + + +int intxx_quad_init(intxx_quad_type* p, int quad) { + // Local type defs + typedef intxx_quad_info_type info_type; + + // Sanity check + if(p == NULL) return INTXX_NULL_QUADPTR; + p->info = NULL; + p->_state_quad = NULL; + p->_state_parm = NULL; + p->npoints = -1; + + if(quad < 0) return INTXX_INVALID_QUAD; + + // Determine quadrature class + int is_prmq = quad & INTXX_PRMQ_MASK; + int is_radq = quad & INTXX_RADQ_MASK; + int is_angq = quad & INTXX_ANGQ_MASK; + int is_sphq = is_radq && is_angq; + + // Passed quadrature had to be something sane + if(!is_prmq && !is_radq && !is_angq) + return INTXX_INVALID_QUAD; + + // Primitive quadratures cannot be mixed with angular or + // radial quadratures + if(is_prmq && is_radq) return INTXX_INVALID_QUAD; + if(is_prmq && is_angq) return INTXX_INVALID_QUAD; + + // Get info + info_type* finfo = (info_type*)malloc(sizeof(info_type)); + intxx_default_quad_info(finfo); // set (invalid) state + p->info = finfo; + + int error; + if(is_prmq) { + // Get primitive quadrature info + error = intxx_get_prmq_info(finfo, quad); + } else if(is_sphq) { + // TODO: Get spherical quadrature info + } else if(is_radq) { + // TODO: Get radial quadrature info + error = intxx_get_radq_info(finfo, quad); + } else { + // Angular by exclusion + error = intxx_get_angq_info(finfo, quad); + } + + if(error) return error; + + if(finfo->ext_params.n && finfo->ext_params.generate) { + error = finfo->ext_params.generate(p); + } + + return error; +} + +void intxx_quad_end(intxx_quad_type* p) { + // Stateless - just bail + if(p == NULL || p->info == NULL) return; + + // Destroy traits if required + if(p->info->ext_params.n && p->info->ext_params.destroy) { + p->info->ext_params.destroy(p); + } + + // Destroy quadrature state if populated + intxx_destroy_quad(p); + + // Free info + free((void*)p->info); + p->info = NULL; +} + + + +int intxx_quad_set_npts(intxx_quad_type* p, int npts) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + // NPTS must be > 0 + if(npts <= 0) return INTXX_INVALID_ARG; + + // TODO: Handle the case when NPTS is derived from params + if(p->info->set_npts) { + // Some quadratures have rules about how to set npts (e.g. angular grids) + return p->info->set_npts(p, npts); + } else { + p->npoints = npts; + return INTXX_SUCCESS; + } +} + +int intxx_quad_get_npts(intxx_quad_type* p, int* npts) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + // Return memory must be valid + if(npts == NULL) return INTXX_INVALID_ARG; + + // TODO: Handle the case when NPTS is derived from params + *npts = p->npoints; + return *npts > 0 ? INTXX_SUCCESS : INTXX_INVALID_OUT; +} + +int intxx_generate_quad(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(p->info->generate) + return p->info->generate(p); + else return INTXX_SUCCESS; +} + +int intxx_destroy_quad(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(p->info->destroy) + return p->info->destroy(p); + else return INTXX_SUCCESS; +} + +int intxx_get_ext_param_name(intxx_quad_type* p, const char* name, double* v) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + if(p->info->ext_params.n <= 0) { + return INTXX_INVALID_OPT; + } + + if(p->info->ext_params.get_name == NULL) { + // Something more descriptive? + return INTXX_INVALID_OPT; + } + + return p->info->ext_params.get_name(p, name, v); +} + +int intxx_set_ext_param_name(intxx_quad_type* p, const char* name, double v) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + if(p->info->ext_params.n <= 0) { + return INTXX_INVALID_OPT; + } + + if(p->info->ext_params.set_name == NULL) { + // Something more descriptive? + return INTXX_INVALID_OPT; + } + + return p->info->ext_params.set_name(p, name, v); +} diff --git a/src/c_api/c_api_util.hpp b/src/c_api/c_api_util.hpp new file mode 100644 index 0000000..e736cfb --- /dev/null +++ b/src/c_api/c_api_util.hpp @@ -0,0 +1,195 @@ +#pragma once +#include "c_internal.h" +#include +#include + +template +auto generate_quadrature(Args&&... args) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + auto ptr = alloc_traits::allocate(alloc, 1); + alloc_traits::construct(alloc, ptr, std::forward(args)...); + + return ptr; +} + +template +int intxx_generate_noparam_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using quad_type = QuadType; + + int npts = p->npoints; + if( npts <= 0 ) { + // Return error code + } + + if(p->_state_quad != NULL) { + // Check if we need to regenerate + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(npts, std::forward(args)... ); + + // Store the pointer + p->_state_quad = (void*)ptr; + + return INTXX_SUCCESS; + +} + + +template +int intxx_generate_param_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm == NULL) { + // Return error code + } + + // Cross your fingers... + const auto _params = reinterpret_cast(p->_state_parm); + return intxx_generate_noparam_impl(p, *_params, std::forward(args)... ); + +} + +template +void destroy_quadrature(void* ptr) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + + // Cross your fingers... + auto quad_ptr = reinterpret_cast(ptr); + + // Destroy and deallocate + alloc_traits::destroy(alloc, quad_ptr); + alloc_traits::deallocate(alloc, quad_ptr, 1); +} + +template +int intxx_destroy_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + using quad_type = QuadType; + + if(p->_state_quad != NULL) { + // Destroy quadrature + destroy_quadrature(p->_state_quad); + + // Null out state + p->_state_quad = NULL; + } + return INTXX_SUCCESS; +} + +template +int intxx_generate_params_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(p->_state_parm != NULL) { + return INTXX_SUCCESS; + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(std::forward(args)... ); + + // Store the pointer + p->_state_parm = (void*)ptr; + + return INTXX_SUCCESS; + +} + +template +int intxx_destroy_params_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm != NULL) { + // Destroy traits object + destroy_quadrature(p->_state_parm); + + // Null out state + p->_state_parm = NULL; + } + return INTXX_SUCCESS; +} + +template +int intxx_get_rad_scal_impl(intxx_quad_type* p, double* R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R == NULL) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Read data + *R = reinterpret_cast(p->_state_parm)->R(); + + return INTXX_SUCCESS; +} + +template +int intxx_set_rad_scal_impl(intxx_quad_type* p, double R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R <= 0.0) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Overwrite + *reinterpret_cast(p->_state_parm) = TraitsType(R); + + // Regenerate if quadrature is populated + if(p->_state_quad) { + intxx_destroy_quad(p); + intxx_generate_quad(p); + } + + return INTXX_SUCCESS; +} + + diff --git a/src/c_api/c_internal.h b/src/c_api/c_internal.h new file mode 100644 index 0000000..748df83 --- /dev/null +++ b/src/c_api/c_internal.h @@ -0,0 +1,30 @@ +#pragma once +#include + +#ifdef __cplusplus +extern "C" { +#endif + +void intxx_default_quad_info(intxx_quad_info_type* p); +int intxx_get_radq_info(intxx_quad_info_type* p, int quad); +int intxx_get_angq_info(intxx_quad_info_type* p, int quad); +int intxx_get_prmq_info(intxx_quad_info_type* p, int quad); + + +int intxx_noparam_info(intxx_quad_info_type* p, + int (*g)(intxx_quad_type*), + int (*d)(intxx_quad_type*), + int (*s)(intxx_quad_type*, int)); + +int intxx_radscal_info(intxx_quad_info_type* p, + int (*sn)(intxx_quad_type*,const char*,double), + int (*gn)(intxx_quad_type*,const char*,double*), + int (*gparm)(intxx_quad_type*), + int (*dparm)(intxx_quad_type*), + int (*gquad)(intxx_quad_type*), + int (*dquad)(intxx_quad_type*)); + +#ifdef __cplusplus +} +#endif + diff --git a/src/c_api/c_misc.cxx b/src/c_api/c_misc.cxx new file mode 100644 index 0000000..3279028 --- /dev/null +++ b/src/c_api/c_misc.cxx @@ -0,0 +1,64 @@ +#include "c_internal.h" +#include + +extern "C" { + +void intxx_default_quad_info(intxx_quad_info_type* p) { + p->number = 0; // Invalid + p->kind = 0; // Invalid + p->dim = 0; // Invalid + p->name = "Default"; // Default state + p->ext_params.n = 0; // To initialize *everything* + p->generate = NULL; // Invalid + p->destroy = NULL; // Invalid +} + +/******************************************************/ +/****** Info generation for quads w/o parameters ******/ +/******************************************************/ + +int intxx_noparam_info(intxx_quad_info_type* p, + int (*g)(intxx_quad_type*), + int (*d)(intxx_quad_type*), + int (*s)(intxx_quad_type*, int)) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + p->ext_params.n = 0; /// No External Parameters + p->generate = g; + p->destroy = d; + p->set_npts = s; + + return INTXX_SUCCESS; +} + +/********************************************************/ +/****** Info generation for quads w rad parameters ******/ +/********************************************************/ + +int intxx_radscal_info(intxx_quad_info_type* p, + int (*sn)(intxx_quad_type*,const char*,double), + int (*gn)(intxx_quad_type*,const char*,double*), + int (*gparm)(intxx_quad_type*), + int (*dparm)(intxx_quad_type*), + int (*gquad)(intxx_quad_type*), + int (*dquad)(intxx_quad_type*)) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + static const char* names[] = {"RAD_SCAL"}; + static const char* desc[] = {"Radial Scaling Factor"}; + p->ext_params = { + 1, names, desc, + sn, /* set_name */ + gn, /* get_name */ + gparm, /* generate */ + dparm /* destroy */ + }; + + p->generate = gquad; + p->destroy = dquad; + p->set_npts = NULL; + + return INTXX_SUCCESS; +} + +} diff --git a/src/c_api/c_primitive.cxx b/src/c_api/c_primitive.cxx new file mode 100644 index 0000000..9cd866e --- /dev/null +++ b/src/c_api/c_primitive.cxx @@ -0,0 +1,116 @@ +#include + +#include "c_internal.h" +#include "c_api_util.hpp" + +extern "C" { + +// C++ types for internal quadrature states +using uniform_quad_type = IntegratorXX::UniformTrapezoid; +using gaussleg_quad_type = IntegratorXX::GaussLegendre; +using gausslob_quad_type = IntegratorXX::GaussLobatto; +using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; +using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; +using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; +using gausscheb2m_quad_type = IntegratorXX::GaussChebyshev2Modified; + +/*************************************************************/ +/****** Forward declare C-API for Primitive Quadratures ******/ +/*************************************************************/ + +#define FWD_PRM(cname) \ +int intxx_get_##cname##_info(intxx_quad_info_type* p); \ +int intxx_generate_##cname(intxx_quad_type* p); \ +int intxx_destroy_##cname(intxx_quad_type* p); + +FWD_PRM(uniform) +FWD_PRM(gaussleg) +FWD_PRM(gausslob) +FWD_PRM(gausscheb1) +FWD_PRM(gausscheb2) +FWD_PRM(gausscheb2m) +FWD_PRM(gausscheb3) + +#undef FWD_PRM + +/*************************************************************/ +/****** Runtime Generator for Primitive Quadrature Info ******/ +/*************************************************************/ + +int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_PRMQ_MASK; + p->number = quad; + p->kind = INTXX_PRM_QUAD; + p->dim = 1; + switch(quad) { + case INTXX_PRMQ_UNIFORM: + p->name = "UNIFORM"; + return intxx_get_uniform_info(p); + case INTXX_PRMQ_GAUSSLEG: + p->name = "GAUSS_LEGENDRE"; + return intxx_get_gaussleg_info(p); + case INTXX_PRMQ_GAUSSCHEB_1: + p->name = "GAUSS_CHEBYSHEV_1"; + return intxx_get_gausscheb1_info(p); + case INTXX_PRMQ_GAUSSCHEB_2: + p->name = "GAUSS_CHEBYSHEV_2"; + return intxx_get_gausscheb2_info(p); + case INTXX_PRMQ_GAUSSCHEB_2MOD: + p->name = "GAUSS_CHEBYSHEV_2MOD"; + return intxx_get_gausscheb2m_info(p); + case INTXX_PRMQ_GAUSSCHEB_3: + p->name = "GAUSS_CHEBYSHEV_3"; + return intxx_get_gausscheb3_info(p); + case INTXX_PRMQ_GAUSSLOB: + p->name = "GAUSS_LOBATTO"; + return intxx_get_gausslob_info(p); + default: + return INTXX_INVALID_QUAD; + } +} + +/*******************************************************/ +/****** Info generation for Primitive Quadratures ******/ +/*******************************************************/ + +#define INTXX_NOPARAM_GET_INFO_IMPL(cname) \ +int intxx_get_##cname##_info(intxx_quad_info_type* p) { \ + return intxx_noparam_info(p, &intxx_generate_##cname, \ + &intxx_destroy_##cname, NULL); \ +} + +INTXX_NOPARAM_GET_INFO_IMPL(uniform); +INTXX_NOPARAM_GET_INFO_IMPL(gaussleg); +INTXX_NOPARAM_GET_INFO_IMPL(gausslob); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb1); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2m); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb3); + +#undef INTXX_NOPARAM_GET_INFO_IMPL + +/*************************************************************/ +/****** Allocate/Deallocate Primitive Quadrature States ******/ +/*************************************************************/ + +#define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ +int intxx_generate_##cname(intxx_quad_type* p) { \ + return intxx_generate_noparam_impl(p, ##__VA_ARGS__); \ +} \ +int intxx_destroy_##cname(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ +} + +INTXX_GD_BASIC_IMPL(uniform, uniform_quad_type, 0.0, 1.0); +INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type ); +INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb1, gausscheb1_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2, gausscheb2_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb3, gausscheb3_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2m, gausscheb2m_quad_type); + +#undef INTXX_GD_BASIC_IMPL + +} // extern C diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx new file mode 100644 index 0000000..f736755 --- /dev/null +++ b/src/c_api/c_radial.cxx @@ -0,0 +1,153 @@ +#include + +#include "c_internal.h" +#include "c_api_util.hpp" + + +extern "C" { + + +// C++ types for internal quadrature states +using becke_traits_type = IntegratorXX::BeckeRadialTraits; +using mhl_traits_type = IntegratorXX::MurrayHandyLamingRadialTraits<2>; +using ta_traits_type = IntegratorXX::TreutlerAhlrichsRadialTraits; +using mk_traits_type = IntegratorXX::MuraKnowlesRadialTraits; + +using becke_quad_type = IntegratorXX::Becke; +using mhl_quad_type = IntegratorXX::MurrayHandyLaming; +using ta_quad_type = IntegratorXX::TreutlerAhlrichs; +using mk_quad_type = IntegratorXX::MuraKnowles; + +/**********************************************************/ +/****** Forward declare C-API for Radial Quadratures ******/ +/**********************************************************/ + +#define FWD_RAD(cname) \ +int intxx_get_##cname##_quad_info(intxx_quad_info_type* p); \ +int intxx_generate_##cname##_quad(intxx_quad_type* p); \ +int intxx_destroy_##cname##_quad(intxx_quad_type* p); \ +int intxx_generate_##cname##_quad_params(intxx_quad_type* p); \ +int intxx_destroy_##cname##_quad_params(intxx_quad_type* p); \ +int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v); \ +int intxx_set_name_##cname##_quad(intxx_quad_type* p, const char* name, double v); + +FWD_RAD(becke) +FWD_RAD(mhl) +FWD_RAD(ta) +FWD_RAD(mk) + +#undef FWD_RAD + +/**********************************************************/ +/****** Runtime Generator for Radial Quadrature Info ******/ +/**********************************************************/ + +int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_RADQ_MASK; + p->number = quad; + p->kind = INTXX_RAD_QUAD; + p->dim = 1; + switch(quad) { + case INTXX_RADQ_BECKE: + p->name = "BECKE"; + return intxx_get_becke_quad_info(p); + case INTXX_RADQ_MHL: + p->name = "MURRAY_HANDY_LAMING"; + return intxx_get_mhl_quad_info(p); + case INTXX_RADQ_TA: + p->name = "TREUTLER_AHLRICHS"; + return intxx_get_ta_quad_info(p); + case INTXX_RADQ_MK: + p->name = "MURA_KNOWLES"; + return intxx_get_mk_quad_info(p); + default: + return INTXX_INVALID_QUAD; + } + return INTXX_SUCCESS; +} + +/****************************************************/ +/****** Info generation for Radial Quadratures ******/ +/****************************************************/ + +#define INTXX_RADSCAL_GET_INFO_IMPL(cname) \ +int intxx_get_##cname##_quad_info(intxx_quad_info_type* p) { \ + return intxx_radscal_info(p, \ + intxx_set_name_##cname##_quad, \ + intxx_get_name_##cname##_quad, \ + intxx_generate_##cname##_quad_params, \ + intxx_destroy_##cname##_quad_params, \ + intxx_generate_##cname##_quad, \ + intxx_destroy_##cname##_quad \ + ); \ +} + +INTXX_RADSCAL_GET_INFO_IMPL(becke) +INTXX_RADSCAL_GET_INFO_IMPL(mhl) +INTXX_RADSCAL_GET_INFO_IMPL(ta) +INTXX_RADSCAL_GET_INFO_IMPL(mk) + +#undef INTXX_RADSCAL_GET_INFO_IMPL + +/***********************************************/ +/****** Allocate/Deallocate Radial Traits ******/ +/***********************************************/ + +#define INTXX_RAD_TRAITS_IMPL(cname, traits_type) \ +int intxx_generate_##cname##_quad_params(intxx_quad_type* p) { \ + return intxx_generate_params_impl(p); \ +} \ +int intxx_destroy_##cname##_quad_params(intxx_quad_type* p) { \ + return intxx_destroy_params_impl(p); \ +} + +INTXX_RAD_TRAITS_IMPL(becke, becke_traits_type); +INTXX_RAD_TRAITS_IMPL(mhl, mhl_traits_type ); +INTXX_RAD_TRAITS_IMPL(ta, ta_traits_type ); +INTXX_RAD_TRAITS_IMPL(mk, mk_traits_type ); + +#undef INTXX_RAD_TRAITS_IMPL + +/**********************************/ +/****** Radial Quad GET_NAME ******/ +/**********************************/ + +#define INTXX_RAD_GETSET_IMPL(cname, traits_type) \ +int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v) { \ + if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ + return intxx_get_rad_scal_impl(p, v);\ +}\ +int intxx_set_name_##cname##_quad(intxx_quad_type* p, const char* name, double v) { \ + if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ + return intxx_set_rad_scal_impl(p, v);\ +} + +INTXX_RAD_GETSET_IMPL(becke, becke_traits_type); +INTXX_RAD_GETSET_IMPL(mhl, mhl_traits_type ); +INTXX_RAD_GETSET_IMPL(ta, ta_traits_type ); +INTXX_RAD_GETSET_IMPL(mk, mk_traits_type ); + +#undef INTXX_RAD_GETSET_IMPL + +/**********************************************************/ +/****** Allocate/Deallocate Radial Quadrature States ******/ +/**********************************************************/ + +#define INTXX_GD_RAD_IMPL(cname, qname, tname) \ +int intxx_generate_##cname##_quad(intxx_quad_type* p) { \ + return intxx_generate_param_impl(p); \ +}\ +int intxx_destroy_##cname##_quad(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ +} + +INTXX_GD_RAD_IMPL(becke, becke_quad_type, becke_traits_type); +INTXX_GD_RAD_IMPL(mhl, mhl_quad_type, mhl_traits_type); +INTXX_GD_RAD_IMPL(ta, ta_quad_type, ta_traits_type); +INTXX_GD_RAD_IMPL(mk, mk_quad_type, mk_traits_type); + +#undef INTXX_GD_RAD_IMPL + +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 865e52e..ea59a35 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -24,6 +24,9 @@ target_link_libraries( integratorxx_common_ut PUBLIC Catch2::Catch2WithMain inte add_executable( quadrature_manipulation quadrature_manipulation.cxx ) target_link_libraries( quadrature_manipulation PUBLIC integratorxx_common_ut ) +add_executable( c_api c_api.cxx ) +target_link_libraries( c_api PUBLIC integratorxx_common_ut ) + add_executable( 1d_quadratures 1d_quadratures.cxx ) target_link_libraries( 1d_quadratures PUBLIC integratorxx_common_ut ) @@ -51,3 +54,12 @@ add_test( NAME QUADRATURES_LEGENDRE COMMAND gausslegendre ) add_test( NAME QUADRATURES_LOBATTO COMMAND gausslobatto ) add_test( NAME QUADRATURES_CHEBYSHEV1 COMMAND gausschebyshev1 ) add_test( NAME QUADRATURES_CHEBYSHEV2 COMMAND gausschebyshev2 ) +add_test( NAME C_API COMMAND c_api ) + +# Memory check +find_program( VALGRIND_EXECUTABLE valgrind ) +if(VALGRIND_EXECUTABLE) + message(STATUS "Found Valgrind: ${VALGRIND_EXECUTABLE}") + message(STATUS "Enabling C_API_MEMORY test") + add_test( NAME C_API_MEMORY COMMAND ${VALGRIND_EXECUTABLE} --leak-check=full --error-exitcode=100 $ ) +endif() diff --git a/test/c_api.cxx b/test/c_api.cxx new file mode 100644 index 0000000..2afe0a0 --- /dev/null +++ b/test/c_api.cxx @@ -0,0 +1,374 @@ +#include "catch2/catch_all.hpp" +#include "quad_matcher.hpp" +#include +#include + + +using namespace IntegratorXX; + +TEST_CASE("C API") { + + intxx_quad_type quad; + int error; + + SECTION("Invalid") { + error = intxx_quad_init(&quad, -1); + REQUIRE(error == INTXX_INVALID_QUAD); + REQUIRE(quad.info == NULL); + REQUIRE(quad.npoints == -1); + } + + + SECTION("Primitive") { + + const char* name; + const int base_npts = 100; + using base_quad_type = QuadratureBase, std::vector>; + std::unique_ptr base_quad = nullptr; + + int quad_num; + SECTION("Uniform") { + using quad_type = UniformTrapezoid; + quad_num = INTXX_PRMQ_UNIFORM; + name = "UNIFORM"; + base_quad = std::make_unique(base_npts, 0.0, 1.0); + } + + SECTION("Gauss-Legendre") { + using quad_type = GaussLegendre; + quad_num = INTXX_PRMQ_GAUSSLEG; + name = "GAUSS_LEGENDRE"; + base_quad = std::make_unique(base_npts); + } + + SECTION("Gauss-Lobatto") { + using quad_type = GaussLobatto; + quad_num = INTXX_PRMQ_GAUSSLOB; + name = "GAUSS_LOBATTO"; + base_quad = std::make_unique(base_npts); + } + + SECTION("Gauss-Chebyshev 1") { + using quad_type = GaussChebyshev1; + quad_num = INTXX_PRMQ_GAUSSCHEB_1; + name = "GAUSS_CHEBYSHEV_1"; + base_quad = std::make_unique(base_npts); + } + + SECTION("Gauss-Chebyshev 2") { + using quad_type = GaussChebyshev2; + quad_num = INTXX_PRMQ_GAUSSCHEB_2; + name = "GAUSS_CHEBYSHEV_2"; + base_quad = std::make_unique(base_npts); + } + + SECTION("Gauss-Chebyshev 2MOD") { + using quad_type = GaussChebyshev2Modified; + quad_num = INTXX_PRMQ_GAUSSCHEB_2MOD; + name = "GAUSS_CHEBYSHEV_2MOD"; + base_quad = std::make_unique(base_npts); + } + + SECTION("Gauss-Chebyshev 3") { + using quad_type = GaussChebyshev3; + quad_num = INTXX_PRMQ_GAUSSCHEB_3; + name = "GAUSS_CHEBYSHEV_3"; + base_quad = std::make_unique(base_npts); + } + + // Initialize + error = intxx_quad_init(&quad, quad_num); + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_PRM_QUAD); + REQUIRE(quad.info->dim == 1); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); + REQUIRE(quad.info->ext_params.n == 0); + REQUIRE(!strcmp(quad.info->name, name)); + + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS + error = intxx_quad_set_npts(&quad, base_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == base_npts); + + // Get NPTS + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == base_npts); + + // Check Quadrature Generation and Destruction + if(base_quad) { + intxx_generate_quad(&quad); + REQUIRE(quad._state_quad != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad->weights()[i], 1e-15)); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state_quad == NULL); + } + } + + + SECTION("Radial") { + const char* name; + const int base_npts = 100; + const double RSCAL = 2.0; + using base_quad_type = QuadratureBase, std::vector>; + std::unique_ptr base_quad_default = nullptr; + std::unique_ptr base_quad_scaled = nullptr; + + int quad_num; + SECTION("Becke") { + using quad_type = Becke; + using traits_type = typename quad_type::traits_type; + quad_num = INTXX_RADQ_BECKE; + name = "BECKE"; + base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); + } + + SECTION("MHL") { + using quad_type = MurrayHandyLaming; + using traits_type = typename quad_type::traits_type; + quad_num = INTXX_RADQ_MHL; + name = "MURRAY_HANDY_LAMING"; + base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); + } + + SECTION("TA") { + using quad_type = TreutlerAhlrichs; + using traits_type = typename quad_type::traits_type; + quad_num = INTXX_RADQ_TA; + name = "TREUTLER_AHLRICHS"; + base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); + } + + SECTION("MK") { + using quad_type = MuraKnowles; + using traits_type = typename quad_type::traits_type; + quad_num = INTXX_RADQ_MK; + name = "MURA_KNOWLES"; + base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); + } + + // Initialize + error = intxx_quad_init(&quad, quad_num); + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_RAD_QUAD); + REQUIRE(quad.info->dim == 1); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); + REQUIRE(!strcmp(quad.info->name, name)); + + + // Check default parameters + REQUIRE(quad.info->ext_params.n == 1); + REQUIRE(!strcmp(quad.info->ext_params.names[0], "RAD_SCAL")); + REQUIRE(!strcmp(quad.info->ext_params.descriptions[0], "Radial Scaling Factor")); + + double R; + intxx_get_ext_param_name(&quad, "RAD_SCAL", &R); + REQUIRE(R == 1.0); + + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS + error = intxx_quad_set_npts(&quad, base_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == base_npts); + + // Get NPTS + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == base_npts); + + // Check Quadrature Generation and Destruction + if(base_quad_default) { + intxx_generate_quad(&quad); + REQUIRE(quad._state_quad != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state_quad == NULL); + } + + // Check Scaling Works + if(base_quad_scaled) { + intxx_set_ext_param_name(&quad, "RAD_SCAL", RSCAL); + + // Check change stuck + intxx_get_ext_param_name(&quad, "RAD_SCAL", &R); + REQUIRE(R == RSCAL); + + // Regenerate and check + intxx_generate_quad(&quad); + + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_scaled->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_scaled->weights()[i], 1e-15)); + } + + + // Set without destroying to make sure quadratures are regerated + intxx_set_ext_param_name(&quad, "RAD_SCAL", 1.0); + state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + } + + + } + } + + SECTION("Angular") { + const char* name; + using base_quad_type = QuadratureBase>, std::vector>; + std::unique_ptr base_quad_default = nullptr; + //std::unique_ptr base_quad_other = nullptr; + + int default_npts, other_npts; + int quad_num; + SECTION("LebedevLaikov") { + using quad_type = LebedevLaikov; + quad_num = INTXX_ANGQ_LEB; + name = "LEBEDEV_LAIKOV"; + default_npts = 302; + other_npts = 974; + } + + SECTION("Delley") { + using quad_type = Delley; + quad_num = INTXX_ANGQ_DEL; + name = "DELLEY"; + default_npts = 302; + other_npts = 974; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); + } + + SECTION("AB") { + using quad_type = AhrensBeylkin; + quad_num = INTXX_ANGQ_AB; + name = "AHRENS_BEYLKIN"; + default_npts = 372; + other_npts = 972; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); + } + + SECTION("WOM") { + using quad_type = Womersley; + quad_num = INTXX_ANGQ_WOM; + name = "WOMERSLEY"; + default_npts = 393; + other_npts = 969; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); + } + + // Initialize + error = intxx_quad_init(&quad, quad_num); + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_ANG_QUAD); + REQUIRE(quad.info->dim == 3); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); + REQUIRE(!strcmp(quad.info->name, name)); + + + // Check default parameters + REQUIRE(quad.info->ext_params.n == 0); + + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS (incorrect) + error = intxx_quad_set_npts(&quad, default_npts+1); + REQUIRE(error == INTXX_INVALID_ARG); + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS (correct) + error = intxx_quad_set_npts(&quad, default_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == default_npts); + + // Get NPTS + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == default_npts); + + // Check Quadrature Generation and Destruction + if(base_quad_default) { + intxx_generate_quad(&quad); + REQUIRE(quad._state_quad != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE(state_as_quad->points()[i] == base_quad_default->points()[i]); + REQUIRE(state_as_quad->weights()[i] == base_quad_default->weights()[i]); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state_quad == NULL); + } + + } + + intxx_quad_end(&quad); + +}