3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
integrate.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_COMMON_INTEGRATE_HH
25#define DUMUX_COMMON_INTEGRATE_HH
26
27#include <cmath>
28#include <type_traits>
29
30#include <dune/geometry/quadraturerules.hh>
31#include <dune/common/concept.hh>
32
33#if HAVE_DUNE_FUNCTIONS
34#include <dune/functions/gridfunctions/gridfunction.hh>
35#endif
36
41
42namespace Dumux {
43
44// implementation details of the integrate functions
45#ifndef DOXYGEN
46namespace Detail {
47
48struct HasLocalFunction
49{
50 template<class F>
51 auto require(F&& f) -> decltype(
52 localFunction(f),
53 localFunction(f).unbind()
54 );
55};
56
57template<class F>
58static constexpr bool hasLocalFunction()
59{ return Dune::models<HasLocalFunction, F>(); }
60
61template<class Error,
62 typename std::enable_if_t<IsIndexable<Error>::value, int> = 0>
63Error sqrtNorm(const Error& error)
64{
65 using std::sqrt;
66 auto e = error;
67 for (int i = 0; i < error.size(); ++i)
68 e[i] = sqrt(error[i]);
69 return e;
70}
71
72template<class Error,
73 typename std::enable_if_t<!IsIndexable<Error>::value, int> = 0>
74Error sqrtNorm(const Error& error)
75{
76 using std::sqrt;
77 return sqrt(error);
78}
79
80template<class T, typename = int>
81struct FieldTypeImpl
82{
83 using type = T;
84};
85
86template<class T>
87struct FieldTypeImpl<T, typename std::enable_if<(sizeof(std::declval<T>()[0]) > 0), int>::type>
88{
89 using type = typename FieldTypeImpl<std::decay_t<decltype(std::declval<T>()[0])>>::type;
90};
91
92template<class T>
93using FieldType = typename FieldTypeImpl<T>::type;
94
95} // end namespace Detail
96#endif
97
104template<class GridGeometry, class SolutionVector,
105 typename std::enable_if_t<!Detail::hasLocalFunction<SolutionVector>(), int> = 0>
106auto integrateGridFunction(const GridGeometry& gg,
107 const SolutionVector& sol,
108 std::size_t order)
109{
110 using GridView = typename GridGeometry::GridView;
111 using Scalar = typename Detail::FieldType< std::decay_t<decltype(sol[0])> >;
112
113 Scalar integral(0.0);
114 for (const auto& element : elements(gg.gridView()))
115 {
116 const auto elemSol = elementSolution(element, sol, gg);
117 const auto geometry = element.geometry();
118 const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
119 for (auto&& qp : quad)
120 {
121 auto value = evalSolution(element, geometry, gg, elemSol, geometry.global(qp.position()));
122 value *= qp.weight()*geometry.integrationElement(qp.position());
123 integral += value;
124 }
125 }
126 return integral;
127}
128
137template<class GridGeometry, class Sol1, class Sol2,
138 typename std::enable_if_t<!Detail::hasLocalFunction<Sol1>(), int> = 0>
139auto integrateL2Error(const GridGeometry& gg,
140 const Sol1& sol1,
141 const Sol2& sol2,
142 std::size_t order)
143{
144 using GridView = typename GridGeometry::GridView;
145 using Scalar = typename Detail::FieldType< std::decay_t<decltype(sol1[0])> >;
146
147 Scalar l2norm(0.0);
148 for (const auto& element : elements(gg.gridView()))
149 {
150 const auto elemSol1 = elementSolution(element, sol1, gg);
151 const auto elemSol2 = elementSolution(element, sol2, gg);
152
153 const auto geometry = element.geometry();
154 const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
155 for (auto&& qp : quad)
156 {
157 const auto& globalPos = geometry.global(qp.position());
158 const auto value1 = evalSolution(element, geometry, gg, elemSol1, globalPos);
159 const auto value2 = evalSolution(element, geometry, gg, elemSol2, globalPos);
160 const auto error = (value1 - value2);
161 l2norm += (error*error)*qp.weight()*geometry.integrationElement(qp.position());
162 }
163 }
164
165 using std::sqrt;
166 return sqrt(l2norm);
167}
168
169#if HAVE_DUNE_FUNCTIONS
170
178template<class GridView, class F,
179 typename std::enable_if_t<Detail::hasLocalFunction<F>(), int> = 0>
180auto integrateGridFunction(const GridView& gv,
181 const F& f,
182 std::size_t order)
183{
184 auto fLocal = localFunction(f);
185
186 using Element = typename GridView::template Codim<0>::Entity;
187 using LocalPosition = typename Element::Geometry::LocalCoordinate;
188 using Scalar = typename Detail::FieldType< std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))> >;
189
190 Scalar integral(0.0);
191 for (const auto& element : elements(gv))
192 {
193 fLocal.bind(element);
194
195 const auto geometry = element.geometry();
196 const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
197 for (auto&& qp : quad)
198 {
199 auto value = fLocal(qp.position());
200 value *= qp.weight()*geometry.integrationElement(qp.position());
201 integral += value;
202 }
203
204 fLocal.unbind();
205 }
206 return integral;
207}
208
218template<class GridView, class F, class G,
219 typename std::enable_if_t<Detail::hasLocalFunction<F>(), int> = 0>
220auto integrateL2Error(const GridView& gv,
221 const F& f,
222 const G& g,
223 std::size_t order)
224{
225 auto fLocal = localFunction(f);
226 auto gLocal = localFunction(g);
227
228 using Element = typename GridView::template Codim<0>::Entity;
229 using LocalPosition = typename Element::Geometry::LocalCoordinate;
230 using Scalar = typename Detail::FieldType< std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))> >;
231
232 Scalar l2norm(0.0);
233 for (const auto& element : elements(gv))
234 {
235 fLocal.bind(element);
236 gLocal.bind(element);
237
238 const auto geometry = element.geometry();
239 const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
240 for (auto&& qp : quad)
241 {
242 const auto error = fLocal(qp.position()) - gLocal(qp.position());
243 l2norm += (error*error)*qp.weight()*geometry.integrationElement(qp.position());
244 }
245
246 gLocal.unbind();
247 fLocal.unbind();
248 }
249
250 using std::sqrt;
251 return sqrt(l2norm);
252}
253#endif
254
263template<class Scalar, class Function,
264 typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
265Scalar integrateScalarFunction(const Function& f,
266 const Scalar lowerBound,
267 const Scalar upperBound,
268 const Scalar targetAbsoluteError = 1e-13)
269{
270 return DoubleExponentialIntegrator<Scalar>::integrate(f, lowerBound, upperBound, targetAbsoluteError);
271}
272
273} // end namespace Dumux
274
275#endif
A double exponential integrator.
Element solution classes and factory functions.
free functions for the evaluation of primary variables inside elements.
PrimaryVariables evalSolution(const Element &element, const typename Element::Geometry &geometry, const typename FVElementGeometry::GridGeometry &gridGeometry, const BoxElementSolution< FVElementGeometry, PrimaryVariables > &elemSol, const typename Element::Geometry::GlobalCoordinate &globalPos, bool ignoreState=false)
Interpolates a given box element solution at a given global position. Uses the finite element cache o...
Definition: evalsolution.hh:94
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:115
Definition: adapt.hh:29
Scalar integrateScalarFunction(const Function &f, const Scalar lowerBound, const Scalar upperBound, const Scalar targetAbsoluteError=1e-13)
Integrate a scalar function.
Definition: integrate.hh:265
auto integrateGridFunction(const GridGeometry &gg, const SolutionVector &sol, std::size_t order)
Integrate a grid function over a grid view.
Definition: integrate.hh:106
auto integrateL2Error(const GridGeometry &gg, const Sol1 &sol1, const Sol2 &sol2, std::size_t order)
Integrate a function over a grid view.
Definition: integrate.hh:139
static Scalar integrate(const Function &f, const Scalar a, const Scalar b, const Scalar targetAbsoluteError, int &numFunctionEvaluations, Scalar &errorEstimate)
Integrate an analytic function over a finite interval.
Definition: doubleexpintegrator.hh:86