Integrate.hpp 3.77 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
8
9

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

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

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

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

    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(FWD(expr), gridView);
60
61

    // test whether the gridFct model `Concepts::HasLocalFunctionOrder`
62
    using GF = TYPEOF(gridFct);
63
    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
    using Rules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
    auto makeQuad = [](auto&& t, auto&& lf) { return Rules::rule(t, order(lf)); };

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

73
74

  /// \brief Integrate expression with quadrature rule provided
75
76
77
78
79
80
81
  /**
   * **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
   * ```
   **/
82
83
84
85
  template <class Expr, class GridView,
    class QuadratureRule = Dune::QuadratureRule<typename GridView::ctype, GridView::dimension>>
  auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad)
  {
86
87
    auto&& gridFct = makeGridFunction(FWD(expr), gridView);
    return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto&&, auto&&) { return quad; });
88
  }
89

90
91

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

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

} // end namespace AMDiS