functionintegrator.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef FUNCTIONINTEGRATOR_HH
#define FUNCTIONINTEGRATOR_HH
 
#include <dune/common/fvector.hh>
 
#include <dune/geometry/quadraturerules.hh>
 
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
class FunctionIntegrator
{
    public:
 
        template <class GridView, class FunctionType>
        static auto integrate(const GridView& gv, const FunctionType& integrand, const int quad_order)
        {
            using GlobalDomain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
            using LocalDomain = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
            using RangeType = std::decay_t<std::invoke_result_t<FunctionType, GlobalDomain>>;
            using FieldType = typename Dune::template FieldTraits<RangeType>::field_type;
 
            RangeType ret(0.0);
            RangeType f_pos(0.0);
 
            auto gridFunctionIntegrand = Dune::Functions::makeGridViewFunction(integrand, gv);
            auto localIntegrand = localFunction(gridFunctionIntegrand);
 
            for (const auto& element : elements(gv))
            {
                auto&& geometry = element.geometry();
 
                localIntegrand.bind(element);
 
                QuadratureRuleKey quadKey(element.type(), quad_order);
 
                const auto& quad = QuadratureRuleCache<FieldType, LocalDomain::dimension>::rule(quadKey);
 
                // loop over quadrature points
                for (size_t pt=0; pt < quad.size(); ++pt)
                {
                    // get quadrature point
                    const auto& quadPos = quad[pt].position();
 
                    // get integration factor
                    const double integrationElement = geometry.integrationElement(quadPos);
 
                    // compute values of function
                    f_pos = localIntegrand(quadPos);
 
                    // add to integral
                    ret.axpy(quad[pt].weight()*integrationElement, f_pos);
                }
            }
 
            return ret;
        }
 
        template <class RangeType, class BasisType, class VectorType>
        static auto integrateDiscreteFunction(const BasisType& basis, const VectorType& integrandCoeff, const int quad_order)
        {
            auto integrand = Dune::Functions::makeDiscreteGlobalBasisFunction<RangeType>(basis, integrandCoeff);
            return integrate(basis.getGridView(), integrand, quad_order);
        }
};
 
 
#endif