[common] Add function to integrate FV-style solutions

dumux/common/integrate.hh

@@ -28,17 +28,117 @@
 #include <type_traits>
 
 #include <dune/geometry/quadraturerules.hh>
+#include <dune/common/concept.hh>
+
+#if HAVE_DUNE_FUNCTIONS
+#include <dune/functions/gridfunctions/gridfunction.hh>
+#endif
+
+#include <dumux/discretization/evalsolution.hh>
+#include <dumux/discretization/elementsolution.hh>
 
 namespace Dumux {
+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>(); }
+
+} // end namespace Impl
+
+/*!
+ * \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 GridGeometry, class SolutionVector,
+         typename std::enable_if_t<!Impl::hasLocalFunction<SolutionVector>(), int> = 0>
+auto integrateGridFunction(const GridGeometry& gg,
+                           SolutionVector&& sol,
+                           std::size_t order)
+{
+    using Scalar = std::decay_t<decltype(sol[0][0])>;
+    using GridView = typename GridGeometry::GridView;
+
+    Scalar integral = 0.0;
+    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
+ * \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 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 Scalar = std::decay_t<decltype(sol1[0][0])>;
+    using GridView = typename GridGeometry::GridView;
+
+    Scalar l2norm = 0.0;
+    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);
+            l2norm += error*error*qp.weight()*geometry.integrationElement(qp.position());
+        }
+    }
+    using std::sqrt;
+    return sqrt(l2norm);
+}
+
+#if HAVE_DUNE_FUNCTIONS
 
 /*!
  * \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
+ * \note overload for a Dune::Funtions::GridFunction
  */
-template<class GridView, class F>
-auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
+template<class GridView, class F,
+         typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
+auto integrateGridFunction(const GridView& gv,
+                           F&& f,
+                           std::size_t order)
 {
     auto fLocal = localFunction(f);
 
@@ -54,7 +154,11 @@ auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
         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());
+        {
+            auto value = fLocal(qp.position());
+            value *= qp.weight()*geometry.integrationElement(qp.position());
+            integral += value;
+        }
 
         fLocal.unbind();
     }
@@ -67,10 +171,15 @@ auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
  * \param f the first function
  * \param g the second function
  * \param order the order of the quadrature rule
+ * \note overload for a Dune::Funtions::GridFunction
  * \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)
+template<class GridView, class F, class G,
+         typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
+auto integrateL2Error(const GridView& gv,
+                      F&& f,
+                      G&& g,
+                      std::size_t order)
 {
     auto fLocal = localFunction(f);
     auto gLocal = localFunction(g);
@@ -102,6 +211,7 @@ auto integrateL2Error(const GridView& gv, F&& f, G&& g, std::size_t order)
     using std::sqrt;
     return sqrt(l2norm);
 }
+#endif
 
 } // end namespace Dumux
 

test/common/functions/test_function_l2norm.cc

@@ -12,7 +12,8 @@
 #include <dumux/common/integrate.hh>
 #include <dumux/multidomain/glue.hh>
 #include <dumux/discretization/projection/projector.hh>
-#include <dumux/discretization/fem/fegridgeometry.hh>
+#include <dumux/discretization/basegridgeometry.hh>
+#include <dumux/discretization/box/fvgridgeometry.hh>
 #include <dumux/io/grid/gridmanager_yasp.hh>
 
 int main (int argc, char *argv[]) try
@@ -35,38 +36,58 @@ int main (int argc, char *argv[]) try
     const auto gvFine = gridManagerFine.grid().leafGridView();
     const auto gvCoarse = gridManagerCoarse.grid().leafGridView();
 
-    // a quadratic function on the grid
-    auto f = [](const auto& pos){ return pos.two_norm2(); };
-
-    // make discrete functions
+    // get some types
     using R = Dune::FieldVector<double, 1>;
     using SolutionVector = Dune::BlockVector<R>;
     using GridView = Grid::LeafGridView;
     using namespace Dune::Functions;
     using P2 = LagrangeBasis<GridView, 2>;
