ISTL solvers and data structured usable

examples/ellipt.cc

@@ -15,7 +15,7 @@ using namespace AMDiS;
 using namespace Dune::Indices;
 
 // 1 component with polynomial degree 1
-using Param   = YaspGridBasis<AMDIS_DIM, 1>;
+using Param   = YaspGridBasis<AMDIS_DIM, 2>;
 using ElliptProblem = ProblemStat<Param>;
 
 int main(int argc, char** argv)

examples/init/ellipt.dat.2d

@@ -3,10 +3,10 @@ elliptMesh->global refinements: 0
 ellipt->mesh:                   elliptMesh
 
 ellipt->solver->name:           bicgstab
-ellipt->solver->max iteration:  1000
-ellipt->solver->relative tolerance: 1.e-8
-ellipt->solver->info:           10
-ellipt->solver->left precon:    diag
+ellipt->solver->max iteration:  10000
+ellipt->solver->relative tolerance: 1.e-10
+ellipt->solver->info:           1
+ellipt->solver->left precon:    ilu
 ellipt->solver->right precon:   no
 
 ellipt->output[0]->filename:    ellipt.2d

examples/stokes0.cc

@@ -80,10 +80,12 @@ int main(int argc, char** argv)
     // assemble and solve system
     prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
     // write matrix to file
     mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
     out << prob.getSystemMatrix().matrix();
     std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
 
     prob.solve(adaptInfo);
 

examples/stokes1.cc

@@ -81,10 +81,12 @@ int main(int argc, char** argv)
     // assemble and solve system
     prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
     // write matrix to file
     mtl::io::matrix_market_ostream out("matrix_stokes1.mtx");
     out << prob.getSystemMatrix().matrix();
     std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
 
     prob.solve(adaptInfo);
 

examples/vecellipt.cc

@@ -42,6 +42,7 @@ int main(int argc, char** argv)
 
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
   // write matrix to file
   if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 4) {
     mtl::io::matrix_market_ostream out("matrix.mtx");
@@ -49,6 +50,7 @@ int main(int argc, char** argv)
 
     std::cout << prob.getSystemMatrix().matrix() << '\n';
   }
+#endif
 
   prob.solve(adaptInfo);
   prob.writeFiles(adaptInfo, true);

src/amdis/Assembler.inc.hpp

