replace mtl types with dune types and added backend parameter to DOFVector and DOFMatrix

src/amdis/Assembler.hpp

@@ -3,11 +3,7 @@
 #include <memory>
 #include <tuple>
 
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-
 #include <amdis/DirichletBC.hpp>
-#include <amdis/LinearAlgebra.hpp>
 #include <amdis/LocalAssemblerList.hpp>
 #include <amdis/common/Mpl.hpp>
 

src/amdis/Assembler.inc.hpp

 #pragma once
 
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
 #include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <amdis/utility/TreePath.hpp>
 #include <amdis/utility/Visitor.hpp>
@@ -24,14 +26,14 @@ void Assembler<Traits>::assemble(
   // 2. create a local matrix and vector
   std::size_t localSize = localView.maxSize();
 
-  mtl::mat::dense2D<double> elementMatrix(localSize, localSize);
-  mtl::vec::dense_vector<double> elementVector(localSize);
+  Dune::DynamicMatrix<double> elementMatrix(localSize, localSize);
+  Dune::DynamicVector<double> elementVector(localSize);
 
   // 3. traverse grid and assemble operators on the elements
   for (auto const& element : elements(globalBasis_.gridView()))
   {
-    set_to_zero(elementMatrix);
-    set_to_zero(elementVector);
+    elementMatrix = 0;
+    elementVector = 0;
 
     localView.bind(element);
     auto geometry = element.geometry();
@@ -75,24 +77,9 @@ void Assembler<Traits>::assemble(
       });
     });
 
-    // add element-matrix to system-matrix
-    for (std::size_t i = 0; i < localView.size(); ++i) {
-      auto const row = localView.index(i);
-      for (std::size_t j = 0; j < localView.size(); ++j) {
-        if (std::abs(elementMatrix(i,j)) > threshold<double>) {
-          auto const col = localView.index(j);
-          matrix(row,col) += elementMatrix(i,j);
-        }
-      }
-    }
-
-    // add element-vector to system-vector
-    for (std::size_t i = 0; i < localView.size(); ++i) {
-      if (std::abs(elementVector[i]) > threshold<double>) {
-        auto const idx = localView.index(i);
-        rhs[idx] += elementVector[i];
-      }
-    }
+    // add element-matrix to system-matrix and element-vector to rhs
+    matrix.insert(localView, localView, elementMatrix);
+    rhs.insert(localView, elementVector);
 
     // unbind all operators
     forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {

src/amdis/GridFunctionOperator.hpp

@@ -6,6 +6,7 @@
 #include <amdis/GridFunctions.hpp>
 #include <amdis/LocalOperator.hpp>
 #include <amdis/QuadratureFactory.hpp>
+#include <amdis/common/Transposed.hpp>
 #include <amdis/common/Utility.hpp>
 #include <amdis/utility/FiniteElementType.hpp>
 
@@ -172,7 +173,7 @@ namespace AMDiS
                           ColNode const& colNode,
                           ElementMatrix& elementMatrix)
     {
-      auto elementMatrixTransposed = trans(elementMatrix);
+      auto elementMatrixTransposed = transposed(elementMatrix);
       transposedOp_.getElementMatrix(
         context, colNode, rowNode, elementMatrixTransposed);
     }

src/amdis/LocalAssemblerBase.hpp

@@ -2,7 +2,8 @@
 
 #include <type_traits>
 
-#include <boost/numeric/mtl/mtl.hpp>
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
 
 #include <amdis/ContextGeometry.hpp>
 
@@ -22,8 +23,8 @@ namespace AMDiS
     static_assert( numNodes == 1 || numNodes == 2,
       "VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
 
-    using ElementMatrix = mtl::mat::dense2D<double>; // TODO: choose correct value_type
-    using ElementVector = mtl::vec::dense_vector<double>;
+    using ElementMatrix = Dune::DynamicMatrix<double>; // TODO: choose correct value_type
+    using ElementVector = Dune::DynamicVector<double>;
 
     /// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
     using ElementMatrixVector = std::conditional_t<

src/amdis/ProblemInstat.hpp

@@ -60,6 +60,8 @@ namespace AMDiS
     /// Returns \ref oldSolution.
     std::unique_ptr<SystemVector> getOldSolutionVector() const
     {
+      test_exit_dbg(oldSolution,
+        "OldSolution need to be created. Call initialize with INIT_UH_OLD.");
       return *oldSolution;
     }
 

src/amdis/ProblemInstat.inc.hpp

@@ -53,7 +53,7 @@ template <class Traits>
 void ProblemInstat<Traits>::initTimestep(AdaptInfo&)
 {
   if (oldSolution)
-    oldSolution->copy(problemStat.getSolutionVector());
+    *oldSolution = problemStat.getSolutionVector();
 }
 
 } // end namespace AMDiS

