VectorBackend.hpp

#pragma once

#include <dune/istl/bvector.hh>
#include <dune/functions/backends/istlvectorbackend.hh>

#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
#include <amdis/linearalgebra/VectorBase.hpp>
#include <amdis/typetree/MultiIndex.hpp>

namespace AMDiS
{
  template <class T, class = void>
  struct BlockVectorType
  {
    using type = Dune::FieldVector<T,1>;
  };

  template <class T>
  struct BlockVectorType<T, typename T::field_type>
  {
    using type = T;
  };

  template <class T>
  class ISTLBlockVector
      : public VectorBase<ISTLBlockVector<T>>
  {
  public:
    using BaseVector = Dune::BlockVector<typename BlockVectorType<T>::type>;

    using block_type = typename BaseVector::block_type;

    using value_type = T;

    using field_type = typename block_type::field_type;

    using size_type  = typename BaseVector::size_type;

  public:
    template <class Basis>
    explicit ISTLBlockVector(Basis const&) {}

    BaseVector const& vector() const
    {
      return vector_;
    }

    BaseVector& vector()
    {
      return vector_;
    }

    std::size_t size() const
    {
      return vector_.size();
    }

    // NOTE: this resize works recursively on the hierarchy of the vector
    template <class SizeInfo>
    void init(SizeInfo const& sizeInfo, bool clear)
    {
      Dune::Functions::istlVectorBackend(vector_).resize(sizeInfo);
      if (clear)
        vector_ = 0;
    }

    template <class MultiIndex>
    value_type& at(MultiIndex const& idx)
    {
      return Dune::Functions::istlVectorBackend(vector_)[idx];
    }

    template <class MultiIndex>
    value_type const& at(MultiIndex const& idx) const
    {
      return Dune::Functions::istlVectorBackend(vector_)[idx];
    }

  private:
    BaseVector vector_;
  };

} // end namespace AMDiS