BilinearForm.h

/*
 * ########################################################################
 * 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_BILINEARFORM_H
#define BSPLINE_INTEGRATION_BILINEARFORM_H
 
#include <bspline/Spline.h>
#include <bspline/operators/GenericOperators.h>
 
namespace bspline::integration {
 
template <typename O1, typename O2,
          std::enable_if_t<operators::are_operators_v<O1, O2>, bool> = true>
class BilinearForm final {
 private:
  O1 _o1;
  O2 _o2;
 
  template <typename T, size_t sizea, size_t sizeb>
  static T evaluateInterval(const std::array<T, sizea> &a,
                            const std::array<T, sizeb> &b, const T &dxhalf) {
    std::array<T, (sizea + sizeb) / 2> coefficients;
    coefficients.fill(static_cast<T>(0));
    for (size_t i = 0; i < sizea; i++) {
      for (size_t j = i % 2; j < sizeb; j += 2) {
        coefficients[(i + j) / 2] += a[i] * b[j];
      }
    }
 
    // Use Horner's scheme to evaluate.
    constexpr int endIndex = static_cast<int>(coefficients.size()) - 1;
    static_assert(endIndex >= 0);
 
    T result = coefficients[endIndex] / static_cast<T>(2 * endIndex + 1);
    const T dxhalf_squared = dxhalf * dxhalf;
 
    for (int i = endIndex - 1; i >= 0; i--) {
      result =
          result * dxhalf_squared + coefficients[i] / static_cast<T>(2 * i + 1);
    }
    return static_cast<T>(2) * dxhalf * result;
  }
 
 public:
  BilinearForm(O1 o1, O2 o2) : _o1(std::move(o1)), _o2(std::move(o2)){};
 
  explicit BilinearForm(O2 o2) : _o1(O1{}), _o2(std::move(o2)){};
 
  BilinearForm() : _o1(O1{}), _o2(O2{}){};
 
  template <typename T, size_t ordera, size_t orderb>
  T evaluate(const Spline<T, ordera> &a, const Spline<T, orderb> &b) const {
    // Will also check whether the two grids are equivalent.
    const support::Support integrandSupport =
        a.getSupport().calcIntersection(b.getSupport());
    const size_t nintervals = integrandSupport.numberOfIntervals();
 
    const auto &grid = integrandSupport.getGrid();
 
    T result = static_cast<T>(0);
 
    for (size_t interv = 0; interv < nintervals; interv++) {
      const auto absIndex = integrandSupport.absoluteFromRelative(interv);
 
      const auto aIndex =
          a.getSupport().intervalIndexFromAbsolute(absIndex).value();
      const auto bIndex =
          b.getSupport().intervalIndexFromAbsolute(absIndex).value();
 
      const T dxhalf = (a.getSupport()[aIndex + 1] - a.getSupport()[aIndex]) /
                       static_cast<T>(2);
 
      result += evaluateInterval(
          _o1.transform(a.getCoefficients()[aIndex], grid, absIndex),
          _o2.transform(b.getCoefficients()[bIndex], grid, absIndex), dxhalf);
    }
    return result;
  }
 
  template <typename T, size_t ordera, size_t orderb>
  T operator()(const Spline<T, ordera> &a, const Spline<T, orderb> &b) const {
    return evaluate(a, b);
  }
};
 
template <typename O2>
BilinearForm(O2 o2) -> BilinearForm<operators::IdentityOperator, O2>;
 
BilinearForm()
    ->BilinearForm<operators::IdentityOperator, operators::IdentityOperator>;
 
using ScalarProduct =
    BilinearForm<operators::IdentityOperator, operators::IdentityOperator>;
 
}  // namespace bspline::integration
#endif  // BSPLINE_INTEGRATION_BILINEARFORM_H