src/amdis/common/CMakeLists.txt

@@ -18,6 +18,7 @@ install(FILES
     ScalarTypes.hpp
     Size.hpp
     Tags.hpp
+    Transposed.hpp
     TupleUtility.hpp
     Utility.hpp
     ValueCategory.hpp

src/amdis/common/Math.hpp

@@ -118,7 +118,7 @@ namespace AMDiS
   }
 
   template <class T>
-  constexpr T threshold = T(1.e-18L); //Math::sqr(std::numeric_limits<T>::epsilon());
+  constexpr T threshold = std::numeric_limits<T>::epsilon();
 
 
   /// Calculates factorial of i

src/amdis/common/Transposed.hpp

+#pragma once
+
+#include <type_traits>
+
+namespace AMDiS
+{
+  template <class Matrix>
+  class TransposedMatrix
+  {
+    using RawMatrix = std::decay_t<Matrix>;
+
+  public:
+    using size_type = typename RawMatrix::size_type;
+    using value_type = typename RawMatrix::value_type;
+
+  private:
+    struct ConstRowProxy
+    {
+      RawMatrix const* mat;
+      size_type row;
+
+      value_type const& operator[](size_type col) const
+      {
+        return (*mat)[col][row];
+      }
+    };
+
+    struct MutableRowProxy
+    {
+      RawMatrix* mat;
+      size_type row;
+
+      value_type& operator[](size_type col)
+      {
+        return (*mat)[col][row];
+      }
+    };
+
+  public:
+    template <class M>
+    TransposedMatrix(M&& matrix)
+      : matrix_(std::forward<M>(matrix))
+    {}
+
+    ConstRowProxy operator[](size_type row) const
+    {
+      return ConstRowProxy{&matrix_, row};
+    }
+
+    template <class SizeType, class M = Matrix,
+      std::enable_if_t<not std::is_const<M>::value, int> = 0>
+    MutableRowProxy operator[](SizeType row)
+    {
+      return MutableRowProxy{&matrix_, row};
+    }
+
+    size_type N() const
+    {
+      return matrix_.M();
+    }
+
+    size_type M() const
+    {
+      return matrix_.N();
+    }
+
+  private:
+    Matrix& matrix_;
+  };
+
+  template <class Matrix>
+  auto transposed(Matrix&& matrix)
+  {
+    using M = std::remove_reference_t<Matrix>;
+    return TransposedMatrix<M>{std::forward<Matrix>(matrix)};
+  }
+
+} // end namespace AMDiS

