VectorBackend.hpp

#pragma once

#include <Eigen/Dense>

#include <dune/common/ftraits.hh>

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

namespace AMDiS
{
  template <class T>
  class EigenVector
      : public VectorBase<EigenVector<T>>
  {
  public:
    using BaseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;

    using value_type = typename BaseVector::Scalar;

    using size_type  = typename BaseVector::Index;

  private:
    using block_type = value_type;

    using field_type = typename Dune::FieldTraits<value_type>::field_type;

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

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

    BaseVector& vector()
    {
      return vector_;
    }

    size_type size() const
    {
      return vector_.size();
    }


    template <class SizeInfo>
    void init(SizeInfo const& size, bool clear)
    {
      vector_.resize(size_type(size));
      if (clear)
        vector_.setZero();
    }


    template <class MultiIndex>
    value_type& at(MultiIndex const& idx)
    {
      const size_type i = flatMultiIndex(idx);
      test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
      return vector_.coeffRef(i);
    }

    template <class MultiIndex>
    value_type const& at(MultiIndex const& idx) const
    {
      const size_type i = flatMultiIndex(idx);
      test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
      return vector_.coeff(i);
    }

  private:
    BaseVector vector_;
  };


  namespace Recursive
  {
    template <class T>
    struct ForEach<EigenVector<T>> : ForEach<VectorBase<EigenVector<T>>> {};

    template <class T>
    struct Transform<EigenVector<T>> : Transform<VectorBase<EigenVector<T>>> {};

    template <class T>
    struct InnerProduct<EigenVector<T>> : InnerProduct<VectorBase<EigenVector<T>>> {};

  } // end namespace Recursive
} // end namespace AMDiS