localfunctionadapter.hh

#ifndef DUNE_FEM_LOCALFUNCTIONADAPTER_HH
#define DUNE_FEM_LOCALFUNCTIONADAPTER_HH
 
#include <functional>
#include <memory>
#include <set>
#include <tuple>
#include <type_traits>
 
#include <dune/common/hybridutilities.hh>
 
#include <dune/fem/function/common/discretefunction.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    template< class LocalFunctionImpl >
    class LocalFunctionAdapter;
 
    template< class LocalFunctionImpl >
    class LocalFunctionAdapterLocalFunction;
 
 
    template< class LocalFunctionImpl >
    struct LocalFunctionAdapterTraits
    {
      typedef typename LocalFunctionImpl::FunctionSpaceType FunctionSpaceType;
      typedef typename LocalFunctionImpl::GridPartType GridPartType;
 
      // use storage as reference if no copy constructor
      template< class T, bool >
      struct LocalFuncType
      {
        typedef T& Type;
      };
      // otherwise use storage as copy
      template< class T >
      struct LocalFuncType< T, true >
      {
        typedef T Type;
      };
 
      // store the local function object by value if it can be copy constructed - otherwise store a reference
      static constexpr bool localFunctionHasCopyConstructor = std::is_copy_constructible<LocalFunctionImpl>::value;
      typedef typename LocalFuncType< LocalFunctionImpl, localFunctionHasCopyConstructor >::Type LocalFuncStorageType;
 
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename GridPartType :: GridType GridType;
      typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
      typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
      typedef typename GridPartType :: IndexSetType IndexSetType;
 
      typedef DiscreteFunctionSpaceAdapter< FunctionSpaceType, GridPartType > DiscreteFunctionSpaceType;
 
      typedef LocalFunctionAdapter< LocalFunctionImpl > DiscreteFunctionType;
      typedef LocalFunctionAdapterLocalFunction< LocalFunctionImpl > LocalFunctionType;
    };
 
 
 
    template< class LocalFunctionImpl >
    class LocalFunctionAdapter
    : public Function< typename LocalFunctionImpl::FunctionSpaceType, LocalFunctionAdapter< LocalFunctionImpl > >,
      public HasLocalFunction
    {
      typedef LocalFunctionAdapter< LocalFunctionImpl > ThisType;
      typedef Function< typename LocalFunctionImpl::FunctionSpaceType, ThisType > BaseType;
 
      friend class LocalFunctionAdapterLocalFunction< LocalFunctionImpl >;
 
    public:
      typedef ThisType  DiscreteFunctionType;
 
      typedef LocalFunctionImpl LocalFunctionImplType;
 
      typedef typename BaseType::FunctionType FunctionType;
 
      typedef LocalFunctionAdapterTraits< LocalFunctionImplType > Traits;
 
      typedef typename Traits::GridPartType GridPartType;
 
      typedef typename Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      typedef typename Traits::GridType GridType;
      typedef typename Traits::DomainFieldType DomainFieldType ;
      typedef typename Traits::RangeFieldType RangeFieldType ;
      typedef typename Traits::DomainType DomainType;
      typedef typename Traits::RangeType RangeType;
      typedef typename Traits::JacobianRangeType JacobianRangeType;
      typedef typename Traits::EntityType EntityType;
 
      typedef typename Traits::LocalFunctionType LocalFunctionType;
 
    protected:
      struct ArgumentIF
      {
        virtual void initialize( LocalFunctionType* ) const = 0;
 
        virtual ~ArgumentIF ()
        {}
      };
 
      template <class ArgType>
      struct ArgumentInitializer : public ArgumentIF
      {
        const ArgType arg_;
        const double time_;
 
        ArgumentInitializer( const ArgType& arg, double time )
        : arg_( arg ), time_( time )
        {}
 
        ~ArgumentInitializer ()
        {}
 
        virtual void initialize( LocalFunctionType* lf ) const
        {
          lf->initialize( arg_, time_ );
        }
      };
 
      typedef typename Traits::LocalFuncStorageType LocalFuncStorageType;
 
    public:
      LocalFunctionAdapter ( const std::string &name,
                             LocalFunctionImplType &localFunctionImpl,
                             const GridPartType &gridPart,
                             unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : space_( gridPart, order ),
        localFunctionImpl_( localFunctionImpl ),
        lfList_(),
        argInitializer_(),
        name_( name ),
        order_( order )
      {}
 
      LocalFunctionAdapter ( const std::string &name,
                             const LocalFunctionImplType &localFunctionImpl,
                             const GridPartType &gridPart,
                             unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : space_( gridPart, order ),
        localFunctionImpl_( localFunctionImpl ),
        lfList_(),
        argInitializer_(),
        name_( name ),
        order_( order )
      {}
 
      LocalFunctionAdapter ( const std::string &name,
                             LocalFunctionImplType &&localFunctionImpl,
                             const GridPartType &gridPart,
                             unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : space_( gridPart, order ),
        localFunctionImpl_( localFunctionImpl ),
        lfList_(),
        argInitializer_(),
        name_( name ),
        order_( order )
      {}
 
      template < class ...Args >
      LocalFunctionAdapter ( const std::string &name,
                             const GridPartType &gridPart,
                             unsigned int order,
                             Args&... args)
      : space_( gridPart, order ),
        localFunctionImpl_( args... ),
        lfList_(),
        argInitializer_(),
        name_( name ),
        order_( order )
      {}
 
      template < class ...Args >
      LocalFunctionAdapter ( const std::string &name,
                             const GridPartType &gridPart,
                             const std::tuple< Args&... > &args,
                             unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : LocalFunctionAdapter( name, gridPart, args, order,
           std::make_index_sequence< std::tuple_size< std::tuple<Args&...> >::value >{} )
      {}
 
      LocalFunctionAdapter( const ThisType &other )
      : space_( other.space_ ),
        localFunctionImpl_( other.localFunctionImpl_ ),
        lfList_(),
        argInitializer_(),
        name_( other.name_ ),
        order_( other.order_ )
      {}
 
      unsigned int order() const
      {
        return order_;
      }
 
      inline bool continuous () const
      {
        return space().continuous();
      }
 
      const LocalFuncStorageType& localFunctionImpl() const
      {
        return localFunctionImpl_;
      }
 
      LocalFuncStorageType& localFunctionImpl()
      {
        return localFunctionImpl_;
      }
 
      void evaluate(const DomainType& global, RangeType& result) const
      {
        DUNE_THROW( NotImplemented, "LocalFunctionAdapter::evaluate is not implemented." );
      }
 
      LocalFunctionType localFunction( const EntityType &entity )
      {
        return LocalFunctionType( entity, *this );
      }
 
      const LocalFunctionType localFunction( const EntityType &entity ) const
      {
        return LocalFunctionType( entity, *this );
      }
 
      const std::string &name() const
      {
        return name_;
      }
 
      const DiscreteFunctionSpaceType &space () const
      {
        return space_;
      }
 
      const GridPartType &gridPart () const
      {
        return space().gridPart();
      }
 
      template <class DFType>
      DiscreteFunctionType &operator+= ( const DFType &g )
      {
        DUNE_THROW( NotImplemented, "LocalFunctionAdapter::operator += is not implemented." );
        return *this;
      }
 
      template <class DFType>
      DiscreteFunctionType& operator -= (const DFType& g)
      {
        DUNE_THROW( NotImplemented, "LocalFunctionAdapter::operator -= is not implemented." );
        return *this;
      }
 
      DiscreteFunctionType &operator*= ( const RangeFieldType &scalar )
      {
        DUNE_THROW( NotImplemented, "LocalFunctionAdapter::operator *= is not implemented." );
        return *this;
      }
 
      DiscreteFunctionType &operator/= ( const RangeFieldType &scalar )
      {
        DUNE_THROW( NotImplemented, "LocalFunctionAdapter::operator /= is not implemented." );
        return *this;
      }
 
      template< class ArgumentType >
      void initialize( const ArgumentType &arg, double time )
      {
        constexpr bool hasCopyConstructor = Traits::localFunctionHasCopyConstructor;
        if( hasCopyConstructor)
        {
          argInitializer_ = std::make_unique< ArgumentInitializer< ArgumentType > >( arg, time );
          for( auto& localFunctionPtr : lfList_ )
            Hybrid::ifElse( hasCopyConstructor, [ & ]( auto&& ){ argInitializer_->initialize( localFunctionPtr ); } );
        }
        else
          DUNE_THROW(NotImplemented,"LocalFunctionAdapter::initialize is not implemented");
      }
 
      void registerLocalFunction( LocalFunctionType* lf ) const
      {
        if( Traits::localFunctionHasCopyConstructor )
        {
          if( argInitializer_ != nullptr )
            argInitializer_->initialize( lf );
          lfList_.insert( lf );
        }
      }
 
      void deleteLocalFunction( LocalFunctionType* lf ) const
      {
        if( Traits::localFunctionHasCopyConstructor )
          lfList_.erase( lf );
      }
 
    protected:
      template <class ...Args, std::size_t ...index>
      LocalFunctionAdapter ( const std::string &name,
                             const GridPartType &gridPart,
                             const std::tuple< Args&... > &args,
                             unsigned int order,
                             std::index_sequence< index... > )
      : space_( gridPart, order ),
        localFunctionImpl_( std::get< index >( args )... ),
        lfList_(),
        argInitializer_(),
        name_( name ),
        order_( order )
      {}
 
      DiscreteFunctionSpaceType space_;
      LocalFuncStorageType localFunctionImpl_;
      mutable std::set< LocalFunctionType * > lfList_;
      std::unique_ptr< ArgumentIF > argInitializer_ ;
      const std::string name_;
      const unsigned int order_;
    };
 
 
 
    template< class LocalFunctionImpl >
    class LocalFunctionAdapterLocalFunction
    {
      typedef LocalFunctionAdapterLocalFunction< LocalFunctionImpl > ThisType;
 
    public:
      typedef LocalFunctionImpl LocalFunctionImplType;
 
      typedef LocalFunctionAdapterTraits< LocalFunctionImplType > Traits;
 
      typedef typename Traits::FunctionSpaceType FunctionSpaceType;
 
      typedef typename Traits::DomainFieldType DomainFieldType;
      typedef typename Traits::RangeFieldType RangeFieldType;
      typedef typename Traits::DomainType DomainType;
      typedef typename Traits::RangeType RangeType;
      typedef typename Traits::JacobianRangeType JacobianRangeType;
      typedef typename Traits::HessianRangeType HessianRangeType;
 
      typedef typename Traits::DiscreteFunctionType DiscreteFunctionType;
      typedef typename Traits::EntityType EntityType;
 
      typedef typename Traits::LocalFuncStorageType LocalFuncStorageType;
 
      LocalFunctionAdapterLocalFunction ( const EntityType &entity, const DiscreteFunctionType &adapter )
      : adapter_( adapter ), localFunctionImpl_( adapter.localFunctionImpl_ )
      {
        // add local function to list
        adapter_.registerLocalFunction( this );
        localFunctionImpl().init( entity );
      }
 
      explicit LocalFunctionAdapterLocalFunction ( const DiscreteFunctionType &adapter )
      : adapter_( adapter ), localFunctionImpl_( adapter.localFunctionImpl_ )
      {
        // add local function to list
        adapter_.registerLocalFunction( this );
      }
 
      LocalFunctionAdapterLocalFunction ( const ThisType &other )
      : adapter_( other.adapter_ ), localFunctionImpl_( other.localFunctionImpl_ )
      {
        // add local function to list
        adapter_.registerLocalFunction( this );
      }
 
      ~LocalFunctionAdapterLocalFunction ()
      {
        // remove local function from list
        adapter_.deleteLocalFunction( this );
      }
 
      unsigned int order() const
      {
        return adapter_.order();
      }
 
      template< class PointType >
      void evaluate ( const PointType &x, RangeType &ret ) const
      {
        localFunctionImpl().evaluate(x,ret);
      }
 
      template< class PointType >
      void jacobian ( const PointType &x, JacobianRangeType &ret ) const
      {
        localFunctionImpl().jacobian( x, ret );
      }
 
      // hessian of local function
      template< class PointType >
      void hessian ( const PointType &x, HessianRangeType &ret ) const
      {
        localFunctionImpl().hessian( x, ret );
      }
 
      template< class QuadratureType, class ... Vectors  >
      void evaluateQuadrature( const QuadratureType &quad, Vectors& ... result ) const
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore = std::make_tuple(
            ( evaluateQuadrature( quad, result ), 1 ) ... );
      }
 
      template< class QuadratureType, class ... Vectors  >
      void jacobianQuadrature( const QuadratureType &quad, Vectors& ... result ) const
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore = std::make_tuple(
            ( evaluateQuadrature( quad, result ), 1 ) ... );
      }
 
      void init( const EntityType& entity )
      {
        localFunctionImpl().init( entity );
        entity_ = &entity;
      }
 
      const EntityType& entity() const
      {
        assert( entity_ );
        return *entity_;
      }
 
      template <class ArgumentType>
      void initialize ( const ArgumentType& arg, double time )
      {
        localFunctionImpl().initialize( arg, time );
      }
 
      template< class QuadratureType, class VectorType  >
      auto evaluateQuadrature( const QuadratureType &quad, VectorType &result ) const
      -> std::enable_if_t< std::is_same< std::decay_t< decltype(result[ 0 ] ) >, RangeType >::value >
      {
        const size_t quadNop = quad.nop();
        for(size_t i = 0; i<quadNop; ++i)
          evaluate( quad[ i ], result[ i ] );
      }
 
      template< class QuadratureType, class VectorType  >
      auto evaluateQuadrature( const QuadratureType &quad, VectorType &result ) const
      -> std::enable_if_t< std::is_same< std::decay_t< decltype(result[ 0 ] ) >, JacobianRangeType >::value >
      {
        const size_t quadNop = quad.nop();
        for(size_t i = 0; i<quadNop; ++i)
          jacobian( quad[ i ], result[ i ] );
      }
 
    protected:
      const LocalFuncStorageType& localFunctionImpl() const
      {
        return localFunctionImpl_;
      }
 
      LocalFuncStorageType& localFunctionImpl()
      {
        return localFunctionImpl_;
      }
 
      EntityType const* entity_ = nullptr;
      const DiscreteFunctionType &adapter_;
      LocalFuncStorageType localFunctionImpl_;
    };
 
 
 
    template<class DiscreteFunctionSpaceImpl>
    class LocalAnalyticalFunctionBinder
    {
    public:
      typedef DiscreteFunctionSpaceImpl DiscreteFunctionSpaceType;
      typedef LocalAnalyticalFunctionBinder<DiscreteFunctionSpaceType> ThisType;
 
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
 
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef std::function<RangeType(const DomainType&,double,const EntityType&)> AnalyticalFunctionType;
      typedef std::function<JacobianRangeType(const DomainType&,double,const EntityType&)> AnalyticalJacobianType;
      typedef std::function<HessianRangeType(const DomainType&,double,const EntityType&)> AnalyticalHessianType;
 
      LocalAnalyticalFunctionBinder(const AnalyticalFunctionType& f=[](const auto& ,auto ,const auto& ){return RangeType(0.0);},
                                    const AnalyticalJacobianType& j=[](const auto& ,auto ,const auto& ){return JacobianRangeType(0.0);},
                                    const AnalyticalHessianType& h=[](const auto& ,auto ,const auto& ){return HessianRangeType(0.0);},
                                    double t=0.0):
        f_(f),j_(j),h_(h),t_(t)
      {}
 
      LocalAnalyticalFunctionBinder(const ThisType& )=default;
      LocalAnalyticalFunctionBinder(ThisType&& )=default;
      ThisType& operator=(const ThisType& )=default;
      ThisType& operator=(ThisType&& )=default;
 
      AnalyticalFunctionType& function()
      {
        return f_;
      }
 
      const AnalyticalFunctionType& function() const
      {
        return f_;
      }
 
      AnalyticalJacobianType& jacobian()
      {
        return j_;
      }
 
      const AnalyticalJacobianType& jacobian() const
      {
        return j_;
      }
 
      AnalyticalHessianType& hessian()
      {
        return h_;
      }
 
      const AnalyticalHessianType& hessian() const
      {
        return h_;
      }
 
      template<class PointType>
      void evaluate(const PointType& x,RangeType& ret) const
      {
        ret=f_(entity().geometry().global(coordinate(x)),t_,entity());
      }
 
      template<class PointType>
      void jacobian(const PointType& x,JacobianRangeType& ret) const
      {
        ret=j_(entity().geometry().global(coordinate(x)),t_,entity());
      }
 
      template<class PointType>
      void hessian(const PointType& x,HessianRangeType& ret ) const
      {
        ret=h_(entity().geometry().global(coordinate(x)),t_,entity());
      }
 
      void init(const EntityType& entity)
      {
        entity_=&entity;
      }
 
      template<typename Arg>
      void initialize(Arg&& ,double time)
      {
        t_=time;
      }
 
      const EntityType& entity() const
      {
        assert(entity_);
        return *entity_;
      }
 
      double time() const
      {
        return t_;
      }
 
    private:
      EntityType const* entity_;
      AnalyticalFunctionType f_;
      AnalyticalJacobianType j_;
      AnalyticalHessianType h_;
      double t_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
 
#endif // #ifndef DUNE_FEM_LOCALFUNCTIONADAPTER_HH