combinedfunction.hh

#ifndef DUNE_FEM_FUNCTION_COMBINEDFUNCTION_COMBINEDFUNCTION_HH
#define DUNE_FEM_FUNCTION_COMBINEDFUNCTION_COMBINEDFUNCTION_HH
 
#include <algorithm>
#include <iostream>
#include <string>
 
#include <dune/fem/common/stackallocator.hh>
#include <dune/fem/function/common/discretefunction.hh>
#include <dune/fem/function/localfunction/mutable.hh>
#include <dune/fem/space/combinedspace/combinedspace.hh>
#include <dune/fem/space/common/dofmanager.hh>
#include <dune/fem/space/common/dofstorage.hh>
#include <dune/fem/storage/referencevector.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    //- Forward declarations
    template< class ContainedDiscreteFunction, int N >
    class CombinedDiscreteFunction;
 
    template< class ContainedDiscreteFunction, int N >
    class CombinedDiscreteFunctionDofIterator;
 
 
    //- Class definitions
    template< class ContainedDiscreteFunction, int N >
    struct DiscreteFunctionTraits< CombinedDiscreteFunction< ContainedDiscreteFunction, N > >
    {
      typedef ContainedDiscreteFunction ContainedDiscreteFunctionType;
      typedef CombinedDiscreteFunction< ContainedDiscreteFunctionType, N >
        DiscreteFunctionType;
 
      typedef typename ContainedDiscreteFunctionType :: DiscreteFunctionSpaceType
        ContainedDiscreteFunctionSpaceType;
      typedef CombinedSpace< ContainedDiscreteFunctionSpaceType, N, VariableBased >
        DiscreteFunctionSpaceType;
 
      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 RangeFieldType DofType;
      typedef typename DiscreteFunctionSpaceType :: MapperType MapperType;
      typedef typename DiscreteFunctionSpaceType :: GridType GridType;
      typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
 
      typedef CombinedDiscreteFunctionDofIterator< ContainedDiscreteFunctionType, N >
        DofIteratorType;
      typedef ConstDofIteratorDefault< DofIteratorType > ConstDofIteratorType;
 
      typedef typename ContainedDiscreteFunctionType :: DofBlockType DofBlockType;
      typedef typename ContainedDiscreteFunctionType :: ConstDofBlockType
        ConstDofBlockType;
      typedef typename ContainedDiscreteFunctionType :: DofBlockPtrType
        DofBlockPtrType;
      typedef typename ContainedDiscreteFunctionType :: ConstDofBlockPtrType
        ConstDofBlockPtrType;
 
      typedef ThreadSafeValue< UninitializedObjectStack > LocalDofVectorStackType;
      typedef StackAllocator< DofType, LocalDofVectorStackType* > LocalDofVectorAllocatorType;
      typedef DynamicReferenceVector< DofType, LocalDofVectorAllocatorType > LocalDofVectorType;
 
      typedef MutableLocalFunction< DiscreteFunctionType > LocalFunctionType;
    };
 
 
    template <class ContainedDiscreteFunctionImp,int N >
    class CombinedDiscreteFunction
    : public DiscreteFunctionDefault< CombinedDiscreteFunction< ContainedDiscreteFunctionImp, N > >
    {
      typedef CombinedDiscreteFunction< ContainedDiscreteFunctionImp, N > ThisType;
      typedef DiscreteFunctionDefault< CombinedDiscreteFunction< ContainedDiscreteFunctionImp, N > >
        BaseType;
 
    public:
      typedef ContainedDiscreteFunctionImp ContainedDiscreteFunctionType;
      typedef ContainedDiscreteFunctionType SubDiscreteFunctionType;
 
      typedef DiscreteFunctionTraits< ThisType > Traits;
 
      typedef typename BaseType::GridType GridType;
 
      typedef typename BaseType::GridPartType GridPartType;
 
      typedef typename BaseType::DiscreteFunctionType
      DiscreteFunctionType;
      typedef typename BaseType::DiscreteFunctionSpaceType
      DiscreteFunctionSpaceType;
      typedef typename Traits::ContainedDiscreteFunctionSpaceType
      ContainedDiscreteFunctionSpaceType;
      typedef ContainedDiscreteFunctionSpaceType
      SubDiscreteFunctionSpaceType;
      typedef typename BaseType::DofType DofType;
      typedef typename BaseType::RangeFieldType RangeFieldType;
      typedef typename BaseType::DomainFieldType DomainFieldType;
      typedef typename BaseType::RangeType RangeType;
      typedef typename BaseType::DomainType DomainType;
      typedef typename Traits::MapperType MapperType;
 
      typedef typename BaseType::DofIteratorType DofIteratorType;
      typedef typename BaseType::ConstDofIteratorType
      ConstDofIteratorType;
 
      typedef typename BaseType :: DofBlockPtrType DofBlockPtrType;
      typedef typename BaseType :: ConstDofBlockPtrType ConstDofBlockPtrType;
 
      typedef typename BaseType :: LocalDofVectorAllocatorType LocalDofVectorAllocatorType;
 
      using BaseType :: assign; // needs DofIterator!
      using BaseType :: axpy;
      using BaseType :: space;
 
      CombinedDiscreteFunction( const ContainedDiscreteFunctionType& func )
        : BaseType( "combined_"+func.name(), createSpace( func.space().gridPart() ), LocalDofVectorAllocatorType( &ldvStack_ ) ),
          ldvStack_( std::max( sizeof( DofType ), sizeof( DofType* ) ) * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize )
      {
        for (int i=0; i<N; ++i)
          func_[i] = new ContainedDiscreteFunctionType(func);
      }
 
      CombinedDiscreteFunction(const std::string &name,
                               const ContainedDiscreteFunctionSpaceType& spc)
        : BaseType( "combined_"+name, createSpace( spc.gridPart() ), LocalDofVectorAllocatorType( &ldvStack_ ) ),
          ldvStack_( std::max( sizeof( DofType ), sizeof( DofType* ) ) * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize )
      {
        for (int i=0; i<N; ++i)
          func_[i] = new ContainedDiscreteFunctionType(name,spc);
      }
 
      CombinedDiscreteFunction(const std::string &name,
                               const DiscreteFunctionSpaceType& spc )
        : BaseType( "combined_"+name, createSpace( spc.gridPart() ), LocalDofVectorAllocatorType( &ldvStack_) ),
          ldvStack_( std::max( sizeof( DofType ), sizeof( DofType* ) ) * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize )
      {
        for (int i=0; i<N; ++i)
        {
          func_[i] = new ContainedDiscreteFunctionType(name,space().containedSpace());
        }
      }
 
      CombinedDiscreteFunction(const ThisType &other)
        : BaseType( other.name()+"_copy", createSpace( other.space().gridPart() ), LocalDofVectorAllocatorType( &ldvStack_ ) ),
          ldvStack_( other.ldvStack_ )
      {
        for (int i=0; i<N; ++i)
          func_[i] = new ContainedDiscreteFunctionType(other.subFunction(i));
      }
 
      ~CombinedDiscreteFunction()
      {
        for (int i=0; i<N; ++i)
          delete func_[i];
 
        delete &space();
      }
 
      CombinedDiscreteFunction() = delete;
      ThisType& operator= ( const ThisType& ) = delete;
      ThisType& operator= ( ThisType&& ) = delete;
 
      void clear()
      {
        for (int i=0; i<N; ++i)
          func_[i]->clear();
      }
 
      void assign( const ThisType &g )
      {
        for( int i=0; i<N; ++i)
          func_[i]->assign( g.subFunction( i ) );
      }
 
      int size() const
      {
        return func_[0]->size()*N;
      }
 
      ThisType &operator+= ( const ThisType &g )
      {
        for (int i=0; i<N; ++i)
          *func_[ i ] +=  g.subFunction( i );
        return *this;
      }
 
      using BaseType::operator-=;
      ThisType &operator-= ( const ThisType &g )
      {
        for( int i = 0; i < N; ++i )
          *func_[ i ] -=  g.subFunction( i );
        return *this;
      }
 
      DiscreteFunctionType& operator *= (const RangeFieldType &scalar)
      {
        for (int i=0; i<N; ++i)
          *func_[i] *= scalar;
        return *this;
      }
 
      DiscreteFunctionType& operator /= (const RangeFieldType &scalar)
      {
        for (int i=0; i<N; ++i)
          *func_[i] /= scalar;
        return *this;
      }
 
      void addScaled( const ThisType &g, const RangeFieldType &s )
      {
        axpy( g, s );
      }
 
      void axpy( const RangeFieldType &s, const ThisType &g )
      {
        for (int i=0; i<N; ++i)
          func_[i]->axpy( s, g.subFunction( i ) );
      }
 
      RangeFieldType scalarProductDofs ( const ThisType &other ) const
      {
        RangeFieldType ret( 0 );
        for( int i = 0; i < N; ++i )
          ret += func_[ i ]->scalarProductDofs( other.subFunction( i ) );
        return ret;
      }
 
      template< class StreamTraits >
      void read ( InStreamInterface< StreamTraits >& in)
      {
        for (int i=0; i<N; ++i)
          func_[i]->read(in);
      }
      template< class StreamTraits >
      void write ( OutStreamInterface< StreamTraits >& out) const
      {
        for (int i=0; i<N; ++i)
          func_[i]->write(out);
      }
 
      void print( std :: ostream &out ) const
      {
        for (int i=0; i<N; ++i)
          func_[i]->print(out);
      }
 
      bool dofsValid () const
      {
        bool ret = func_[0]->dofsValid();
        for (int i=1;i<N;i++)
          ret |= func_[i]->dofsValid();
        return ret;
      }
 
      ConstDofBlockPtrType block ( unsigned int index ) const
      {
        // This is wrong with the current implementation of CombinedSpace
        const int containedSize = func_[ 0 ]->space().blockMapper().size();
        const int component = index / containedSize;
        const int containedIndex = index % containedSize;
        const ContainedDiscreteFunctionType& func = *(func_[ component ]);
        return func.block( containedIndex );
      }
 
      DofBlockPtrType block ( unsigned int index )
      {
        // This is wrong with the current implementation of CombinedSpace
        const int containedSize = func_[ 0 ]->space().blockMapper().size();
        const int component = index / containedSize;
        const int containedIndex = index % containedSize;
        return func_[ component ]->block( containedIndex );
      }
 
      const RangeFieldType &dof(unsigned int index) const
      {
        int variable = index / func_[0]->size();
        int point    = index % func_[0]->size();
        return func_[variable]->dof(point);
      }
 
      RangeFieldType &dof ( unsigned int index )
      {
        int variable = index / func_[0]->size();
        int point    = index % func_[0]->size();
        return func_[variable]->dof(point);
      }
 
      ConstDofIteratorType dbegin () const
      {
        return ConstDofIteratorType(DofIteratorType(*this));
      }
      ConstDofIteratorType dend () const
      {
        return ConstDofIteratorType(DofIteratorType(false,*this));
      }
      DofIteratorType dbegin ()
      {
        return DofIteratorType(*this);
      }
      DofIteratorType dend ()
      {
        return DofIteratorType(false,*this);
      }
 
      ContainedDiscreteFunctionType& subFunction( const int i )
      {
        return *(func_[i]);
      }
 
      const ContainedDiscreteFunctionType& subFunction( const int i ) const
      {
        return *(func_[i]);
      }
 
      ContainedDiscreteFunctionSpaceType& subSpace()
      {
        return space().containedSpace();
      }
 
    private:
      const ThisType& interface() const
      {
        return *this;
      }
 
      DiscreteFunctionSpaceType& createSpace( GridPartType& gp )
      {
        // we need to delete the space in the destructor
        return *(new DiscreteFunctionSpaceType( gp ));
      }
 
      typename Traits :: LocalDofVectorStackType ldvStack_;
      ContainedDiscreteFunctionType* func_[N];
      friend class CombinedDiscreteFunctionDofIterator<ContainedDiscreteFunctionType,N>;
    };
 
    template <class ContainedDiscreteFunctionImp,int N>
    class CombinedDiscreteFunctionDofIterator
    : public DofIteratorDefault <
      typename ContainedDiscreteFunctionImp::DofType ,
      CombinedDiscreteFunctionDofIterator<ContainedDiscreteFunctionImp,N> >
    {
    public:
      typedef CombinedDiscreteFunctionDofIterator<ContainedDiscreteFunctionImp,N> ThisType;
      typedef DiscreteFunctionTraits<CombinedDiscreteFunction< ContainedDiscreteFunctionImp,N> > Traits;
      typedef typename Traits::DiscreteFunctionType DiscreteFunctionType;
      typedef typename Traits::ContainedDiscreteFunctionType ContainedDiscreteFunctionType;
      typedef typename ContainedDiscreteFunctionType::DofIteratorType ContainedDofIteratorType;
      typedef typename ContainedDiscreteFunctionType::ConstDofIteratorType ContainedConstDofIteratorType;
      typedef typename Traits::DofType DofType;
 
      CombinedDiscreteFunctionDofIterator(bool end,const DiscreteFunctionType& df) :
        df_(const_cast<DiscreteFunctionType&>(df)),
        comp_(N-1),
        iter_(df.func_[N-1]->dend()),
        endIter_(df.func_[N-1]->dend())
      {}
      CombinedDiscreteFunctionDofIterator(const DiscreteFunctionType& df) :
        df_(const_cast<DiscreteFunctionType&>(df)),
        comp_(0),
        iter_(df.func_[0]->dbegin()),
        endIter_(df.func_[0]->dend())
      {}
      CombinedDiscreteFunctionDofIterator(bool end,DiscreteFunctionType& df) :
        df_(df),
        comp_(N-1),
        iter_(df.func_[N-1]->dend()),
        endIter_(df.func_[N-1]->dend())
      {}
      CombinedDiscreteFunctionDofIterator(DiscreteFunctionType& df) :
        df_(df),
        comp_(0),
        iter_(df.func_[0]->dbegin()),
        endIter_(df.func_[0]->dend())
      {}
      CombinedDiscreteFunctionDofIterator(const ThisType& other):
        df_(other.df_),
        comp_(other.comp_),
        iter_(other.iter_),
        endIter_(other.endIter_)
      {}
 
      ThisType& operator=(const ThisType& other)
      {
        df_ = other.df_;
        comp_ = other.comp_;
        iter_ = other.iter_;
        endIter_ = other.endIter_;
        return *this;
      }
 
      DofType& operator *()
      {
        return *iter_;
      }
 
      const DofType& operator * () const
      {
        return *iter_;
      }
 
      ThisType& operator++ ()
      {
        ++iter_;
        if (iter_==endIter_ && comp_<N-1)
        {
          ++comp_;
          iter_ = df_.func_[comp_]->dbegin();
          endIter_ = df_.func_[comp_]->dend();
        }
        return *this;
      }
 
      bool operator == (const ThisType & I ) const
      {
        return (comp_ == I.comp_) && (iter_ == I.iter_);
      }
 
      bool operator != (const ThisType & I ) const
      {
        return !((*this) == I);
      }
 
      void reset()
      {
        comp_    = 0;
        iter_    = df_.func_[ 0 ]->dbegin();
        endIter_ = df_.func_[ 0 ]->dend();
      }
 
  private:
      DiscreteFunctionType& df_;
      mutable int comp_;
      mutable ContainedDofIteratorType iter_,endIter_;
 
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_FUNCTION_COMBINEDFUNCTION_COMBINEDFUNCTION_HH