numerical.h

/*
 * This file contains an additional numerical integration routine for the
 * Splines based on the Gauss-Legendre routines in boost/math/quadrature.
 *
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_INTEGRATION_NUMERICAL_H
#define BSPLINE_INTEGRATION_NUMERICAL_H
 
#include <bspline/Spline.h>
#include <bspline/exceptions/BSplineException.h>
#include <bspline/internal/misc.h>
 
#include <boost/math/quadrature/gauss.hpp>
 
namespace bspline::integration {
using namespace boost::math::quadrature;
using bspline::support::Support;
using namespace bspline::exceptions;
 
template <size_t ordergl, typename T, typename F, size_t order1, size_t order2>
T integrate(const F &f, const bspline::Spline<T, order1> &m1,
            const bspline::Spline<T, order2> &m2) {
  // Will also check whether the two grids are equivalent.
  const Support newSupport = m1.getSupport().calcIntersection(m2.getSupport());
  const size_t nintervals = newSupport.numberOfIntervals();
 
  T result = static_cast<T>(0);
  for (size_t interv = 0; interv < nintervals; interv++) {
    const auto ai = newSupport.absoluteFromRelative(interv);
    const auto m1Index = m1.getSupport().intervalIndexFromAbsolute(ai).value();
    const auto m2Index = m2.getSupport().intervalIndexFromAbsolute(ai).value();
 
    const T &xstart = m1.getSupport().at(m1Index);
    const T &xend = m1.getSupport().at(m1Index + 1);
    const T xm = (xstart + xend) / static_cast<T>(2);
    const auto &c1 = m1.getCoefficients().at(m1Index);
    const auto &c2 = m2.getCoefficients().at(m2Index);
    result += gauss<T, ordergl>::integrate(
        [&c1, &c2, &xm, &f](const T &x) {
          using namespace bspline::internal;
          return f(x) * evaluateInterval(x, c1, xm) *
                 evaluateInterval(x, c2, xm);
        },
        xstart, xend);
  }
  return result;
}
}  // namespace bspline::integration
#endif  // BSPLINE_INTEGRATION_NUMERICAL_H