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

3
4
#include <type_traits>
#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 GV, class QP>
    auto integrateImpl(GF&& gf, GV const& gridView, QP makeQuad)
13
    {
14
      auto localFct = localFunction(FWD(gf));
15

16
      using GridFct = remove_cvref_t<GF>;
17
      using Range = typename GridFct::Range;
18
19
20
21
22
23
24
25
26
27
28
29
30
      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();
      }

Praetorius, Simon's avatar
Praetorius, Simon committed
31
      return gridView.comm().sum(result);
32
33
    }

34
35
    template <class GF, class GV, class QP>
    auto integrateImpl(GF&& gf, GV const& gv, QP makeQuad, std::true_type)
36
    {
37
      return integrateImpl(FWD(gf), gv, makeQuad);
38
39
40
41
42
    }

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

43
44
45
  } // end namespace Impl


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

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

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

69
    return Impl::integrateImpl(FWD(gridFct), gridView, makeQuad, bool_<expr_has_order>);
70
  }
71

72
73

  /// \brief Integrate expression with quadrature rule provided
74
75
76
77
78
79
80
  /**
   * **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
   * ```
   **/
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)
  {
85
86
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
    return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto&&, auto&&) { return quad; });
87
  }
88

89
90

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

103
104
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
    return Impl::integrateImpl(FWD(gridFct), gridView,
105
      [&](auto const& type, auto&&) { return Rules::rule(type, degree, qt); });
106
107
108
  }

} // end namespace AMDiS