integrate.hh 9.39 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup Common
 * \brief Define helper functions for integration
 */
#ifndef DUMUX_COMMON_INTEGRATE_HH
#define DUMUX_COMMON_INTEGRATE_HH

#include <cmath>
#include <type_traits>

#include <dune/geometry/quadraturerules.hh>
31
#include <dune/common/concept.hh>
32
#include <dune/common/std/type_traits.hh>
33
34
35
36
37

#if HAVE_DUNE_FUNCTIONS
#include <dune/functions/gridfunctions/gridfunction.hh>
#endif

38
#include <dumux/common/doubleexpintegrator.hh>
39
40
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/elementsolution.hh>
41
#include <dumux/common/typetraits/typetraits.hh>
42
43

namespace Dumux {
44
45
46

// implementation details of the integrate functions
#ifndef DOXYGEN
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
namespace Impl {

struct HasLocalFunction
{
    template<class F>
    auto require(F&& f) -> decltype(
      localFunction(f),
      localFunction(f).unbind()
    );
};

template<class F>
static constexpr bool hasLocalFunction()
{ return Dune::models<HasLocalFunction, F>(); }

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
template<class Error,
         typename std::enable_if_t<IsIndexable<Error>::value, int> = 0>
Error sqrtNorm(const Error& error)
{
    using std::sqrt;
    auto e = error;
    for (int i = 0; i < error.size(); ++i)
        e[i] = sqrt(error[i]);
    return e;
}

template<class Error,
         typename std::enable_if_t<!IsIndexable<Error>::value, int> = 0>
Error sqrtNorm(const Error& error)
{
    using std::sqrt;
    return sqrt(error);
}

template<class T, typename = int>
struct FieldTypeImpl
{
    using type = T;
};

template<class T>
struct FieldTypeImpl<T, typename std::enable_if<(sizeof(std::declval<T>()[0]) > 0), int>::type>
{
    using type = typename FieldTypeImpl<std::decay_t<decltype(std::declval<T>()[0])>>::type;
};

template<class T>
using FieldType = typename FieldTypeImpl<T>::type;

96
} // end namespace Impl
97
#endif
98
99
100

/*!
 * \brief Integrate a grid function over a grid view
Dennis Gläser's avatar
Dennis Gläser committed
101
 * \param gg the grid geometry
102
 * \param sol the solution vector
103
104
105
106
107
 * \param order the order of the quadrature rule
 */
template<class GridGeometry, class SolutionVector,
         typename std::enable_if_t<!Impl::hasLocalFunction<SolutionVector>(), int> = 0>
auto integrateGridFunction(const GridGeometry& gg,
108
                           const SolutionVector& sol,
109
110
111
                           std::size_t order)
{
    using GridView = typename GridGeometry::GridView;
112
    using Scalar = typename Impl::FieldType< std::decay_t<decltype(sol[0])> >;
113

114
    Scalar integral(0.0);
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
    for (const auto& element : elements(gg.gridView()))
    {
        const auto elemSol = elementSolution(element, sol, gg);
        const auto geometry = element.geometry();
        const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
        {
            auto value = evalSolution(element, geometry, gg, elemSol, geometry.global(qp.position()));
            value *= qp.weight()*geometry.integrationElement(qp.position());
            integral += value;
        }
    }
    return integral;
}

/*!
 * \brief Integrate a function over a grid view
132
133
134
 * \param gg the grid geometry
 * \param sol1 the first function
 * \param sol2 the second function
135
136
137
138
139
140
141
142
143
144
145
 * \param order the order of the quadrature rule
 * \note dune functions currently doesn't support composing two functions
 */
template<class GridGeometry, class Sol1, class Sol2,
         typename std::enable_if_t<!Impl::hasLocalFunction<Sol1>(), int> = 0>
auto integrateL2Error(const GridGeometry& gg,
                      const Sol1& sol1,
                      const Sol2& sol2,
                      std::size_t order)
{
    using GridView = typename GridGeometry::GridView;
146
    using Scalar = typename Impl::FieldType< std::decay_t<decltype(sol1[0])> >;
147

148
    Scalar l2norm(0.0);
149
150
151
152
153
154
155
156
157
158
159
160
161
    for (const auto& element : elements(gg.gridView()))
    {
        const auto elemSol1 = elementSolution(element, sol1, gg);
        const auto elemSol2 = elementSolution(element, sol2, gg);

        const auto geometry = element.geometry();
        const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
        {
            const auto& globalPos = geometry.global(qp.position());
            const auto value1 = evalSolution(element, geometry, gg, elemSol1, globalPos);
            const auto value2 = evalSolution(element, geometry, gg, elemSol2, globalPos);
            const auto error = (value1 - value2);
162
            l2norm += (error*error)*qp.weight()*geometry.integrationElement(qp.position());
163
164
        }
    }
165

166
167
    using std::sqrt;
    return sqrt(l2norm);
168
169
170
}

#if HAVE_DUNE_FUNCTIONS
171
172
173
174
175
176

/*!
 * \brief Integrate a grid function over a grid view
 * \param gv the grid view
 * \param f the grid function
 * \param order the order of the quadrature rule
177
 * \note overload for a Dune::Funtions::GridFunction
178
 */
179
180
181
template<class GridView, class F,
         typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
auto integrateGridFunction(const GridView& gv,
182
                           const F& f,
183
                           std::size_t order)
184
185
186
187
188
{
    auto fLocal = localFunction(f);

    using Element = typename GridView::template Codim<0>::Entity;
    using LocalPosition = typename Element::Geometry::LocalCoordinate;
189
    using Scalar = typename Impl::FieldType< std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))> >;
