From e744545fd873a66cf4461fc5b32ffe023d346a15 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Sat, 29 Sep 2018 23:26:00 -0400
Subject: [PATCH] replace mtl types with dune types and added backend parameter
to DOFVector and DOFMatrix
---
src/amdis/Assembler.hpp | 4 -
src/amdis/Assembler.inc.hpp | 31 ++---
src/amdis/GridFunctionOperator.hpp | 3 +-
src/amdis/LocalAssemblerBase.hpp | 7 +-
src/amdis/ProblemInstat.hpp | 2 +
src/amdis/ProblemInstat.inc.hpp | 2 +-
src/amdis/common/CMakeLists.txt | 1 +
src/amdis/common/Math.hpp | 2 +-
src/amdis/common/Transposed.hpp | 78 +++++++++++++
src/amdis/linear_algebra/CMakeLists.txt | 4 +
src/amdis/linear_algebra/DOFMatrixBase.hpp | 70 +++++++++++-
src/amdis/linear_algebra/DOFVectorBase.hpp | 91 +++++++++++++--
.../linear_algebra/HierarchicWrapper.hpp | 5 +
src/amdis/linear_algebra/istl/CMakeLists.txt | 4 +-
src/amdis/linear_algebra/istl/DOFMatrix.hpp | 71 +++++-------
src/amdis/linear_algebra/istl/DOFVector.hpp | 101 ++++------------
src/amdis/linear_algebra/mtl/CMakeLists.txt | 1 -
src/amdis/linear_algebra/mtl/DOFMatrix.hpp | 68 +++++------
src/amdis/linear_algebra/mtl/DOFVector.hpp | 108 +++++-------------
19 files changed, 369 insertions(+), 284 deletions(-)
create mode 100644 src/amdis/common/Transposed.hpp
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 8b43adaf..e487c99f 100644
--- a/src/amdis/Assembler.hpp
+++ b/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>
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index ac0ccb05..2854493e 100644
--- a/src/amdis/Assembler.inc.hpp
+++ b/src/amdis/Assembler.inc.hpp
@@ -1,5 +1,7 @@
#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&&) {
diff --git a/src/amdis/GridFunctionOperator.hpp b/src/amdis/GridFunctionOperator.hpp
index 318506fa..167892ff 100644
--- a/src/amdis/GridFunctionOperator.hpp
+++ b/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);
}
diff --git a/src/amdis/LocalAssemblerBase.hpp b/src/amdis/LocalAssemblerBase.hpp
index 44e756a2..2504c42c 100644
--- a/src/amdis/LocalAssemblerBase.hpp
+++ b/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<
diff --git a/src/amdis/ProblemInstat.hpp b/src/amdis/ProblemInstat.hpp
index 2a676802..a61b9b3b 100644
--- a/src/amdis/ProblemInstat.hpp
+++ b/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;
}
diff --git a/src/amdis/ProblemInstat.inc.hpp b/src/amdis/ProblemInstat.inc.hpp
index 502f9d9a..5e93f0e8 100644
--- a/src/amdis/ProblemInstat.inc.hpp
+++ b/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
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index 92e85985..7b4b4ac2 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/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
diff --git a/src/amdis/common/Math.hpp b/src/amdis/common/Math.hpp
index 555d1585..9897ad0c 100644
--- a/src/amdis/common/Math.hpp
+++ b/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
diff --git a/src/amdis/common/Transposed.hpp b/src/amdis/common/Transposed.hpp
new file mode 100644
index 00000000..b590e3ed
--- /dev/null
+++ b/src/amdis/common/Transposed.hpp
@@ -0,0 +1,78 @@
+#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
diff --git a/src/amdis/linear_algebra/CMakeLists.txt b/src/amdis/linear_algebra/CMakeLists.txt
index ed57a338..ff79554f 100644
--- a/src/amdis/linear_algebra/CMakeLists.txt
+++ b/src/amdis/linear_algebra/CMakeLists.txt
@@ -1,6 +1,10 @@
#install headers
install(FILES
+ Common.hpp
+ DOFMatrixBase.hpp
+ DOFVectorBase.hpp
+ DOFVectorInterface.hpp
HierarchicWrapper.hpp
LinearSolverInterface.hpp
PreconditionerInterface.hpp
diff --git a/src/amdis/linear_algebra/DOFMatrixBase.hpp b/src/amdis/linear_algebra/DOFMatrixBase.hpp
index a37148e0..9080a167 100644
--- a/src/amdis/linear_algebra/DOFMatrixBase.hpp
+++ b/src/amdis/linear_algebra/DOFMatrixBase.hpp
@@ -1,8 +1,11 @@
#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
diff --git a/src/amdis/linear_algebra/DOFVectorBase.hpp b/src/amdis/linear_algebra/DOFVectorBase.hpp
index 68cb1c99..33fa8ea8 100644
--- a/src/amdis/linear_algebra/DOFVectorBase.hpp
+++ b/src/amdis/linear_algebra/DOFVectorBase.hpp
@@ -1,13 +1,10 @@
#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
diff --git a/src/amdis/linear_algebra/HierarchicWrapper.hpp b/src/amdis/linear_algebra/HierarchicWrapper.hpp
index 609cf475..79f1f4c1 100644
--- a/src/amdis/linear_algebra/HierarchicWrapper.hpp
+++ b/src/amdis/linear_algebra/HierarchicWrapper.hpp
@@ -1,6 +1,9 @@
#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());
diff --git a/src/amdis/linear_algebra/istl/CMakeLists.txt b/src/amdis/linear_algebra/istl/CMakeLists.txt
index 26238d75..f1e73c05 100644
--- a/src/amdis/linear_algebra/istl/CMakeLists.txt
+++ b/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)
diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
index 22568854..61837bf6 100644
--- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp
+++ b/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
{
- return matrix_->nonzeroes();
+ return matrix_.nonzeroes();
}
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 = (i == cIt.index() ? 1 : 0);
@@ -133,10 +120,12 @@ namespace AMDiS
}
private:
- ClonablePtr<BaseMatrix> matrix_;
- Dune::MatrixIndexSet occupationPattern_;
+ BaseMatrix matrix_;
bool initialized_ = false;
};
+ template <class RowBasisType, class ColBasisType, class ValueType = double>
+ using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, IstlMatrix<ValueType>>;
+
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DOFVector.hpp b/src/amdis/linear_algebra/istl/DOFVector.hpp
index af5690fb..a1234090 100644
--- a/src/amdis/linear_algebra/istl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/istl/DOFVector.hpp
@@ -3,9 +3,7 @@
#include <dune/istl/bvector.hh>
#include <amdis/Output.hpp>
-#include <amdis/common/ClonablePtr.hpp>
#include <amdis/linear_algebra/DOFVectorBase.hpp>
-#include <amdis/utility/MultiIndex.hpp>
namespace AMDiS
{
@@ -21,17 +19,10 @@ namespace AMDiS
using type = T;
};
- template <class BasisType, class ValueType = double>
- class DOFVector
- : public DOFVectorBase<BasisType>
+ template <class ValueType>
+ class IstlVector
{
- using Self = DOFVector;
- using Super = DOFVectorBase<BasisType>;
-
public:
- /// The type of the finite element space / basis
- using Basis = BasisType;
-
/// The type of the elements of the DOFVector
using value_type = typename BlockVectorType<ValueType>::type;
@@ -39,97 +30,61 @@ namespace AMDiS
using BaseVector = Dune::BlockVector<value_type>;
/// The index/size - type
- using size_type = typename Basis::size_type;
+ using size_type = typename BaseVector::size_type;
+ public:
/// Constructor. Constructs new BaseVector.
- DOFVector(BasisType const& basis)
- : Super(basis)
- , vector_(ClonablePtr<BaseVector>::make())
- {
- compress();
- }
-
- /// Constructor. Takes pointer of data-vector.
- DOFVector(BasisType const& basis, Dune::BlockVector<ValueType>& vector)
- : Super(basis)
- , vector_(vector)
- {}
+ IstlVector() = default;
/// Return the data-vector \ref vector
BaseVector const& vector() const
{
- return *vector_;
+ return vector_;
}
/// Return the data-vector \ref vector
BaseVector& vector()
{
- return *vector_;
+ return vector_;
}
- /// Resize the \ref vector to the size of the \ref basis.
- virtual void compress() override
+ /// Return the current size of the \ref vector_
+ size_type size() const
{
- if (vector_->size() != this->size()) {
- resize(this->size());
- *vector_ = value_type(0);
- }
+ return vector_.size();
}
- template <class SizeInfo>
- void resize(SizeInfo s)
+ /// Resize the \ref vector_ to the size \p s
+ void resize(size_type s)
{
- vector_->resize(size_type(s));
+ vector_.resize(s);
}
/// Access the entry \p i of the \ref vector with read-access.
- template <class Index>
- value_type const& operator[](Index idx) const
+ value_type const& operator[](size_type i) const
{
- size_type i = flatMultiIndex(idx);
- test_exit_dbg(i < vector_->size(),
- "Index ", i, " out of range [0,", vector_->size(), ")" );
- return (*vector_)[i];
+ test_exit_dbg(i < vector_.size(),
+ "Index ", i, " out of range [0,", vector_.size(), ")" );
+ return vector_[i];
}
/// Access the entry \p i of the \ref vector with write-access.
- template <class Index>
- value_type& operator[](Index idx)
- {
- size_type i = flatMultiIndex(idx);
- test_exit_dbg(i < vector_->size(),
- "Index ", i, " out of range [0,", vector_->size(), ")" );
- return (*vector_)[i];
- }
-
- /// Sets each DOFVector to the scalar \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator=(Scalar s)
- {
- (*vector_) = s;
- return *this;
- }
-
- /// Calls the copy assignment operator of the BaseVector \ref vector
- void copy(Self const& that)
+ value_type& operator[](size_type i)
{
- *vector_ = that.vector();
+ test_exit_dbg(i < vector_.size(),
+ "Index ", i, " out of range [0,", vector_.size(), ")" );
+ return vector_[i];
}
private:
- ClonablePtr<BaseVector> vector_;
+ BaseVector vector_;
};
-#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- // Deduction rule
- template <class Basis, class T>
- DOFVector(Basis const& basis, Dune::BlockVector<T>& coefficients)
- -> DOFVector<Basis, typename T::field_type>;
-#endif
+ template <class BasisType, class VectorType = double>
+ using DOFVector = DOFVectorBase<BasisType, IstlVector<VectorType>>;
/// Constructor a dofvector from given basis and name
template <class ValueType = double, class Basis>
@@ -139,12 +94,4 @@ namespace AMDiS
return {basis};
}
- /// Constructor a dofvector from given basis, name, and coefficient vector
- template <class Basis, class T>
- DOFVector<Basis, typename T::field_type>
- makeDOFVector(Basis const& basis, Dune::BlockVector<T>& coefficients)
- {
- return {basis, coefficients};
- }
-
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/CMakeLists.txt b/src/amdis/linear_algebra/mtl/CMakeLists.txt
index e5d7c4d7..e295f73f 100644
--- a/src/amdis/linear_algebra/mtl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/mtl/CMakeLists.txt
@@ -4,7 +4,6 @@ install(FILES
ITL_Preconditioner.hpp
ITL_Solver.hpp
KrylovRunner.hpp
- LinearSolver.hpp
Preconditioner.hpp
UmfpackRunner.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/mtl)
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index 6edf7b90..5d3824c1 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -11,7 +11,6 @@
#include <boost/numeric/mtl/utility/range_wrapper.hpp>
#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>
@@ -20,13 +19,9 @@ namespace AMDiS
{
/// \brief The basic container that stores a base matrix and a corresponding
/// row/column feSpace.
- template <class RowBasisType, class ColBasisType,
- class ValueType = double>
- class DOFMatrix
- : public DOFMatrixBase<RowBasisType, ColBasisType>
+ template <class ValueType>
+ class MtlMatrix
{
- using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
-
public:
/// The matrix type of the underlying base matrix
using BaseMatrix = mtl::compressed2D<ValueType>;
@@ -34,63 +29,55 @@ namespace AMDiS
/// The type of the elements of the DOFMatrix
using value_type = ValueType;
+ /// The index/size - type
+ using size_type = typename BaseMatrix::size_type;
+
/// The type of the inserter used when filling the matrix. NOTE: need not be public
using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
public:
/// Constructor. Constructs new BaseMatrix.
- DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
- : Super(rowBasis, colBasis)
- , matrix_(ClonablePtr<BaseMatrix>::make())
- {}
-
- /// Constructor. Takes a reference to a base matrix
- DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
- : Super(rowBasis, colBasis)
- , matrix_(matrix)
- {}
+ MtlMatrix() = default;
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix()
{
assert( !inserter_ );
- return *matrix_;
+ return matrix_;
}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix const& matrix() const
{
assert( !inserter_ );
- return *matrix_;
+ return matrix_;
}
/// \brief Returns an update-proxy of the inserter, to inserte/update a value at
/// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
/// be created by \ref init and must be closed by \ref finish after insertion.
- template <class Index>
- auto operator()(Index row, Index col)
+ void insert(size_type r, size_type c, value_type const& value)
{
- auto r = flatMultiIndex(row), c = flatMultiIndex(col);
test_exit_dbg(inserter_, "Inserter not initilized!");
- test_exit_dbg(r < this->rows() && c < this->cols(),
- "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
- return (*inserter_)[r][c];
+ test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
+ "Indices out of range [0,", num_rows(matrix_), ")x[0,", num_cols(matrix_), ")" );
+ (*inserter_)[r][c] += value;
}
/// Create inserter. Assumes that no inserter is currently active on this matrix.
- void init(bool prepareForInsertion)
+ template <class RowBasis, class ColBasis>
+ void init(RowBasis const& rowBasis, ColBasis const& colBasis,
+ bool prepareForInsertion)
{
test_exit(!inserter_, "Matrix already in insertion mode!");
calculateSlotSize();
- matrix_->change_dim(this->rows(), this->cols());
- test_exit(num_rows(*matrix_) == this->rows() && num_cols(*matrix_) == this->cols(),
- "Wrong dimensions in matrix");
+ matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
if (prepareForInsertion) {
- set_to_zero(*matrix_);
- inserter_ = new Inserter(*matrix_, slotSize_);
+ set_to_zero(matrix_);
+ inserter_ = new Inserter(matrix_, slotSize_);
}
}
@@ -105,7 +92,7 @@ namespace AMDiS
/// Return the number of nonzeros in the matrix
std::size_t nnz() const
{
- return matrix_->nnz();
+ return matrix_.nnz();
}
/// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
@@ -113,12 +100,12 @@ namespace AMDiS
auto applyDirichletBC(std::vector<char> const& dirichletNodes)
{
// Define the property maps
- auto row = mtl::mat::row_map(*matrix_);
- auto col = mtl::mat::col_map(*matrix_);
- auto value = mtl::mat::value_map(*matrix_);
+ auto row = mtl::mat::row_map(matrix_);
+ auto col = mtl::mat::col_map(matrix_);
+ auto value = mtl::mat::value_map(matrix_);
// iterate over the matrix
- for (auto r : mtl::rows_of(*matrix_)) { // rows of the matrix
+ for (auto r : mtl::rows_of(matrix_)) { // rows of the matrix
if (dirichletNodes[r.value()]) {
for (auto i : mtl::nz_of(r)) { // non-zeros within
// set identity row
@@ -136,8 +123,8 @@ namespace AMDiS
{
slotSize_ = 0;
- if (num_rows(*matrix_) != 0)
- slotSize_ = int(double(matrix_->nnz()) / num_rows(*matrix_) * 1.2);
+ if (num_rows(matrix_) != 0)
+ slotSize_ = int(double(matrix_.nnz()) / num_rows(matrix_) * 1.2);
if (slotSize_ < MIN_NNZ_PER_ROW)
slotSize_ = MIN_NNZ_PER_ROW;
}
@@ -147,7 +134,7 @@ namespace AMDiS
static constexpr int MIN_NNZ_PER_ROW = 50;
/// The data-matrix (can hold a new BaseMatrix or a pointer to external data
- ClonablePtr<BaseMatrix> matrix_;
+ BaseMatrix matrix_;
/// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
Inserter* inserter_ = nullptr;
@@ -156,4 +143,7 @@ namespace AMDiS
int slotSize_ = MIN_NNZ_PER_ROW;
};
+ template <class RowBasisType, class ColBasisType, class ValueType = double>
+ using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, MtlMatrix<ValueType>>;
+
} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp
index 212685d9..01842741 100644
--- a/src/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -3,125 +3,79 @@
#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/DOFVectorBase.hpp>
-#include <amdis/utility/MultiIndex.hpp>
namespace AMDiS
{
/// The basic container that stores a base vector and a corresponding basis
- template <class BasisType, class ValueType = double>
- class DOFVector
- : public DOFVectorBase<BasisType>
+ template <class ValueType>
+ class MtlVector
{
- using Self = DOFVector;
- using Super = DOFVectorBase<BasisType>;
-
public:
- /// The type of the \ref basis
- using Basis = BasisType;
-
/// The type of the base vector
using BaseVector = mtl::vec::dense_vector<ValueType>;
/// The index/size - type
- using size_type = typename BasisType::size_type;
+ using size_type = typename BaseVector::size_type;
/// The type of the elements of the DOFVector
using value_type = ValueType;
+ public:
/// Constructor. Constructs new BaseVector.
- DOFVector(BasisType const& basis)
- : Super(basis)
- , vector_(ClonablePtr<BaseVector>::make())
- {
- compress();
- }
-
- /// Constructor. Takes reference to existing BaseVector
- DOFVector(BasisType const& basis, BaseVector& vector)
- : Super(basis)
- , vector_(vector)
- {}
+ MtlVector() = default;
/// Return the data-vector \ref vector_
BaseVector const& vector() const
{
- return *vector_;
+ return vector_;
}
/// Return the data-vector \ref vector_
BaseVector& vector()
{
- return *vector_;
+ return vector_;
}
- /// Resize the \ref vector to the size of the \ref basis.
- virtual void compress() override
+ /// Return the current size of the \ref vector_
+ size_type size() const
{
- if (mtl::vec::size(*vector_) != this->size()) {
- resize(this->size());
- *vector_ = value_type(0);
- }
+ return mtl::vec::size(vector_);
}
- template <class SizeInfo>
- void resize(SizeInfo s)
+ /// Resize the \ref vector_ to the size \p s
+ void resize(size_type s)
{
- vector_->change_dim(size_type(s));
+ vector_.change_dim(s);
}
/// Access the entry \p i of the \ref vector with read-access.
- template <class Index>
- value_type const& operator[](Index idx) const
+ value_type const& operator[](size_type i) const
{
- size_type i = flatMultiIndex(idx);
- test_exit_dbg(i < mtl::vec::size(*vector_),
- "Index ", i, " out of range [0,", mtl::vec::size(*vector_), ")" );
- return (*vector_)[i];
+ test_exit_dbg(i < mtl::vec::size(vector_),
+ "Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
+ return vector_[i];
}
/// Access the entry \p i of the \ref vector with write-access.
- template <class Index>
- value_type& operator[](Index idx)
+ value_type& operator[](size_type i)
{
- size_type i = flatMultiIndex(idx);
- test_exit_dbg(i < mtl::vec::size(*vector_),
- "Index ", i, " out of range [0,", mtl::vec::size(*vector_), ")" );
- return (*vector_)[i];
- }
-
- // /// Scale each DOFVector by the factor \p s.
- // template <class Scalar>
- // std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- // operator*=(Scalar s)
- // {
- // (*vector_) *= s;
- // return *this;
- // }
-
- /// Sets each DOFVector to the scalar \p s.
- template <class Scalar>
- std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
- operator=(Scalar s)
- {
- (*vector_) = s;
- return *this;
- }
-
- /// Calls the copy assignment operator of the BaseVector \ref vector
- void copy(Self const& that)
- {
- *vector_ = *that.vector_;
+ test_exit_dbg(i < mtl::vec::size(vector_),
+ "Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
+ return vector_[i];
}
private:
/// The data-vector (can hold a new BaseVector or a pointer to external data
- ClonablePtr<BaseVector> vector_;
+ BaseVector vector_;
};
+
+ template <class BasisType, class VectorType = double>
+ using DOFVector = DOFVectorBase<BasisType, MtlVector<VectorType>>;
+
+
/// Constructor a dofvector from given basis and name
template <class ValueType = double, class Basis>
DOFVector<Basis, ValueType>
@@ -130,12 +84,4 @@ namespace AMDiS
return {basis};
}
- /// Constructor a dofvector from given basis, name, and coefficient vector
- template <class Basis, class ValueType>
- DOFVector<Basis, ValueType>
- makeDOFVector(Basis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
- {
- return {basis, coefficients};
- }
-
} // end namespace AMDiS
--
GitLab