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 (localFctOrder >= 0)
        return localFctOrder;
      else 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