Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

Integrate.hpp 3.43 KB
Newer Older
1 2
#pragma once

3 4
#include <type_traits>
#include <dune/geometry/quadraturerules.hh>
5
#include <dune/grid/common/partitionset.hh>
6
#include <amdis/GridFunctions.hpp>
7
#include <amdis/common/Order.hpp>
8 9 10

namespace AMDiS
{
11 12
  namespace Impl
  {
13 14
    template <class GF, class GV, class QP>
    auto integrateImpl(GF&& gf, GV const& gridView, QP makeQuad)
15
    {
16
      auto localFct = localFunction(FWD(gf));
17

18
      using GridFct = remove_cvref_t<GF>;
19
      using Range = typename GridFct::Range;
20 21
      Range result(0);

22
      for (auto const& element : elements(gridView, Dune::Partitions::interior)) {
23 24 25 26 27 28 29 30 31 32
        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();
      }

Praetorius, Simon's avatar
Praetorius, Simon committed
33
      return gridView.comm().sum(result);
34 35 36 37 38
    }

  } // end namespace Impl


39 40 41 42 43 44 45 46 47
  /// \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
   * ```
   **/
48 49 50
  template <class Expr, class GridView>
  auto integrate(Expr&& expr, GridView const& gridView)
  {
51
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
52

53 54 55
    using LocalFct = TYPEOF(localFunction(gridFct));
    static_assert(Concepts::Polynomial<LocalFct>, "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()`.");

56
    using Rules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
57 58 59 60 61
    if constexpr(Concepts::Polynomial<LocalFct>)
      return Impl::integrateImpl(FWD(gridFct), gridView,
        [](auto&& t, auto&& lf) { return Rules::rule(t, order(lf)); });
    else
      return 0.0;
62
  }
63

64 65

  /// \brief Integrate expression with quadrature rule provided
66 67 68 69 70 71 72
  /**
   * **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
   * ```
   **/
73 74 75 76
  template <class Expr, class GridView,
    class QuadratureRule = Dune::QuadratureRule<typename GridView::ctype, GridView::dimension>>
  auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad)
  {
77
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
78 79
    return Impl::integrateImpl(FWD(gridFct), gridView,
      [&](auto&&, auto&&) { return quad; });
80
  }
81

82 83

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

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

} // end namespace AMDiS