defaultblockvectors.hh

#ifndef DUNE_FEM_SIMPLEBLOCKVECTOR_HH
#define DUNE_FEM_SIMPLEBLOCKVECTOR_HH
 
#include <algorithm>
#include <cassert>
#include <memory>
 
#include <dune/common/densevector.hh>
#include <dune/common/dynvector.hh>
#include <dune/common/exceptions.hh>
 
#include <dune/fem/common/hybrid.hh>
#include <dune/fem/misc/debug.hh>
#include <dune/fem/storage/dynamicarray.hh>
#include <dune/fem/storage/envelope.hh>
#include <dune/fem/storage/subvector.hh>
 
#if HAVE_DUNE_ISTL
#include <dune/istl/bvector.hh>
#endif
 
namespace Dune {
 
  namespace Fem {
 
  // tag for block vectors
  class IsBlockVector {};
 
 
 
  template< class Imp, class Field >
  class BlockVectorInterface
  : public IsBlockVector
  {
  protected:
    typedef DebugCounter<size_t> CounterType;
 
    BlockVectorInterface () {}
 
    typedef Imp    ThisType;
  public:
    typedef Field  FieldType;
    typedef Field  DofType;
 
    const ThisType& operator= ( const ThisType &other )
    {
      if( &asImp() != &other )
      {
        assign( other );
        sequence_ = other.sequence_;
      }
      return asImp();
    }
 
    const ThisType &operator+= ( const ThisType &other )
    {
      assert( asImp().size() == other.size() );
      const auto endit = asImp().end();
      auto oit = other.begin();
      for( auto it = asImp().begin(); it != endit; ++it, ++oit )
        *it += *oit;
 
      ++sequence_;
      return asImp();
    }
 
    const ThisType &operator-= ( const ThisType &other )
    {
      assert( asImp().size() == other.size() );
      const auto endit = asImp().end();
      auto oit = other.begin();
      for( auto it = asImp().begin(); it != endit; ++it, ++oit )
        *it -= *oit;
 
      ++sequence_;
      return asImp();
    }
 
    FieldType operator* ( const ThisType &other ) const
    {
      assert( asImp().size() == other.size() );
      FieldType sum( 0 );
      const auto endit = asImp().end();
      auto oit = other.asImp().begin();
      for( auto it = asImp().begin(); it != endit; ++it, ++oit )
        sum += (*it * *oit);
 
      return sum;
    }
 
    const ThisType &operator*= ( const FieldType &scalar )
    {
      const auto endit = asImp().end();
      for( auto it = asImp().begin(); it != endit; ++it )
        *it *= scalar;
 
      ++sequence_;
      return asImp();
    }
 
    void axpy ( const FieldType &scalar, const ThisType &other )
    {
      assert( asImp().size() == other.size() );
      const auto endit = asImp().end();
      auto oit = other.begin();
      for( auto it = asImp().begin(); it != endit; ++it, ++oit )
        *it += scalar * (*oit);
 
      ++sequence_;
    }
 
    void clear ()
    {
      std::fill( asImp().begin(), asImp().end(), FieldType( 0 ) );
      ++sequence_;
    }
 
    std::size_t usedMemorySize() const
    {
      return asImp().numDofs() * sizeof( FieldType ) ;
    }
 
    void copyContent( const size_t newIndex, const size_t oldIndex )
    {
      asImp()[ newIndex ] = asImp()[ oldIndex ];
    }
 
    void memMoveBackward(const size_t length, const size_t oldStartIdx, const size_t newStartIdx)
    {
      assert( newStartIdx >= oldStartIdx );
      // get new end of block which is offSet + (length of block - 1)
      size_t newIdx = newStartIdx + length - 1;
      assert( newIdx < asImp().size() );
      // copy all entries backwards
      for(size_t oldIdx = oldStartIdx + length-1; oldIdx >= oldStartIdx; --oldIdx, --newIdx )
      {
        assert( oldIdx < asImp().size() );
        // copy to new location
        copyContent( newIdx, oldIdx );
      }
    }
 
    void memMoveForward(const size_t length, const size_t oldStartIdx, const size_t newStartIdx)
    {
      assert( newStartIdx <= oldStartIdx );
      const size_t upperBound = oldStartIdx + length;
      // get new off set that should be smaller then old one
      size_t newIdx = newStartIdx;
      for(size_t oldIdx = oldStartIdx; oldIdx<upperBound; ++oldIdx, ++newIdx )
      {
        // copy to new location
        copyContent( newIdx, oldIdx );
      }
    }
 
    void setMemoryFactor( const double memFactor ) {}
 
  protected:
    // Copy block vectors.
    //    Note: No '++sequence_' here, sequence_ is only changed in public methods
    void assign ( const ThisType &other )
    {
      assert( asImp().size() == other.size() );
      std::copy( other.begin(), other.end(), asImp().begin() );
    }
 
    ThisType& asImp() { return static_cast< ThisType& > (*this); }
    const ThisType& asImp() const { return static_cast< const ThisType& > (*this); }
 
    mutable CounterType sequence_; // for consistency checks...
  };
 
 
 