@@ -24,8 +24,8 @@ void Assembler<Traits>::assemble(
   // 2. create a local matrix and vector
   std::size_t localSize = localView.maxSize();
 
-  mtl::mat::dense2D<typename SystemMatrixType::value_type> elementMatrix(localSize, localSize);
-  mtl::vec::dense_vector<typename SystemVectorType::value_type> elementVector(localSize);
+  mtl::mat::dense2D<double> elementMatrix(localSize, localSize);
+  mtl::vec::dense_vector<double> elementVector(localSize);
 
   // 3. traverse grid and assemble operators on the elements
   for (auto const& element : elements(globalBasis_.gridView()))

src/amdis/LinearAlgebra.hpp

 #pragma once
 
-#include <amdis/linear_algebra/LinearSolverInterface.hpp>
+#include <amdis/linear_algebra/LinearSolver.hpp>
 #include <amdis/linear_algebra/SolverInfo.hpp>
+
+#if HAVE_MTL
 #include <amdis/linear_algebra/mtl/DOFVector.hpp>
 #include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
-#include <amdis/linear_algebra/mtl/LinearSolver.hpp>
 #include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
+#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
+#else
+#include <amdis/linear_algebra/istl/DOFVector.hpp>
+#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
+#endif
\ No newline at end of file

src/amdis/linear_algebra/DOFMatrixBase.hpp

+#pragma once
+
+#include <cstddef>
+
+namespace AMDiS
+{
+  template <class T>
+  struct Triplet
+  {
+    std::size_t row, cols;
+    T value;
+  };
+
+} // end namespace AMDiS

src/amdis/linear_algebra/DOFVectorBase.hpp

+#pragma once
+
+#include <algorithm>
+#include <memory>
+#include <string>
+
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+#include <amdis/Output.hpp>
+#include <amdis/common/ClonablePtr.hpp>
+#include <amdis/common/ScalarTypes.hpp>
+#include <amdis/linear_algebra/DOFVectorInterface.hpp>
+#include <amdis/utility/MultiIndex.hpp>
+
+namespace AMDiS
+{
+  /// The basic container that stores a base vector and a corresponding basis
+  template <class BasisType>
+  class DOFVectorBase
+      : public DOFVectorInterface
+  {
+    using Self = DOFVectorBase;
+
+  public:
+    /// The type of the \ref basis
+    using Basis = BasisType;
+
+    /// The index/size - type
+    using size_type  = typename BasisType::size_type;
+
+    /// Constructor. Constructs new BaseVector.
+    DOFVectorBase(BasisType const& basis)
+      : basis_(&basis)
+    {}
+
+    /// Return the basis \ref basis_ of the vector
+    Basis const& basis() const
+    {
+        return *basis_;
+    }
+
+    /// Return the size of the \ref basis
+    size_type size() const
+    {
+      return basis_->dimension();
+    }
+
+  private:
+    /// The finite element space / basis associated with the data vector
+    Basis const* basis_;
+  };
+
+} // end namespace AMDiS

src/amdis/linear_algebra/mtl/LinearSolver.hpp → src/amdis/linear_algebra/LinearSolver.hpp

@@ -14,18 +14,18 @@ namespace AMDiS
   /** \ingroup Solver
    *
    * \brief
-   * Wrapper class for various MTL4 solvers. These solvers
+   * Wrapper class for various solvers. These solvers
    * are parametrized by MatrixType and VectorType. The different
    * algorithms, like krylov subspace methods or other external
-   * solvers where MTL4 provides an interface, can be assigned
+   * solvers where the backend provides an interface, can be assigned
    * by different Runner objects.
    **/
-  template <class Matrix, class Vector, class Runner>
+  template <class Matrix, class VectorX, class VectorB, class Runner>
   class LinearSolver
-      : public LinearSolverInterface<Matrix, Vector, Vector>
+      : public LinearSolverInterface<Matrix, VectorX, VectorB>
   {
     using Self = LinearSolver;
-    using Super = LinearSolverInterface<Matrix, Vector, Vector>;
+    using Super = LinearSolverInterface<Matrix, VectorX, VectorB>;
     using RunnerBase = typename Super::RunnerBase;
 
   public:
@@ -44,15 +44,15 @@ namespace AMDiS
       : runner_(std::make_shared<Runner>(prefix))
     {}
 
-    /// Implements \ref LinearSolverInterface::getRunner()
-    virtual std::shared_ptr<RunnerBase> getRunner() override
+    /// Implements \ref LinearSolverInterface::runner()
+    virtual std::shared_ptr<RunnerBase> runner() override
     {
       return runner_;
     }
 
   private:
     /// Implements \ref LinearSolverInterface::solveSystemImpl()
-    virtual void solveImpl(Matrix const& A, Vector& x, Vector const& b,
+    virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
                            SolverInfo& solverInfo) override
     {
       Dune::Timer t;
@@ -61,9 +61,8 @@ namespace AMDiS
         runner_->init(A);
       }
 
-      if (solverInfo.getInfo() > 0) {
-        msg("fill MTL4 matrix needed ", t.elapsed(), " seconds");
-      }
+      if (solverInfo.getInfo() > 0)
+        msg("fill matrix needed ", t.elapsed(), " seconds");
 
       int error = runner_->solve(A, x, b, solverInfo);
       solverInfo.setError(error);

src/amdis/linear_algebra/LinearSolverInterface.hpp

@@ -48,7 +48,7 @@ namespace AMDiS
     }
 
     // return the runner/worker corresponding to this solver (optional)
-    virtual std::shared_ptr<RunnerBase> getRunner() { return {}; };
+    virtual std::shared_ptr<RunnerBase> runner() { return {}; };
 
   private:
     /// main methods that all solvers must implement

src/amdis/linear_algebra/istl/CMakeLists.txt

@@ -11,7 +11,4 @@ install(FILES
     ISTLRunner.hpp
     LinearSolver.hpp
     Preconditioner.hpp
-    SystemMatrix.hpp
-    SystemVector.hpp
-    UmfpackRunner.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/istl)

src/amdis/linear_algebra/istl/DOFMatrix.hpp

@@ -8,111 +8,126 @@
 
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
-  template <class RowFeSpaceType, class ColFeSpaceType,
-            class ValueType = Dune::FieldMatrix<double,1,1>>
+  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 RowBasisType, class ColBasisType,
+            class ValueType = double>
   class DOFMatrix
   {
   public:
-    using RowFeSpace = RowFeSpaceType;
-    using ColFeSpace = ColFeSpaceType;
-    using BaseMatrix = Dune::BCRSMatrix<ValueType>;
+    /// The type of the finite element space / basis of the row
+    using RowBasis = RowBasisType;
+
+    /// The type of the finite element space / basis of the column
+    using ColBasis = ColBasisType;
 
-    using size_type  = typename RowFeSpaceType::size_type;
-    using field_type = typename ValueType::field_type;
+    /// The type of the elements of the DOFMatrix
+    using value_type = typename BlockMatrixType<ValueType>::type;
 
-    using value_type = ValueType;
+    /// 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(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
-      : rowFeSpace(rowFeSpace)
-      , colFeSpace(colFeSpace)
-      , name(name)
-      , matrix(ClonablePtr<BaseMatrix>::make())
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis)
+      : rowBasis_(&rowBasis)
+      , colBasis_(&colBasis)
+      , matrix_(ClonablePtr<BaseMatrix>::make())
     {}
 
-    /// Constructor. Takes pointer of data-vector.
-    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name,
-              BaseMatrix& matrix_ref)
-      : rowFeSpace(rowFeSpace)
-      , colFeSpace(colFeSpace)
-      , name(name)
-      , matrix(matrix_ref)
+    /// Constructor. Takes reference to a base matrix.
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix)
+      : rowBasis_(&rowBasis)
+      , colBasis_(&colBasis)
+      , matrix_(matrix)
     {}
 
