gridfunctionadapter.hh

#ifndef DUNE_FEM_GRIDFUNCTIONADAPTER_HH
#define DUNE_FEM_GRIDFUNCTIONADAPTER_HH
 
#include <type_traits>
#include <utility>
 
#include <dune/common/exceptions.hh>
 
//- local includes
#include <dune/fem/version.hh>
#include <dune/fem/common/coordinate.hh>
#include <dune/fem/function/common/discretefunction.hh>
#include <dune/fem/storage/entitygeometry.hh>
 
// for compatibility
#include <dune/fem/function/common/localfunctionadapter.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    //- forward declaration
    template <class FunctionImp, class GridPartImp>
    class BasicGridFunctionAdapter;
 
    template <class FunctionImp, class GridPartImp>
    struct BasicGridFunctionAdapterTraits
    {
      typedef typename std::decay_t< FunctionImp >::FunctionSpaceType FunctionSpaceType;
 
      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 DiscreteFunctionSpaceAdapter<FunctionSpaceType,GridPartImp> DiscreteFunctionSpaceType;
 
      typedef GridPartImp GridPartType;
      typedef typename GridPartType :: GridType GridType;
      typedef typename GridPartType :: template Codim<0> :: EntityType EntityType;
      typedef typename GridPartType :: IntersectionType IntersectionType;
      typedef typename GridPartType :: template Codim<0> :: IteratorType IteratorType;
      typedef typename GridPartType :: IndexSetType IndexSetType;
 
      typedef BasicGridFunctionAdapter<FunctionImp,GridPartImp> DiscreteFunctionType;
    };
 
 
 
    // BasicGridFunctionAdapter
    // ------------------------
 
    template< class FunctionImp, class GridPartImp >
    class BasicGridFunctionAdapter
    : public Function< typename std::decay_t< FunctionImp >::FunctionSpaceType,
                       BasicGridFunctionAdapter< FunctionImp, GridPartImp > >,
      public HasLocalFunction
    {
      typedef BasicGridFunctionAdapter< FunctionImp, GridPartImp > ThisType;
      typedef Function< typename std::decay_t< FunctionImp >::FunctionSpaceType, ThisType > BaseType;
 
      // Make sure the function is not a discrete functon
      static_assert( !(std::is_convertible< FunctionImp, HasLocalFunction >::value),
                     "FunctionType may not be a discrete function type." );
 
    public:
      typedef BasicGridFunctionAdapterTraits< FunctionImp, GridPartImp > Traits;
 
      typedef std::decay_t< FunctionImp > FunctionType;
 
      typedef GridPartImp GridPartType;
 
      typedef typename Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      // type of discrete function space
      typedef typename Traits::FunctionSpaceType FunctionSpaceType;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType ;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType ;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType ;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType ;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
 
      typedef typename Traits :: EntityType  EntityType;
      typedef typename Traits :: IntersectionType IntersectionType;
 
    private:
      class LocalFunction;
 
    public:
      typedef LocalFunction LocalFunctionType;
 
      // reference to function this local belongs to
      BasicGridFunctionAdapter ( std::string name, FunctionImp f, const GridPartType &gridPart, unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : space_( gridPart, order ),
        function_( std::move( f ) ),
        name_( std::move( name ) )
      {}
 
      // reference to function this local belongs to
      BasicGridFunctionAdapter ( const ThisType &other )
      : space_( other.space_ ),
        function_( other.function_ ),
        name_( other.name_ )
      {}
 
      void evaluate ( const DomainType &global, RangeType &result ) const
      {
        function_.evaluate( global, result );
      }
 
      void jacobian ( const DomainType &global, JacobianRangeType &result ) const
      {
        function_.jacobian(global,result);
      }
 
      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();
      }
 
      inline int order () const
      {
        return space().order();
      }
 
      inline bool continuous () const
      {
        return space().continuous();
      }
 
    protected:
      DiscreteFunctionSpaceType space_;
      FunctionImp function_;
      const std::string name_;
    };
 
 
 
    // BasicGridFunctionAdapter::LocalFunction
    // ---------------------------------------
 
    template< class Function, class GridPart >
    class BasicGridFunctionAdapter< Function, GridPart >::LocalFunction
      : public EntityGeometryStorage< typename Traits::EntityType >
    {
      typedef EntityGeometryStorage< typename Traits::EntityType > BaseType;
      typedef LocalFunction ThisType;
      typedef BasicGridFunctionAdapter< Function, GridPart > DiscreteFunctionType;
 
    public:
      typedef typename Traits::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
      static const int dimDomain = GridPart::GridType::dimensionworld;
      static const int dimRange = FunctionSpaceType::dimRange;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename Traits::EntityType EntityType;
      typedef typename EntityType::Geometry Geometry;
      typedef typename Traits::IntersectionType IntersectionType;
      typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
      static const int dimLocal = LocalCoordinateType::dimension;
 
      LocalFunction ( const EntityType &entity, const DiscreteFunctionType &df )
      : BaseType( entity ),
        function_( &df.function_ ),
        order_( df.space().order() )
      {}
 
      LocalFunction ( const DiscreteFunctionType &df )
      : BaseType(),
        function_( &df.function_ ),
        order_( df.space().order() )
      {}
 
      using BaseType :: entity;
      using BaseType :: geometry;
      using BaseType :: bind;
      using BaseType :: unbind;
 
      LocalFunction ( const LocalFunction &other ) = default;
 
      template< class PointType >
      void evaluate ( const PointType &x, RangeType &ret ) const
      {
        auto global = geometry().global( coordinate( x ) );
        function().evaluate( global, ret );
      }
      template< class PointType >
      RangeType operator() ( const PointType &x ) const
      {
        RangeType ret;
        evaluate(x,ret);
        return ret;
      }
 
      template< class PointType >
      void jacobian ( const PointType &x, JacobianRangeType &ret ) const
      {
        const auto &geom = geometry();
        auto global = geom.global( coordinate( x ) );
        function().jacobian( global, ret );
 
        if( dimLocal != dimDomain )
        {
          // This computes the projection to the tangential space
          // (i.e. the hyperplane this entity is contained in). This
          // is done in a generic way by first projecting to the local
          // tangential space of the reference elment, and then
          // projecting back to the ambient space.
 
          const auto gjt = geom.jacobianTransposed( coordinate( x ) );
          const auto gjit = geom.jacobianInverseTransposed( coordinate( x ) );
 
          FieldVector< RangeFieldType, dimLocal > tmp;
          for( auto i = 0; i < dimRange; ++i )
          {
            gjit.mtv( ret[ i ], tmp );
            gjt.mtv( tmp, ret[ i ] );
          }
        }
      }
 
      template< class PointType >
      void hessian ( const PointType &x, HessianRangeType &ret ) const
      {
        DUNE_THROW( NotImplemented, "Method hessian() not implemented yet" );
      }
 
      template < class QuadratureType, class ... Vectors >
      void evaluateQuadrature( const QuadratureType& quadrature, Vectors& ... values ) const
      {
        static_assert( sizeof...( Vectors ) > 0, "evaluateQuadrature needs to be called with at least one vector." );
        std::ignore = std::make_tuple( ( evaluateQuadratureImp( quadrature, values ), 1 ) ... );
      }
 
      int order () const { return order_; }
 
      void init ( const EntityType &entity )
      {
        BaseType::bind( entity );
      }
 
    protected:
      template < class QuadratureType, class VectorType >
      auto evaluateQuadratureImp( const QuadratureType& quadrature, VectorType& values ) const
      -> std::enable_if_t< std::is_same< std::decay_t< decltype(values[ 0 ] ) >, RangeType >::value >
      {
        for( auto qp : quadrature )
          evaluate( qp, values[ qp.index() ] );
      }
 
      template < class QuadratureType, class VectorType >
      auto evaluateQuadratureImp( const QuadratureType& quadrature, VectorType& values ) const
      -> std::enable_if_t< std::is_same< std::decay_t< decltype(values[ 0 ] ) >, JacobianRangeType >::value >
      {
        for( auto qp : quadrature )
          jacobian( qp, values[ qp.index() ] );
      }
 
      const FunctionType &function () const
      {
        return *function_;
      }
 
      const FunctionType *function_;
      int order_;
    };
 
 
 
    // GridFunctionAdapter
    // -------------------
 
    template< class Function, class GridPart >
    using GridFunctionAdapter = BasicGridFunctionAdapter< const Function &, GridPart >;
 
 
 
    // gridFunctionAdapter
    // -------------------
 
    template< class Function, class GridPart >
    inline static GridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( std::string name, const Function &function, const GridPart &gridPart, unsigned int order )
    {
      return GridFunctionAdapter< Function, GridPart >( std::move( name ), function, gridPart, order );
    }
 
    template< class Function, class GridPart >
    inline static GridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( const Function &function, const GridPart &gridPart, unsigned int order )
    {
      return GridFunctionAdapter< Function, GridPart >( std::string(), function, gridPart, order );
    }
 
    template< class Function, class GridPart >
    inline static GridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( std::string name, Function &function, const GridPart &gridPart, unsigned int order )
    {
      const Function& cf = function;
      return gridFunctionAdapter( name, cf, gridPart, order );
    }
 
    template< class Function, class GridPart >
    inline static GridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( Function &function, const GridPart &gridPart, unsigned int order )
    {
      const Function& cf = function;
      return gridFunctionAdapter( cf, gridPart, order );
    }
 
    template< class Function, class GridPart >
    inline static BasicGridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( std::string name, Function &&function, const GridPart &gridPart, unsigned int order )
    {
      return BasicGridFunctionAdapter< Function, GridPart >( std::move( name ), std::move( function ), gridPart, order );
    }
 
    template< class Function, class GridPart >
    inline static BasicGridFunctionAdapter< Function, GridPart >
    gridFunctionAdapter ( Function &&function, const GridPart &gridPart, unsigned int order )
    {
      return BasicGridFunctionAdapter< Function, GridPart >( std::string(), std::move( function ), gridPart, order );
    }
 
 
 
    namespace
    {
      template <class FunctionImp,class GridPartType,bool>
      struct ConvertDFTypeHelper;
 
      template <class FunctionImp,class GridPartType>
      struct ConvertDFTypeHelper<FunctionImp,GridPartType,true>
      {
        typedef ConvertDFTypeHelper<FunctionImp,GridPartType,true> ThisType;
        enum {compatible = std::is_convertible<GridPartType,typename FunctionImp::DiscreteFunctionSpaceType::GridPartType>::value};
        typedef FunctionImp FunctionType;
        typedef typename FunctionType::DiscreteFunctionSpaceType DFSType;
        ConvertDFTypeHelper(const std::string& name,const FunctionImp& func,const GridPartType& gp) :
          func_(func)
        {}
        ConvertDFTypeHelper(const ConvertDFTypeHelper& other) :
          func_(other.func_)
        {}
        const FunctionType& function() const
        {
          return func_;
        }
        const DFSType& space() const
        {
          return func_.space();
        }
        private:
        const FunctionImp& func_;
      };
 
      template <class FunctionImp,class GridPartType>
      struct ConvertDFTypeHelper<FunctionImp,GridPartType,false>
        : GridFunctionAdapter<FunctionImp,GridPartType>
      {
        typedef ConvertDFTypeHelper<FunctionImp,GridPartType,false> ThisType;
        typedef GridFunctionAdapter<FunctionImp,GridPartType> BaseType;
        typedef BaseType FunctionType;
        typedef typename FunctionType::DiscreteFunctionSpaceType DFSType;
        ConvertDFTypeHelper(const std::string& name,const FunctionImp& func,const GridPartType& gp) :
          BaseType(name,func,gp)
        {}
        ConvertDFTypeHelper(const ConvertDFTypeHelper& other) :
          BaseType(other)
        {}
        const FunctionType& function() const
        {
          return *this;
        }
        const DFSType& space() const
        {
          return BaseType::space();
        }
      };
    }
 
    template< class FunctionImp, class GridPartImp >
    class ConvertToGridFunction
    : public Function< typename FunctionImp::FunctionSpaceType,
                       ConvertToGridFunction< FunctionImp, GridPartImp > >,
      public HasLocalFunction
    {
      typedef ConvertToGridFunction< FunctionImp, GridPartImp > ThisType;
      typedef Function< typename FunctionImp::FunctionSpaceType, ThisType > BaseType;
      static const bool hasLocalFunction = std::is_convertible< FunctionImp, HasLocalFunction >::value;
      typedef ConvertDFTypeHelper< FunctionImp, GridPartImp, hasLocalFunction > Helper;
      typedef typename Helper::FunctionType ConvertedType;
 
    public:
      typedef FunctionImp FunctionType;
      typedef GridPartImp GridPartType;
 
      typedef typename ConvertedType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      // type of discrete function space
      typedef typename ConvertedType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename DiscreteFunctionSpaceType::GridType GridType;
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType ;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType ;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType ;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType ;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
 
      typedef typename GridPartType :: template Codim<0> :: EntityType  EntityType;
 
      typedef typename ConvertedType::LocalFunctionType LocalFunctionType;
 
      ConvertToGridFunction ( const std::string &name,
                              const FunctionImp &function,
                              const GridPartType &gridPart )
      : name_( name ),
        helper_( name, function, gridPart )
      {}
 
      ConvertToGridFunction ( const ThisType &other )
      : name_( other.name_ ),
        helper_( other.helper_ )
      {}
 
      void evaluate ( const DomainType &global, RangeType &result ) const
      {
        helper_.function().evaluate(global,result);
      }
 
      const LocalFunctionType localFunction( const EntityType &entity ) const
      {
        return helper_.function().localFunction(entity);
      }
 
      LocalFunctionType localFunction( const EntityType &entity )
      {
        return helper_.function().localFunction(entity);
      }
 
      const std::string &name() const
      {
        return name_;
      }
 
      const DiscreteFunctionSpaceType &space() const
      {
        return helper_.space();
      }
 
    private:
      const std::string name_;
      Helper helper_;
    };
 
    template< class Function, class GridPart >
    inline ConvertToGridFunction< Function, GridPart >
    convertToGridFunction ( const std::string &name,
                            const Function &function,
                            const GridPart &gridPart )
    {
      return ConvertToGridFunction< Function, GridPart >( name, function, gridPart );
    }
 
  } // namespace Fem
 
} // namespace Dune
 
 
#endif // #ifndef DUNE_DISCRETEFUNCTIONADAPTER_HH