localfunction.hh

#ifndef DUNE_FEM_FUNCTION_LOCALFUNCTION_LOCALFUNCTION_HH
#define DUNE_FEM_FUNCTION_LOCALFUNCTION_LOCALFUNCTION_HH
 
#include <utility>
#include <type_traits>
#include <tuple>
 
#include <dune/common/deprecated.hh>
#include <dune/common/fvector.hh>
 
 
namespace Dune
{
  namespace Fem {
    // forward declaration for DenseMatVecTraits
    template< class BasisFunctionSet, class LocalDofVector >
    class LocalFunction;
  }
 
 
  // DenseMatVecTraits for LocalFunction for Discrete Functions
  // ----------------------------------------------------------
 
  template< class BasisFunctionSet, class LocalDofVector >
  struct DenseMatVecTraits< Fem::LocalFunction< BasisFunctionSet, LocalDofVector > >
  {
    typedef Fem::LocalFunction< BasisFunctionSet, LocalDofVector >  derived_type;
    typedef typename LocalDofVector::value_type                     value_type;
    //typedef LocalDofVector                                          container_type;
    typedef typename LocalDofVector::size_type                      size_type;
  };
 
 
 
  // FieldTraits for LocalContribution for Discrete Functions
  // --------------------------------------------------------
 
  template< class BasisFunctionSet, class LocalDofVector >
  struct FieldTraits< Fem::LocalFunction< BasisFunctionSet, LocalDofVector > >
  {
    typedef typename FieldTraits< typename LocalDofVector::value_type >::field_type field_type;
    typedef typename FieldTraits< typename LocalDofVector::value_type >::real_type  real_type;
  };
 
 
  namespace Fem
  {
 
    // LocalFunction
    // -------------
 
    template< class BasisFunctionSet, class LocalDofVector >
    class LocalFunction
      : public Dune::DenseVector< LocalFunction< BasisFunctionSet, LocalDofVector > >
    {
      typedef LocalFunction< BasisFunctionSet, LocalDofVector > ThisType;
      typedef Dune::DenseVector< ThisType > BaseType;
    public:
 
      typedef BasisFunctionSet BasisFunctionSetType;
 
      typedef LocalDofVector LocalDofVectorType;
 
      typedef typename LocalDofVectorType::value_type DofType;
 
      typedef typename LocalDofVectorType::size_type SizeType;
 
      typedef typename BasisFunctionSetType::EntityType EntityType;
 
      typedef typename EntityType::Geometry  Geometry;
 
      typedef typename BasisFunctionSetType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename Geometry::LocalCoordinate LocalCoordinateType;
 
      static const int dimDomain = FunctionSpaceType::dimDomain;
      static const int dimRange = FunctionSpaceType::dimRange;
 
      LocalFunction () {}
 
      explicit LocalFunction ( const BasisFunctionSetType &basisFunctionSet )
      : basisFunctionSet_( basisFunctionSet )
      {
        localDofVector_.resize( basisFunctionSet.size() );
      }
 
      explicit LocalFunction ( const LocalDofVectorType &localDofVector ) : localDofVector_( localDofVector ) {}
 
      LocalFunction ( const BasisFunctionSetType &basisFunctionSet, const LocalDofVector &localDofVector )
      : basisFunctionSet_( basisFunctionSet ),
        localDofVector_( localDofVector )
      {
        localDofVector_.resize( basisFunctionSet.size() );
      }
 
      explicit LocalFunction ( LocalDofVectorType &&localDofVector ) : localDofVector_( localDofVector ) {}
 
      LocalFunction ( const BasisFunctionSetType &basisFunctionSet, LocalDofVector &&localDofVector )
      : basisFunctionSet_( basisFunctionSet ),
        localDofVector_( localDofVector )
      {
        localDofVector_.resize( basisFunctionSet.size() );
      }
 
      LocalFunction ( ThisType && other )
      : basisFunctionSet_( std::move( other.basisFunctionSet_ ) ),
        localDofVector_( std::move( other.localDofVector_ ) )
      {}
 
      LocalFunction ( const ThisType & other )
      : basisFunctionSet_( other.basisFunctionSet_ ),
        localDofVector_( other.localDofVector_ )
      {}
 
      const DofType &operator[] ( SizeType num ) const { return localDofVector_[ num ]; }
 
      DofType &operator[] ( SizeType num ) { return localDofVector_[ num ]; }
 
      using BaseType :: operator +=;
      using BaseType :: operator -=;
      using BaseType :: operator *=;
      using BaseType :: operator /=;
 
      template< class T >
      void assign ( const LocalFunction< BasisFunctionSet, T > &other )
      {
        localDofVector() = other.localDofVector();
      }
 
      void clear ()
      {
        std::fill( localDofVector().begin(), localDofVector().end(), DofType( 0 ) );
      }
 
#if 0
      template< class T >
      ThisType &axpy ( const RangeFieldType s, const LocalFunction< BasisFunctionSet, T > &other )
      {
        localDofVector().axpy( s, other.localDofVector() );
        return *this;
      }
#endif
      using BaseType :: axpy;
 
      template< class PointType >
      void axpy ( const PointType &x, const RangeType &factor )
      {
        basisFunctionSet().axpy( x, factor, localDofVector() );
      }
 
      template< class PointType >
      void axpy ( const PointType &x, const JacobianRangeType &factor)
      {
        basisFunctionSet().axpy( x, factor, localDofVector() );
      }
      template< class PointType >
      void axpy ( const PointType &x, const HessianRangeType &factor)
      {
        basisFunctionSet().axpy( x, factor, localDofVector() );
      }
 
      template< class PointType >
      void axpy ( const PointType &x, const RangeType &factor1, const JacobianRangeType &factor2 )
      {
        basisFunctionSet().axpy( x, factor1, factor2, localDofVector() );
      }
 
      int order () const { return basisFunctionSet().order(); }
 
      const BasisFunctionSetType &basisFunctionSet () const { return basisFunctionSet_; }
 
      const EntityType &entity () const { return basisFunctionSet().entity(); }
 
      const Geometry &geometry () const { return basisFunctionSet().geometry(); }
 
 
      template< class PointType >
      void evaluate ( const PointType &x, RangeType &ret ) const
      {
        basisFunctionSet().evaluateAll( x, localDofVector(), ret );
      }
 
      template< class PointType >
      void jacobian ( const PointType &x, JacobianRangeType &ret ) const
      {
        basisFunctionSet().jacobianAll( x, localDofVector(), ret );
      }
 
      template< class PointType >
      void hessian ( const PointType &x, HessianRangeType &ret ) const
      {
        basisFunctionSet().hessianAll( x, localDofVector(), ret );
      }
 
      int numDofs () const { return localDofVector().size(); }
 
      SizeType size () const { return localDofVector().size(); }
 
      template< class QuadratureType, class ... Vectors >
      void axpyQuadrature ( const QuadratureType &quad, const Vectors& ... values )
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore =
          std::make_tuple( ( basisFunctionSet().axpy( quad, values, localDofVector() ), 1 )...);
      }
 