-    /// Return the row-basis \ref rowFeSpace of the matrix
-    RowFeSpace const& getRowFeSpace() const
+    /// Return the row-basis \ref rowBasis of the matrix
+    RowBasis const& rowBasis() const
     {
-      return rowFeSpace;
+      return *rowBasis_;
     }
 
-    /// Return the col-basis \ref colFeSpace of the matrix
-    ColFeSpace const& getColFeSpace() const
+    /// Return the col-basis \ref colBasis of the matrix
+    ColBasis const& colBasis() const
     {
-      return colFeSpace;
+      return *colBasis_;
     }
 
     /// Return the data-vector \ref vector
-    BaseMatrix const& getMatrix() const
+    BaseMatrix const& matrix() const
     {
-      return *matrix;
+      return *matrix_;
     }
 
     /// Return the data-vector \ref vector
-    BaseMatrix& getMatrix()
+    BaseMatrix& matrix()
     {
-      return *matrix;
+      return *matrix_;
     }
 
     /// Return the size of the \ref feSpace
-    size_type N() const
-    {
-      return rowFeSpace.size();
-    }
-
-    size_type M() const
+    size_type rows() const
     {
-      return colFeSpace.size();
+      return rowBasis_->size();
     }
 
-    /// Return the \ref name of this vector
-    std::string getName() const
+    size_type cols() const
     {
-      return name;
+      return colBasis_->size();
     }
 
 
     /// Access the entry \p i of the \ref vector with read-access.
