diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 8b43adaffeff0a9d9723d661d76e42973ffb62dd..e487c99f6b2814ea3195cd9de0ca5da46fa843f1 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 ac0ccb050f6ac3ace69290cae99331a663c0b948..2854493e03077196120f629a32ecb0fccbca0308 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 318506fa12473ce9c561258b8a9f781d4cbb2975..167892ffa6c2c05ef08667453d05cb9ed154e319 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 44e756a240510764a02121d99b93cbb63a1aa5e1..2504c42c1edb863fa8992202e055829d2b0a6c73 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 2a67680284ac7f7943224caaec4a2c9b13a12aae..a61b9b3b55b5e88736fb289a74e5b8b4817c6c83 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 502f9d9a9bc36bd98eabad835a1b9f8b53e2d778..5e93f0e889caf94fec46c4bd1c18675e98b71f5b 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 92e85985864553ed516303d2dad36a3edec12d01..7b4b4ac25b8cc8d50d016821c0f0c2bfb9565e00 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 555d1585c0cbb0cb97bd83e727f7afa681845834..9897ad0c7898ec30295d730648b2ad83a7d420af 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 0000000000000000000000000000000000000000..b590e3ed2f9c68d786c3aaed71b8a5c82fa76568
--- /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 ed57a3387c17abd2a34f6750aecf0e353d02806f..ff79554f15af37891a18244704545418050e4a33 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 a37148e092fd989c879145d3058d9733d1ab7a6f..9080a167cf70a842aaef3dbad509b17db985e851 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 68cb1c990c2cb8555db5b40571101ade7df8eb38..33fa8ea8277796806c2c4dbebd00de08fff58c0c 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 609cf4750388f4ed95b3df11e06ebcb08ffbef73..79f1f4c1b576b65bcffe2787fd35cf0deb536a6e 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 26238d75de2f02d584150a6e059b496557369d57..f1e73c05ad38ceff6076ef8155c52301e49c5f6b 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 225688548a2ff2ce26bcdabb72daf420ada07fe3..61837bf6083cdb3830c306245e1edcff9880edfd 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 af5690fb9f9c5f6ff81a2a5a7d2d88f6623cbcf0..a1234090740cb542c806b6e13b54bf5817ff1c2a 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 e5d7c4d7641ac63672f3b2068cebe024e5901cb6..e295f73f7767c0b95987825dd993f89720bc30f5 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 6edf7b909ee2220e470c249c6c6287f022046f15..5d3824c10d697d0d442f7714904e57f3c0b6a69c 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 212685d95d65e4356b092b0c0be850892cc806c4..01842741a5ffae416f7fe4866204bab847b3072a 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