  template< class Container, int BlockSize >
  class SimpleBlockVector
  : public BlockVectorInterface< SimpleBlockVector< Container, BlockSize>, typename Container::value_type >
  {
    typedef BlockVectorInterface< SimpleBlockVector< Container, BlockSize>, typename Container::value_type > BaseType;
    typedef SimpleBlockVector< Container, BlockSize > ThisType;
    typedef Container                                 ArrayType;
 
    typedef DebugCounter<size_t>                      CounterType;
 
  public:
    typedef ArrayType DofContainerType;
 
    typedef typename ArrayType::value_type                  FieldType;
    typedef typename ArrayType::iterator                    IteratorType;
    typedef typename ArrayType::const_iterator              ConstIteratorType;
    typedef typename ArrayType::size_type                   SizeType;
 
    typedef FieldType           value_type;
    typedef SizeType            size_type;
 
    typedef SubVector< DofContainerType, StaticOffsetSubMapper< BlockSize > > DofBlockType;
    typedef DofBlockType ConstDofBlockType;
 
    typedef Fem::Envelope< DofBlockType >       DofBlockPtrType;
    typedef Fem::Envelope< ConstDofBlockType >  ConstDofBlockPtrType;
 
    enum { blockSize = BlockSize };
 
    typedef Hybrid::IndexRange< int, blockSize > BlockIndices;
 
  protected:
    template <class Array>
    struct ContainerAccess
    {
      static FieldType* data( Array& array ) { return array.data(); }
      static const FieldType* data( const Array& array ) { return array.data(); }
    };
 
    template <class K>
    struct ContainerAccess< Dune::DynamicVector< K > >
    {
      typedef Dune::DynamicVector< K >  Array;
      static FieldType* data( Array& array ) { return array.container().data(); }
      static const FieldType* data( const Array& array ) { return array.container().data(); }
    };
 
  public:
    explicit SimpleBlockVector ( ArrayType& array )
    : array_( array )
    {}
 
    const ThisType& operator= ( const ThisType& other )
    {
      BaseType::operator=( other );
      return *this;
    }
 
    ConstDofBlockType operator[] ( const unsigned int i ) const
    {
      assert( i < size() );
      return ConstDofBlockType( array_, StaticOffsetSubMapper< blockSize >( i*blockSize ) );
    }
 
    DofBlockType operator[] ( const unsigned int i )
    {
      assert( i < size() );
      return DofBlockType( array_, StaticOffsetSubMapper< blockSize >( i*blockSize ) );
    }
 
    ConstDofBlockPtrType blockPtr( const unsigned int i ) const
    {
      return ConstDofBlockPtrType( this->operator[] ( i ) );
    }
 
    DofBlockPtrType blockPtr( const unsigned int i )
    {
      return DofBlockPtrType( this->operator[] ( i ) );
    }
 
    IteratorType begin() { return array().begin(); }
 
    ConstIteratorType begin() const { return array().begin(); }
 
    IteratorType end() { return array().end(); }
 
    ConstIteratorType end() const { return array().end(); }
 
    IteratorType beforeEnd() { return array().beforeEnd(); }
 
    ConstIteratorType beforeEnd() const { return array().beforeEnd(); }
 
    IteratorType beforeBegin() { return array().beforeBegin(); }
 
    ConstIteratorType beforeBegin() const { return array().beforeBegin(); }
 
    IteratorType find( const SizeType dof ) { return array().find( dof ); }
 
    ConstIteratorType find( const SizeType dof ) const { return array().find( dof ); }
 
    SizeType size () const { return array().size() / blockSize; }
 
    SizeType numDofs() const { return array().size(); }
 
    FieldType* data() { return ContainerAccess< ArrayType >::data( array() ); }
    const FieldType* data() const { return ContainerAccess< ArrayType >::data( array() ); }
 
    const ArrayType &array () const { return array_; }
    ArrayType &array () { return array_; }
 
  protected:
    ArrayType& array_;
  };
 
 
 
