diff --git a/src/amdis/linear_algebra/Common.hpp b/src/amdis/linear_algebra/Common.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e7bf16d923776bf233d6114d4706986d9929d9b3 --- /dev/null +++ b/src/amdis/linear_algebra/Common.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include <cstddef> + +namespace AMDiS +{ + template <class T> + struct Triplet + { + std::size_t row, cols; + T value; + }; + +} // end namespace AMDiS diff --git a/src/amdis/linear_algebra/DOFMatrixBase.hpp b/src/amdis/linear_algebra/DOFMatrixBase.hpp index e7bf16d923776bf233d6114d4706986d9929d9b3..a37148e092fd989c879145d3058d9733d1ab7a6f 100644 --- a/src/amdis/linear_algebra/DOFMatrixBase.hpp +++ b/src/amdis/linear_algebra/DOFMatrixBase.hpp @@ -1,14 +1,54 @@ #pragma once -#include <cstddef> - namespace AMDiS { - template <class T> - struct Triplet + template <class RowBasisType, class ColBasisType> + class DOFMatrixBase { - std::size_t row, cols; - T value; + public: + /// 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; + + /// The index/size - type + using size_type = typename RowBasis::size_type; + + public: + /// Constructor. Constructs new BaseVector. + DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis) + : rowBasis_(&rowBasis) + , colBasis_(&colBasis) + {} + + /// Return the row-basis \ref rowBasis of the matrix + RowBasis const& rowBasis() const + { + return *rowBasis_; + } + + /// Return the col-basis \ref colBasis of the matrix + ColBasis const& colBasis() const + { + return *colBasis_; + } + + /// Return the size of the \ref rowBasis_ + size_type rows() const + { + return rowBasis_->dimension(); + } + + /// Return the size of the \ref colBasis_ + size_type cols() const + { + return colBasis_->dimension(); + } + + protected: + RowBasis const* rowBasis_; + ColBasis const* colBasis_; }; } // end namespace AMDiS diff --git a/src/amdis/linear_algebra/RunnerInterface.hpp b/src/amdis/linear_algebra/RunnerInterface.hpp index 082c023452f2b7a58ebf4a69018045b67ab65fda..ae7fc03da6e7f36ff8de4cd2a474f6ad93e1b889 100644 --- a/src/amdis/linear_algebra/RunnerInterface.hpp +++ b/src/amdis/linear_algebra/RunnerInterface.hpp @@ -34,8 +34,8 @@ namespace AMDiS return 0; } - virtual std::shared_ptr<PreconBase> getLeftPrecon() { return {}; } - virtual std::shared_ptr<PreconBase> getRightPrecon() { return {}; } + virtual std::shared_ptr<PreconBase> leftPrecon() { return {}; } + virtual std::shared_ptr<PreconBase> rightPrecon() { return {}; } }; } // end namespace AMDiS diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp index 09739652d3b2649e12946e69876bc6ce5b64428e..225688548a2ff2ce26bcdabb72daf420ada07fe3 100644 --- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp +++ b/src/amdis/linear_algebra/istl/DOFMatrix.hpp @@ -1,5 +1,6 @@ #pragma once +#include <list> #include <string> #include <memory> @@ -8,6 +9,7 @@ #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> @@ -28,50 +30,30 @@ namespace AMDiS template <class RowBasisType, class ColBasisType, class ValueType = double> class DOFMatrix + : public DOFMatrixBase<RowBasisType, ColBasisType> { - public: - /// 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 Super = DOFMatrixBase<RowBasisType, ColBasisType>; + public: /// The type of the elements of the DOFMatrix using value_type = typename BlockMatrixType<ValueType>::type; /// The matrix type of the underlying base matrix using BaseMatrix = Dune::BCRSMatrix<value_type>; - /// The index/size - type - using size_type = typename BaseMatrix::size_type; - public: /// Constructor. Constructs new BaseVector. - DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis) - : rowBasis_(&rowBasis) - , colBasis_(&colBasis) + DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis) + : Super(rowBasis, colBasis) , matrix_(ClonablePtr<BaseMatrix>::make()) {} /// Constructor. Takes reference to a base matrix. - DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix) - : rowBasis_(&rowBasis) - , colBasis_(&colBasis) + DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix) + : Super(rowBasis, colBasis) , matrix_(matrix) {} - /// Return the row-basis \ref rowBasis of the matrix - RowBasis const& rowBasis() const - { - return *rowBasis_; - } - - /// Return the col-basis \ref colBasis of the matrix - ColBasis const& colBasis() const - { - return *colBasis_; - } - /// Return the data-vector \ref vector BaseMatrix const& matrix() const { @@ -84,50 +66,28 @@ namespace AMDiS return *matrix_; } - /// Return the size of the \ref feSpace - size_type rows() const - { - return rowBasis_->size(); - } - - size_type cols() const - { - return colBasis_->size(); - } - - - /// Access the entry \p i of the \ref vector with read-access. - template <class Index> - value_type const& operator()(Index row, Index col) const - { - 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. template <class Index> value_type& operator()(Index row, Index col) { - size_type r = flatMultiIndex(row), c = flatMultiIndex(col); + auto 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(), ")" ); + test_exit_dbg( r < this->rows() && c < this->cols() , + "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" ); return (*matrix_)[r][c]; } /// create occupation pattern and apply it to the matrix void init(bool prepareForInsertion) { - occupationPattern_.resize(rowBasis_->size(), colBasis_->size()); + occupationPattern_.resize(this->rows(), this->cols()); // A loop over all elements of the grid - auto rowLocalView = rowBasis_->localView(); - auto colLocalView = colBasis_->localView(); + auto rowLocalView = this->rowBasis_->localView(); + auto colLocalView = this->colBasis_->localView(); - for (const auto& element : elements(rowBasis_->gridView())) { + for (const auto& element : elements(this->rowBasis_->gridView())) { rowLocalView.bind(element); colLocalView.bind(element); @@ -169,14 +129,10 @@ namespace AMDiS } } - std::list<Triplet<value_type>> columns; // symmetric dbc not yet implemented - return columns; + return std::list<Triplet<value_type>>{}; } private: - RowBasis const* rowBasis_; - ColBasis const* colBasis_; - ClonablePtr<BaseMatrix> matrix_; Dune::MatrixIndexSet occupationPattern_; diff --git a/src/amdis/linear_algebra/istl/DOFVector.hpp b/src/amdis/linear_algebra/istl/DOFVector.hpp index ddcfd1516bb746e92feb275bf58c0b9217950933..af5690fb9f9c5f6ff81a2a5a7d2d88f6623cbcf0 100644 --- a/src/amdis/linear_algebra/istl/DOFVector.hpp +++ b/src/amdis/linear_algebra/istl/DOFVector.hpp @@ -123,4 +123,28 @@ namespace AMDiS ClonablePtr<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 + + + /// Constructor a dofvector from given basis and name + template <class ValueType = double, class Basis> + DOFVector<Basis, ValueType> + makeDOFVector(Basis const& basis) + { + 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/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp index 67c702b32c0ed743f16fda44e65c5ffaaf59a35c..6edf7b909ee2220e470c249c6c6287f022046f15 100644 --- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp +++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp @@ -12,6 +12,7 @@ #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> @@ -22,20 +23,14 @@ namespace AMDiS template <class RowBasisType, class ColBasisType, class ValueType = double> class DOFMatrix + : public DOFMatrixBase<RowBasisType, ColBasisType> { - public: - /// 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 Super = DOFMatrixBase<RowBasisType, ColBasisType>; + public: /// The matrix type of the underlying base matrix using BaseMatrix = mtl::compressed2D<ValueType>; - /// The index/size - type - using size_type = typename BaseMatrix::size_type; - /// The type of the elements of the DOFMatrix using value_type = ValueType; @@ -44,31 +39,17 @@ namespace AMDiS public: /// Constructor. Constructs new BaseMatrix. - DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis) - : rowBasis_(&rowBasis) - , colBasis_(&colBasis) + DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis) + : Super(rowBasis, colBasis) , matrix_(ClonablePtr<BaseMatrix>::make()) {} /// Constructor. Takes a reference to a base matrix - DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix) - : rowBasis_(&rowBasis) - , colBasis_(&colBasis) + DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix) + : Super(rowBasis, colBasis) , matrix_(matrix) {} - /// Return the row-basis \ref rowBasis_ of the matrix - RowBasis const& rowBasis() const - { - return *rowBasis_; - } - - /// Return the col-basis \ref colBasis_ of the matrix - ColBasis const& colBasis() const - { - return *colBasis_; - } - /// Return a reference to the data-matrix \ref matrix BaseMatrix& matrix() { @@ -83,18 +64,6 @@ namespace AMDiS return *matrix_; } - /// Return the size of the row \ref feSpace - size_type rows() const - { - return rowBasis_->size(); - } - - /// Return the size of the column \ref feSpace - size_type cols() const - { - return colBasis_->size(); - } - /// \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 @@ -102,10 +71,10 @@ namespace AMDiS template <class Index> auto operator()(Index row, Index col) { - size_type r = flatMultiIndex(row), c = flatMultiIndex(col); + auto r = flatMultiIndex(row), c = flatMultiIndex(col); test_exit_dbg(inserter_, "Inserter not initilized!"); - test_exit_dbg(r < rows() && c < cols(), - "Indices out of range [0,", rows(), ")x[0,", cols(), ")" ); + test_exit_dbg(r < this->rows() && c < this->cols(), + "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" ); return (*inserter_)[r][c]; } @@ -116,8 +85,8 @@ namespace AMDiS test_exit(!inserter_, "Matrix already in insertion mode!"); calculateSlotSize(); - matrix_->change_dim(rowBasis_->dimension(), colBasis_->dimension()); - test_exit(num_rows(*matrix_) == rowBasis_->dimension() && num_cols(*matrix_) == colBasis_->dimension(), + matrix_->change_dim(this->rows(), this->cols()); + test_exit(num_rows(*matrix_) == this->rows() && num_cols(*matrix_) == this->cols(), "Wrong dimensions in matrix"); if (prepareForInsertion) { set_to_zero(*matrix_); @@ -143,30 +112,13 @@ namespace AMDiS /// a one on the diagonal of the matrix. auto applyDirichletBC(std::vector<char> const& dirichletNodes) { - std::list<Triplet<value_type>> columns; - // 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_); // iterate over the matrix -#if 0 - for (auto r : mtl::major_of(*matrix)) { // rows or columns - for (auto i : mtl::nz_of(r)) { // non-zeros within - if (dirichletNodes[row(i)]) { - // set identity row - value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) ); - } - else if (dirichletNodes[col(i)]) { - columns.push_back({row(i), col(i), value(i)}); - value(i, value_type(0)); - } - } - } -#endif - - for (auto r : mtl::rows_of(*matrix_)) { // rows or columns + 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 @@ -175,67 +127,8 @@ namespace AMDiS } } - return columns; - } - -#if 0 - void applyPeriodicBC(std::vector<char> const& periodicNodes, - std::map<size_t, size_t> const& association, bool applyRow, bool applyCol) - { - bool apply = applyRow && applyCol; - - // 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); - - std::vector<std::map<size_t, std::list<value_type>>> row_values(num_cols(*matrix)); - std::vector<std::map<size_t, std::list<value_type>>> col_values(num_rows(*matrix)); - std::vector<char> dualNodes(periodicNodes.size(), 0); - - // iterate over the matrix to collect the values and erase the source-row/col - for (auto r : mtl::major_of(*matrix)) { // rows or columns - for (auto i : mtl::nz_of(r)) { // non-zeros within - if (applyRow && periodicNodes[row(i)]) { - row_values[col(i)][association[row(i)]].push_back(value(i)); - value(i, value_type(0)); - dualNodes[association[row(i)]] = 1; - } else if (applyCol && periodicNodes[col(i)]) { - col_values[row(i)][association[col(i)]].push_back(value(i)); - value(i, value_type(0)); - } - } - } - - // TODO: use inserter for assignment of values!!! - // iterate over the matrix to assign the values - for (auto r : mtl::major_of(*matrix)) { // rows or columns - for (auto i : mtl::nz_of(r)) { // non-zeros within - if (applyRow && dualNodes[row(i)]) { - value(i, std::accumulate(row_values[col(i)][row(i)].begin(), - row_values[col(i)][row(i)].end(), - value(i)) ); - } - if (applyCol && dualNodes[col(i)]) { - value(i, std::accumulate(col_values[row(i)][col(i)].begin(), - col_values[row(i)][col(i)].end(), - value(i)) ); - } - } - } - - // assign values 1, -1 to diagonal and associated entries - if (apply) { - Inserter ins(*matrix); - for (auto const& a : association) { - ins[a.first][a.first] << value_type( 1); - ins[a.second][a.second] << value_type( 1); - ins[a.first][a.second] << value_type(-1); - ins[a.second][a.first] << value_type(-1); - } - } + return std::list<Triplet<value_type>>{}; } -#endif protected: // Estimates the slot-size used in the inserter to reserve enough space per row. @@ -253,12 +146,6 @@ namespace AMDiS /// The minimal number of nonzeros per row static constexpr int MIN_NNZ_PER_ROW = 50; - /// The finite element space / basis of the row - RowBasis const* rowBasis_; - - /// The finite element space / basis of the column - ColBasis const* colBasis_; - /// The data-matrix (can hold a new BaseMatrix or a pointer to external data ClonablePtr<BaseMatrix> matrix_; diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp index 7287b807d33fe77f7770839cad744ffe95c22f48..212685d95d65e4356b092b0c0be850892cc806c4 100644 --- a/src/amdis/linear_algebra/mtl/DOFVector.hpp +++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp @@ -123,7 +123,7 @@ namespace AMDiS }; /// Constructor a dofvector from given basis and name - template <class Basis, class ValueType = double> + template <class ValueType = double, class Basis> DOFVector<Basis, ValueType> makeDOFVector(Basis const& basis) { diff --git a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp index a226241f895e2eba3597577c6e53a1843884936f..2927909c207baba7066ab66cba64476561ef54c2 100644 --- a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp +++ b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp @@ -39,14 +39,14 @@ namespace AMDiS Parameters::get(prefix + "->print cycle", printCycle_); } - /// Implements \ref RunnerInterface::getLeftPrecon(), Returns \ref L_ - virtual std::shared_ptr<PreconBase> getLeftPrecon() override + /// Implements \ref RunnerInterface::lLeftPrecon(), Returns \ref L_ + virtual std::shared_ptr<PreconBase> leftPrecon() override { return L_; } - /// Implements \ref RunnerInterface::getRightPrecon(), Returns \ref R_ - virtual std::shared_ptr<PreconBase> getRightPrecon() override + /// Implements \ref RunnerInterface::rightPrecon(), Returns \ref R_ + virtual std::shared_ptr<PreconBase> rightPrecon() override { return R_; } @@ -107,14 +107,14 @@ namespace AMDiS void initPrecon(std::string prefix) { // Creator for the left preconditioner - std::string leftPrecon = "", rightPrecon = ""; - Parameters::get(prefix + "->left precon", leftPrecon); - Parameters::get(prefix + "->right precon", rightPrecon); + std::string leftPreconName = "", rightPreconName = ""; + Parameters::get(prefix + "->left precon", leftPreconName); + Parameters::get(prefix + "->right precon", rightPreconName); auto leftCreator - = CreatorMap<PreconBase>::getCreator(leftPrecon, prefix + "->left precon"); + = CreatorMap<PreconBase>::getCreator(leftPreconName, prefix + "->left precon"); auto rightCreator - = CreatorMap<PreconBase>::getCreator(rightPrecon, prefix + "->right precon"); + = CreatorMap<PreconBase>::getCreator(rightPreconName, prefix + "->right precon"); L_ = leftCreator->create(); R_ = rightCreator->create(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b5eb3aaba4cf03f3d7803773612b3456642963a3..afb4a3ccb701712555aa9c954ca56d6c346907e1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -8,6 +8,9 @@ dune_add_test(SOURCES ClonablePtrTest.cpp dune_add_test(SOURCES ConceptsTest.cpp LINK_LIBRARIES amdis) +dune_add_test(SOURCES DOFVectorTest.cpp + LINK_LIBRARIES amdis) + dune_add_test(SOURCES ExpressionsTest.cpp LINK_LIBRARIES amdis CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d") diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a1706a5b90537cb71dbd8bb2ec952e1fc5f5285e --- /dev/null +++ b/test/DOFVectorTest.cpp @@ -0,0 +1,56 @@ +#include <dune/common/filledarray.hh> +#include <dune/grid/yaspgrid.hh> +#include <dune/functions/functionspacebases/compositebasis.hh> +#include <dune/functions/functionspacebases/powerbasis.hh> +#include <dune/functions/functionspacebases/lagrangebasis.hh> + +#include <amdis/LinearAlgebra.hpp> + +#include "Tests.hpp" + +using namespace AMDiS; +using namespace Dune::Functions::BasisFactory; + +template <class B, class T> +void test_dofvector(B const& basis, DOFVector<B,T>& vec) +{ + AMDIS_TEST( vec.size() == basis.dimension() ); + vec.compress(); + + vec = T(0); + vec[0] = T(1); + auto m0 = std::min_element(std::begin(vec.vector()), std::end(vec.vector())); + auto m1 = std::max_element(std::begin(vec.vector()), std::end(vec.vector())); + AMDIS_TEST( *m0 == T(0) ); + AMDIS_TEST( *m1 == T(1) ); +} + +int main() +{ + // create grid + Dune::FieldVector<double, 2> L; L = 1.0; + auto s = Dune::filledArray<2>(1); + Dune::YaspGrid<2> grid(L, s); + + // create basis + auto basis = makeBasis(grid.leafGridView(), + composite(power<2>(lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic())); + + using Basis = decltype(basis); + DOFVector<Basis> vec1(basis); + test_dofvector(basis, vec1); + + DOFVector<Basis, float> vec2(basis); + test_dofvector(basis, vec2); + + auto vec3 = makeDOFVector(basis); + test_dofvector(basis, vec3); + + auto vec4 = makeDOFVector<float>(basis); + test_dofvector(basis, vec4); + +#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + DOFVector vec5(basis); + test_dofvector(basis, vec5); +#endif +} \ No newline at end of file