l2projection.hh

#ifndef DUNE_FEM_OPERATOR_PROJECTION_LOCAL_L2PROJECTION_HH
#define DUNE_FEM_OPERATOR_PROJECTION_LOCAL_L2PROJECTION_HH
 
#include <cstddef>
 
#include <type_traits>
#include <utility>
 
#include <dune/common/dynvector.hh>
#include <dune/common/bartonnackmanifcheck.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    // LocalL2Projection
    // -----------------
 
    template< class BasisFunctionSet, class Implementation >
    class LocalL2Projection
    {
   public:
      typedef BasisFunctionSet BasisFunctionSetType;
 
    protected:
#ifndef DOXYGEN
 
      LocalL2Projection () {}
 
#endif // #ifndef DOXYGEN
 
      LocalL2Projection ( const LocalL2Projection & ) = default;
 
      LocalL2Projection ( LocalL2Projection && ) {}
 
      LocalL2Projection &operator= ( const LocalL2Projection & ) = default;
 
      LocalL2Projection &operator= ( LocalL2Projection && ) {}
 
    public:
      BasisFunctionSet basisFunctionSet () const
      {
        CHECK_INTERFACE_IMPLEMENTATION( impl().basisFunctionSet() );
        return impl().basisFunctionSet();
      }
 
      template< class LocalFunction, class LocalDofVector >
      void operator() ( const LocalFunction &localFunction, LocalDofVector &localDofVector ) const
      {
        apply( localFunction, localDofVector );
      }
 
      template< class LocalFunction, class LocalDofVector >
      void apply ( const LocalFunction &localFunction, LocalDofVector &localDofVector ) const
      {
        CHECK_INTERFACE_IMPLEMENTATION( impl().apply( localFunction, localDofVector ) );
        impl().apply( localFunction, localDofVector );
      }
 
    protected:
      const Implementation &impl () const
      {
        return static_cast< const Implementation & >( *this );
      }
    };
 
 
 
    // DefaultLocalL2Projection
    // ------------------------
 
    template< class LocalRieszProjection, class Quadrature >
    class DefaultLocalL2Projection
    : public LocalL2Projection< typename LocalRieszProjection::BasisFunctionSetType, DefaultLocalL2Projection< LocalRieszProjection, Quadrature > >
    {
      typedef DefaultLocalL2Projection< LocalRieszProjection, Quadrature > ThisType;
      typedef LocalL2Projection< typename LocalRieszProjection::BasisFunctionSetType, DefaultLocalL2Projection< LocalRieszProjection, Quadrature > > BaseType;
 
    public:
      typedef LocalRieszProjection LocalRieszProjectionType;
      typedef typename BaseType::BasisFunctionSetType BasisFunctionSetType;
 
    private:
      typedef typename BasisFunctionSetType::RangeType RangeType;
      typedef typename RangeType::value_type RangeFieldType;
 
    public:
      template< class... Args >
      explicit DefaultLocalL2Projection ( Args &&... args )
        : rieszProjection_( std::forward< Args >( args )... )
      {}
 
      DefaultLocalL2Projection ( const ThisType & ) = default;
 
#ifndef DOXYGEN
 
      DefaultLocalL2Projection ( ThisType &other )
        : DefaultLocalL2Projection( static_cast< const ThisType & >( other ) )
      {}
 
#endif // #ifndef DOXYGEN
 
      DefaultLocalL2Projection ( ThisType &&other )
        : rieszProjection_( std::move( other.rieszProjection_ ) )
      {}
 
      DefaultLocalL2Projection &operator= ( const ThisType & ) = default;
 
      DefaultLocalL2Projection &operator= ( ThisType &&other )
      {
        rieszProjection_ = std::move( other.rieszProjection_ );
        return *this;
      }
 
      void unbind() {}
 
      BasisFunctionSetType basisFunctionSet () const
      {
        return rieszProjection_.basisFunctionSet();
      }
 
      template< class LocalFunction, class LocalDofVector >
      void apply ( const LocalFunction &localFunction, LocalDofVector &dofs ) const
      {
        static_assert( std::is_same< RangeType, typename LocalFunction::RangeType >::value,
                       "RangeType and LocalFunction::RangeType have to be the same type" );
 
        const BasisFunctionSetType basisFunctionSet = this->basisFunctionSet();
        const auto geometry = basisFunctionSet.entity().geometry();
 
        f_.resize( basisFunctionSet.size() );
        f_ = static_cast< RangeFieldType >( 0 );
 
        Quadrature quadrature( geometry.type(), localFunction.order() + basisFunctionSet.order() );
        const std::size_t nop = quadrature.nop();
        for( std::size_t qp = 0; qp < nop; ++qp )
        {
          RangeType value;
          localFunction.evaluate( quadrature[ qp ], value );
          value *= quadrature.weight( qp )*geometry.integrationElement( quadrature.point( qp ) );
          basisFunctionSet.axpy( quadrature[ qp ], value, f_ );
        }
 
        rieszProjection_.apply( f_, dofs );
      }
 
    private:
      LocalRieszProjectionType rieszProjection_;
      mutable Dune::DynamicVector< RangeFieldType > f_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_OPERATOR_PROJECTION_LOCAL_L2PROJECTION_HH