  template< class Container, unsigned int BlockSize >
  class MutableBlockVector
  : public SimpleBlockVector< Container, BlockSize >
  {
    typedef MutableBlockVector< Container, BlockSize > ThisType;
    typedef SimpleBlockVector < Container, BlockSize > BaseType;
 
    typedef Container                                  ArrayType;
    using BaseType :: array_;
    using BaseType :: sequence_;
  public:
 
    using BaseType :: array;
    using BaseType :: blockSize ;
    typedef typename BaseType :: SizeType SizeType;
 
    explicit MutableBlockVector ( SizeType size )
    : BaseType( *(new Container( size*blockSize ) ) )
    {}
 
    MutableBlockVector ( const ThisType &other )
    : BaseType( *(new Container( other.array().size() ) ) )
    {
      assign( other );
    }
 
    ~MutableBlockVector()
    {
      delete &array_;
    }
 
    void setMemoryFactor( const double memFactor )
    {
      doSetMemoryFactor( array_, memFactor );
    }
 
    void reserve ( const int size )
    {
      array().reserve( size*blockSize );
    }
 
    void resize ( SizeType size )
    {
      array().resize( size*blockSize );
      ++sequence_;
    }
 
  private:
    template <class T, class Allocator>
    void doSetMemoryFactor( Dune::Fem::DynamicArray< T, Allocator >& array, const double memFactor )
    {
      array_.setMemoryFactor( memFactor );
    }
 
    template <class Array>
    void doSetMemoryFactor( Array& , const double memFactor ) {}
 
  };
 
 
 
  template< class Field, unsigned int BlockSize >
  class MutableBlockVector< DynamicArray< Field >, BlockSize >
  : public SimpleBlockVector< StaticArray< Field >, BlockSize >
  {
    typedef StaticArray< Field >          StaticContainer ;
    typedef DynamicArray< Field >         MutableContainer ;
    typedef SimpleBlockVector< StaticContainer, BlockSize >   BaseType;
    typedef MutableBlockVector< MutableContainer, BlockSize > ThisType;
 
  protected:
    using BaseType :: array_;
    using BaseType :: sequence_;
 
    std::unique_ptr< MutableContainer > container_;
  public:
    using BaseType :: blockSize ;
    typedef typename BaseType :: SizeType SizeType;
 
    explicit MutableBlockVector ( SizeType size )
    : BaseType( allocateContainer( size*blockSize ) ),
      container_( static_cast< MutableContainer* > (&array_) )
    {}
 
    MutableBlockVector ( const ThisType &other )
    : BaseType( allocateContainer( other.array().size() ) ),
      container_( static_cast< MutableContainer* > (&array_) )
    {
      assign( other );
    }
 
    void reserve ( const int size )
    {
      assert( container_ );
      container_->reserve( size*blockSize );
    }
 
    void resize ( SizeType size )
    {
      assert( container_ );
      container_->resize( size*blockSize );
      ++sequence_;
    }
 
  protected:
    StaticContainer& allocateContainer( const SizeType size )
    {
      MutableContainer* container = new MutableContainer( size );
      return *container;
    }
  };
 
 
 