190

191
    Scalar integral(0.0);
192
193
194
195
196
197
198
    for (const auto& element : elements(gv))
    {
        fLocal.bind(element);

        const auto geometry = element.geometry();
        const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
199
200
201
202
203
        {
            auto value = fLocal(qp.position());
            value *= qp.weight()*geometry.integrationElement(qp.position());
            integral += value;
        }
204
205
206
207
208
209
210
211
212
213
214
215

        fLocal.unbind();
    }
    return integral;
}

/*!
 * \brief Integrate a function over a grid view
 * \param gv the grid view
 * \param f the first function
 * \param g the second function
 * \param order the order of the quadrature rule
216
 * \note overload for a Dune::Funtions::GridFunction
217
218
 * \note dune functions currently doesn't support composing two functions
 */
219
220
221
template<class GridView, class F, class G,
         typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
auto integrateL2Error(const GridView& gv,
222
223
                      const F& f,
                      const G& g,
224
                      std::size_t order)
225
226
227
228
229
230
{
    auto fLocal = localFunction(f);
    auto gLocal = localFunction(g);

    using Element = typename GridView::template Codim<0>::Entity;
    using LocalPosition = typename Element::Geometry::LocalCoordinate;
231
    using Scalar = typename Impl::FieldType< std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))> >;
232

233
    Scalar l2norm(0.0);
234
235
236
237
238
239
240
241
242
243
    for (const auto& element : elements(gv))
    {
        fLocal.bind(element);
        gLocal.bind(element);

        const auto geometry = element.geometry();
        const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
        {
            const auto error = fLocal(qp.position()) - gLocal(qp.position());
244
            l2norm += (error*error)*qp.weight()*geometry.integrationElement(qp.position());
245
246
247
248
249
        }

        gLocal.unbind();
        fLocal.unbind();
    }
250

251
252
    using std::sqrt;
    return sqrt(l2norm);
253
}
254
#endif
255

256
257
258
/*!
 * \brief Integrate a scalar function
 * \param f the integrand (invocable with a single scalar)
Melanie Lipp's avatar
Melanie Lipp committed
259
260
 * \param lowerBound lower integral bound
 * \param upperBound upper integral bound
261
262
263
264
265
266
267
268
269
270
271
272
273
 * \param targetAbsoluteError desired absolute error in the result
 * \return The value of the integral
 */
template<class Scalar, class Function,
          typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::value>...>
Scalar integrateScalarFunction(const Function& f,
                               const Scalar lowerBound,
                               const Scalar upperBound,
                               const Scalar targetAbsoluteError = 1e-13)
{
    return DoubleExponentialIntegrator<Scalar>::integrate(f, lowerBound, upperBound, targetAbsoluteError);
}

274
275
276
} // end namespace Dumux

#endif