integrator.hh

#ifndef DUNE_FEM_INTEGRATOR_HH
#define DUNE_FEM_INTEGRATOR_HH
 
#include <dune/fem/quadrature/quadrature.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    template< class Quadrature >
    class Integrator
    {
    public:
      typedef Quadrature QuadratureType;
 
      typedef typename QuadratureType :: EntityType EntityType;
 
    protected:
      const int order_;
 
    public:
      explicit Integrator ( unsigned int order )
      : order_( order )
      {}
 
      template< class Function >
      void integrateAdd ( const EntityType &entity,
                          const Function &function,
                          typename Function :: RangeType &ret ) const
      {
        typedef typename Function :: RangeType RangeType;
        typedef typename Function :: RangeFieldType RangeFieldType;
        typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
 
        const QuadratureType quadrature( entity, order_ );
        const auto& geometry = entity.geometry();
        for( const auto& qp : quadrature )
        {
          // evaluate function in quadrature point
          RangeType phi;
          function.evaluate( qp, phi );
 
          // calculate the weight of the quadrature point
          const RealType weight = geometry.integrationElement( qp.position() ) * qp.weight();
 
          ret.axpy( weight, phi );
        }
      }
 
      template< class Function >
      void integrate ( const EntityType &entity,
                       const Function &function,
                       typename Function :: RangeType &ret ) const
      {
        ret = 0;
        integrateAdd( entity, function, ret );
      }
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif  // #ifndef DUNE_FEM_INTEGRATOR_HH