Commit 7bd6ff1c authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/cleanup_dofvector' into 'develop'

Issue/cleanup dofvector

See merge request !29
parents 9d6a1e4c 55a910fd
Pipeline #1231 failed with stage
in 4 seconds
#pragma once
#include <algorithm>
#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
#include <amdis/linear_algebra/mtl/SystemVector.hpp>
#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
#include <amdis/linear_algebra/mtl/Mapper.hpp>
namespace AMDiS
{
// requires MatrixType to have an inserter
template <class Value, std::size_t _N, std::size_t _M, class TargetType>
void copy(BlockMTLMatrix<Value, _N, _M> const& source, TargetType& target)
{
using Mapper = std::conditional_t<(_N == _M), BlockMapper, RectangularMapper>;
using value_type = typename TargetType::value_type;
using BaseInserter = mtl::mat::inserter<TargetType, mtl::operations::update_plus<value_type>>;
using MappedInserter = typename Mapper::template Inserter<BaseInserter>::type;
std::size_t nnz = 0;
for (std::size_t rb = 0; rb < _N; ++rb)
for (std::size_t cb = 0; cb < _M; ++cb)
nnz += source[rb][cb].nnz();
{
target.change_dim(num_rows(source), num_cols(source));
set_to_zero(target);
Mapper mapper(source);
MappedInserter ins(target, mapper, int(1.2 * nnz / target.dim1()));
for (std::size_t rb = 0; rb < _N; ++rb) {
mapper.setRow(rb);
for (std::size_t cb = 0; cb < _M; ++cb) {
mapper.setCol(cb);
ins << source[rb][cb];
}
}
}
}
template <class FeSpaces, class Value, class TargetType>
void copy(SystemMatrix<FeSpaces, Value> const& source, TargetType& target)
{
copy(source.getMatrix(), target);
}
} // end namespace AMDiS
\ No newline at end of file
#pragma once
#include <list>
#include <map>
#include <memory>
#include <string>
#include <vector>
......@@ -11,148 +10,74 @@
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>
#include <dune/common/std/optional.hh>
#include <dune/common/reservedvector.hh>
#include <dune/functions/functionspacebases/flatmultiindex.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>
namespace AMDiS
{
/// \brief The basic container that stores a base matrix and a corresponding
/// row/column feSpace.
template <class RowFeSpaceType,
class ColFeSpaceType,
class ValueType = double>
class DOFMatrix
template <class ValueType>
class MtlMatrix
{
public:
/// The type of the finite element space / basis of the row
using RowFeSpace = RowFeSpaceType;
/// The type of the finite element space / basis of the column
using ColFeSpace = ColFeSpaceType;
/// 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;
/// The underlying field-type (typically the same as \ref value_type)
using field_type = typename BaseMatrix::value_type;
/// 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(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
std::string name)
: rowFeSpace(rowFeSpace)
, colFeSpace(colFeSpace)
, name(name)
, matrix(ClonablePtr<BaseMatrix>::make())
{}
/// Constructor. Takes pointer of data-matrix.
DOFMatrix(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
std::string name,
BaseMatrix& matrix_ref)
: rowFeSpace(rowFeSpace)
, colFeSpace(colFeSpace)
, name(name)
, matrix(matrix_ref)
{}
/// Return the row-basis \ref rowFeSpace of the matrix
RowFeSpace const& getRowFeSpace() const
{
return rowFeSpace;
}
/// Return the col-basis \ref colFeSpace of the matrix
ColFeSpace const& getColFeSpace() const
{
return colFeSpace;
}
MtlMatrix() = default;
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& getMatrix()
BaseMatrix& matrix()
{
assert( !insertionMode );
return *matrix;
assert( !inserter_ );
return matrix_;
}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix const& getMatrix() const
BaseMatrix const& matrix() const
{
assert( !insertionMode );
return *matrix;
assert( !inserter_ );
return matrix_;
}
/// Return the size of the row \ref feSpace
size_type N() const
{
return rowFeSpace.size();
}
/// Return the size of the column \ref feSpace
size_type M() const
{
return colFeSpace.size();
}
/// Return the \ref name of this matrix
std::string getName() const
{
return name;
}
using FlatIndex = Dune::Functions::FlatMultiIndex<size_type>;
using BlockIndex = Dune::ReservedVector<size_type,1>;
/// \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.
auto operator()(FlatIndex row, FlatIndex col)
void insert(size_type r, size_type c, value_type const& value)
{
size_type r = row[0], c = col[0];
test_exit_dbg( insertionMode, "Inserter not initilized!");
test_exit_dbg( r < N() && c < M() ,
"Indices out of range [0,", N(), ")x[0,", M(), ")" );
return (*inserter)[r][c];
test_exit_dbg(inserter_, "Inserter not initilized!");
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;
}
auto operator()(BlockIndex row, BlockIndex col)
{
size_type r = row[0], c = col[0];
test_exit_dbg( insertionMode, "Inserter not initilized!");
test_exit_dbg( r < N() && c < M() ,
"Indices out of range [0,", N(), ")x[0,", M(), ")" );
return (*inserter)[r][c];
}
/// 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(!insertionMode && !inserter, "Matrix already in insertion mode!");
test_exit(!inserter_, "Matrix already in insertion mode!");
calculateSlotSize();
matrix->change_dim(rowFeSpace.dimension(), colFeSpace.dimension());
test_exit(num_rows(*matrix) == rowFeSpace.dimension() && num_cols(*matrix) == colFeSpace.dimension(),
"Wrong dimensions in matrix");
matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
if (prepareForInsertion) {
set_to_zero(*matrix);
inserter = new Inserter(*matrix, slot_size);
insertionMode = true;
set_to_zero(matrix_);
inserter_ = new Inserter(matrix_, slotSize_);
}
}
......@@ -160,45 +85,27 @@ namespace AMDiS
/// final construction of the matrix.
void finish()
{
delete inserter;
inserter = NULL;
insertionMode = false;
delete inserter_;
inserter_ = nullptr;
}
/// Return the number of nonzeros in the matrix
std::size_t getNnz() const
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
/// 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);
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
......@@ -207,107 +114,36 @@ 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.
void calculateSlotSize()
{
slot_size = 0;
slotSize_ = 0;
if (num_rows(*matrix) != 0)
slot_size = int(double(matrix->nnz()) / num_rows(*matrix) * 1.2);
if (slot_size < MIN_NNZ_PER_ROW)
slot_size = MIN_NNZ_PER_ROW;
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;
}
private:
/// The minimal number of nonzeros per row
static constexpr int MIN_NNZ_PER_ROW = 50;
/// The finite element space / basis of the row
RowFeSpace const& rowFeSpace;
/// The finite element space / basis of the column
ColFeSpace const& colFeSpace;
/// The name of the DOFMatrix
std::string name;
/// 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 = NULL;
/// A flag that indicates whether the matrix is in insertion mode
bool insertionMode = false;
Inserter* inserter_ = nullptr;
/// The size of the slots used to insert values per row
int slot_size = MIN_NNZ_PER_ROW;
// friend class declarations
template <class, class> friend class SystemMatrix;
int slotSize_ = MIN_NNZ_PER_ROW;
};
template <class RowBasisType, class ColBasisType, class ValueType = double>
using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, MtlMatrix<ValueType>>;
} // end namespace AMDiS
#pragma once
#include <algorithm>
#include <memory>
#include <string>
#include <dune/common/reservedvector.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#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/mtl/MTLDenseVector.hpp>
#include <amdis/utility/MultiIndex.hpp>
#include <amdis/linear_algebra/DOFVectorBase.hpp>
namespace AMDiS
{
/// The basic container that stores a base vector and a corresponding feSpace
template <class Traits, class ValueType = double>
class DOFVector
/// The basic container that stores a base vector and a corresponding basis
template <class ValueType>
class MtlVector
{
using Self = DOFVector;
public:
/// The type of the \ref feSpace
using GlobalBasis = typename Traits::GlobalBasis;
/// The type of the base vector
using BaseVector = MTLDenseVector<ValueType>;
using BaseVector = mtl::vec::dense_vector<ValueType>;
/// The index/size - type
using size_type = typename GlobalBasis::size_type;
using size_type = typename BaseVector::size_type;
/// The type of the elements of the DOFVector
using value_type = ValueType;
/// The underlying field-type (typically the same as \ref value_type)
using field_type = typename BaseVector::value_type;
public:
/// Constructor. Constructs new BaseVector.
DOFVector(GlobalBasis const& feSpace, std::string name)
: feSpace(feSpace)
, name(name)
, vector(ClonablePtr<BaseVector>::make())
{
compress();
}
/// Constructor. Takes reference to existing BaseVector
DOFVector(GlobalBasis const& feSpace, std::string name,
BaseVector& vector_ref)
: feSpace(feSpace)
, name(name)
, vector(vector_ref)
{}
/// Return the basis \ref feSpace of the vector
GlobalBasis const& getFeSpace() const
{
return feSpace;
}
/// Return the data-vector \ref vector
BaseVector const& getVector() const
{
return *vector;
}
/// Return the data-vector \ref vector
BaseVector& getVector()
{
return *vector;
}
MtlVector() = default;
/// Return the size of the \ref feSpace
size_type getSize() const
/// Return the data-vector \ref vector_
BaseVector const& vector() const
{
return feSpace.dimension();
return vector_;
}
/// Return the \ref name of this vector
std::string getName() const
/// Return the data-vector \ref vector_
BaseVector& vector()
{
return name;
return vector_;
}
/// Resize the \ref vector to the size of the \ref feSpace.
void compress()
/// Return the current size of the \ref vector_
size_type size() const
{
if (num_rows(*vector) != getSize()) {
resize(getSize());
*vector = ValueType(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
{
auto const i = flatMultiIndex(idx);
test_exit_dbg( i < num_rows(*vector) ,
"Index ", i, " out of range [0,", num_rows(*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)
{
auto const i = flatMultiIndex(idx);
test_exit_dbg( i < num_rows(*vector) ,
"Index ", i, " out of range [0,", num_rows(*vector), ")" );
return (*vector)[i];
}
/// \brief interpolate a function \p f to the basis \ref feSpace and store the
/// coefficients in \ref vector.
template <class F>
void interpol(F&& f)
{
Dune::Functions::interpolate(feSpace, *this, std::forward<F>(f));
}
/// 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.getVector();
test_exit_dbg(i < mtl::vec::size(vector_),
"Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
return vector_[i];
}
private:
/// The finite element space / basis associated with the data vector
GlobalBasis const& feSpace;
/// The name of the DOFVector (used in file-writing)