Added ability to generate and destroy quadratures, tested with Unifor…

…mTrapezoid

include/integratorxx/c_api.h

@@ -66,12 +66,14 @@ typedef struct {
 
   intxx_quad_params_type ext_params; ///< External params
 
-  //void (*init)(struct intxx_quad_type* p);
+  int (*generate)(struct intxx_quad_type* p);
+  int (*destroy) (struct intxx_quad_type* p);
 } intxx_quad_info_type;
 
 struct intxx_quad_type {
   int npoints;
   const intxx_quad_info_type* info;
+  void* _state; ///< Internal state, NOT FOR PUBLIC USE
 };
 
 typedef struct intxx_quad_type intxx_quad_type;
@@ -135,6 +137,10 @@ int intxx_quad_set_npts(intxx_quad_type* p, int npts);
  */ 
 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);
+
 #ifdef __cplusplus
 }
 #endif

src/c_api/c_api.c

@@ -11,6 +11,7 @@ int intxx_quad_init(intxx_quad_type* p, int quad) {
   // Sanity check
   if(p == NULL) return INTXX_NULL_QUADPTR;
   p->info    = NULL;
+  p->_state  = NULL;
   p->npoints = -1;
 
   if(quad < 0) return INTXX_INVALID_QUAD;
@@ -84,4 +85,17 @@ int intxx_quad_get_npts(intxx_quad_type* p, int* npts) {
   *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;
+
+  return p->info->generate(p);
+}
 
+int intxx_destroy_quad(intxx_quad_type* p) {
+  if(p == NULL) return INTXX_NULL_QUADPTR;
+  if(p->info == NULL) return INTXX_NULL_INFOPTR;
+
+  return p->info->destroy(p);
+}

src/c_api/c_radial.cxx

@@ -1,3 +1,5 @@
+#include <integratorxx/quadratures/primitive/uniform.hpp>
+
 #include "c_internal.h"
 #include <cstddef>
 
@@ -11,11 +13,15 @@ int intxx_get_gausscheb2_info(intxx_quad_info_type* p);
 int intxx_get_gausscheb2m_info(intxx_quad_info_type* p);
 int intxx_get_gausscheb3_info(intxx_quad_info_type* p);
 
+int intxx_generate_uniform(intxx_quad_type* p);
+int intxx_destroy_uniform(intxx_quad_type* p);
+
 void intxx_default_quad_info(intxx_quad_info_type* p) {
   p->number = 0; // Invalid
   p->kind   = 0; // Invalid
   p->name   = "Default"; // Default state
-  //p->init   = NULL;
+  p->generate = NULL;
+  p->destroy  = NULL;
 }
 
 int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) {
@@ -51,34 +57,98 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) {
   }
 }
 
-int intxx_noparam_info(intxx_quad_info_type* p) {
+int intxx_noparam_info(intxx_quad_info_type* p, 
+  int (*g)(intxx_quad_type*),
+  int (*d)(intxx_quad_type*)) {
   if(p == NULL) return INTXX_NULL_INFOPTR;
 
   p->ext_params.n = 0; /// No External Parameters
+  p->generate = g;
+  p->destroy  = d;
 
   return INTXX_SUCCESS;
 }
 
 int intxx_get_uniform_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, &intxx_generate_uniform,
+    &intxx_destroy_uniform);
 }
 int intxx_get_gausslob_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
 }
 int intxx_get_gaussleg_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
 }
 int intxx_get_gausscheb1_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
 }
 int intxx_get_gausscheb2_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
 }
 int intxx_get_gausscheb2m_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
 }
 int intxx_get_gausscheb3_info(intxx_quad_info_type* p) {
-  return intxx_noparam_info(p);
+  return intxx_noparam_info(p, NULL, NULL);
+}
+
+
+
+int intxx_generate_uniform(intxx_quad_type* p) {
+  if(p == NULL) return INTXX_NULL_QUADPTR;
+  if(p->info == NULL) return INTXX_NULL_INFOPTR;
+
+  using namespace IntegratorXX;
+  using quad_type = UniformTrapezoid<double,double>;
+  using alloc_type = std::allocator<quad_type>;
+  using alloc_traits = std::allocator_traits<alloc_type>;
+
+  int npts = p->npoints;
+  if( npts <= 0 ) {
+    // Return error code
+  }
+
+  if(p->_state != NULL) {
+    // Check if we need to regenerate
+  }
+
+
+  // Allocate and construct the quadrature instance
+  alloc_type alloc;
+  auto ptr = alloc_traits::allocate(alloc, 1);
+  alloc_traits::construct(alloc, ptr, npts, 0.0, 1.0);
+
+  // Store the pointer
+  p->_state = (void*)ptr;
+
+  return INTXX_SUCCESS;
+
+} 
+
+int intxx_destroy_uniform(intxx_quad_type* p) {
+  if(p == NULL) return INTXX_NULL_QUADPTR;
+  if(p->info == NULL) return INTXX_NULL_INFOPTR;
+
+  using namespace IntegratorXX;
+  using quad_type = UniformTrapezoid<double,double>;
+  using alloc_type = std::allocator<quad_type>;
+  using alloc_traits = std::allocator_traits<alloc_type>;
+
+  if(p->_state != NULL) { 
+    alloc_type alloc;
+
+    // Cross your fingers...
+    auto ptr = p->_state;
+    auto quad_ptr = reinterpret_cast<quad_type*>(ptr);
+
+    // Destroy and deallocate
+    alloc_traits::destroy(alloc, quad_ptr);
+    alloc_traits::deallocate(alloc, quad_ptr, 1);
+
+    // Null out state
+    p->_state = NULL;
+  }
+  return INTXX_SUCCESS;
 }
 
 }

test/c_api.cxx

@@ -1,7 +1,11 @@
 #include "catch2/catch_all.hpp"
+#include "quad_matcher.hpp"
 #include <integratorxx/quadratures/all.hpp>
 #include <integratorxx/c_api.h>
 
+
+using namespace IntegratorXX;
+
 TEST_CASE("C API") {
 
   intxx_quad_type quad;
@@ -18,10 +22,15 @@ TEST_CASE("C API") {
 
     const char* name;
     const int base_npts = 100;
+    using base_quad_type = QuadratureBase<std::vector<double>, std::vector<double>>;
+
+    std::unique_ptr<base_quad_type> base_quad = nullptr;
 
     SECTION("Uniform") {
+      using quad_type = UniformTrapezoid<double,double>;
       error = intxx_quad_init(&quad, INTXX_PRMQ_UNIFORM);
       name = "UNIFORM";
+      base_quad = std::make_unique<quad_type>(base_npts, 0.0, 1.0);
     }
 
     SECTION("Gauss-Legendre") {
@@ -59,6 +68,7 @@ TEST_CASE("C API") {
     // Meta data
     REQUIRE(quad.info != NULL);
     REQUIRE(quad.npoints == -1);
+    REQUIRE(quad._state == NULL);
     REQUIRE(quad.info->ext_params.n == 0);
     REQUIRE(!strcmp(quad.info->name, name));
 
@@ -77,6 +87,23 @@ TEST_CASE("C API") {
     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 != NULL);
+
+      // Check validity of the state
+      auto state_as_quad = reinterpret_cast<base_quad_type*>(quad._state);
+      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 == NULL);
+    }
   }
 
   intxx_quad_end(&quad);