From 343b3c1ae22bb9082b7fde19037003baab98f4d8 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 31 May 2019 20:37:20 +0200
Subject: [PATCH] [common] Implement function to integrate grid functions
We assume the grid function to model a Dune::Functions::GridFunction.
---
dumux/common/integrate.hh | 108 ++++++++++++++++++++++++++++++++++++++
1 file changed, 108 insertions(+)
create mode 100644 dumux/common/integrate.hh
diff --git a/dumux/common/integrate.hh b/dumux/common/integrate.hh
new file mode 100644
index 0000000000..9ea524a6da
--- /dev/null
+++ b/dumux/common/integrate.hh
@@ -0,0 +1,108 @@
+// -*- 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>
+
+namespace Dumux {
+
+/*!
+ * \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
+ */
+template<class GridView, class F>
+auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
+{
+ auto fLocal = localFunction(f);
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using LocalPosition = typename Element::Geometry::LocalCoordinate;
+ using Scalar = std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))>;
+
+ Scalar integral = 0.0;
+ 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)
+ integral += fLocal(qp.position())*qp.weight()*geometry.integrationElement(qp.position());
+
+ 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
+ * \note dune functions currently doesn't support composing two functions
+ */
+template<class GridView, class F, class G>
+auto integrateL2Error(const GridView& gv, F&& f, G&& g, std::size_t order)
+{
+ auto fLocal = localFunction(f);
+ auto gLocal = localFunction(g);
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using LocalPosition = typename Element::Geometry::LocalCoordinate;
+ using FScalar = std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))>;
+ using GScalar = std::decay_t<decltype(gLocal(std::declval<LocalPosition>()))>;
+
+ using Scalar = std::decay_t<decltype((std::declval<FScalar>()-std::declval<GScalar>())
+ *(std::declval<FScalar>()-std::declval<GScalar>()))>;
+ Scalar l2norm = 0.0;
+ 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());
+ l2norm += error*error*qp.weight()*geometry.integrationElement(qp.position());
+ }
+
+ gLocal.unbind();
+ fLocal.unbind();
+ }
+ using std::sqrt;
+ return sqrt(l2norm);
+}
+
+} // end namespace Dumux
+
+#endif
--
GitLab