#pragma once

#include <list>
#include <string>
#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>

#include <amdis/Output.hpp>
#include <amdis/linearalgebra/Common.hpp>
#include <amdis/linearalgebra/DOFMatrixBase.hpp>

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

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

  template <class ValueType>
  class IstlMatrix
  {
  public:
    /// The type of the elements of the DOFMatrix
    using value_type = typename BlockMatrixType<ValueType>::type;

    /// The matrix type of the underlying base matrix
    using BaseMatrix = Dune::BCRSMatrix<value_type>;

    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

  public:
    /// Constructor. Constructs new BaseVector.
    IstlMatrix() = default;

    /// Return the data-vector \ref vector
    BaseMatrix const& matrix() const
    {
      return matrix_;
    }

    /// Return the data-vector \ref vector
    BaseMatrix& matrix()
    {
      return matrix_;
    }


    /// Insert a single value into the matrix (add to existing value)
    void insert(size_type r, size_type c, value_type const& value)
    {
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
      test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
          "Indices out of range [0,{})x[0,{})", matrix_.N(), matrix_.M() );
      matrix_[r][c] += value;
    }

    /// create occupation pattern and apply it to the matrix
    template <class RowBasis, class ColBasis>
    void init(RowBasis const& rowBasis, ColBasis const& colBasis,
              bool prepareForInsertion)
    {
      auto occupationPattern = Dune::MatrixIndexSet{rowBasis.dimension(), colBasis.dimension()};

      auto rowLocalView = rowBasis.localView();
      auto colLocalView = colBasis.localView();
      for (const auto& element : elements(rowBasis.gridView())) {
        rowLocalView.bind(element);
        colLocalView.bind(element);

        for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
          // The global index of the i-th vertex of the element
          auto row = rowLocalView.index(i);
          for (std::size_t j = 0; j < colLocalView.size(); ++j) {
            // The global index of the j-th vertex of the element
            auto col = colLocalView.index(j);
            occupationPattern.add(row, col);
          }
        }
      }
      occupationPattern.exportIdx(matrix_);

      initialized_ = true;
    }

    void finish()
    {
      initialized_ = false;
    }

    std::size_t nnz() const
    {
      return matrix_.nonzeroes();
    }

  private:
    BaseMatrix matrix_;

    bool initialized_ = false;
  };

  template <class RowBasisType, class ColBasisType, class ValueType = double>
  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, IstlMatrix<ValueType>>;

} // end namespace AMDiS