gridfunctionadapter.hh

#ifndef __DUNE_ACFEM_GRIDFUNCTIONADAPTER_HH__
#define __DUNE_ACFEM_GRIDFUNCTIONADAPTER_HH__
 
//- local includes
#include <dune/fem/version.hh>
#include <dune/fem/function/common/discretefunction.hh>
#include <dune/fem/quadrature/quadrature.hh>
 
// for compatibility
#include <dune/fem/function/common/localfunctionadapter.hh>
 
#include "../expressions/expressionoperations.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
#if __DUNE_ACFEM_MAKE_CHECK__
    // ugly pre-processor stuff, but so what.
    static inline
    size_t& globalToLocalCount()
    {
      static size_t count;
 
      return count;
    }
#endif
 
    //- forward declaration
    template <class FunctionImp, class GridPartImp>
    class GridFunctionAdapter;
 
    // Also ACFem::GridFunctionAdapter constructs should be stored as copy, as
    // they are already kind of an expression.
    template<class Function, class GridPart>
    struct ExpressionStorage<GridFunctionAdapter<Function, GridPart> >
      : public ExpressionCopyStorage<GridFunctionAdapter<Function, GridPart> >
    {
      typedef GridFunctionAdapter<Function, GridPart> ObjectType;
      typedef ExpressionCopyStorage<ObjectType> BaseType;
 
      using BaseType::BaseType;
    };
 
    template <class FunctionImp, class GridPartImp>
    struct GridFunctionAdapterTraits
    {
      typedef typename 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 Fem::DiscreteFunctionSpaceAdapter<FunctionSpaceType,GridPartImp> DiscreteFunctionSpaceType;
 
      typedef GridPartImp GridPartType;
      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 GridFunctionAdapter<FunctionImp,GridPartImp> DiscreteFunctionType;
    };
 
 
 
    // GridFunctionAdapter
    // -------------------
 
    template< class FunctionImp, class GridPartImp >
    class GridFunctionAdapter
      : public Fem::Function< typename FunctionImp::FunctionSpaceType,
                       GridFunctionAdapter< FunctionImp, GridPartImp > >,
        public Fem::HasLocalFunction
    {
      typedef GridFunctionAdapter< FunctionImp, GridPartImp > ThisType;
      typedef Fem::Function< typename FunctionImp::FunctionSpaceType, ThisType > BaseType;
 
#if 0
      // Make sure the function is not a discrete functon
      static_assert( !(std::is_base_of< HasLocalFunction,  FunctionImp >::value),
                     "FunctionType may not be a discrete function type." );
#endif
 
    public:
      typedef GridFunctionAdapterTraits< FunctionImp, GridPartImp > Traits;
 
      typedef FunctionImp FunctionType;
 
      typedef GridPartImp GridPartType;
 
      typedef ThisType DiscreteFunctionType;
 
      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 DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename Traits :: EntityType  EntityType;
 
    private:
      class LocalFunction;
      class LocalFunctionStorage;
 
      public:
      typedef LocalFunction LocalFunctionType;
      typedef LocalFunctionStorage LocalFunctionStorageType;
 
      // reference to function this local belongs to
      GridFunctionAdapter ( const std::string &name,
                            const FunctionType &f,
                            const GridPartType &gridPart,
                            unsigned int order = DiscreteFunctionSpaceType::polynomialOrder )
      : space_( gridPart, order ),
        localFunctionStorage_( *this ),
        function_( f ),
        name_( name )
      {}
 
      // reference to function this local belongs to
      GridFunctionAdapter ( const ThisType &other )
      : space_( other.space_ ),
        localFunctionStorage_( *this ),
        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);
      }
 
      void hessian ( const DomainType &global, HessianRangeType &result ) const
      {
        function_.hessian(global,result);
      }
 
      const LocalFunctionType localFunction ( const EntityType &entity ) const
      {
        return localFunctionStorage().localFunction( entity );
      }
 
      LocalFunctionType localFunction ( const EntityType &entity )
      {
        return localFunctionStorage().localFunction( entity );
      }
 
      LocalFunctionStorageType &localFunctionStorage () const
      {
        return localFunctionStorage_;
      }
 
      const std::string &name () const
      {
        return name_;
      }
 
      const DiscreteFunctionSpaceType &space () const
      {
        return space_;
      }
 
      const GridPartType &gridPart () const { return space().gridPart(); }
 
    protected:
      DiscreteFunctionSpaceType space_;
      mutable LocalFunctionStorageType localFunctionStorage_;
      const FunctionType& function_;
      const std::string name_;
    };
 
    // GridFunctionAdapter::LocalFunction
    // ----------------------------------
 
    template< class Function, class GridPart >
    class GridFunctionAdapter< Function, GridPart >::LocalFunction
    {
      typedef LocalFunction ThisType;
    public:
      typedef GridFunctionAdapter< Function, GridPart > DiscreteFunctionType;
 
      typedef typename Traits::EntityType EntityType;
 
      typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      static const int dimRange = DiscreteFunctionSpaceType::dimRange;
      static const int dimDomain = GridPart::GridType::dimensionworld;
      static const int dimLocal = GridPart::GridType::dimension;
 
      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 DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
 
      LocalFunction ( const EntityType &entity, const DiscreteFunctionType &df )
      : function_( &df.function_ ),
        entity_( &entity ),
        order_( df.space().order() )
      {}
 
      LocalFunction ( const DiscreteFunctionType &df )
      : function_( &df.function_ ),
        entity_( 0 ),
        order_( df.space().order() )
      {}
 
      LocalFunction ( LocalFunctionStorage &storage )
      : function_( &storage.function().function_ ),
        entity_( 0 ),
        order_( storage.function().space().order() )
      {}
 
      template< class PointType >
      void evaluate ( const PointType &x, RangeType &ret ) const
      {
        typedef typename EntityType::Geometry GeometryType;
        const GeometryType geometry = entity().geometry();
        DomainType global = geometry.global( Fem::coordinate( x ) );
        function().evaluate( global, ret );
#if __DUNE_ACFEM_MAKE_CHECK__
        ++globalToLocalCount();
#endif
      }
 
      template< class PointType >
      void jacobian ( const PointType &x, JacobianRangeType &ret ) const
      {
        typedef typename EntityType::Geometry GeometryType;
        const GeometryType geometry = entity().geometry();
 
        const typename GeometryType::LocalCoordinate cx = Fem::coordinate( x );
        DomainType global = geometry.global( cx );
        function().jacobian( global, ret );
 
        if( dimLocal != dimDomain )
        {
          const typename GeometryType::JacobianTransposed gjt = geometry.jacobianTransposed( cx );
          const typename GeometryType::JacobianInverseTransposed gjit = geometry.jacobianInverseTransposed( cx );
 
          FieldVector< RangeFieldType, dimLocal > tmp;
          for( int 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
      {
        typedef typename EntityType::Geometry GeometryType;
        const GeometryType geometry = entity().geometry();
 
        const typename GeometryType::LocalCoordinate cx = Fem::coordinate( x );
        DomainType global = geometry.global( cx );
        function().hessian( global, ret );
 
        static_assert(dimLocal == dimDomain, "dimLocal != dimDomain, please FIXME");
      }
 
      template < class QuadratureType, class VectorType >
      void evaluateQuadrature( const QuadratureType& quadrature, VectorType& values ) const
      {
        assert( values.size() == quadrature.nop() );
        evaluateQuadratureImp( quadrature, values, values[ 0 ] );
      }
 
      int order () const { return order_; }
 
      void init ( const EntityType &entity )
      {
        entity_ = &entity;
      }
 
      const EntityType &entity () const
      {
        assert( entity_ );
        return *entity_;
      }
 
    protected:
      template < class QuadratureType, class VectorType >
      void evaluateQuadratureImp( const QuadratureType& quadrature, VectorType& values, const RangeType& ) const
      {
        const unsigned int nop = quadrature.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
        {
          evaluate( quadrature[ qp ], values[ qp ] );
        }
      }
 
      template < class QuadratureType, class VectorType >
      void evaluateQuadratureImp( const QuadratureType& quadrature, VectorType& values, const JacobianRangeType& ) const
      {
        const unsigned int nop = quadrature.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
        {
          jacobian( quadrature[ qp ], values[ qp ] );
        }
      }
 
      const FunctionType &function () const
      {
        return *function_;
      }
 
      const FunctionType *function_;
      const EntityType *entity_;
      int order_;
    };
 
 
 
    // GridFunctionAdapter::LocalFunctionStorage
    // ---------------------------------------------
 
    template< class Function, class GridPart >
    class GridFunctionAdapter< Function, GridPart >::LocalFunctionStorage
    {
      typedef LocalFunctionStorage ThisType;
      typedef GridFunctionAdapter< Function, GridPart > DiscreteFunctionType;
 
    public:
      typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
 
      explicit LocalFunctionStorage ( DiscreteFunctionType &discreteFunction )
      : discreteFunction_( discreteFunction )
      {}
 
      LocalFunctionType localFunction ()
      {
        return LocalFunctionType( discreteFunction_ );
      }
 
      template< class Entity >
      const LocalFunctionType localFunction ( const Entity &entity ) const
      {
        return LocalFunctionType( entity, discreteFunction_ );
      }
 
      template< class Entity >
      LocalFunctionType localFunction ( const Entity &entity )
      {
        return LocalFunctionType( entity, discreteFunction_ );
      }
 
      DiscreteFunctionType &function ()
      {
        return discreteFunction_;
      }
 
    private:
      LocalFunctionStorage ( const ThisType & );
      ThisType operator= ( const ThisType & );
 
      DiscreteFunctionType &discreteFunction_;
    };
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
 
#endif // #ifndef __DUNE_ACFEM_DISCRETEFUNCTIONADAPTER_HH__