-    auto basisCoarse = std::make_shared<P2>(gvCoarse), basisFine = std::make_shared<P2>(gvFine);
 
+    // construct glue, any grid geometry should be fine here
+    Dumux::BaseGridGeometry<GridView, Dumux::DefaultMapperTraits<GridView>> ggCoarse(gvCoarse), ggFine(gvFine);
+    auto glue = makeGlue(ggFine, ggCoarse);
+
+    // a quadratic function on the grid
+    auto f = [](const auto& pos){ return pos.two_norm2(); };
+
+    // make bases
+    P2 basisCoarse(gvCoarse), basisFine(gvFine);
     SolutionVector solCoarse, solFine;
-    interpolate(*basisCoarse, solCoarse, f);
-    interpolate(*basisFine, solFine, f);
-    auto gfCoarse = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse);
-    auto gfFine = makeDiscreteGlobalBasisFunction<R>(*basisFine, solFine);
+    interpolate(basisCoarse, solCoarse, f);
+    interpolate(basisFine, solFine, f);
+    auto gfCoarse = makeDiscreteGlobalBasisFunction<R>(basisCoarse, solCoarse);
+    auto gfFine = makeDiscreteGlobalBasisFunction<R>(basisFine, solFine);
 
-    // project fine discrete function onto coarse grid
-    Dumux::FEGridGeometry<P2> ggCoarse(basisCoarse), ggFine(basisFine);
-    const auto projector = Dumux::makeProjector(*basisFine, *basisCoarse, makeGlue(ggFine, ggCoarse));
+    // project fine discrete function onto coarse grid and convert back to grid function
+    const auto projector = Dumux::makeProjector(basisFine, basisCoarse, glue);
     SolutionVector solCoarse2;
     projector.project(solFine, solCoarse2);
-    auto gfCoarse2 = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse);
+    auto gfCoarse2 = makeDiscreteGlobalBasisFunction<R>(basisCoarse, solCoarse);
 
     // check that these functions are identical and that they exactly represent the analytical function
     auto l2norm = Dumux::integrateL2Error(gvCoarse, gfCoarse, gfCoarse2, 3);
     if (l2norm > 1e-14)
-        DUNE_THROW(Dune::MathError, "L2 norm of projection and interpolation is non-zero!");
+        DUNE_THROW(Dune::MathError, "L2 norm of projection and interpolation is non-zero (" << l2norm << ")");
 
     l2norm = Dumux::integrateL2Error(gvCoarse, makeAnalyticGridViewFunction(f, gvCoarse), gfCoarse2, 3);
     if (l2norm > 1e-14)
-        DUNE_THROW(Dune::MathError, "L2 norm of analytical solution and interpolation is non-zero!");
+        DUNE_THROW(Dune::MathError, "L2 norm of analytical solution and interpolation is non-zero (" << l2norm << ")");
+
+    // check the FV-style interface of integrateL2Error
+    Dumux::BoxFVGridGeometry<double, GridView> ggBoxCoarse(gvCoarse), ggBoxFine(gvFine);
+    auto p1BasisCoarse = Dumux::getFunctionSpaceBasis(ggBoxCoarse);
+    auto p1BasisFine = Dumux::getFunctionSpaceBasis(ggBoxFine);
+    auto f2 = [](const auto& pos){ return pos[0]; };
+    interpolate(p1BasisFine, solFine, f2);
+    interpolate(p1BasisCoarse, solCoarse, f2);
+
+    // project fine discrete function onto coarse grid
+    const auto projector2 = Dumux::makeProjector(p1BasisFine, p1BasisCoarse, glue);
+    projector2.project(solFine, solCoarse2);
+
+    l2norm = Dumux::integrateL2Error(ggBoxCoarse, solCoarse, solCoarse2, 3);
+    if (l2norm > 1e-14)
+        DUNE_THROW(Dune::MathError, "L2 norm of FV-style projected solutions is non-zero (" << l2norm << ")");
 
     std::cout << "All tests passed!" << std::endl;
     return 0;