src/amdis/linear_algebra/CMakeLists.txt

 #install headers
 
 install(FILES
+    Common.hpp
+    DOFMatrixBase.hpp
+    DOFVectorBase.hpp
+    DOFVectorInterface.hpp
     HierarchicWrapper.hpp
     LinearSolverInterface.hpp
     PreconditionerInterface.hpp

src/amdis/linear_algebra/DOFMatrixBase.hpp

 #pragma once
 
+#include <cmath>
+#include <amdis/common/Math.hpp>
+
 namespace AMDiS
 {
-  template <class RowBasisType, class ColBasisType>
+  template <class RowBasisType, class ColBasisType, class Backend>
   class DOFMatrixBase
   {
   public:
@@ -15,6 +18,8 @@ namespace AMDiS
     /// The index/size - type
     using size_type = typename RowBasis::size_type;
 
+    using BaseMatrix = typename Backend::BaseMatrix;
+
   public:
     /// Constructor. Constructs new BaseVector.
     DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis)
@@ -34,6 +39,18 @@ namespace AMDiS
       return *colBasis_;
     }
 
+    /// Return the data-matrix
+    BaseMatrix const& matrix() const
+    {
+      return backend_.matrix();
+    }
+
+    /// Return the data-matrix
+    BaseMatrix& matrix()
+    {
+      return backend_.matrix();
+    }
+
     /// Return the size of the \ref rowBasis_
     size_type rows() const
     {
@@ -46,9 +63,60 @@ namespace AMDiS
       return colBasis_->dimension();
     }
 
+    /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern
+    /// If \p setToZero is true, the matrix is set to 0
+    void init(bool setToZero)
+    {
+      backend_.init(*rowBasis_, *colBasis_, setToZero);
+    }
+
+    /// Finish the matrix insertion, e.g. cleanup or final insertion
+    void finish()
+    {
+      backend_.finish();
+    }
+
+    /// Insert a block of values into the matrix (add to existing values)
+    template <class ElementMatrix>
+    void insert(typename RowBasis::LocalView const& rowLocalView,
+                typename ColBasis::LocalView const& colLocalView,
+                ElementMatrix const& elementMatrix)
+    {
+      for (size_type i = 0; i < rowLocalView.size(); ++i) {
+        size_type const row = flatMultiIndex(rowLocalView.index(i));
+        for (size_type j = 0; j < colLocalView.size(); ++j) {
+          if (std::abs(elementMatrix[i][j]) > threshold<double>) {
+            size_type const col = flatMultiIndex(colLocalView.index(j));
+            backend_.insert(row, col, elementMatrix[i][j]);
+          }
+        }
+      }
+    }
+
+    /// Insert a single value into the matrix (add to existing value)
+    template <class RowIndex, class ColIndex>
+    void insert(RowIndex row, ColIndex col, typename Backend::value_type const& value)
+    {
+      backend_.insert(flatMultiIndex(row), flatMultiIndex(col), value);
+    }
+
+    /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
+    /// a one on the diagonal of the matrix.
+    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
+    {
+      return backend_.applyDirichletBC(dirichletNodes);
+    }
+
+    size_type nnz() const
+    {
+      return backend_.nnz();
+    }
+
   protected:
     RowBasis const*	rowBasis_;
     ColBasis const*	colBasis_;
+
+    Backend backend_;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/DOFVectorBase.hpp

 #pragma once
 
-#include <algorithm>
-#include <memory>
-#include <string>
+#include <cmath>
 
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <dune/functions/functionspacebases/sizeinfo.hh>
 
-#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
+#include <amdis/common/Math.hpp>
 #include <amdis/common/ScalarTypes.hpp>
 #include <amdis/linear_algebra/DOFVectorInterface.hpp>
 #include <amdis/utility/MultiIndex.hpp>
@@ -15,7 +12,7 @@
 namespace AMDiS
 {
   /// The basic container that stores a base vector and a corresponding basis
-  template <class BasisType>
+  template <class BasisType, class Backend>
   class DOFVectorBase
       : public DOFVectorInterface
   {
@@ -28,15 +25,35 @@ namespace AMDiS
     /// The index/size - type
     using size_type  = typename BasisType::size_type;
 
+    /// The type of the elements of the DOFVector
+    using value_type = typename Backend::value_type;
+
+    using BaseVector = typename Backend::BaseVector;
+
+  public:
     /// Constructor. Constructs new BaseVector.
     DOFVectorBase(BasisType const& basis)
       : basis_(&basis)
-    {}
+    {
+      compress();
+    }
 
-    /// Return the basis \ref basis_ of the vector
+    /// Return the basis \ref basis_ associated to the vector
     Basis const& basis() const
     {
-        return *basis_;
+      return *basis_;
+    }
+
+    /// Return the data-vector
+    BaseVector const& vector() const
+    {
+      return backend_.vector();
+    }
+
+    /// Return the data-vector
+    BaseVector& vector()
+    {
+      return backend_.vector();
     }
 
     /// Return the size of the \ref basis
@@ -45,9 +62,63 @@ namespace AMDiS
       return basis_->dimension();
     }
 
+    void resize(Dune::Functions::SizeInfo<Basis> const& s)
+    {
+      backend_.resize(size_type(s));
+    }
+
+    /// Resize the \ref vector to the size of the \ref basis.
+    virtual void compress() override
+    {
+      if (backend_.size() != size()) {
+        backend_.resize(size());
+        vector() = 0;
+      }
+    }
+
+    /// Access the entry \p idx of the \ref vector with read-access.
+    template <class Index>
+    auto const& operator[](Index idx) const
+    {
+      size_type i = flatMultiIndex(idx);
+      return backend_[i];
+    }
+
+    /// Access the entry \p idx of the \ref vector with write-access.
+    template <class Index>
+    auto& operator[](Index idx)
+    {
+      size_type i = flatMultiIndex(idx);
+      return backend_[i];
+    }
+
+    /// Insert a block of values into the matrix (add to existing values)
+    template <class ElementVector>
+    void insert(typename Basis::LocalView const& localView,
+                ElementVector const& elementVector)
+    {
+      for (size_type i = 0; i < localView.size(); ++i) {
+        if (std::abs(elementVector[i]) > threshold<double>) {
+          size_type const idx = flatMultiIndex(localView.index(i));
+          backend_[idx] += elementVector[i];
+        }
+      }
+    }
+
+    /// Sets each DOFVector to the scalar \p value.
+    template <class Scalar,
+      std::enable_if_t<Concepts::Arithmetic<Scalar>, int> = 0>
+    Self& operator=(Scalar value)
+    {
+      vector() = value;
+      return *this;
+    }
+
   private:
     /// The finite element space / basis associated with the data vector
     Basis const* basis_;
+
+    Backend backend_;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/HierarchicWrapper.hpp

 #pragma once
 
+#if HAVE_MTL
 #include <boost/numeric/mtl/mtl_fwd.hpp>
+#endif
+
 #include <dune/functions/functionspacebases/sizeinfo.hh>
 
 #include <amdis/utility/MultiIndex.hpp>
@@ -86,11 +89,13 @@ namespace AMDiS
       vec.change_dim(std::size_t(s));
     }
 
+#if HAVE_MTL
     template <class... Params>
     std::size_t sizeImpl(mtl::vec::dense_vector<Params...> const& v) const
     {
       return mtl::size(v);
     }
+#endif
 
     template <class V>
     using HasSizeMember = decltype(std::declval<V>().size());

src/amdis/linear_algebra/istl/CMakeLists.txt

@@ -4,11 +4,11 @@
 #)
 
 install(FILES
+    DirectRunner.hpp
     DOFMatrix.hpp
     DOFVector.hpp
+    Fwd.hpp
     ISTL_Preconditioner.hpp
     ISTL_Solver.hpp
     ISTLRunner.hpp
-    LinearSolver.hpp
-    Preconditioner.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/istl)

src/amdis/linear_algebra/istl/DOFMatrix.hpp

@@ -8,7 +8,6 @@
 #include <dune/istl/matrixindexset.hh>
 
 #include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
 #include <amdis/linear_algebra/Common.hpp>
 #include <amdis/linear_algebra/DOFMatrixBase.hpp>
 #include <amdis/utility/MultiIndex.hpp>
@@ -27,13 +26,9 @@ namespace AMDiS
     using type = T;
   };
 
