Integrate.hpp 3.79 KB
Newer Older
1
#pragma once
2
#include "config.h"
3
4
#include <dune/geometry/quadraturerules.hh>

5
#include <amdis/GridFunctions.hpp>
6
7
8

namespace AMDiS
{
9
10
  namespace Impl
  {
11
12
    template <class GF, class GridView, class QuadProvider>
    auto integrateImpl(GF&& gf, GridView const& gridView, QuadProvider makeQuad)
13
    {
14
      auto localFct = localFunction(std::forward<GF>(gf));
15

16
17
      using GridFct = std::decay_t<GF>;
      using Range = typename GridFct::Range;
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
      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 result;
    }

  } // end namespace Impl


37
38
39
40
41
42
43
44
45
  /// \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
   * ```
   **/
46
47
48
  template <class Expr, class GridView>
  auto integrate(Expr&& expr, GridView const& gridView)
  {
49
50
51
52
    auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), gridView);
    static const bool expr_has_order = Concepts::HasLocalFunctionOrder<std::decay_t<decltype(gridFct)>>;
    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()`.");
53
54

    using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
55
    auto makeQuad = [](auto&& type, auto&& localFct) { return QuadratureRules::rule(type, order(localFct)); };
56

57
#if AMDIS_HAS_CXX_CONSTEXPR_IF
58
    if constexpr(expr_has_order)
59
      return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView, makeQuad);
60
61
62
63
    else
      return 0.0;
#else
    return Dune::Hybrid::ifElse(bool_<expr_has_order>,
64
      [&](auto) { return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView, makeQuad); },
65
66
67
      [ ](auto) { return 0.0; });
#endif
  }
68

69
  /// Integrate expression with quadrature rule provided
70
71
72
73
74
75
76
  /**
   * **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
   * ```
   **/
77
78
79
80
81
82
83
84
  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(std::forward<Expr>(expr), gridView);
    return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView,
      [&](auto&&, auto&&) { return quad; });
  }
85

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

99
100
    auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), gridView);
    return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView,
101
      [&](auto const& type, auto&&) { return QuadratureRules::rule(type, degree, qt); });
102
103
104
  }

} // end namespace AMDiS