      template< class QuadratureType, class RangeVectorType, class JacobianRangeVectorType >
      void axpyQuadrature ( const QuadratureType &quad,
                            const RangeVectorType& rangeVector,
                            const JacobianRangeVectorType& jacobianVector )
      {
        basisFunctionSet().axpy( quad, rangeVector, jacobianVector, localDofVector() );
      }
 
      template< class QuadratureType, class ... Vectors >
      void evaluateQuadrature( const QuadratureType &quad, Vectors& ... vec ) const
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore =
          std::make_tuple( ( evaluateQuadrature( quad, vec, typename std::decay< decltype( vec[ 0 ] ) >::type() ), 1 )
            ... );
      }
 
      template< class QuadratureType, class ... Vectors >
      void jacobianQuadrature( const QuadratureType &quad, Vectors& ... vec ) const
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore =
          std::make_tuple( ( evaluateQuadrature( quad, vec, typename std::decay< decltype( vec[ 0 ] ) >::type() ), 1 )
            ... );
      }
 
      template< class QuadratureType, class ... Vectors >
      void hessianQuadrature( const QuadratureType &quad, Vectors& ... vec ) const
      {
        // make sure vec size if large enough
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore =
           std::make_tuple( ( evaluateQuadrature( quad, vec, typename std::decay< decltype( vec[ 0 ] ) >::type() ), 1 )
           ... );
      }
 
      const LocalDofVectorType &localDofVector () const { return localDofVector_; }
 
      LocalDofVectorType &localDofVector () { return localDofVector_; }
 
 
      bool valid () const
      {
        return basisFunctionSet_.valid();
      }
 
    protected:
      void init ( const EntityType& entity )
      {
        DUNE_THROW(NotImplemented,"LocalFunction::init( entity ) must be overloaded on derived class!");
      }
 
      void bind ( const EntityType& entity )
      {
        DUNE_THROW(NotImplemented,"LocalFunction::bind( entity ) must be overloaded on derived class!");
      }
 
      void unbind ()
      {
        // basically sets entity pointer inside basis function set to nullptr
        basisFunctionSet_ = BasisFunctionSetType();
      }
 
      void init ( const BasisFunctionSetType &basisFunctionSet )
      {
        basisFunctionSet_ = basisFunctionSet;
        localDofVector_.resize( basisFunctionSet.size() );
      }
 
      // evaluate local function and store results in vector of RangeTypes
      // this method only helps to identify the correct method on
      // the basis function set
      template< class QuadratureType, class VectorType >
      void evaluateQuadrature( const QuadratureType &quad, VectorType &result, const RangeType & ) const
      {
        basisFunctionSet().evaluateAll( quad, localDofVector(), result );
      }
 
      // evaluate jacobian of local function and store result in vector of
      // JacobianRangeTypes, this method only helps to identify the correct method on
      // the basis function set
      template< class QuadratureType, class VectorType >
      void evaluateQuadrature( const QuadratureType &quad, VectorType &result, const JacobianRangeType & ) const
      {
        basisFunctionSet().jacobianAll( quad, localDofVector(), result );
      }
 
      // evaluate jacobian of local function and store result in vector of
      // JacobianRangeTypes, this method only helps to identify the correct method on
      // the basis function set
      template< class QuadratureType, class VectorType >
      void evaluateQuadrature( const QuadratureType &quad, VectorType &result, const HessianRangeType & ) const
      {
        basisFunctionSet().hessianAll( quad, localDofVector(), result );
      }
 
      BasisFunctionSetType basisFunctionSet_;
      LocalDofVectorType localDofVector_;
    };
 
  } // end namespace Fem
 
} // end namespace Dune
 
#endif // #ifndef DUNE_FEM_FUNCTION_LOCALFUNCTION_LOCALFUNCTION_HH