  template< class DofBlock >
  class ISTLBlockVector
  : public BlockVectorInterface< ISTLBlockVector< DofBlock >, typename DofBlock :: value_type >
  {
    typedef ISTLBlockVector< DofBlock>                            ThisType;
#if HAVE_DUNE_ISTL
    typedef BlockVector< DofBlock >                               ArrayType;
#else
    // fallback in case dune-istl is not present
    typedef Dune::DynamicVector< DofBlock >                       ArrayType;
#endif
    typedef BlockVectorInterface< ISTLBlockVector< DofBlock >, typename DofBlock :: value_type  >  BaseType;
 
 
    using BaseType :: sequence_;
 
  public:
    ISTLBlockVector ( const ThisType& ) = default;
 
    typedef ArrayType DofContainerType;
 
    enum { blockSize = DofBlock :: dimension };
    typedef Hybrid::IndexRange< int, blockSize > BlockIndices;
 
    typedef typename DofBlock :: value_type FieldType;
 
  protected:
    template <class EmbeddedIterator, class V>
    class Iterator
      : public ForwardIteratorFacade< Iterator< EmbeddedIterator,V >, V >
    {
    public:
      typedef V FieldType;
    protected:
      mutable EmbeddedIterator it_;
#ifndef NDEBUG
      EmbeddedIterator end_;
#endif
      int index_;
    public:
      Iterator( const EmbeddedIterator& it
#ifndef NDEBUG
              , const EmbeddedIterator& end = EmbeddedIterator()
#endif
              )
        : Iterator( it, 0
#ifndef NDEBUG
                  , end
#endif
                  )
      {}
 
      Iterator( const EmbeddedIterator& it, const int index
#ifndef NDEBUG
              , const EmbeddedIterator& end = EmbeddedIterator()
#endif
              )
        : it_( it ),
#ifndef NDEBUG
          end_( end ),
#endif
          index_( index )
      {}
 
      FieldType& dereference () const
      {
        assert( it_ != end_ );
        assert( index_ < blockSize );
        return (*it_)[ index_ ];
      }
 
      void increment ()
      {
        ++index_;
        if( index_ >= blockSize )
        {
          index_ = 0;
          ++it_;
        }
      }
 
      bool equals ( const Iterator &other ) const
      {
        return (it_ == other.it_) && (index_ == other.index_);
      }
 
    }; // end DofIteratorBlockVectorDiscreteFunction
 
  public:
    typedef Iterator< typename ArrayType::Iterator, FieldType >            IteratorType;
    typedef Iterator< typename ArrayType::ConstIterator, const FieldType > ConstIteratorType;
 
    typedef DofBlock DofBlockType;
    typedef const DofBlock ConstDofBlockType;
 
    typedef DofBlockType*  DofBlockPtrType;
    typedef ConstDofBlockType*  ConstDofBlockPtrType;
 
    typedef typename ArrayType::size_type   SizeType;
    typedef typename ArrayType::value_type  value_type;
 
    explicit ISTLBlockVector ( ArrayType* array )
    : array_( array )
    {}
 
    ISTLBlockVector () = default;
 
    const ThisType& operator= ( const ThisType& other )
    {
      if( this != &other )
      {
        array() = other.array();
      }
      return *this;
    }
 
    DofBlockPtrType blockPtr(const unsigned int i ) { return &array()[ i ]; }
    ConstDofBlockPtrType blockPtr(const unsigned int i ) const { return &array()[ i ]; }
 
    DofBlockType& operator[] (const unsigned int i ) { return array()[ i ]; }
    ConstDofBlockType& operator[] (const unsigned int i ) const { return array()[ i ]; }
 
    IteratorType begin() { return IteratorType( array().begin()
#ifndef NDEBUG
                                              , array().end()
#endif
                                              ); }
    ConstIteratorType begin() const
    {
      return ConstIteratorType( array().begin()
#ifndef NDEBUG
                              , array().end()
#endif
                              ); }
 
    IteratorType end() { return IteratorType( array().end() ); }
    ConstIteratorType end() const  { return ConstIteratorType( array().end() ); }
 
    IteratorType beforeEnd()
    {
      DUNE_THROW(NotImplemented,"ISTLBlockVector::beforeEnd not implemented yet");
      return array().end();
    }
 
    ConstIteratorType beforeEnd() const
    {
      DUNE_THROW(NotImplemented,"ISTLBlockVector::beforeEnd not implemented yet");
      return array().end();
    }
 
    IteratorType beforeBegin()
    {
      DUNE_THROW(NotImplemented,"ISTLBlockVector::beforeBegin not implemented yet");
      return array().end();
    }
 
    ConstIteratorType beforeBegin() const
    {
      DUNE_THROW(NotImplemented,"ISTLBlockVector::beforeBegin not implemented yet");
      return array().end();
    }
 
    IteratorType find( const SizeType dof )
    {
      const SizeType block = dof / blockSize;
      const SizeType index = dof % blockSize;
      return IteratorType( array().find( block ), index
#ifndef NDEBUG
                          , array().end()
#endif
                         );
    }
 
    ConstIteratorType find( const SizeType dof ) const
    {
      const SizeType block = dof / blockSize;
      const SizeType index = dof % blockSize;
      return ConstIteratorType( array().find( block ), index
#ifndef NDEBUG
                              , array().end()
#endif
                              );
    }
 
    SizeType size() const { return array().size(); }
 
    SizeType numDofs() const { return array().size() * DofBlock::dimension; }
 
    void reserve ( const int size )
    {
      array().reserve( size );
    }
 
    void resize ( SizeType size )
    {
      array().resize( size );
      ++sequence_;
    }
 
    ArrayType& array() { assert( array_ ); return *array_; }
    const ArrayType& array() const { assert( array_ ); return *array_; }
 
  protected:
    // ISTL BlockVector
    ArrayType* array_;
  };
 
} // namespace Fem
} // namespace Dune
 
#endif // DUNE_FEM_REFERENCEBLOCKVECTOR_HH