Commit 8521711f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

added DOFVector test

parent 4a878787
#pragma once
#include <cstddef>
namespace AMDiS
{
template <class T>
struct Triplet
{
std::size_t row, cols;
T value;
};
} // end namespace AMDiS
#pragma once #pragma once
#include <cstddef>
namespace AMDiS namespace AMDiS
{ {
template <class T> template <class RowBasisType, class ColBasisType>
struct Triplet class DOFMatrixBase
{ {
std::size_t row, cols; public:
T value; /// 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 } // end namespace AMDiS
...@@ -34,8 +34,8 @@ namespace AMDiS ...@@ -34,8 +34,8 @@ namespace AMDiS
return 0; return 0;
} }
virtual std::shared_ptr<PreconBase> getLeftPrecon() { return {}; } virtual std::shared_ptr<PreconBase> leftPrecon() { return {}; }
virtual std::shared_ptr<PreconBase> getRightPrecon() { return {}; } virtual std::shared_ptr<PreconBase> rightPrecon() { return {}; }
}; };
} // end namespace AMDiS } // end namespace AMDiS
#pragma once #pragma once
#include <list>
#include <string> #include <string>
#include <memory> #include <memory>
...@@ -8,6 +9,7 @@ ...@@ -8,6 +9,7 @@
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp> #include <amdis/common/ClonablePtr.hpp>
#include <amdis/linear_algebra/Common.hpp>
#include <amdis/linear_algebra/DOFMatrixBase.hpp> #include <amdis/linear_algebra/DOFMatrixBase.hpp>
#include <amdis/utility/MultiIndex.hpp> #include <amdis/utility/MultiIndex.hpp>
...@@ -28,50 +30,30 @@ namespace AMDiS ...@@ -28,50 +30,30 @@ namespace AMDiS
template <class RowBasisType, class ColBasisType, template <class RowBasisType, class ColBasisType,
class ValueType = double> class ValueType = double>
class DOFMatrix class DOFMatrix
: public DOFMatrixBase<RowBasisType, ColBasisType>
{ {
public: using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
/// 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;
public:
/// The type of the elements of the DOFMatrix /// The type of the elements of the DOFMatrix
using value_type = typename BlockMatrixType<ValueType>::type; using value_type = typename BlockMatrixType<ValueType>::type;
/// The matrix type of the underlying base matrix /// The matrix type of the underlying base matrix
using BaseMatrix = Dune::BCRSMatrix<value_type>; using BaseMatrix = Dune::BCRSMatrix<value_type>;
/// The index/size - type
using size_type = typename BaseMatrix::size_type;
public: public:
/// Constructor. Constructs new BaseVector. /// Constructor. Constructs new BaseVector.
DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis) DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
: rowBasis_(&rowBasis) : Super(rowBasis, colBasis)
, colBasis_(&colBasis)
, matrix_(ClonablePtr<BaseMatrix>::make()) , matrix_(ClonablePtr<BaseMatrix>::make())
{} {}
/// Constructor. Takes reference to a base matrix. /// Constructor. Takes reference to a base matrix.
DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix) DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
: rowBasis_(&rowBasis) : Super(rowBasis, colBasis)
, colBasis_(&colBasis)
, matrix_(matrix) , 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 /// Return the data-vector \ref vector
BaseMatrix const& matrix() const BaseMatrix const& matrix() const
{ {
...@@ -84,50 +66,28 @@ namespace AMDiS ...@@ -84,50 +66,28 @@ namespace AMDiS
return *matrix_; 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. /// Access the entry \p i of the \ref vector with write-access.
template <class Index> template <class Index>
value_type& operator()(Index row, Index col) 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( initialized_, "Occupation pattern not initialized!");
test_exit_dbg( r < rows() && c < cols() , test_exit_dbg( r < this->rows() && c < this->cols() ,
"Indices out of range [0,", rows(), ")x[0,", cols(), ")" ); "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
return (*matrix_)[r][c]; return (*matrix_)[r][c];
} }
/// create occupation pattern and apply it to the matrix /// create occupation pattern and apply it to the matrix
void init(bool prepareForInsertion) void init(bool prepareForInsertion)
{ {
occupationPattern_.resize(rowBasis_->size(), colBasis_->size()); occupationPattern_.resize(this->rows(), this->cols());
// A loop over all elements of the grid // A loop over all elements of the grid
auto rowLocalView = rowBasis_->localView(); auto rowLocalView = this->rowBasis_->localView();
auto colLocalView = colBasis_->localView(); auto colLocalView = this->colBasis_->localView();
for (const auto& element : elements(rowBasis_->gridView())) { for (const auto& element : elements(this->rowBasis_->gridView())) {
rowLocalView.bind(element); rowLocalView.bind(element);
colLocalView.bind(element); colLocalView.bind(element);
...@@ -169,14 +129,10 @@ namespace AMDiS ...@@ -169,14 +129,10 @@ namespace AMDiS
} }
} }
std::list<Triplet<value_type>> columns; // symmetric dbc not yet implemented return std::list<Triplet<value_type>>{};
return columns;
} }
private: private:
RowBasis const* rowBasis_;
ColBasis const* colBasis_;
ClonablePtr<BaseMatrix> matrix_; ClonablePtr<BaseMatrix> matrix_;
Dune::MatrixIndexSet occupationPattern_; Dune::MatrixIndexSet occupationPattern_;
......
...@@ -123,4 +123,28 @@ namespace AMDiS ...@@ -123,4 +123,28 @@ namespace AMDiS
ClonablePtr<BaseVector> vector_; 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 } // end namespace AMDiS
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ClonablePtr.hpp> #include <amdis/common/ClonablePtr.hpp>
#include <amdis/linear_algebra/Common.hpp>
#include <amdis/linear_algebra/DOFMatrixBase.hpp> #include <amdis/linear_algebra/DOFMatrixBase.hpp>
#include <amdis/utility/MultiIndex.hpp> #include <amdis/utility/MultiIndex.hpp>
...@@ -22,20 +23,14 @@ namespace AMDiS ...@@ -22,20 +23,14 @@ namespace AMDiS
template <class RowBasisType, class ColBasisType, template <class RowBasisType, class ColBasisType,
class ValueType = double> class ValueType = double>
class DOFMatrix class DOFMatrix
: public DOFMatrixBase<RowBasisType, ColBasisType>
{ {
public: using Super = DOFMatrixBase<RowBasisType, ColBasisType>;
/// 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;
public:
/// The matrix type of the underlying base matrix /// The matrix type of the underlying base matrix
using BaseMatrix = mtl::compressed2D<ValueType>; using BaseMatrix = mtl::compressed2D<ValueType>;
/// The index/size - type
using size_type = typename BaseMatrix::size_type;
/// The type of the elements of the DOFMatrix /// The type of the elements of the DOFMatrix
using value_type = ValueType; using value_type = ValueType;
...@@ -44,31 +39,17 @@ namespace AMDiS ...@@ -44,31 +39,17 @@ namespace AMDiS
public: public:
/// Constructor. Constructs new BaseMatrix. /// Constructor. Constructs new BaseMatrix.
DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis) DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis)
: rowBasis_(&rowBasis) : Super(rowBasis, colBasis)
, colBasis_(&colBasis)
, matrix_(ClonablePtr<BaseMatrix>::make()) , matrix_(ClonablePtr<BaseMatrix>::make())
{} {}
/// Constructor. Takes a reference to a base matrix /// Constructor. Takes a reference to a base matrix
DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix) DOFMatrix(RowBasisType const& rowBasis, ColBasisType const& colBasis, BaseMatrix& matrix)
: rowBasis_(&rowBasis) : Super(rowBasis, colBasis)
, colBasis_(&colBasis)
, matrix_(matrix) , 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 /// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix() BaseMatrix& matrix()
{ {
...@@ -83,18 +64,6 @@ namespace AMDiS ...@@ -83,18 +64,6 @@ namespace AMDiS
return *matrix_; 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 /// \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 /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
...@@ -102,10 +71,10 @@ namespace AMDiS ...@@ -102,10 +71,10 @@ namespace AMDiS
template <class Index> template <class Index>
auto operator()(Index row, Index col) 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(inserter_, "Inserter not initilized!");
test_exit_dbg(r < rows() && c < cols(), test_exit_dbg(r < this->rows() && c < this->cols(),
"Indices out of range [0,", rows(), ")x[0,", cols(), ")" ); "Indices out of range [0,", this->rows(), ")x[0,", this->cols(), ")" );
return (*inserter_)[r][c]; return (*inserter_)[r][c];
} }
...@@ -116,8 +85,8 @@ namespace AMDiS ...@@ -116,8 +85,8 @@ namespace AMDiS
test_exit(!inserter_, "Matrix already in insertion mode!"); test_exit(!inserter_, "Matrix already in insertion mode!");
calculateSlotSize(); calculateSlotSize();
matrix_->change_dim(rowBasis_->dimension(), colBasis_->dimension()); matrix_->change_dim(this->rows(), this->cols());
test_exit(num_rows(*matrix_) == rowBasis_->dimension() && num_cols(*matrix_) == colBasis_->dimension(), test_exit(num_rows(*matrix_) == this->rows() && num_cols(*matrix_) == this->cols(),
"Wrong dimensions in matrix"); "Wrong dimensions in matrix");
if (prepareForInsertion) { if (prepareForInsertion) {
set_to_zero(*matrix_); set_to_zero(*matrix_);
...@@ -143,30 +112,13 @@ namespace AMDiS ...@@ -143,30 +112,13 @@ namespace AMDiS
/// a one on the diagonal of the matrix. /// a one on the diagonal of the matrix.
auto applyDirichletBC(std::vector<char> const& dirichletNodes) auto applyDirichletBC(std::vector<char> const& dirichletNodes)
{ {
std::list<Triplet<value_type>> columns;
// Define the property maps // Define the property maps
auto row = mtl::mat::row_map(*matrix_); auto row = mtl::mat::row_map(*matrix_);
auto col = mtl::mat::col_map(*matrix_); auto col = mtl::mat::col_map(*matrix_);
auto value = mtl::mat::value_map(*matrix_); auto value = mtl::mat::value_map(*matrix_);
// iterate over the matrix // iterate over the matrix
#if 0 for (auto r : mtl::rows_of(*matrix_)) { // rows of the matrix
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
if (dirichletNodes[r.value()]) { if (dirichletNodes[r.value()]) {
for (auto i : mtl::nz_of(r)) { // non-zeros within for (auto i : mtl::nz_of(r)) { // non-zeros within
// set identity row // set identity row
...@@ -175,67 +127,8 @@ namespace AMDiS ...@@ -175,67 +127,8 @@ namespace AMDiS
} }
} }
return columns; return std::list<Triplet<value_type>>{};
}
#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);
}
}
} }
#endif
protected: protected:
// Estimates the slot-size used in the inserter to reserve enough space per row. // Estimates the slot-size used in the inserter to reserve enough space per row.
...@@ -253,12 +146,6 @@ namespace AMDiS ...@@ -253,12 +146,6 @@ namespace AMDiS
/// The minimal number of nonzeros per row /// The minimal number of nonzeros per row
static constexpr int MIN_NNZ_PER_ROW = 50; 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 /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
ClonablePtr<BaseMatrix> matrix_; ClonablePtr<BaseMatrix> matrix_;
......
...@@ -123,7 +123,7 @@ namespace AMDiS ...@@ -123,7 +123,7 @@ namespace AMDiS
}; };
/// Constructor a dofvector from given basis and name /// Constructor a dofvector from given basis and name
template <class Basis, class ValueType = double> template <class ValueType = double, class Basis>
DOFVector<Basis, ValueType>