Integrate.hpp

#pragma once

#include <type_traits>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/common/partitionset.hh>

#include <amdis/GridFunctions.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Order.hpp>

namespace AMDiS
{

  template <class Expr, class GV>
  auto integrate(Expr&& expr, GV const& gridView, int localFctOrder = -1)
  {
    auto gridFct = makeGridFunction(FWD(expr), gridView);
    auto localFct = localFunction(gridFct);

    using LF = TYPEOF(localFct);
    auto quadOrder = [localFctOrder](auto const& lf) {
      if constexpr(Concepts::Polynomial<LF>)
        return order(lf);
      else
        return localFctOrder;
    };

    test_exit(Concepts::Polynomial<LF> || localFctOrder >= 0, "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 Range = typename decltype(localFct)::Range;
    Range result(0);

    using Rules = Dune::QuadratureRules<typename GV::ctype, GV::dimension>;
    for (auto const& element : elements(gridView, Dune::Partitions::interior)) {
      auto geometry = element.geometry();

      localFct.bind(element);
      auto const& quad = Rules::rule(element.type(), quadOrder(localFct));

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

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

} // end namespace AMDiS