functionintegrator.hh

#ifndef FUNCTIONINTEGRATOR_HH
#define FUNCTIONINTEGRATOR_HH

#include <dune/common/fvector.hh>

#include <dune/geometry/quadraturerules.hh>

#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>

/**  \brief provides static methods for numerical integration
  *
  */
class FunctionIntegrator
{
    public:

        /**  \brief numerically integrates a given function
          *
          *  Numerically integrates the given \a integrand over the domain covered by the
          *  given GridView using a Dune::QuadratureRule of order \a quad_order on the grid specified by the GridView.
          *
          *  \tparam GridView the type of the GridView provided
          *  \tparam FunctionType the type of the function provided, needs to provide method <tt>void evalall(...)</tt> and the FunctionType's rangetype has to be RT
          *  \param gv the GridView describing the domain of integration and the Grid for quadrature
          *  \param integrand the function to be integrated
          *  \param quad_order the order of the quadrature rule employed
          */
        template <class GridView, class FunctionType>
        static typename FunctionType::RangeType integrate(const GridView& gv, const FunctionType& integrand, const int quad_order)
        {
            typedef typename FunctionType::DomainType DomainType;
            typedef typename GridView::template Codim<0>::Geometry::LocalCoordinate LocalDomainType;
            typedef typename DomainType::field_type DFT;
            typedef typename FunctionType::RangeType RangeType;

            RangeType ret(0.0);
            RangeType f_pos(0.0);

            typedef typename GridView::template Codim<0>::Iterator ElementIterator;
            typedef typename GridView::template Codim<0>::Entity GridElement;
            typedef typename GridElement::Geometry Geometry;

            ElementIterator eltIt = gv.template begin<0>();
            ElementIterator eltEnd = gv.template end<0>();

            for (; eltIt!=eltEnd; ++eltIt)
            {
//                const Dune::template QuadratureRule<DT, GridView::dimension>& quad = Dune::QuadratureRules<DT, GridView::dimension>::rule(eltIt->type(), quad_order);
                const Dune::template QuadratureRule<DFT, DomainType::dimension>& quad = QuadratureRuleCache<DFT, DomainType::dimension>::rule(eltIt->type(), quad_order, IsRefinedLocalFiniteElementInFunction<FunctionType, GridElement>::value(integrand, *eltIt) );
//                std::cout << IsRefinedLocalFiniteElementInFunction<FunctionType, GridElement>::value(integrand, *eltIt) << std::endl;

                // loop over quadrature points
                for (size_t pt=0; pt < quad.size(); ++pt)
                {
                    // get quadrature point
                    const LocalDomainType& quadPos = quad[pt].position();

                    const Geometry geometry = eltIt->geometry();

                    // get integration factor
                    const double integrationElement = geometry.integrationElement(quadPos);

                    // compute values of function
                    if (dynamic_cast<const VirtualGridViewFunction<GridView,RangeType>* >(&integrand))
                        dynamic_cast<const VirtualGridViewFunction<GridView,RangeType>* >(&integrand)->evaluateLocal(*eltIt, quadPos, f_pos);
                    else
                        integrand.evaluate(geometry.global(quadPos), f_pos);

                    // add to integral
                    ret.axpy(quad[pt].weight()*integrationElement, f_pos);
                }
            }

            return ret;
        }

        /**  \brief numerically integrates a gridfunction given by a coefficient vector
          *
          *  Numerically integrates the gridfunction represented by a coefficient vector in a given basis
          *  using a Dune::QuadratureRule of order quad_order.
          *
          *  \tparam BasisType the functionspace basis for the coefficient representation
          *  \tparam VectorType the type of the coefficient vector
          *  \param basis the functionspace basis of the gridfunction
          *  \param integrand the coefficient vector
          *  \param quad_order the order of the quadrature rule employed
          */
        template <class BasisType, class VectorType>
        static typename VectorType::block_type integrateDiscreteFunction(const BasisType& basis, const VectorType& integrand, const int quad_order)
        {
            BasisGridFunction<BasisType, VectorType> integrand_fn(basis,integrand);

            return integrate(basis.getGridView(), integrand_fn, quad_order);
        }
};


#endif