Integrate.hpp

#pragma once

#include <type_traits>
#include <dune/geometry/quadraturerules.hh>
#include <amdis/GridFunctions.hpp>

namespace AMDiS
{
  namespace Impl
  {
    template <class GF, class GV, class QP>
    auto integrateImpl(GF&& gf, GV const& gridView, QP makeQuad)
    {
      auto localFct = localFunction(FWD(gf));

      using GridFct = remove_cvref_t<GF>;
      using Range = typename GridFct::Range;
      Range result(0);

      for (auto const& element : elements(gridView)) {
        auto geometry = element.geometry();

        localFct.bind(element);
        auto const& quad = makeQuad(element.type(), localFct);

        for (auto const qp : quad)
          result += localFct(qp.position()) * geometry.integrationElement(qp.position()) * qp.weight();
        localFct.unbind();
      }

      return gridView.comm().sum(result);
    }

    template <class GF, class GV, class QP>
    auto integrateImpl(GF&& gf, GV const& gv, QP makeQuad, std::true_type)
    {
      return integrateImpl(FWD(gf), gv, makeQuad);
    }

    template <class GF, class GV, class QP>
    auto integrateImpl(GF&&, GV&&, QP&&, std::false_type) { return 0.0; }

  } // end namespace Impl


  /// \brief Integrate expression with quadrature rule determined by polynomial order of expression
  /**
   * **Example:**
   * ```
   * double i1 = integrate(prob.solution(0), prob.gridView());
   * double i2 = integrate(X(0) + X(1), prob.gridView());
   * double i3 = integrate([](auto const& x) { return x[0] + x[1]; }, prob.gridView()); // ERROR
   * ```
   **/
  template <class Expr, class GridView>
  auto integrate(Expr&& expr, GridView const& gridView)
  {
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);

    // test whether the gridFct model `Concepts::HasLocalFunctionOrder`
    using GF = TYPEOF(gridFct);
    static const bool expr_has_order = Concepts::HasLocalFunctionOrder<GF>;
    static_assert(expr_has_order,
      "Polynomial degree of expression can not be deduced. You need to provide an explicit value for the quadrature degree or a quadrature rule in `integrate()`.");

    using Rules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
    auto makeQuad = [](auto&& t, auto&& lf) { return Rules::rule(t, order(lf)); };

    return Impl::integrateImpl(FWD(gridFct), gridView, makeQuad, bool_<expr_has_order>);
  }


  /// \brief Integrate expression with quadrature rule provided
  /**
   * **Example:**
   * ```
   * auto quad = Dune::QuadratureRules<double,DIM>::rule(Dune::GeometryTypes::triangle, 1);
   * double i4 = integrate([](auto const& x) { return x[0]; }, prob.gridView(), quad); // OK
   * ```
   **/
  template <class Expr, class GridView,
    class QuadratureRule = Dune::QuadratureRule<typename GridView::ctype, GridView::dimension>>
  auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad)
  {
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
    return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto&&, auto&&) { return quad; });
  }


  /// \brief Integrate expression with quadrature rule determined by provided polynomial `degree`
  /**
   * **Example:**
   * ```
   * double i5 = integrate([](auto const& x) { return x[0]; }, prob.gridView(), 1); // OK
   * ```
   **/
  template <class Expr, class GridView>
  auto integrate(Expr&& expr, GridView const& gridView, int degree,
                 Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
  {
    using Rules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;

    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
    return Impl::integrateImpl(FWD(gridFct), gridView,
      [&](auto const& type, auto&&) { return Rules::rule(type, degree, qt); });
  }

} // end namespace AMDiS