integrate.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_COMMON_INTEGRATE_HH
#define DUMUX_COMMON_INTEGRATE_HH
 
#include <cmath>
#include <type_traits>
 
#include <dune/common/typetraits.hh>
#include <dune/common/ftraits.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/common/concept.hh>
 
#if HAVE_DUNE_FUNCTIONS
#include <dune/functions/gridfunctions/gridfunction.hh>
#endif
 
#include <dumux/common/doubleexpintegrator.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/elementsolution.hh>
 
namespace Dumux {
 
// implementation details of the integrate functions
#ifndef DOXYGEN
namespace Detail {
 
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>(); }
 
template<class Error,
         typename std::enable_if_t<Dune::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<!Dune::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;
 
// a scalar integrand, which may be given as a vector with a single entry
template<class T>
auto scalarIntegrand(const T& t)
{
    if constexpr (Dune::IsNumber<T>::value)
        return t;
    else
    {
        static_assert(T::size() == 1, "The integrand has to be a scalar or a vector with a single entry");
        return t[0];
    }
}
 
// squared Euclidean norm, real-valued also for complex-valued entries
template<class T>
auto squaredNorm(const T& t)
{
    if constexpr (Dune::IsNumber<T>::value)
    {
        using std::abs;
        const auto a = abs(t);
        return a*a;
    }
    else
        return t.two_norm2();
}
 
} // end namespace Detail
#endif

template<class GridGeometry, class SolutionVector,
         typename std::enable_if_t<!Detail::hasLocalFunction<SolutionVector>(), int> = 0>
auto integrateGridFunction(const GridGeometry& gg,
                           const SolutionVector& sol,
                           std::size_t order)
{
    using GridView = typename GridGeometry::GridView;
    using Scalar = typename Detail::FieldType< std::decay_t<decltype(sol[0])> >;
 
    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<typename GridView::ctype, 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 += Detail::scalarIntegrand(value);
        }
    }
    return integral;
}

template<class GridGeometry, class Sol1, class Sol2,
         typename std::enable_if_t<!Detail::hasLocalFunction<Sol1>(), int> = 0>
auto integrateL2Error(const GridGeometry& gg,
                      const Sol1& sol1,
                      const Sol2& sol2,
                      std::size_t order)
{
    using GridView = typename GridGeometry::GridView;
    using Scalar = typename Dune::FieldTraits<Detail::FieldType< std::decay_t<decltype(sol1[0])> >>::real_type;
 
    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<typename GridView::ctype, 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 += Detail::squaredNorm(error)*qp.weight()*geometry.integrationElement(qp.position());
        }
    }
 
    using std::sqrt;
    return sqrt(l2norm);
}
 
#if HAVE_DUNE_FUNCTIONS

template<class GridView, class F,
         typename std::enable_if_t<Detail::hasLocalFunction<F>(), int> = 0>
auto integrateGridFunction(const GridView& gv,
                           const 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 = typename Detail::FieldType< 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<typename GridView::ctype, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
        {
            auto value = fLocal(qp.position());
            value *= qp.weight()*geometry.integrationElement(qp.position());
            integral += Detail::scalarIntegrand(value);
        }
 
        fLocal.unbind();
    }
    return integral;
}

template<class GridView, class F, class G,
         typename std::enable_if_t<Detail::hasLocalFunction<F>(), int> = 0>
auto integrateL2Error(const GridView& gv,
                      const F& f,
                      const 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 Scalar = typename Dune::FieldTraits<Detail::FieldType< std::decay_t<decltype(fLocal(std::declval<LocalPosition>()))> >>::real_type;
 
    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<typename GridView::ctype, GridView::dimension>::rule(geometry.type(), order);
        for (auto&& qp : quad)
        {
            const auto error = fLocal(qp.position()) - gLocal(qp.position());
            l2norm += Detail::squaredNorm(error)*qp.weight()*geometry.integrationElement(qp.position());
        }
 
        gLocal.unbind();
        fLocal.unbind();
    }
 
    using std::sqrt;
    return sqrt(l2norm);
}
#endif

template<class Scalar, class Function,
          typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
Scalar integrateScalarFunction(const Function& f,
                               const Scalar lowerBound,
                               const Scalar upperBound,
                               const Scalar targetAbsoluteError = 1e-13)
{
    return DoubleExponentialIntegrator<Scalar>::integrate(f, lowerBound, upperBound, targetAbsoluteError);
}
 
} // end namespace Dumux
 
#endif