-  template <class RowBasisType, class ColBasisType,
-            class ValueType = double>
-  class DOFMatrix
-      : public DOFMatrixBase<RowBasisType, ColBasisType>
+  template <class ValueType>
+  class IstlMatrix
   {
-    using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
-
   public:
     /// The type of the elements of the DOFMatrix
     using value_type = typename BlockMatrixType<ValueType>::type;
@@ -41,53 +36,45 @@ namespace AMDiS
     /// 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.
-    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
-      : Super(rowBasis, colBasis)
-      , matrix_(ClonablePtr<BaseMatrix>::make())
-    {}
-
-    /// Constructor. Takes reference to a base matrix.
-    DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
-      : Super(rowBasis, colBasis)
-      , matrix_(matrix)
-    {}
+    IstlMatrix() = default;
 
     /// Return the data-vector \ref vector
     BaseMatrix const& matrix() const
     {
-      return *matrix_;
+      return matrix_;
     }
 
     /// Return the data-vector \ref vector
     BaseMatrix& matrix()
     {
-      return *matrix_;
+      return matrix_;
     }
 
 
-    /// Access the entry \p i of the \ref vector with write-access.
-    template <class Index>
-    value_type& operator()(Index row, Index col)
+    /// Insert a single value into the matrix (add to existing value)
+    void insert(size_type r, size_type c, value_type const& value)
     {
-      auto r = flatMultiIndex(row), c = flatMultiIndex(col);
       test_exit_dbg( initialized_, "Occupation pattern not initialized!");
-      test_exit_dbg( r < this->rows() && c < this->cols() ,
-          "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
-      return (*matrix_)[r][c];
+      test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
+          "Indices out of range [0,", matrix_.N(), ")x[0,", matrix_.M(), ")" );
+      matrix_[r][c] += value;
     }
 
     /// create occupation pattern and apply it to the matrix
-    void init(bool prepareForInsertion)
+    template <class RowBasis, class ColBasis>
+    void init(RowBasis const& rowBasis, ColBasis const& colBasis,
+              bool prepareForInsertion)
     {
-      occupationPattern_.resize(this->rows(), this->cols());
+      auto occupationPattern = Dune::MatrixIndexSet{rowBasis.dimension(), colBasis.dimension()};
 
-      // A loop over all elements of the grid
-      auto rowLocalView = this->rowBasis_->localView();
-      auto colLocalView = this->colBasis_->localView();
-
-      for (const auto& element : elements(this->rowBasis_->gridView())) {
+      auto rowLocalView = rowBasis.localView();
+      auto colLocalView = colBasis.localView();
+      for (const auto& element : elements(rowBasis.gridView())) {
         rowLocalView.bind(element);
         colLocalView.bind(element);
 
@@ -97,11 +84,11 @@ namespace AMDiS
           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.add(row, col);
           }
         }
       }
-      occupationPattern_.exportIdx(*matrix_);
+      occupationPattern.exportIdx(matrix_);
 
       initialized_ = true;
     }
@@ -113,16 +100,16 @@ namespace AMDiS
 
     std::size_t nnz() const
     {