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

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

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

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

35
36
37
38
39
40
41
42
43
    template <class GF, class GridView, class QuadProvider>
    auto integrateImpl(GF&& gf, GridView const& gv, QuadProvider makeQuad, std::true_type)
    {
      return integrateImpl(std::forward<GF>(gf), gv, makeQuad);
    }

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

44
45
46
  } // end namespace Impl


47
48
49
50
51
52
53
54
55
  /// \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
   * ```
   **/
56
57
58
  template <class Expr, class GridView>
  auto integrate(Expr&& expr, GridView const& gridView)
  {
59
    auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), gridView);
60
61
62
63

    // test whether the gridFct model `Concepts::HasLocalFunctionOrder`
    using GF = std::decay_t<decltype(gridFct)>;
    static const bool expr_has_order = Concepts::HasLocalFunctionOrder<GF>;
64
65
    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()`.");
66

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

    return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView, makeQuad,
      bool_<expr_has_order>);
72
  }
73

74
75

  /// \brief Integrate expression with quadrature rule provided
76
77
78
79
80
81
82
  /**
   * **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
   * ```
   **/
83
84
85
86
87
88
89
90
  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; });
  }
91

92
93

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

106
107
    auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), gridView);
    return Impl::integrateImpl(std::forward<decltype(gridFct)>(gridFct), gridView,
108
      [&](auto const& type, auto&&) { return Rules::rule(type, degree, qt); });
109
110
111
  }

} // end namespace AMDiS