-    value_type const& operator()(size_type r, size_type c) const
+    template <class Index>
+    value_type const& operator()(Index row, Index col) const
     {
-      test_exit_dbg( initialized, "Occupation pattern not initialized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
-      return (*matrix)[r][c];
+      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
+      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
+      test_exit_dbg( r < rows() && c < cols() ,
+          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
+      return (*matrix_)[r][c];
     }
 
     /// Access the entry \p i of the \ref vector with write-access.
-    value_type& operator()(size_type r, size_type c)
+    template <class Index>
+    value_type& operator()(Index row, Index col)
     {
-      test_exit_dbg( initialized, "Occupation pattern not initialized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
-      return (*matrix)[r][c];
+      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
+      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
+      test_exit_dbg( r < rows() && c < cols() ,
+          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
+      return (*matrix_)[r][c];
     }
 
     /// create occupation pattern and apply it to the matrix
-    void init()
+    void init(bool prepareForInsertion)
     {
-      occupationPattern.resize(rowFeSpace.size(), colFeSpace.size());
-      auto meshView = rowFeSpace.gridView();
+      occupationPattern_.resize(rowBasis_->size(), colBasis_->size());
 
       // A loop over all elements of the grid
-      auto rowLocalView = rowFeSpace.localView();
-      auto colLocalView = colFeSpace.localView();
+      auto rowLocalView = rowBasis_->localView();
+      auto colLocalView = colBasis_->localView();
 
-      for (const auto& element : elements(meshView)) {
+      for (const auto& element : elements(rowBasis_->gridView())) {
         rowLocalView.bind(element);
         colLocalView.bind(element);
 
@@ -122,54 +137,50 @@ 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;
+      initialized_ = true;
     }
 
     void finish()
     {
-      initialized = false;
+      initialized_ = false;
     }
 
-    std::size_t nnz()
+    std::size_t nnz() const
     {
-      return matrix->nonzeroes();
+      return matrix_->nonzeroes();
     }
 
-    auto clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
+    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
     {
       // loop over the matrix rows
-      for (std::size_t i = 0; i < matrix->N(); ++i) {
+      for (std::size_t i = 0; i < matrix_->N(); ++i) {
         if (dirichletNodes[i]) {
-            auto cIt = (*matrix)[i].begin();
-            auto cEndIt = (*matrix)[i].end();
+            auto cIt = (*matrix_)[i].begin();
+            auto cEndIt = (*matrix_)[i].end();
             // loop over nonzero matrix entries in current row
             for (; cIt != cEndIt; ++cIt)
-                *cIt = (apply && i == cIt.index() ? 1 : 0);
+                *cIt = (i == cIt.index() ? 1 : 0);
         }
       }
 
-      std::vector<std::map<size_t,value_type>> columns; // symmetric dbc not yet implemented
+      std::list<Triplet<value_type>> columns; // symmetric dbc not yet implemented
       return columns;
     }
 
   private:
-    RowFeSpace const&	rowFeSpace;
-    ColFeSpace const&	colFeSpace;
-    std::string	        name;
-
-    ClonablePtr<BaseMatrix> matrix;
-    Dune::MatrixIndexSet    occupationPattern;
+    RowBasis const*	rowBasis_;
+    ColBasis const*	colBasis_;
 
-    bool initialized = false;
+    ClonablePtr<BaseMatrix> matrix_;
+    Dune::MatrixIndexSet occupationPattern_;
 
-    // friend class declarations
-    template <class, class, class> friend class SystemMatrix;
+    bool initialized_ = false;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/istl/DOFVector.hpp

 #pragma once
 
-#include <string>
-#include <memory>
-
 #include <dune/istl/bvector.hh>
 
-#include <dune/functions/functionspacebases/interpolate.hh>
-
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
-    template <class FeSpaceType, class ValueType = Dune::FieldVector<double,1>>
-    class DOFVector
-    {
-        using Self = DOFVector;
+  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 BasisType, class ValueType = double>
+	class DOFVector
+      : public DOFVectorBase<BasisType>
+	{
+	  using Self = DOFVector;
+    using Super = DOFVectorBase<BasisType>;
 
-    public:
-	using FeSpace    = FeSpaceType;
-	using BaseVector = Dune::BlockVector<ValueType>;
+	public:
+    /// The type of the finite element space / basis
+	  using Basis = BasisType;
 
-	using size_type  = typename FeSpace::size_type;
-	using field_type = typename ValueType::field_type;
+    /// The type of the elements of the DOFVector
+	  using value_type = typename BlockVectorType<ValueType>::type;
 
-	using value_type = ValueType;
+    /// The vector type of the underlying base vector
+	  using BaseVector = Dune::BlockVector<value_type>;
 
-	/// Constructor. Constructs new BaseVector.
-	DOFVector(FeSpace const& feSpace, std::string name)
-	  : feSpace(feSpace)
-	  , name(name)
-	  , vector(ClonablePtr<BaseVector>::make())
-	{
-	    compress();
-	}
-
-	/// Constructor. Takes pointer of data-vector.
-	DOFVector(FeSpace const& feSpace, std::string name,
-		  BaseVector& vector_ref)
-	  : feSpace(feSpace)
-	  , name(name)
-	  , vector(vector_ref)
-	{}
-
-	/// Return the basis \ref feSpace of the vector
-	FeSpace const& getFeSpace() const
-	{
-	    return feSpace;
-	}
+    /// The index/size - type
+	  using size_type  = typename Basis::size_type;
 
-	/// Return the data-vector \ref vector
-	BaseVector const& vector() const
-	{
-	    return *vector;
-	}
-
-	/// Return the data-vector \ref vector
-	BaseVector& vector()
-	{
-	    return *vector;
-	}
-
-	/// Return the size of the \ref feSpace
-	size_type getSize() const
-	{
-	    return feSpace.size();
-	}
+    /// Constructor. Constructs new BaseVector.