From bd9cc7c0eb25d2f3c527cb15c882cbc53d186a2d Mon Sep 17 00:00:00 2001
From: "Praetorius, Simon" <simon.praetorius@tu-dresden.de>
Date: Mon, 30 Dec 2019 23:32:03 +0100
Subject: [PATCH] Matrix-Vector facade class
---
doc/doxygen/Doxylocal | 5 +-
src/amdis/BiLinearForm.hpp | 33 +-
src/amdis/BiLinearForm.inc.hpp | 16 +-
src/amdis/DOFVector.hpp | 53 +-
src/amdis/DOFVector.inc.hpp | 8 +-
src/amdis/DataTransfer.inc.hpp | 4 +-
src/amdis/LinearAlgebra.hpp | 4 +-
src/amdis/LinearForm.hpp | 34 +-
src/amdis/LinearForm.inc.hpp | 16 +-
src/amdis/Observer.hpp | 75 +--
src/amdis/ProblemStat.inc.hpp | 4 +-
src/amdis/functions/CMakeLists.txt | 1 +
src/amdis/functions/Interpolate.hpp | 7 +-
src/amdis/functions/NodeIndices.hpp | 8 +-
src/amdis/functions/ParallelGlobalBasis.hpp | 6 +-
src/amdis/functions/SizeInfo.hpp | 89 ++++
src/amdis/gridfunctions/DiscreteFunction.hpp | 255 +++++-----
.../gridfunctions/DiscreteFunction.inc.hpp | 26 +-
.../DiscreteLocalFunction.inc.hpp | 32 +-
src/amdis/io/FileWriterCreator.hpp | 28 +-
src/amdis/linearalgebra/CMakeLists.txt | 6 +-
src/amdis/linearalgebra/Constraints.hpp | 12 +-
src/amdis/linearalgebra/DOFMapping.hpp | 2 +
src/amdis/linearalgebra/MatrixBase.hpp | 133 -----
src/amdis/linearalgebra/MatrixFacade.hpp | 107 +++++
src/amdis/linearalgebra/SparsityPattern.hpp | 155 ++++++
src/amdis/linearalgebra/VectorBase.hpp | 454 +++---------------
src/amdis/linearalgebra/VectorFacade.hpp | 370 ++++++++++++++
src/amdis/linearalgebra/eigen/CMakeLists.txt | 1 +
src/amdis/linearalgebra/eigen/Constraints.hpp | 12 +-
.../linearalgebra/eigen/MatrixBackend.hpp | 18 +-
src/amdis/linearalgebra/eigen/MatrixSize.hpp | 60 +++
src/amdis/linearalgebra/eigen/Traits.hpp | 19 +-
.../linearalgebra/eigen/VectorBackend.hpp | 85 +---
src/amdis/linearalgebra/istl/Constraints.hpp | 10 +-
.../linearalgebra/istl/MatrixBackend.hpp | 47 +-
src/amdis/linearalgebra/istl/Traits.hpp | 44 +-
.../linearalgebra/istl/VectorBackend.hpp | 98 ++--
src/amdis/linearalgebra/mtl/CMakeLists.txt | 1 +
src/amdis/linearalgebra/mtl/Constraints.hpp | 12 +-
src/amdis/linearalgebra/mtl/MatrixBackend.hpp | 45 +-
src/amdis/linearalgebra/mtl/SlotSize.hpp | 90 ++++
src/amdis/linearalgebra/mtl/Traits.hpp | 20 +-
src/amdis/linearalgebra/mtl/VectorBackend.hpp | 88 +---
.../linearalgebra/petsc/Communication.hpp | 10 +-
src/amdis/linearalgebra/petsc/Constraints.hpp | 12 +-
.../linearalgebra/petsc/MatrixBackend.hpp | 64 +--
.../petsc/MatrixNnzStructure.hpp | 43 +-
.../petsc/MatrixNnzStructure.inc.hpp | 16 +-
src/amdis/linearalgebra/petsc/Traits.hpp | 13 +-
.../linearalgebra/petsc/VectorBackend.hpp | 87 ++--
src/amdis/typetree/MultiIndex.hpp | 9 +
test/DOFVectorTest.cpp | 57 ++-
test/DiscreteFunctionTest.cpp | 4 +-
test/ExpressionsTest.cpp | 2 +-
test/IntegrateTest.cpp | 2 +-
test/ObserverTest.cpp | 12 +-
test/PETScCommTest.cpp | 2 +-
test/ProblemStatTest.cpp | 13 +-
59 files changed, 1651 insertions(+), 1288 deletions(-)
create mode 100644 src/amdis/functions/SizeInfo.hpp
delete mode 100644 src/amdis/linearalgebra/MatrixBase.hpp
create mode 100644 src/amdis/linearalgebra/MatrixFacade.hpp
create mode 100644 src/amdis/linearalgebra/SparsityPattern.hpp
create mode 100644 src/amdis/linearalgebra/VectorFacade.hpp
create mode 100644 src/amdis/linearalgebra/eigen/MatrixSize.hpp
create mode 100644 src/amdis/linearalgebra/mtl/SlotSize.hpp
diff --git a/doc/doxygen/Doxylocal b/doc/doxygen/Doxylocal
index 5c018991..220d9233 100644
--- a/doc/doxygen/Doxylocal
+++ b/doc/doxygen/Doxylocal
@@ -10,6 +10,7 @@ HIDE_SCOPE_NAMES = YES
HIDE_UNDOC_CLASSES = NO
INTERNAL_DOCS = NO
MARKDOWN_SUPPORT = YES
+GENERATE_XML = YES
EXCLUDE_SYMBOLS = AMDiS::Impl \
AMDiS::Math::Impl_ \
@@ -30,6 +31,7 @@ PREDEFINED += HAVE_UMFPACK \
INPUT += @top_srcdir@/src/amdis \
@top_srcdir@/src/amdis/common \
+ @top_srcdir@/src/amdis/functions \
@top_srcdir@/src/amdis/gridfunctions \
@top_srcdir@/src/amdis/linearalgebra \
@top_srcdir@/src/amdis/linearalgebra/mtl \
@@ -47,7 +49,8 @@ INPUT += @top_srcdir@/src/amdis \
# subdirectory from a directory tree whose root is specified with the INPUT tag.
EXCLUDE += @top_srcdir@/src/amdis/linearalgebra/eigen \
- @top_srcdir@/src/amdis/linearalgebra/istl
+ @top_srcdir@/src/amdis/linearalgebra/istl \
+ @top_srcdir@/src/amdis/linearalgebra/petsc
# The EXAMPLE_PATH tag can be used to specify one or more files or
# directories that contain example code fragments that are included (see
diff --git a/src/amdis/BiLinearForm.hpp b/src/amdis/BiLinearForm.hpp
index d3864d6d..f7d3b8e2 100644
--- a/src/amdis/BiLinearForm.hpp
+++ b/src/amdis/BiLinearForm.hpp
@@ -20,14 +20,12 @@ namespace AMDiS
* \tparam CB Basis of matrix columns
* \tparam T Coefficient type to store in the matrix
**/
- template <class RB, class CB, class T = double>
+ template <class RB, class CB, class T = double, class Traits = BackendTraits<RB,T>>
class BiLinearForm
- : public MatrixBase<RB,CB,MatrixBackend<BackendTraits<RB,T>>>
+ : public MatrixFacade<T, typename Traits::SparsityPattern, Traits::template MatrixImpl>
{
using Self = BiLinearForm;
- using Traits = BackendTraits<RB,T>;
- using Backend = MatrixBackend<Traits>;
- using Super = MatrixBase<RB,CB,Backend>;
+ using Super = MatrixFacade<T, typename Traits::SparsityPattern, Traits::template MatrixImpl>;
public:
/// The type of the finite element space / basis of the row
@@ -47,23 +45,21 @@ namespace AMDiS
public:
/// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into
/// a new shared_ptr.
- template <class RB_, class CB_>
- BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
- : Super(FWD(rowBasis), FWD(colBasis))
+ BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
+ : Super(*rowBasis, *colBasis)
+ , rowBasis_(rowBasis)
+ , colBasis_(colBasis)
{
- operators_.init(*Super::rowBasis(), *Super::colBasis());
- }
-
- /// Prepare the matrix for insertion. Clears all the entries.
- void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
- {
- Super::init(symmetry);
+ operators_.init(*rowBasis, *colBasis);
- auto const rowSize = this->rowBasis()->localView().maxSize();
- auto const colSize = this->colBasis()->localView().maxSize();
+ auto const rowSize = rowBasis->localView().maxSize();
+ auto const colSize = colBasis->localView().maxSize();
elementMatrix_.resize(rowSize, colSize);
}
+ std::shared_ptr<RowBasis const> const& rowBasis() const { return rowBasis_; }
+ std::shared_ptr<ColBasis const> const& colBasis() const { return colBasis_; }
+
/// \brief Associate a local operator with this BiLinearForm
/**
* Stores an operator in a list that gets assembled during a call to \ref assemble().
@@ -118,6 +114,9 @@ namespace AMDiS
/// List of operators associated to row/col node
MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;
+
+ std::shared_ptr<RowBasis const> rowBasis_;
+ std::shared_ptr<ColBasis const> colBasis_;
};
diff --git a/src/amdis/BiLinearForm.inc.hpp b/src/amdis/BiLinearForm.inc.hpp
index 3e66984c..83cd7fb3 100644
--- a/src/amdis/BiLinearForm.inc.hpp
+++ b/src/amdis/BiLinearForm.inc.hpp
@@ -9,9 +9,9 @@
namespace AMDiS {
-template <class RB, class CB, class T>
+template <class RB, class CB, class T, class Traits>
template <class ContextTag, class Expr, class RowTreePath, class ColTreePath>
-void BiLinearForm<RB,CB,T>::
+void BiLinearForm<RB,CB,T,Traits>::
addOperator(ContextTag contextTag, Expr const& expr,
RowTreePath row, ColTreePath col)
{
@@ -24,16 +24,16 @@ addOperator(ContextTag contextTag, Expr const& expr,
auto j = child(this->colBasis()->localView().tree(), makeTreePath(col));
using LocalContext = typename ContextTag::type;
- using Traits = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
+ using Tr = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis()->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i, j));
+ auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), i, j));
operators_[i][j].push(contextTag, std::move(localAssembler));
}
-template <class RB, class CB, class T>
-void BiLinearForm<RB,CB,T>::
+template <class RB, class CB, class T, class Traits>
+void BiLinearForm<RB,CB,T,Traits>::
assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
{
elementMatrix_.resize(rowLocalView.size(), colLocalView.size());
@@ -58,8 +58,8 @@ assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
}
-template <class RB, class CB, class T>
-void BiLinearForm<RB,CB,T>::
+template <class RB, class CB, class T, class Traits>
+void BiLinearForm<RB,CB,T,Traits>::
assemble(SymmetryStructure symmetry)
{
auto rowLocalView = this->rowBasis()->localView();
diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp
index f70b4d91..1dc708b9 100644
--- a/src/amdis/DOFVector.hpp
+++ b/src/amdis/DOFVector.hpp
@@ -33,38 +33,40 @@ namespace AMDiS
/**
* \tparam GB Basis of the vector
* \tparam T Type of the coefficients
+ * \tparam TraitsType Collection of parameter for the linear-algebra backend
**/
- template <class GB, class T = double>
+ template <class GB, class T = double, class TraitsType = BackendTraits<GB,T>>
class DOFVector
- : public VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>
+ : public VectorFacade<T, TraitsType::template VectorImpl>
, private Observer<event::preAdapt>
, private Observer<event::adapt>
, private Observer<event::postAdapt>
{
using Self = DOFVector;
- using Super = VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>;
+ using Coefficients = VectorFacade<T, TraitsType::template VectorImpl>;
public:
- using Backend = VectorBackend<BackendTraits<GB,T>>;
-
using GlobalBasis = GB;
/// The index/size - type
using size_type = typename GlobalBasis::size_type;
/// The type of the elements of the DOFVector
- using value_type = typename Backend::value_type;
+ using value_type = T;
/// Wrapper for the implementation of the transfer of data during grid adaption
using DataTransfer = DataTransferWrapper<Self>;
+ /// Collection of parameter for the linear-algebra backend
+ using Traits = TraitsType;
+
public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
DOFVector(std::shared_ptr<GB> const& basis,
DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- : Super(basis)
+ : Coefficients(*basis)
, Observer<event::preAdapt>(basis->gridView().grid())
- , Observer<event::adapt>(*basis)
+ , Observer<event::adapt>(basis)
, Observer<event::postAdapt>(basis->gridView().grid())
, dataTransfer_(op, basis)
, basis_(basis)
@@ -83,19 +85,20 @@ namespace AMDiS
: DOFVector(Dune::wrap_or_move(FWD(basis)), op)
{}
+ std::shared_ptr<GlobalBasis const> const& basis() const { return basis_; }
template <class TreePath = RootTreePath>
auto child(TreePath const& path = {})
{
auto&& tp = makeTreePath(path);
- return makeDiscreteFunction(*this, tp);
+ return makeDiscreteFunction(coefficients(), *basis_, tp);
}
template <class TreePath = RootTreePath>
auto child(TreePath const& path = {}) const
{
auto&& tp = makeTreePath(path);
- return makeDiscreteFunction(*this, tp);
+ return makeDiscreteFunction(coefficients(), *basis_, tp);
}
@@ -103,7 +106,7 @@ namespace AMDiS
/// reference to this DOFVector in the expression.
/// See \ref DiscreteFunction::interpolate_noalias
template <class Expr, class Tag = tag::average>
- void interpolate_noalias(Expr&& expr, Tag strategy)
+ void interpolate_noalias(Expr&& expr, Tag strategy = {})
{
child().interpolate_noalias(FWD(expr), strategy);
}
@@ -111,7 +114,7 @@ namespace AMDiS
/// Interpolation of GridFunction to DOFVector.
/// See \ref DiscreteFunction::interpolate
template <class Expr, class Tag = tag::average>
- void interpolate(Expr&& expr, Tag strategy)
+ void interpolate(Expr&& expr, Tag strategy = {})
{
child().interpolate(FWD(expr), strategy);
}
@@ -126,6 +129,16 @@ namespace AMDiS
}
+ Coefficients const& coefficients() const
+ {
+ return static_cast<Coefficients const&>(*this);
+ }
+
+ Coefficients& coefficients()
+ {
+ return static_cast<Coefficients&>(*this);
+ }
+
/// Write DOFVector to file
void backup(std::string const& filename);
@@ -192,7 +205,7 @@ namespace AMDiS
/// to \ref DataTransferOperation::INTERPOLATE.
DataTransfer dataTransfer_;
- std::shared_ptr<GB> basis_;
+ std::shared_ptr<GlobalBasis const> basis_;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
@@ -218,6 +231,20 @@ namespace AMDiS
return {FWD(basis), op};
}
+ /// A Generator for a mutable \ref DiscreteFunction
+ template <class GB, class T, class Path = RootTreePath>
+ auto makeDiscreteFunction(DOFVector<GB,T>& dofVec, Path const& path = {})
+ {
+ return dofVec.child(path);
+ }
+
+ /// A Generator for a mutable \ref DiscreteFunction
+ template <class GB, class T, class Path = RootTreePath>
+ auto makeDiscreteFunction(DOFVector<GB,T> const& dofVec, Path const& path = {})
+ {
+ return dofVec.child(path);
+ }
+
} // end namespace AMDiS
#include <amdis/DOFVector.inc.hpp>
diff --git a/src/amdis/DOFVector.inc.hpp b/src/amdis/DOFVector.inc.hpp
index fd940199..eb02fe3f 100644
--- a/src/amdis/DOFVector.inc.hpp
+++ b/src/amdis/DOFVector.inc.hpp
@@ -10,8 +10,8 @@
namespace AMDiS {
-template <class GB, class B>
-void DOFVector<GB,B>::
+template <class GB, class T, class Traits>
+void DOFVector<GB,T,Traits>::
backup(std::string const& filename)
{
std::ofstream out(filename, std::ios::binary);
@@ -35,8 +35,8 @@ backup(std::string const& filename)
}
-template <class GB, class B>
-void DOFVector<GB,B>::
+template <class GB, class T, class Traits>
+void DOFVector<GB,T,Traits>::
restore(std::string const& filename)
{
std::ifstream in(filename, std::ios::binary);
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index bf84ccfc..f9434fb8 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -82,7 +82,7 @@ preAdapt(C const& coeff, bool mightCoarsen)
auto maxLevel = gv.grid().maxLevel();
using std::sqrt;
typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
- for (auto const& e : elements(gv, typename C::Backend::Traits::PartitionSet{}))
+ for (auto const& e : elements(gv, typename C::Traits::PartitionSet{}))
{
auto father = e;
while (father.mightVanish() && father.hasFather())
@@ -162,7 +162,7 @@ void DataTransfer<C,B>::adapt(C& coeff)
mapper_.update();
std::vector<bool> finished(mapper_.size(), false);
- for (const auto& e : elements(gv, typename C::Backend::Traits::PartitionSet{}))
+ for (const auto& e : elements(gv, typename C::Traits::PartitionSet{}))
{
auto index = mapper_.index(e);
if (finished[index])
diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp
index f8c5d4de..d1dda9f1 100644
--- a/src/amdis/LinearAlgebra.hpp
+++ b/src/amdis/LinearAlgebra.hpp
@@ -37,7 +37,7 @@
#endif
-#include <amdis/linearalgebra/MatrixBase.hpp>
-#include <amdis/linearalgebra/VectorBase.hpp>
+#include <amdis/linearalgebra/MatrixFacade.hpp>
+#include <amdis/linearalgebra/VectorFacade.hpp>
#include <amdis/linearalgebra/LinearSolver.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
diff --git a/src/amdis/LinearForm.hpp b/src/amdis/LinearForm.hpp
index 8601d904..78d04718 100644
--- a/src/amdis/LinearForm.hpp
+++ b/src/amdis/LinearForm.hpp
@@ -17,24 +17,20 @@ namespace AMDiS
*
* \tparam GB Basis of the vector
* \tparam T Coefficient type to store in the vector
+ * \tparam Traits Collection of parameter for the linear-algebra backend
**/
- template <class GB, class T = double>
+ template <class GB, class T = double, class Traits = BackendTraits<GB,T>>
class LinearForm
- : public VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>
+ : public VectorFacade<T, Traits::template VectorImpl>
{
- using Super = VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>;
+ using Super = VectorFacade<T, Traits::template VectorImpl>;
using Self = LinearForm;
public:
- using Backend = VectorBackend<BackendTraits<GB,T>>;
-
/// The type of the functionspace basis associated to this linearform
using GlobalBasis = GB;
using LocalView = typename GlobalBasis::LocalView;
- /// A traits class collecting several parameters of the linear algebra backend
- using Traits = typename Backend::Traits;
-
/// The type of the elements of the DOFVector
using CoefficientType = typename Traits::CoefficientType;
@@ -43,24 +39,18 @@ namespace AMDiS
public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
- template <class GB_,
- Dune::disableCopyMove<Self, GB_> = 0>
- explicit LinearForm(GB_&& basis)
- : Super(FWD(basis))
+ explicit LinearForm(std::shared_ptr<GB> const& basis)
+ : Super(*basis)
+ , basis_(basis)
{
- operators_.init(*Super::basis());
- }
+ operators_.init(*basis);
- /// Prepare the LinearForm for insertion of values, finish the insertion with
- /// \ref finish().
- void init(bool clear)
- {
- Super::init(clear);
-
- auto const localSize = this->basis()->localView().maxSize();
+ auto const localSize = basis->localView().maxSize();
elementVector_.resize(localSize);
}
+ std::shared_ptr<GlobalBasis const> const& basis() const { return basis_; }
+
/// \brief Associate a local operator with this LinearForm
/**
* Stores an operator in a list that gets assembled during a call to \ref assemble().
@@ -112,6 +102,8 @@ namespace AMDiS
/// List of operators associated to nodes, filled in \ref addOperator().
VectorOperators<GlobalBasis,ElementVector> operators_;
+
+ std::shared_ptr<GlobalBasis const> basis_;
};
diff --git a/src/amdis/LinearForm.inc.hpp b/src/amdis/LinearForm.inc.hpp
index cd39cb80..5277f410 100644
--- a/src/amdis/LinearForm.inc.hpp
+++ b/src/amdis/LinearForm.inc.hpp
@@ -10,9 +10,9 @@
namespace AMDiS {
-template <class GB, class B>
+template <class GB, class T, class Traits>
template <class ContextTag, class Expr, class TreePath>
-void LinearForm<GB,B>::
+void LinearForm<GB,T,Traits>::
addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
{
static_assert( Concepts::PreTreePath<TreePath>,
@@ -21,16 +21,16 @@ addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
auto i = child(this->basis()->localView().tree(), makeTreePath(path));
using LocalContext = typename ContextTag::type;
- using Traits = DefaultAssemblerTraits<LocalContext, ElementVector>;
+ using Tr = DefaultAssemblerTraits<LocalContext, ElementVector>;
auto op = makeLocalOperator<LocalContext>(expr, this->basis()->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i));
+ auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), i));
operators_[i].push(contextTag, std::move(localAssembler));
}
-template <class GB, class B>
-void LinearForm<GB,B>::
+template <class GB, class T, class Traits>
+void LinearForm<GB,T,Traits>::
assemble(typename GB::LocalView const& localView)
{
elementVector_.resize(localView.size());
@@ -53,8 +53,8 @@ assemble(typename GB::LocalView const& localView)
}
-template <class GB, class B>
-void LinearForm<GB,B>::
+template <class GB, class T, class Traits>
+void LinearForm<GB,T,Traits>::
assemble()
{
auto localView = this->basis()->localView();
diff --git a/src/amdis/Observer.hpp b/src/amdis/Observer.hpp
index b031e431..865e5936 100644
--- a/src/amdis/Observer.hpp
+++ b/src/amdis/Observer.hpp
@@ -1,14 +1,15 @@
#pragma once
#include <algorithm>
-#include <list>
#include <memory>
+#include <set>
#include <type_traits>
+#include <dune/common/shared_ptr.hh>
#include <dune/common/typeutilities.hh>
#include <amdis/common/ConceptsBase.hpp>
-#include <amdis/Output.hpp>
+#include <amdis/common/Index.hpp>
namespace AMDiS
{
@@ -75,19 +76,17 @@ namespace AMDiS
/// Attach a new observer that gets called on \ref notify
void attach(ObserverInterface<Event>* o)
{
- observers_.push_back(o);
+ observers_.insert(o);
}
/// Detaches the passed observer from the list, if stored.
void detach(ObserverInterface<Event>* o)
{
- auto it = std::find(observers_.begin(), observers_.end(), o);
- if (it != observers_.end())
- observers_.erase(it);
+ observers_.erase(o);
}
private:
- std::list<ObserverInterface<Event>*> observers_;
+ std::set<ObserverInterface<Event>*> observers_;
};
@@ -98,17 +97,23 @@ namespace AMDiS
{
public:
template <class Notifier>
- Observer(Notifier const& notifier)
- : notifier_(const_cast<Notifier*>(¬ifier))
+ Observer(std::shared_ptr<Notifier> notifier)
+ : notifier_(std::const_pointer_cast<std::remove_const_t<Notifier>>(std::move(notifier)))
{
notifier_->attach(this);
}
+ template <class Notifier,
+ class = void_t<decltype(std::declval<std::remove_const_t<Notifier>>().notify(std::declval<Event>()))> >
+ Observer(Notifier& notifier)
+ : Observer(Dune::stackobject_to_shared_ptr(notifier))
+ {}
+
/// Destructor, detaches from the notifier
virtual ~Observer()
{
- assert(notifier_);
- notifier_->detach(this);
+ if (notifier_)
+ notifier_->detach(this);
}
/// Copy constructor. Attaches this to the copied notifier
@@ -121,23 +126,11 @@ namespace AMDiS
/// Copy-assignment operator, copies the notifier and attaches this.
Observer& operator=(Observer const& other)
{
- notifier_ = other.notifier_;
- notifier_->attach(this);
- return *this;
- }
-
- /// Move assignment operator, implemented in terms of the move
- /// assignment operator
- Observer(Observer&& other)
- {
- *this = std::move(other);
- }
-
- /// Move-assignment operator, copies the notifier and attaches this.
- Observer& operator=(Observer&& other)
- {
- notifier_ = other.notifier_;
- notifier_->attach(this);
+ if (&other != this) {
+ notifier_->detach(this);
+ notifier_ = other.notifier_;
+ notifier_->attach(this);
+ }
return *this;
}
@@ -155,7 +148,31 @@ namespace AMDiS
virtual void updateImpl(Event e, Tags...) = 0;
private:
- Notifier<Event>* notifier_ = nullptr;
+ std::shared_ptr<Notifier<Event>> notifier_ = nullptr;
};
+
+ namespace Impl
+ {
+ template <class Event, class Tags>
+ class ObserverSequenceImpl;
+
+ /// Combination of multiple observers of the same event but with different tags
+ template <class Event, std::size_t... Is>
+ class ObserverSequenceImpl<Event, std::index_sequence<Is...>>
+ : private Observer<Event,index_t<Is>>...
+ {
+ public:
+ template <class... Notifiers,
+ REQUIRES(sizeof...(Notifiers) == sizeof...(Is))>
+ ObserverSequenceImpl(Notifiers&&... notifiers)
+ : Observer<Event,index_t<Is>>(FWD(notifiers))...
+ {}
+ };
+
+ } // end namespace Impl
+
+ template <class Event, std::size_t N>
+ using ObserverSequence = Impl::ObserverSequenceImpl<Event, std::make_index_sequence<N>>;
+
} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index ab4699a6..b83915b6 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -382,8 +382,8 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solverInfo.setStoreMatrixData(storeMatrixData);
solution_->resize();
- linearSolver_->solve(systemMatrix_->backend().matrix(), solution_->backend().vector(),
- rhs_->backend().vector(), globalBasis_->communicator(), solverInfo);
+ linearSolver_->solve(systemMatrix_->impl().matrix(), solution_->impl().vector(),
+ rhs_->impl().vector(), globalBasis_->communicator(), solverInfo);
if (solverInfo.info() > 0) {
msg("solution of discrete system needed {} seconds", t.elapsed());
diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt
index 73697b63..c8ef3913 100644
--- a/src/amdis/functions/CMakeLists.txt
+++ b/src/amdis/functions/CMakeLists.txt
@@ -6,4 +6,5 @@ install(FILES
NodeIndices.hpp
Nodes.hpp
ParallelGlobalBasis.hpp
+ SizeInfo.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)
diff --git a/src/amdis/functions/Interpolate.hpp b/src/amdis/functions/Interpolate.hpp
index dbec178c..e7f23a8f 100644
--- a/src/amdis/functions/Interpolate.hpp
+++ b/src/amdis/functions/Interpolate.hpp
@@ -16,6 +16,7 @@
#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
+#include <amdis/linearalgebra/Traits.hpp>
#include <amdis/operations/Assigner.hpp>
#include <amdis/typetree/Traversal.hpp>
@@ -39,7 +40,7 @@ namespace AMDiS
// set vector to zero at subtree
if (! std::is_same<Assign, Assigner::assign>::value) {
- for (const auto& e : elements(basis.gridView(), typename Vec::Backend::Traits::PartitionSet{}))
+ for (const auto& e : elements(basis.gridView(), typename BackendTraits<B>::PartitionSet{}))
{
localView.bind(e);
auto&& subTree = Dune::TypeTree::child(localView.tree(), treePath);
@@ -54,7 +55,7 @@ namespace AMDiS
vector.init(false);
counter.init(true); // set to zero
- for (const auto& e : elements(basis.gridView(), typename Vec::Backend::Traits::PartitionSet{}))
+ for (const auto& e : elements(basis.gridView(), typename BackendTraits<B>::PartitionSet{}))
{
localView.bind(e);
lf.bind(e);
@@ -68,7 +69,7 @@ namespace AMDiS
auto&& fe = node.finiteElement();
std::size_t feSize = fe.localBasis().size();
- auto bitVecRange = mappedRangeView(Dune::range(feSize), [&](auto i) -> bool {
+ auto bitVecRange = mappedRangeView(Dune::range(feSize), [&](std::size_t i) -> bool {
return bitVec[localView.index(node.localIndex(i))];
});
diff --git a/src/amdis/functions/NodeIndices.hpp b/src/amdis/functions/NodeIndices.hpp
index 98fd9463..3c17182e 100644
--- a/src/amdis/functions/NodeIndices.hpp
+++ b/src/amdis/functions/NodeIndices.hpp
@@ -16,8 +16,8 @@ namespace AMDiS
static_assert(Dune::models<Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>(), "");
static_assert(Dune::models<Concept::BasisTree<typename LocalView::GridView>, Node>(), "");
- return mappedRangeView(Dune::range(node.size()), [&](std::size_t j) -> std::size_t {
- return flatMultiIndex(localView.index(node.localIndex(j)));
+ return mappedRangeView(Dune::range(node.size()), [&](std::size_t j) -> typename LocalView::MultiIndex {
+ return localView.index(node.localIndex(j));
});
}
@@ -28,8 +28,8 @@ namespace AMDiS
using namespace Dune::Functions;
static_assert(Dune::models<Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>(), "");
- return mappedRangeView(Dune::range(localView.size()), [&](std::size_t i) -> std::size_t {
- return flatMultiIndex(localView.index(i));
+ return mappedRangeView(Dune::range(localView.size()), [&](std::size_t i) -> typename LocalView::MultiIndex {
+ return localView.index(i);
});
}
diff --git a/src/amdis/functions/ParallelGlobalBasis.hpp b/src/amdis/functions/ParallelGlobalBasis.hpp
index 0837be09..8e219428 100644
--- a/src/amdis/functions/ParallelGlobalBasis.hpp
+++ b/src/amdis/functions/ParallelGlobalBasis.hpp
@@ -101,11 +101,15 @@ namespace AMDiS
* This will create a new ParallellGlobalBasis. The pre-basis is copied from the constructor
* argument and a new communication object is built.
*/
- // TODO(FM): Replace explicit type with concept
ParallelGlobalBasis(Dune::Functions::DefaultGlobalBasis<PB>&& from)
: ParallelGlobalBasis(std::string(""), from.gridView().grid(), from.preBasis())
{}
+ /// Copy constructor
+ ParallelGlobalBasis(ParallelGlobalBasis const&) = delete;
+
+ /// Move constructor
+ ParallelGlobalBasis(ParallelGlobalBasis&&) = default;
public:
diff --git a/src/amdis/functions/SizeInfo.hpp b/src/amdis/functions/SizeInfo.hpp
new file mode 100644
index 00000000..5f1d2777
--- /dev/null
+++ b/src/amdis/functions/SizeInfo.hpp
@@ -0,0 +1,89 @@
+#pragma once
+
+#include <algorithm>
+#include <functional>
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <vector>
+
+#include <dune/common/typetraits.hh>
+
+namespace AMDiS
+{
+ /// \brief Type-erased wrapper around size providers
+ class SizeInfo
+ {
+ public:
+ using size_type = std::size_t;
+
+ // The SizePrefix can be created from any container with begin-end iterators.
+ // It just stores a copy of the given prefix in a `std::vector` and copies its
+ // content back on request of the `SizeProvider::size` method, where the actual
+ // `SizeProvider::SizePrefix` is known.
+ class SizePrefix
+ : public std::vector<size_type>
+ {
+ using Super = std::vector<size_type>;
+
+ public:
+ SizePrefix() = default;
+
+ // Constructor copies the content of the passed container to the internal storage
+ template <class Container,
+ Dune::disableCopyMove<SizePrefix, Container> = 0,
+ Dune::disableCopyMove<Super, Container> = 0>
+ SizePrefix(Container const& container)
+ : Super(container.begin(), container.end())
+ {}
+
+ using Super::Super;
+
+ friend std::ostream& operator<<(std::ostream& out, SizePrefix const& prefix)
+ {
+ std::copy(prefix.begin(), prefix.end(), std::ostream_iterator<size_type>(out, " "));
+ return out;
+ }
+ };
+
+ public:
+ /// Store the passed `Dune::SizeInfo` object in a type-erased manner
+ /// by providing an own `SizePrefix` and `size_type` type, that can be
+ /// filled from the actual objects.
+ template <class SizeProvider,
+ class = void_t<decltype(std::declval<SizeProvider>()(std::declval<typename SizeProvider::SizePrefix>()))> >
+ explicit SizeInfo(SizeProvider const& provider)
+ : size_([provider](SizePrefix const& prefix) -> size_type {
+ typename SizeProvider::SizePrefix sizePrefix;
+ sizePrefix.resize(prefix.size());
+ std::copy(prefix.begin(), prefix.end(), sizePrefix.begin());
+ return provider(sizePrefix);
+ })
+ , dimension_([provider]() -> size_type {
+ return size_type(provider);
+ })
+ {}
+
+ /// Return number possible values for next position in multi index
+ size_type size(SizePrefix prefix) const
+ {
+ return size_(prefix);
+ }
+
+ /// Return number possible values for next position in multi index
+ size_type operator()(SizePrefix prefix) const
+ {
+ return size_(prefix);
+ }
+
+ /// Return the total dimension of the associated size-provider (basis)
+ operator size_type() const
+ {
+ return dimension_();
+ }
+
+ private:
+ std::function<size_type(SizePrefix const&)> size_;
+ std::function<size_type()> dimension_;
+ };
+}
\ No newline at end of file
diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
index ea97c780..99298277 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -16,71 +16,137 @@
namespace AMDiS
{
- template <class GB, class VT>
- class DOFVector;
-
/// \class DiscreteFunction
/// \brief A view on a subspace of a \ref DOFVector
/**
* \ingroup GridFunctions
*
- * \tparam GB Type of the global basis
- * \tparam VT Coefficient type of the DOFVector
- * \tparam TP A realization of \ref Dune::TypeTree::HybridTreePath
- * \tparam is_const Specifies whether a const or mutable view is implemented.
+ * \tparam Coeff Const or mutable coefficient type of the DOFVector
+ * \tparam GB Thy type of the global basis
+ * \tparam TreePath A realization of \ref Dune::TypeTree::HybridTreePath
*
* **Requirements:**
* - GB models \ref Dune::Functions::Concept::GlobalBasis
**/
- template <class GB, class VT, class TP, bool is_const>
+ template <class Coeff, class GB, class TreePath = Dune::TypeTree::HybridTreePath<>>
class DiscreteFunction;
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- // Deduction guide for DiscreteFunction class
- template <class GB, class VT>
- DiscreteFunction(DOFVector<GB, VT> const& dofVector)
- -> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, true>;
-
- template <class GB, class VT>
- DiscreteFunction(DOFVector<GB, VT>& dofVector)
- -> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, false>;
+ template <class Coeff, class GB, class TreePath = Dune::TypeTree::HybridTreePath<>>
+ DiscreteFunction(Coeff&, GB const&, TreePath = {})
+ -> DiscreteFunction<Coeff, GB, TreePath>;
#endif
-
- /// A Generator for a const \ref DiscreteFunction
- template <class GlobalBasis, class ValueType,
- class PreTreePath = Dune::TypeTree::HybridTreePath<>>
- auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector,
- PreTreePath const& preTreePath = {})
+ /// A Generator for a mutable \ref DiscreteFunction
+ template <class Coeff, class GB, class Path,
+ class = void_t<decltype(std::declval<GB>().localView())> >
+ auto makeDiscreteFunction(Coeff& coefficients, GB const& basis, Path const& path)
{
- auto treePath = makeTreePath(preTreePath);
- return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), true>{dofVector, treePath};
+ auto treePath = makeTreePath(path);
+ return DiscreteFunction<Coeff, GB, decltype(treePath)>{coefficients, basis, treePath};
}
- /// A Generator for a mutable \ref DiscreteFunction
- template <class GlobalBasis, class ValueType,
- class PreTreePath = Dune::TypeTree::HybridTreePath<>>
- auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType>& dofVector,
- PreTreePath const& preTreePath = {})
+
+ /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
+ template <class Coeff, class GB, class TreePath>
+ class DiscreteFunction
+ : public DiscreteFunction<Coeff const,GB,TreePath>
{
- auto treePath = makeTreePath(preTreePath);
- return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), false>{dofVector, treePath};
- }
+ using Self = DiscreteFunction<Coeff,GB,TreePath>;
+ using Super = DiscreteFunction<Coeff const,GB,TreePath>;
+
+ using Coefficients = Coeff;
+ using GlobalBasis = GB;
+ using ValueType = typename Coeff::value_type;
+ public:
+ /// Constructor. Stores a pointer to the mutable `dofvector`.
+ DiscreteFunction(Coefficients& dofVector, GlobalBasis const& basis, TreePath const& treePath = {})
+ : Super(dofVector, basis, treePath)
+ , mutableCoeff_(&dofVector)
+ {}
+
+ public:
+ /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
+ /// reference to this DOFVector in the expression.
+ /**
+ * **Example:**
+ * ```
+ * auto v = makeDiscreteFunction(prob.solutionVector(),0);
+ * v.interpolate_noalias([](auto const& x) { return x[0]; });
+ * ```
+ **/
+ template <class Expr, class Tag = tag::average>
+ void interpolate_noalias(Expr&& expr, Tag strategy = {});
+ /// \brief Interpolation of GridFunction to DOFVector
+ /**
+ * **Example:**
+ * ```
+ * auto v = makeDiscreteFunction(prob.solutionVector(),0);
+ * v.interpolate(v + [](auto const& x) { return x[0]; });
+ * ```
+ * Allows to have a reference to the DOFVector in the expression, e.g. as
+ * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
+ **/
+ template <class Expr, class Tag = tag::average>
+ void interpolate(Expr&& expr, Tag strategy = {});
+
+ /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
+ template <class Expr>
+ Self& operator<<(Expr&& expr)
+ {
+ interpolate(FWD(expr));
+ return *this;
+ }
+
+ /// \brief interpolate `(*this) + expr` to DOFVector
+ template <class Expr>
+ Self& operator+=(Expr&& expr)
+ {
+ interpolate((*this) + expr);
+ return *this;
+ }
+
+ /// \brief interpolate `(*this) - expr` to DOFVector
+ template <class Expr>
+ Self& operator-=(Expr&& expr)
+ {
+ interpolate((*this) - expr);
+ return *this;
+ }
+
+ /// Return the mutable DOFVector
+ Coefficients& coefficients() { return *mutableCoeff_; }
+
+ /// Return the const DOFVector
+ using Super::coefficients;
+
+ template <class Path = RootTreePath>
+ auto child(Path const& path = {})
+ {
+ auto&& tp = makeTreePath(path);
+ return makeDiscreteFunction(*mutableCoeff_, this->basis(), cat(this->treePath_,tp));
+ }
+
+ using Super::child;
+
+ protected:
+ Coefficients* mutableCoeff_;
+ };
/// A Const DiscreteFunction
- template <class GB, class VT, class TP>
- class DiscreteFunction<GB,VT,TP,true>
+ template <class Coeff, class GB, class TreePath>
+ class DiscreteFunction<Coeff const,GB,TreePath>
{
private:
+ using Coefficients = std::remove_const_t<Coeff>;
using GlobalBasis = GB;
- using TreePath = TP;
+ using ValueType = typename Coeff::value_type;
using Tree = typename GlobalBasis::LocalView::Tree;
using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
-
using GridView = typename GlobalBasis::GridView;
public:
@@ -91,7 +157,7 @@ namespace AMDiS
using Domain = typename EntitySet::GlobalCoordinate;
/// Range type of this DiscreteFunction
- using Range = RangeType_t<SubTree,VT>;
+ using Range = RangeType_t<SubTree,ValueType>;
/// \brief This GridFunction has no derivative function, it can be created
/// by \ref DiscreteGridFunction.
@@ -111,11 +177,12 @@ namespace AMDiS
public:
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
- DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath = {})
- : dofVector_(&dofVector)
+ DiscreteFunction(Coefficients const& coefficients, GlobalBasis const& basis, TreePath const& treePath = {})
+ : coefficients_(&coefficients)
+ , basis_(&basis)
, treePath_(treePath)
- , entitySet_(dofVector.basis()->gridView())
- , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*dofVector.basis(), treePath))
+ , entitySet_(basis_->gridView())
+ , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*basis_, treePath))
{}
/// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
@@ -134,9 +201,9 @@ namespace AMDiS
}
/// \brief Return global basis bound to the DOFVector
- std::shared_ptr<GlobalBasis const> basis() const
+ GlobalBasis const& basis() const
{
- return dofVector_->basis();
+ return *basis_;
}
/// \brief Return treePath associated with this discrete function
@@ -146,114 +213,26 @@ namespace AMDiS
}
/// \brief Return const coefficient vector
- DOFVector<GB,VT> const& coefficients() const
+ Coefficients const& coefficients() const
{
- return *dofVector_;
+ return *coefficients_;
}
- template <class TreePath = RootTreePath>
- auto child(TreePath const& path = {}) const
+ template <class Path = RootTreePath>
+ auto child(Path const& path = {}) const
{
auto&& tp = makeTreePath(path);
- return makeDiscreteFunction(*dofVector_, cat(this->treePath_,tp));
+ return makeDiscreteFunction(*coefficients_, *basis_, cat(this->treePath_,tp));
}
protected:
- DOFVector<GB,VT> const* dofVector_;
+ Coefficients const* coefficients_;
+ GlobalBasis const* basis_;
TreePath treePath_;
EntitySet entitySet_;
NodeToRangeEntry nodeToRangeEntry_;
};
-
-
- /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
- template <class GB, class VT, class TP>
- class DiscreteFunction<GB,VT,TP,false>
- : public DiscreteFunction<GB, VT, TP, true>
- {
- using Self = DiscreteFunction<GB, VT, TP, false>;
- using Super = DiscreteFunction<GB, VT, TP, true>;
-
- using GlobalBasis = GB;
- using TreePath = TP;
-
- public:
- /// Constructor. Stores a pointer to the mutable `dofvector`.
- DiscreteFunction(DOFVector<GB,VT>& dofVector, TP const& treePath = {})
- : Super(dofVector, treePath)
- , mutableDofVector_(&dofVector)
- {}
-
- public:
- /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
- /// reference to this DOFVector in the expression.
- /**
- * **Example:**
- * ```
- * auto v = makeDiscreteFunction(prob.solutionVector(),0);
- * v.interpolate_noalias([](auto const& x) { return x[0]; });
- * ```
- **/
- template <class Expr, class Tag = tag::average>
- void interpolate_noalias(Expr&& expr, Tag strategy = {});
-
- /// \brief Interpolation of GridFunction to DOFVector
- /**
- * **Example:**
- * ```
- * auto v = makeDiscreteFunction(prob.solutionVector(),0);
- * v.interpolate(v + [](auto const& x) { return x[0]; });
- * ```
- * Allows to have a reference to the DOFVector in the expression, e.g. as
- * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
- **/
- template <class Expr, class Tag = tag::average>
- void interpolate(Expr&& expr, Tag strategy = {});
-
- /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
- template <class Expr>
- Self& operator<<(Expr&& expr)
- {
- interpolate(FWD(expr));
- return *this;
- }
-
- /// \brief interpolate `(*this) + expr` to DOFVector
- template <class Expr>
- Self& operator+=(Expr&& expr)
- {
- interpolate((*this) + expr);
- return *this;
- }
-
- /// \brief interpolate `(*this) - expr` to DOFVector
- template <class Expr>
- Self& operator-=(Expr&& expr)
- {
- interpolate((*this) - expr);
- return *this;
- }
-
- /// Return the mutable DOFVector
- DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
-
- /// Return the const DOFVector
- using Super::coefficients;
-
- template <class TreePath = RootTreePath>
- auto child(TreePath const& path = {})
- {
- auto&& tp = makeTreePath(path);
- return makeDiscreteFunction(*mutableDofVector_, cat(this->treePath_,tp));
- }
-
- using Super::child;
-
- protected:
- DOFVector<GB,VT>* mutableDofVector_;
- };
-
} // end namespace AMDiS
#include "DiscreteLocalFunction.inc.hpp"
diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
index 3759909c..83f013a2 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
@@ -10,14 +10,14 @@
namespace AMDiS {
// Evaluate DiscreteFunction in global coordinates
-template <class GB, class VT, class TP>
-typename DiscreteFunction<GB,VT,TP,true>::Range DiscreteFunction<GB,VT,TP,true>::
+template <class Coeff, class GB, class TP>
+typename DiscreteFunction<Coeff const,GB,TP>::Range DiscreteFunction<Coeff const,GB,TP>::
operator()(Domain const& x) const
{
using Grid = typename GlobalBasis::GridView::Grid;
using IS = typename GlobalBasis::GridView::IndexSet;
- auto const& gv = this->basis()->gridView();
+ auto const& gv = this->basis().gridView();
Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
auto element = hsearch.findEntity(x);
@@ -29,21 +29,21 @@ typename DiscreteFunction<GB,VT,TP,true>::Range DiscreteFunction<GB,VT,TP,true>:
// Interpolation of GridFunction to DOFVector
-template <class GB, class VT, class TP>
+template <class Coeff, class GB, class TP>
template <class Expr, class Tag>
-void DiscreteFunction<GB,VT,TP,false>::
+void DiscreteFunction<Coeff,GB,TP>::
interpolate_noalias(Expr&& expr, Tag strategy)
{
- auto const& basis = *this->basis();
+ auto const& basis = this->basis();
auto const& treePath = this->treePath();
auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
if (std::is_same<Tag, tag::average>::value) {
- auto counter = coefficients();
+ VectorType_t<short,Coeff> counter(basis);
AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
- coefficients().forEach([&counter](std::size_t dof, auto& coeff)
+ coefficients().forEach([&counter](auto dof, auto& coeff)
{
coeff /= std::max(double(counter.at(dof)), 1.0);
});
@@ -54,19 +54,19 @@ void DiscreteFunction<GB,VT,TP,false>::
// Interpolation of GridFunction to DOFVector
-template <class GB, class VT, class TP>
+template <class Coeff, class GB, class TP>
template <class Expr, class Tag>
-void DiscreteFunction<GB,VT,TP,false>::
+void DiscreteFunction<Coeff,GB,TP>::
interpolate(Expr&& expr, Tag strategy)
{
// create temporary copy of data
- DOFVector<GB,VT> tmp(coefficients());
+ Coeff tmp(coefficients());
- Self tmpView{tmp, this->treePath()};
+ Self tmpView{tmp, this->basis(), this->treePath()};
tmpView.interpolate_noalias(FWD(expr), strategy);
// move data from temporary vector into stored DOFVector
- coefficients().backend() = std::move(tmp.backend());
+ coefficients() = std::move(tmp);
}
} // end namespace AMDiS
diff --git a/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp
index 899818c6..3380bb41 100644
--- a/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp
+++ b/src/amdis/gridfunctions/DiscreteLocalFunction.inc.hpp
@@ -9,8 +9,8 @@
namespace AMDiS {
-template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP,true>::LocalFunction
+template <class Coeff, class GB, class TP>
+class DiscreteFunction<Coeff const,GB,TP>::LocalFunction
{
public:
using Domain = typename EntitySet::LocalCoordinate;
@@ -27,14 +27,14 @@ public:
/// Constructor. Stores a copy of the DiscreteFunction.
LocalFunction(DiscreteFunction const& globalFunction)
: globalFunction_(globalFunction)
- , localView_(globalFunction_.basis()->localView())
+ , localView_(globalFunction_.basis().localView())
, subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
LocalFunction(LocalFunction const& other)
: globalFunction_(other.globalFunction_)
- , localView_(globalFunction_.basis()->localView())
+ , localView_(globalFunction_.basis().localView())
, subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath()))
{}
@@ -126,14 +126,14 @@ private:
LocalView localView_;
SubTree const* subTree_;
- std::vector<VT> localCoefficients_;
+ std::vector<ValueType> localCoefficients_;
bool bound_ = false;
};
-template <class GB, class VT, class TP>
+template <class Coeff, class GB, class TP>
template <class Type>
-class DiscreteFunction<GB,VT,TP,true>::DerivativeLocalFunctionBase
+class DiscreteFunction<Coeff const,GB,TP>::DerivativeLocalFunctionBase
{
using R = typename DiscreteFunction::Range;
using D = typename DiscreteFunction::Domain;
@@ -156,7 +156,7 @@ public:
DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
: globalFunction_(globalFunction)
, type_(type)
- , localView_(globalFunction_.basis()->localView())
+ , localView_(globalFunction_.basis().localView())
, subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath()))
{}
@@ -164,7 +164,7 @@ public:
DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
: globalFunction_(other.globalFunction_)
, type_(other.type_)
- , localView_(globalFunction_.basis()->localView())
+ , localView_(globalFunction_.basis().localView())
, subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_.treePath()))
{}
@@ -214,13 +214,13 @@ protected:
LocalView localView_;
SubTree const* subTree_;
Dune::Std::optional<Geometry> geometry_;
- std::vector<VT> localCoefficients_;
+ std::vector<ValueType> localCoefficients_;
bool bound_ = false;
};
-template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP,true>::GradientLocalFunction
+template <class Coeff, class GB, class TP>
+class DiscreteFunction<Coeff const,GB,TP>::GradientLocalFunction
: public DerivativeLocalFunctionBase<tag::gradient>
{
using Super = DerivativeLocalFunctionBase<tag::gradient>;
@@ -271,8 +271,8 @@ public:
};
-template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP,true>::DivergenceLocalFunction
+template <class Coeff, class GB, class TP>
+class DiscreteFunction<Coeff const,GB,TP>::DivergenceLocalFunction
: public DerivativeLocalFunctionBase<tag::divergence>
{
using Super = DerivativeLocalFunctionBase<tag::divergence>;
@@ -326,8 +326,8 @@ private:
};
-template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP,true>::PartialLocalFunction
+template <class Coeff, class GB, class TP>
+class DiscreteFunction<Coeff const,GB,TP>::PartialLocalFunction
: public DerivativeLocalFunctionBase<tag::partial>
{
using Super = DerivativeLocalFunctionBase<tag::partial>;
diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp
index 2c74bb15..77bf0154 100644
--- a/src/amdis/io/FileWriterCreator.hpp
+++ b/src/amdis/io/FileWriterCreator.hpp
@@ -3,6 +3,8 @@
#include <memory>
#include <string>
+#include <dune/common/typeutilities.hh>
+
#include <amdis/BoundaryManager.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/Output.hpp>
@@ -45,15 +47,15 @@ namespace AMDiS
create(std::string type, std::string prefix, TreePath treePath = {}) const
{
auto data = makeDiscreteFunction(*systemVector_, treePath);
- using Range = typename TYPEOF(data)::Range;
-
- return create_impl(std::move(type), std::move(prefix), data, ValueCategory_t<Range>{});
+ return create_impl(std::move(type), std::move(prefix), data, Dune::PriorityTag<42>{});
}
private:
- template <class Data, class ValCat>
+ template <class Data,
+ REQUIRES(not std::is_same<ValueCategory_t<typename Data::Range>,tag::unknown>::value),
+ REQUIRES(std::is_arithmetic<typename Dune::FieldTraits<typename Data::Range>::field_type>::value)>
std::unique_ptr<FileWriterInterface>
- create_impl(std::string type, std::string prefix, Data const& data, ValCat) const
+ create_impl(std::string type, std::string prefix, Data const& data, Dune::PriorityTag<2>) const
{
GridView const& gridView = systemVector_->basis()->gridView();
@@ -90,10 +92,13 @@ namespace AMDiS
return {};
}
+
// The value-category is unknown, like a composite/hierarchic vector or any unknown type.
- template <class Data>
+ template <class Data,
+ REQUIRES(std::is_same<ValueCategory_t<typename Data::Range>,tag::unknown>::value),
+ REQUIRES(std::is_arithmetic<typename Dune::FieldTraits<typename Data::Range>::field_type>::value)>
std::unique_ptr<FileWriterInterface>
- create_impl(std::string type, std::string prefix, Data const& /*data*/, tag::unknown) const
+ create_impl(std::string type, std::string prefix, Data const& /*data*/, Dune::PriorityTag<1>) const
{
// Backup writer, writing the grid and the solution vector
if (type == "backup")
@@ -107,6 +112,15 @@ namespace AMDiS
return {};
}
+ // fallback implementation
+ template <class Data>
+ std::unique_ptr<FileWriterInterface>
+ create_impl(std::string, std::string, Data const&, Dune::PriorityTag<0>) const
+ {
+ error_exit("Filewriter cannot be used with unsupported Range and field_type");
+ return {};
+ }
+
private:
std::shared_ptr<SystemVector> systemVector_;
std::shared_ptr<BoundaryManager<Grid>> boundaryManager_ = nullptr;
diff --git a/src/amdis/linearalgebra/CMakeLists.txt b/src/amdis/linearalgebra/CMakeLists.txt
index fb2d7db9..d89ab399 100644
--- a/src/amdis/linearalgebra/CMakeLists.txt
+++ b/src/amdis/linearalgebra/CMakeLists.txt
@@ -6,14 +6,16 @@ install(FILES
DOFMapping.inc.hpp
LinearSolver.hpp
LinearSolverInterface.hpp
- MatrixBase.hpp
+ MatrixFacade.hpp
ParallelIndexSet.hpp
PreconditionerInterface.hpp
RunnerInterface.hpp
SolverInfo.hpp
+ SparsityPattern.hpp
SymmetryStructure.hpp
Traits.hpp
- VectorBase.hpp
+ VectorDefaultImplementation.hpp
+ VectorFacade.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra)
if (BACKEND STREQUAL "MTL")
diff --git a/src/amdis/linearalgebra/Constraints.hpp b/src/amdis/linearalgebra/Constraints.hpp
index 03498fd1..b03dee88 100644
--- a/src/amdis/linearalgebra/Constraints.hpp
+++ b/src/amdis/linearalgebra/Constraints.hpp
@@ -5,7 +5,7 @@
namespace AMDiS
{
// forward declaration
- template <class RB, class CB, class T>
+ template <class RB, class CB, class T, class Traits>
class BiLinearForm;
@@ -40,21 +40,21 @@ namespace AMDiS
}
- template <class RB, class CB, class T>
- struct Constraints<BiLinearForm<RB,CB,T>>
+ template <class RB, class CB, class T, class Traits>
+ struct Constraints<BiLinearForm<RB,CB,T,Traits>>
{
- using Matrix = BiLinearForm<RB,CB,T>;
+ using Matrix = BiLinearForm<RB,CB,T,Traits>;
template <class Sol, class Rhs, class BitVec>
static void dirichletBC(Matrix& matrix, Sol& solution, Rhs& rhs, BitVec const& nodes, bool setDiagonal = true)
{
- AMDiS::dirichletBC(matrix.backend(), solution.backend(), rhs.backend(), nodes, setDiagonal);
+ AMDiS::dirichletBC(matrix.impl(), solution.impl(), rhs.impl(), nodes, setDiagonal);
}
template <class Sol, class Rhs, class BitVec, class Assoc>
static void periodicBC(Matrix& matrix, Sol& solution, Rhs& rhs, BitVec const& left, Assoc const& association, bool setDiagonal = true)
{
- AMDiS::periodicBC(matrix.backend(), solution.backend(), rhs.backend(), left, association, setDiagonal);
+ AMDiS::periodicBC(matrix.impl(), solution.impl(), rhs.impl(), left, association, setDiagonal);
}
};
diff --git a/src/amdis/linearalgebra/DOFMapping.hpp b/src/amdis/linearalgebra/DOFMapping.hpp
index 38153fff..9192e4a5 100644
--- a/src/amdis/linearalgebra/DOFMapping.hpp
+++ b/src/amdis/linearalgebra/DOFMapping.hpp
@@ -16,6 +16,7 @@
namespace AMDiS
{
+ /// \brief Fallback for \ref ParallelDofMapping in case there is only one mpi core.
template <class IS, class GI = std::size_t>
class SequentialDofMapping
{
@@ -152,6 +153,7 @@ namespace AMDiS
#if HAVE_MPI
+ /// \brief Mapping of local to global indices
template <class PIS, class GI = std::size_t>
class ParallelDofMapping
{
diff --git a/src/amdis/linearalgebra/MatrixBase.hpp b/src/amdis/linearalgebra/MatrixBase.hpp
deleted file mode 100644
index 38481bda..00000000
--- a/src/amdis/linearalgebra/MatrixBase.hpp
+++ /dev/null
@@ -1,133 +0,0 @@
-#pragma once
-
-#include <memory>
-#include <type_traits>
-
-#include <dune/common/shared_ptr.hh>
-
-#include <amdis/common/Concepts.hpp>
-#include <amdis/common/TypeTraits.hpp>
-#include <amdis/functions/NodeIndices.hpp>
-#include <amdis/linearalgebra/SymmetryStructure.hpp>
-#include <amdis/typetree/MultiIndex.hpp>
-
-namespace AMDiS
-{
- /**
- * Basis implementation of DOFMatrix, i.e. a sparse matrix storing all the
- * assembled Operators indexed with DOF indices. The matrix data is associated
- * to a row and column global basis.
- *
- * \tparam RB Basis of the matrix rows
- * \tparam CB Basis of matrix columns
- * \tparam B A linear-algebra backend for the matrix storage
- **/
- template <class RB, class CB, class B>
- class MatrixBase
- {
- using Self = MatrixBase;
-
- public:
- /// The type of the finite element space / basis of the row
- using RowBasis = RB;
- using RowLocalView = typename RowBasis::LocalView;
-
- /// The type of the finite element space / basis of the column
- using ColBasis = CB;
- using ColLocalView = typename ColBasis::LocalView;
-
- /// The Linear-Algebra backend used to store the assembled coefficients
- using Backend = B;
-
- public:
- /// (1) Constructor. Stores the shared_ptr to the bases and the communication object.
- MatrixBase(std::shared_ptr<RB> rowBasis, std::shared_ptr<CB> colBasis)
- : rowBasis_(std::move(rowBasis))
- , colBasis_(std::move(colBasis))
- , backend_(*rowBasis_, *colBasis_)
- {}
-
- /// (2) Constructor. Forwards to (1) by wrapping reference into non-destroying
- /// shared_ptr, see \ref Dune::wrap_or_move.
- template <class RB_, class CB_,
- REQUIRES(Concepts::Similar<Types<RB_,CB_>, Types<RB,CB>>)>
- MatrixBase(RB_&& rowBasis, CB_&& colBasis)
- : MatrixBase(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
- {}
-
- /// Return the row-basis \ref rowBasis of the matrix
- std::shared_ptr<RowBasis> const& rowBasis() const
- {
- return rowBasis_;
- }
-
- /// Return the col-basis \ref colBasis of the matrix
- std::shared_ptr<ColBasis> const& colBasis() const
- {
- return colBasis_;
- }
-
- /// Return the underlying linear algebra backend
- Backend const& backend() const { return backend_; }
- Backend& backend() { return backend_; }
-
- /// \brief Initialize the matrix for insertion, i.e. allocate the non-zero pattern
- /**
- * With the optional parameter \p symmetry some additional information about the
- * structure of the values or the sparsity pattern can be provided.
- * See \ref SymmetryStructure.
- **/
- void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
- {
- backend_.init(*rowBasis_, *colBasis_, symmetry);
- }
-
- /// Finish the matrix insertion, e.g. cleanup or final insertion
- void finish()
- {
- backend_.finish();
- }
-
- /// Insert a single value into the matrix (add to existing value)
- template <class RowIndex, class ColIndex>
- void insert(RowIndex const& row, ColIndex const& col, typename Backend::value_type const& value)
- {
- backend_.insert(row, col, value);
- }
-
- /// Insert a block of values into the sparse matrix (add to existing values)
- /// The global matrix indices are determined by the corresponding localviews.
- template <class LocalMatrix>
- void scatter(RowLocalView const& r, ColLocalView const& c, LocalMatrix const& localMatrix)
- {
- assert(r.size() * c.size() == localMatrix.size());
- assert(r.size() == localMatrix.rows());
- assert(c.size() == localMatrix.cols());
-
- const bool optimized = std::is_same<RowLocalView,ColLocalView>::value
- && r.tree().treeIndex() == c.tree().treeIndex();
-
- if (optimized)
- backend_.scatter(nodeIndices(r), localMatrix);
- else
- backend_.scatter(nodeIndices(r), nodeIndices(c), localMatrix);
- }
-
- /// Number of nonzeros in the matrix
- std::size_t nnz() const
- {
- return backend_.nnz();
- }
-
- protected:
- /// The finite element space / basis associated with the rows
- std::shared_ptr<RowBasis> rowBasis_;
-
- /// The finite element space / basis associated with the columns
- std::shared_ptr<ColBasis> colBasis_;
-
- /// Data backend
- Backend backend_;
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/MatrixFacade.hpp b/src/amdis/linearalgebra/MatrixFacade.hpp
new file mode 100644
index 00000000..49adf77c
--- /dev/null
+++ b/src/amdis/linearalgebra/MatrixFacade.hpp
@@ -0,0 +1,107 @@
+#pragma once
+
+#include <type_traits>
+
+#include <dune/common/shared_ptr.hh>
+#include <dune/functions/functionspacebases/concepts.hh>
+
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/functions/NodeIndices.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
+
+namespace AMDiS
+{
+ /**
+ * Basis implementation of DOFMatrix, i.e. a sparse matrix storing all the
+ * assembled Operators indexed with DOF indices. The matrix data is associated
+ * to a row and column global basis.
+ *
+ * \tparam T The coefficient type of the matrix
+ * \tparam Pattern The type of the sparsity pattern
+ * \tparam MatrixImpl A linear-algebra backend for the matrix storage
+ **/
+ template <class T, class Pattern, template <class> class MatrixImpl>
+ class MatrixFacade
+ {
+ using Self = MatrixFacade;
+ using Impl = MatrixImpl<T>;
+
+ /// Type of the sparsity pattern of the backend
+ using SparsityPattern = Pattern;
+
+ public:
+ /// Constructor. Forwards the bases to the implementation class and
+ /// constructs a matrix sparsity pattern.
+ template <class RowBasis, class ColBasis>
+ MatrixFacade(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : impl_(rowBasis, colBasis)
+ , pattern_(rowBasis, colBasis)
+ {}
+
+ /// Return the underlying matrix backend
+ Impl const& impl() const { return impl_; }
+ Impl& impl() { return impl_; }
+
+ /// \brief Initialize the matrix for insertion, i.e. allocate the non-zero pattern
+ /**
+ * With the optional parameter \p symmetry some additional information about the
+ * structure of the values or the sparsity pattern can be provided.
+ * See \ref SymmetryStructure.
+ **/
+ void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
+ {
+ impl_.init(pattern_, symmetry);
+ }
+
+ /// Finish the matrix insertion, e.g. cleanup or final insertion
+ void finish()
+ {
+ impl_.finish();
+ }
+
+ /// Insert a single value into the matrix (add to existing value)
+ template <class RowIndex, class ColIndex,
+ REQUIRES(Concepts::MultiIndex<RowIndex>),
+ REQUIRES(Concepts::MultiIndex<ColIndex>)>
+ void insert(RowIndex const& row, ColIndex const& col, typename Impl::value_type const& value)
+ {
+ impl_.insert(row, col, value);
+ }
+
+ /// Insert a block of values into the sparse matrix (add to existing values)
+ /// The global matrix indices are determined by the corresponding localviews.
+ template <class RowLocalView, class ColLocalView, class LocalMatrix,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename RowLocalView::GlobalBasis>, RowLocalView>()),
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename ColLocalView::GlobalBasis>, ColLocalView>())>
+ void scatter(RowLocalView const& r, ColLocalView const& c, LocalMatrix const& localMatrix)
+ {
+ assert(r.size() * c.size() == localMatrix.size());
+ assert(r.size() == localMatrix.rows());
+ assert(c.size() == localMatrix.cols());
+
+ const bool optimized = std::is_same<RowLocalView,ColLocalView>::value
+ && r.tree().treeIndex() == c.tree().treeIndex();
+
+ if (optimized)
+ impl_.scatter(nodeIndices(r), localMatrix);
+ else
+ impl_.scatter(nodeIndices(r), nodeIndices(c), localMatrix);
+ }
+
+ /// Number of nonzeros in the matrix
+ std::size_t nnz() const
+ {
+ return impl_.nnz();
+ }
+
+ protected:
+ /// The matrix backend
+ Impl impl_;
+
+ /// The structure of the non-zeros in the matrix
+ SparsityPattern pattern_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/SparsityPattern.hpp b/src/amdis/linearalgebra/SparsityPattern.hpp
new file mode 100644
index 00000000..93efde60
--- /dev/null
+++ b/src/amdis/linearalgebra/SparsityPattern.hpp
@@ -0,0 +1,155 @@
+#pragma once
+
+#include <memory>
+
+#include <dune/istl/matrixindexset.hh>
+
+#include <amdis/Observer.hpp>
+#include <amdis/common/Index.hpp>
+
+namespace AMDiS
+{
+ /// \brief A general sparsity pattern implementation using the full pattern of the
+ /// basis by adding all local indices
+ template <class RowBasis, class ColBasis>
+ class SparsityPattern
+ : private ObserverSequence<event::adapt,2>
+ {
+ public:
+ /// Constructor. Stores pointers to the passed bases.
+ SparsityPattern(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
+ , rowBasis_(&rowBasis)
+ , colBasis_(&colBasis)
+ {
+ updateImpl3();
+ }
+
+ /// Number of rows in the matrix
+ std::size_t rows() const
+ {
+ return rows_;
+ }
+
+ /// Number of columns in the matrix
+ std::size_t cols() const
+ {
+ return cols_;
+ }
+
+ /// Average number of non-zeros per row
+ std::size_t avgRowSize() const
+ {
+ assert(rows_ == pattern_.rows());
+ std::size_t sum = 0;
+ for (std::size_t r = 0; r < rows_; ++r)
+ sum += rowSize(r);
+ return (sum + rows_ - 1) / rows_;
+ }
+
+ /// Number of non-zeros in row r
+ std::size_t rowSize(std::size_t r) const
+ {
+ return pattern_.rowsize(r);
+ }
+
+ /// Total number of non-zeros
+ std::size_t nnz() const
+ {
+ return pattern_.size();
+ }
+
+ /// Initialize a matrix with the non-zero structure
+ template <class Matrix>
+ void applyTo(Matrix& matrix) const
+ {
+ pattern_.exportIdx(matrix);
+ }
+
+ protected:
+
+ void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
+ void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }
+
+ private:
+
+ // Track for each basis wether updated and if both are, call updateImpl3()
+ template <std::size_t I>
+ void updateImpl2(event::adapt, index_t<I>)
+ {
+ assert(!updateCounter_.test(I));
+ updateCounter_.set(I);
+
+ if (updateCounter_.all()) {
+ updateImpl3();
+ updateCounter_.reset();
+ }
+ }
+
+ void updateImpl3()
+ {
+ rowBasis_ == colBasis_
+ ? updateSameBasis()
+ : updateDifferentBasis();
+ }
+
+ // Update pattern when basis is updated. This method is called if rowBasis == colBasis.
+ void updateSameBasis()
+ {
+ rows_ = rowBasis_->dimension();
+ cols_ = rows_;
+ pattern_.resize(rows_, cols_);
+
+ auto localView = rowBasis_->localView();
+ for (const auto& element : elements(rowBasis_->gridView())) {
+ localView.bind(element);
+
+ for (std::size_t i = 0, size = localView.size(); i < size; ++i) {
+ // The global index of the i-th vertex of the element
+ auto row = localView.index(i);
+ for (std::size_t j = 0; j < size; ++j) {
+ // The global index of the j-th vertex of the element
+ auto col = localView.index(j);
+ pattern_.add(row, col);
+ }
+ }
+ }
+ }
+
+ // Update pattern when basis is updated. This method is called if rowBasis != colBasis.
+ void updateDifferentBasis()
+ {
+ rows_ = rowBasis_->dimension();
+ cols_ = colBasis_->dimension();
+ pattern_.resize(rows_, cols_);
+
+ auto rowLocalView = rowBasis_->localView();
+ auto colLocalView = colBasis_->localView();
+ for (const auto& element : elements(rowBasis_->gridView())) {
+ rowLocalView.bind(element);
+ colLocalView.bind(element);
+
+ for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
+ // The global index of the i-th vertex of the element
+ auto row = rowLocalView.index(i);
+ 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);
+ pattern_.add(row, col);
+ }
+ }
+ }
+ }
+
+ private:
+ RowBasis const* rowBasis_;
+ ColBasis const* colBasis_;
+
+ std::bitset<2> updateCounter_ = 0;
+
+ std::size_t rows_;
+ std::size_t cols_;
+ Dune::MatrixIndexSet pattern_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/VectorBase.hpp b/src/amdis/linearalgebra/VectorBase.hpp
index 25fc10ed..a0668a9a 100644
--- a/src/amdis/linearalgebra/VectorBase.hpp
+++ b/src/amdis/linearalgebra/VectorBase.hpp
@@ -1,445 +1,111 @@
#pragma once
-#include <memory>
-#include <string>
-#include <type_traits>
-#include <utility>
+#include <iterator>
-#include <dune/common/classname.hh>
-#include <dune/common/deprecated.hh>
-#include <dune/common/hybridutilities.hh>
#include <dune/common/rangeutilities.hh>
-#include <dune/common/shared_ptr.hh>
-#include <dune/common/std/type_traits.hh>
-#include <dune/functions/functionspacebases/concepts.hh>
-#include <dune/functions/functionspacebases/sizeinfo.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
-#include <amdis/Output.hpp>
-#include <amdis/common/Concepts.hpp>
#include <amdis/common/FakeContainer.hpp>
-#include <amdis/common/TypeTraits.hpp>
-#include <amdis/functions/NodeIndices.hpp>
-#include <amdis/operations/Assigner.hpp>
-#include <amdis/typetree/MultiIndex.hpp>
+#include <amdis/utility/MappedRangeView.hpp>
namespace AMDiS
{
- /// \brief The basic container that stores a base vector and a corresponding basis
- /**
- * A vector storing all the assembled Operators indexed with DOF indices. The
- * vector data is associated to a global basis.
- *
- * \tparam GB Basis of the vector
- * \tparam B A linear algebra backend implementing the storage and operations.
- **/
- template <class GB, class B>
+ /// \brief CRTP base class for flat vector backends
+ template <class Derived>
class VectorBase
{
- using Self = VectorBase;
-
- public:
- /// The type of the functionspace basis associated to this linearform
- using GlobalBasis = GB;
- using LocalView = typename GlobalBasis::LocalView;
-
- /// The Linear-Algebra backend used to store the assembled coefficients
- using Backend = B;
-
private:
- // proxy class redirecting the mutable element access to the insertValue method
- template <class Index>
- struct AccessProxy;
-
- enum class VectorState
- {
- unknown = 0,
- synchronized = 1,
- insert_values = 2,
- add_values = 3
- };
-
- template <class A>
- using VectorStateOf_t = std::conditional_t<std::is_same<A,Assigner::plus_assign>::value,
- std::integral_constant<VectorState, VectorState::add_values>,
- std::integral_constant<VectorState, VectorState::insert_values>>;
-
- public:
- /// (1) Constructor. Stores the shared_ptr of the basis.
- VectorBase(std::shared_ptr<GB> basis)
- : basis_(std::move(basis))
- , backend_(*basis_)
- {
- resizeZero();
- }
-
- /// (2) Constructor. Forwards to (1) by wrapping reference into non-destroying shared_ptr.
- template <class GB_,
- REQUIRES(Concepts::Similar<GB_,GB>)>
- VectorBase(GB_&& basis)
- : VectorBase(Dune::wrap_or_move(FWD(basis)))
- {}
-
- /// Return the basis \ref basis_ associated to the vector
- std::shared_ptr<GlobalBasis> const& basis() const { return basis_; }
- std::shared_ptr<GlobalBasis> const& basis() { return basis_; }
-
- /// Return the underlying linear algebra backend
- Backend const& backend() const { return backend_; }
- Backend& backend() { return backend_; }
-
-
- template <class B_>
- using HasLocalSize = decltype(std::declval<B_>().localSize());
-
- template <class B_>
- using HasGlobalSize = decltype(std::declval<B_>().globalSize());
-
- /// Return the number of entries in the local part of the vector
- std::size_t localSize() const
- {
- return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasLocalSize,Backend>{},
- [&](auto id) { return id(backend_).localSize(); },
- [&](auto id) { return id(backend_).size(); });
- }
-
- /// Return the number of entries in the global vector
- std::size_t globalSize() const
- {
- return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalSize,Backend>{},
- [&](auto id) { return id(backend_).globalSize(); },
- [&](auto id) { return id(backend_).size(); });
- }
-
- /// Resize the \ref vector to the size of the \ref basis
- void resize()
+ Derived const& asDerived() const
{
- backend_.init(sizeInfo(*basis_), false);
- state_ = VectorState::unknown;
+ return static_cast<Derived const&>(*this);
}
- /// Resize the \ref vector to the size of the \ref basis and set to zero
- void resizeZero()
+ Derived& asDerived()
{
- backend_.init(sizeInfo(*basis_), true);
- state_ = VectorState::synchronized;
+ return static_cast<Derived&>(*this);
}
- /// Prepare the vector for insertion of values, finish the insertion with
- /// \ref finish().
- void init(bool clear)
- {
- backend_.init(sizeInfo(*basis_), clear);
- state_ = clear ? VectorState::synchronized : VectorState::unknown;
- }
-
- /// Finish the insertion of values, see \ref init()
- void finish()
- {
- backend_.finish();
- state_ = VectorState::unknown;
- }
-
- /// Access the entry \p i of the \ref vector with read-only-access.
- template <class Index>
- typename Backend::value_type DUNE_DEPRECATED_MSG("Do not use element-wise vector access")
- operator[](Index const& idx) const
- {
- assert(state_ == VectorState::synchronized || state_ == VectorState::unknown);
- return at(idx);
- }
-
- /// Access the entry \p i of the \ref vector with write-access. Uses a proxy
- /// class the redirects to \ref insertValue.
- template <class Index>
- AccessProxy<Index> DUNE_DEPRECATED_MSG("Do not use element-wise vector access")
- operator[](Index const& idx)
- {
- return AccessProxy<Index>{this, idx};
- }
-
- /// Return the value of the vector at the local index \ref idx
- template <class Index>
- typename Backend::value_type at(Index const& idx) const
- {
- const_cast<Self*>(this)->synchronize();
-
- return backend_.at(flatMultiIndex(idx));
- }
-
- /// \brief Insert a single value into the matrix (add or overwrite to existing value)
- /**
- * Inserts or adds a value into a certain location \p dof (given as dof multi-index)
- * of a vector. The insertion mode is determined by the \p assign functor. Use
- * \ref Assigner::plus_assign for adding values (default) or \ref Assigner::assign
- * for overwriting (setting) values. Different insertion modes can not be mixed!
- *
- * Insertion must be closed with a call to \ref finish().
- *
- * [[not collective]]
- */
- template <class Index, class Assign = Assigner::plus_assign>
- void insert(Index const& dof, typename Backend::value_type const& value, Assign assign = {})
- {
- test_exit_dbg(state_ == VectorStateOf_t<Assign>::value ||
- state_ == VectorState::unknown ||
- state_ == VectorState::synchronized,
- "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
-
- backend_.insert(flatMultiIndex(dof), value, assign);
+ protected:
+ /// Protected constructor prevents from direct construction of this class.
+ /// Only use it as base class in the sense of CRTP.
+ VectorBase() = default;
- // set the state to insert_values or add_values
- state_ = VectorStateOf_t<Assign>::value;
- }
-
- /// See \ref insert for assignment operation \ref Assigner::assign
- template <class Index>
- void set(Index const& dof, typename Backend::value_type const& value)
- {
- insert(dof, value, Assigner::assign{});
- }
+ public:
+ void finish() { /* do nothing */ }
+ void synchronize() { /* do nothing */ }
- /// See \ref insert for assignment operation \ref Assigner::plus_assign
- template <class Index>
- void add(Index const& dof, typename Backend::value_type const& value)
+ /// Insert or add \p value at position \p idx, using the passed \p assign functor.
+ template <class MultiIndex, class ValueType, class Assign>
+ void insert(MultiIndex const& idx, ValueType const& value, Assign assign)
{
- insert(dof, value, Assigner::plus_assign{});
+ assign(value, asDerived().at(idx));
}
-
- /// \brief Extract values from the vector referring to the given local indices
- /// and store it into a buffer
- /**
- * Collect value of indices and store them into a buffer. The buffer must be
- * a vector-like container with `buffer.resize()` and `buffer.begin()`. The
- * indices must be stored in an iterable container.
- *
- * If the vector is not in synchronized state, collects all necessary values
- * possibly from neighbouring processors.
- *
- * [[expects: localView is bound to an element]]
- * [[expects: node is in localView.tree()]]
- * [[possibly neighbor-wise collective operation]]
- */
- template <class Node, class Buffer,
- REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
- void gather(LocalView const& localView, Node const& node, Buffer& buffer) const
+ /// Gather values from the vector into the output range \p buffer decribed as an output iterator
+ template <class IndexRange, class OutputIterator>
+ void gather(IndexRange const& localInd, OutputIterator buffer) const
{
- test_exit(state_ == VectorState::unknown ||
- state_ == VectorState::synchronized,
- "Vector in invalid state {} for gather operations.", to_string(state_));
-
- const_cast<Self*>(this)->synchronize();
-
- buffer.resize(node.size());
- backend_.gather(nodeIndices(localView, node), buffer.begin());
+ for (auto const& idx : localInd)
+ *buffer++ = asDerived().at(idx);
}
- // [[expects: localView is bound to an element]]
- template <class Buffer>
- void gather(LocalView const& localView, Buffer& buffer) const
- {
- gather(localView, localView.tree(), buffer);
- }
- /// Insert a block of values into the vector (add or overwrite to existing values)
- /**
- * Inserts or adds values into certain locations of a vector. Insertion indices
- * are extracted from the given \p localView. The insertion mode is determined
- * by the \p assign functor. Use \ref Assigner::plus_assign for adding values
- * (default) or \ref Assigner::assign for overwriting (setting) values. Different
- * insertion modes can not be mixed! The \p localVector is assumed to be a continuous
- * memory container with a `data()` method to get a pointer to the beginning.
- *
- * The \p mask models a boolean range with at least a `begin()` method. Must
- * be forward iterable for at least `localVector.size()` elements. Does not
- * need an `end()` method. See, e.g. \ref FakeContainer.
- *
- * Insertion must be closed with a call to \ref finish(). It is not allowed
- * to switch insertion mode before calling `finish()`.
- *
- * [[expects: localView is bound to an element]]
- * [[expects: node is in localView.tree()]]
- * [[not collective]]
- */
- template <class Node, class NodeVector, class MaskRange, class Assign,
- REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
- void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector,
- MaskRange const& mask, Assign assign)
+ /// Scatter the values from the local \p vec into the container using the passed \p assign functor.
+ template <class IndexRange, class LocalVec, class Assign>
+ void scatter(IndexRange const& localInd, LocalVec const& vec, FakeContainer<bool,true>, Assign assign)
{
- test_exit(state_ == VectorStateOf_t<Assign>::value ||
- state_ == VectorState::unknown ||
- state_ == VectorState::synchronized,
- "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
-
- assert(localVector.size() == node.size());
-
- // create a range of DOF indices on the node
- backend_.scatter(nodeIndices(localView, node), localVector, mask, assign);
-
- // set the state to insert_values or add_values
- state_ = VectorStateOf_t<Assign>::value;
+ auto vec_it = std::begin(vec);
+ for (auto const& idx : localInd)
+ assign(*vec_it++, asDerived().at(idx));
}
- /// Call \ref scatter with `MaskRange` set to \ref FakeContainer.
- // [[expects: localView is bound to an element]]
- // [[expects: node is in localView.tree()]]
- template <class Node, class NodeVector, class Assign,
- REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
- void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector, Assign assign)
+ /// Scatter some values from the local \p vec into the container using the passed \p assign functor.
+ /// The values to scatter are selected using the \p mask range that needs to be forward iterable.
+ template <class IndexRange, class LocalVec, class MaskRange, class Assign>
+ void scatter(IndexRange const& localInd, LocalVec const& vec, MaskRange const& mask, Assign assign)
{
- scatter(localView, node, localVector, FakeContainer<bool,true>{}, assign);
+ auto vec_it = std::begin(vec);
+ auto mask_it = std::begin(mask);
+ auto ind_it = std::begin(localInd);
+ for (; vec_it != std::end(vec); ++vec_it, ++mask_it, ++ind_it) {
+ if (*mask_it)
+ assign(*vec_it, asDerived().at(*ind_it));
+ }
}
- /// Call \ref scatter with `Node` given by the tree of the \p localView.
- // [[expects: localView is bound to an element]]
- template <class LocalVector, class Assign>
- void scatter(LocalView const& localView, LocalVector const& localVector, Assign assign)
- {
- scatter(localView, localView.tree(), localVector, assign);
- }
- /// Copies a block of values from \p vector to this
- /**
- * The values to copy might be masked with the range \p mask.
- * The assignment mode \p assign must be one of \ref Assigner::assign, or
- * \ref Assigner::plus_Assign.
- *
- * Requires the (local) size of the vector to be equal to the basis dimension.
- *
- * [[not collective]]
- **/
- template <class Vector, class MaskRange, class Assign>
- void copy(Vector const& vector, MaskRange const& mask, Assign assign)
+ /// Apply the functor \p f to all entries of the list of multi-indices \p localInd.
+ template <class IndexRange, class Func>
+ void forEach(IndexRange const& localInd, Func&& f) const
{
- test_exit(state_ == VectorStateOf_t<Assign>::value ||
- state_ == VectorState::unknown ||
- state_ == VectorState::synchronized,
- "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
-
- using size_type = typename Backend::size_type;
- size_type size = sizeInfo(*basis_);
- backend_.scatter(Dune::range(size), vector, mask, assign);
-
- // set the state to insert_values or add_values
- state_ = VectorStateOf_t<Assign>::value;
+ for (auto const& idx : localInd)
+ f(idx, asDerived().at(idx));
}
- /// Call \ref copy with `MaskRange` set to \ref FakeContainer.
- template <class Vector, class Assign>
- void copy(Vector const& vector, Assign assign)
+ /// Apply the functor \p f to all entries of the list of multi-indices \p localInd.
+ template <class IndexRange, class Func>
+ void forEach(IndexRange const& localInd, Func&& f)
{
- copy(vector, FakeContainer<bool,true>{}, assign);
+ for (auto const& idx : localInd)
+ f(idx, asDerived().at(idx));
}
- /// Apply \p func to each value at given indices \p localInd
- /**
- * First, synchronizes the values of the vector, then applies the functor to
- * each local value associated to the local indices \p localInd.
- *
- * [[mutable]]
- **/
- template <class LocalInd, class Func>
- void forEach(LocalInd const& localInd, Func&& func)
- {
- synchronize();
- backend_.forEach(localInd, FWD(func));
- }
- /// Call \ref forEach for all entries in the vector.
+ /// Apply the functor \p f to all entries of the container
template <class Func>
- void forEach(Func&& func)
- {
- using size_type = typename Backend::size_type;
- size_type size = sizeInfo(*basis_);
- forEach(Dune::range(size), FWD(func));
- }
-
- /// Apply \p func to each value at given indices \p localInd
- /**
- * First, synchronizes the values of the vector, then applies the functor to
- * each local value associated to the local indices \p localInd.
- *
- * [[const]]
- **/
- template <class LocalInd, class Func>
- void forEach(LocalInd const& localInd, Func&& func) const
+ void forEach(Func&& f)
{
- const_cast<Self*>(this)->synchronize();
- backend_.forEach(localInd, FWD(func));
+ forEach(mappedRangeView(Dune::range(asDerived().size()),
+ [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
}
- /// Call \ref forEach for all entries in the vector.
+ /// Apply the functor \p f to all entries of the container
template <class Func>
- void forEach(Func&& func) const
- {
- using size_type = typename Backend::size_type;
- size_type size = sizeInfo(*basis_);
- forEach(Dune::range(size), FWD(func));
- }
-
-
- private:
- // synchronizes ghost values. Does not touch owned values
- void synchronize()
- {
- if (state_ != VectorState::synchronized)
- backend_.synchronize();
-
- state_ = VectorState::synchronized;
- }
-
- // print the vector state to string for debugging purposes
- static std::string to_string(VectorState state)
- {
- switch (state) {
- case VectorState::synchronized: return "synchronized";
- case VectorState::insert_values: return "insert_values";
- case VectorState::add_values: return "add_values";
- default: return "unknown";
- }
- }
-
- private:
- /// The finite element space / basis associated with the data vector
- std::shared_ptr<GlobalBasis> basis_;
-
- /// Data backend
- Backend backend_;
-
- /// The current state of the vector, one of {synchronized, insert_values,
- /// add_values, unknown}
- VectorState state_ = VectorState::unknown;
- };
-
-
- template <class GB, class B>
- template <class Index>
- struct VectorBase<GB,B>::AccessProxy
- {
- using value_type = typename Backend::value_type;
-
- void operator=(value_type const& value)
- {
- self->set(i, value);
- }
-
- void operator+=(value_type const& value)
- {
- self->add(i, value);
- }
-
- void operator-=(value_type const& value)
- {
- self->add(i, -value);
- }
-
- operator value_type() const
+ void forEach(Func&& f) const
{
- return self->at(i);
+ forEach(mappedRangeView(Dune::range(asDerived().size()),
+ [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
}
-
- VectorBase* self;
- Index i;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/VectorFacade.hpp b/src/amdis/linearalgebra/VectorFacade.hpp
new file mode 100644
index 00000000..a60d20e5
--- /dev/null
+++ b/src/amdis/linearalgebra/VectorFacade.hpp
@@ -0,0 +1,370 @@
+#pragma once
+
+#include <string>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/classname.hh>
+#include <dune/common/deprecated.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/common/std/type_traits.hh>
+#include <dune/functions/functionspacebases/concepts.hh>
+#include <dune/functions/functionspacebases/sizeinfo.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/FakeContainer.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/functions/NodeIndices.hpp>
+#include <amdis/functions/SizeInfo.hpp>
+#include <amdis/operations/Assigner.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
+
+namespace AMDiS
+{
+ /// \brief The basic container that stores a base vector and a corresponding basis
+ /**
+ * A vector storing all the assembled Operators indexed with DOF indices. The
+ * vector data is associated to a global basis.
+ *
+ * \tparam T The coefficient type of the vector
+ * \tparam VectorImpl A linear-algebra backend for the vector storage
+ **/
+ template <class T, template <class> class VectorImpl>
+ class VectorFacade
+ {
+ using Self = VectorFacade;
+ using Impl = VectorImpl<T>;
+
+ private:
+ // States the vector can be in. Is changed on various insertion or gathering methods.
+ enum class VectorState
+ {
+ unknown = 0,
+ synchronized = 1,
+ insert_values = 2,
+ add_values = 3
+ };
+
+ template <class A>
+ using VectorStateOf_t = std::conditional_t<std::is_same<A,Assigner::plus_assign>::value,
+ std::integral_constant<VectorState, VectorState::add_values>,
+ std::integral_constant<VectorState, VectorState::insert_values>>;
+
+ public:
+ using size_type = typename Impl::size_type;
+ using value_type = typename Impl::value_type;
+
+ /// Constructor. Forwards the basis to the implementation class and
+ /// constructs a (type-erased) size-info object.
+ template <class GlobalBasis,
+ class = void_t<decltype(std::declval<GlobalBasis>().dimension())> >
+ VectorFacade(GlobalBasis const& basis)
+ : impl_(basis)
+ , sizeInfo_(sizeInfo(basis))
+ {
+ resizeZero();
+ }
+
+ /// Return the underlying linear algebra backend
+ Impl const& impl() const { return impl_; }
+ Impl& impl() { return impl_; }
+
+
+ template <class V>
+ using HasLocalSize = decltype(std::declval<V>().localSize());
+
+ template <class V>
+ using HasGlobalSize = decltype(std::declval<V>().globalSize());
+
+ /// Return the number of entries in the local part of the vector
+ std::size_t localSize() const
+ {
+ return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasLocalSize,Impl>{},
+ [&](auto id) { return id(impl_).localSize(); },
+ [&](auto id) { return id(impl_).size(); });
+ }
+
+ /// Return the number of entries in the global vector
+ std::size_t globalSize() const
+ {
+ return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalSize,Impl>{},
+ [&](auto id) { return id(impl_).globalSize(); },
+ [&](auto id) { return id(impl_).size(); });
+ }
+
+ /// Resize the \ref vector to the size of the \ref basis
+ void resize()
+ {
+ init(false);
+ }
+
+ /// Resize the \ref vector to the size of the \ref basis and set to zero
+ void resizeZero()
+ {
+ init(true);
+ }
+
+ /// Prepare the vector for insertion of values, finish the insertion with
+ /// \ref finish().
+ void init(bool clear)
+ {
+ impl_.init(sizeInfo_, clear);
+ state_ = clear ? VectorState::synchronized : VectorState::unknown;
+ }
+
+ /// Finish the insertion of values, see \ref init()
+ void finish()
+ {
+ impl_.finish();
+ state_ = VectorState::unknown;
+ }
+
+ /// Return the value of the vector at the local index \ref idx
+ template <class Index,
+ REQUIRES(Concepts::MultiIndex<Index>)>
+ typename Impl::value_type at(Index const& idx) const
+ {
+ const_cast<Self*>(this)->synchronize();
+ return impl_.at(idx);
+ }
+
+ /// \brief Insert a single value into the matrix (add or overwrite to existing value)
+ /**
+ * Inserts or adds a value into a certain location \p dof (given as dof multi-index)
+ * of a vector. The insertion mode is determined by the \p assign functor. Use
+ * \ref Assigner::plus_assign for adding values (default) or \ref Assigner::assign
+ * for overwriting (setting) values. Different insertion modes can not be mixed!
+ *
+ * Insertion must be closed with a call to \ref finish().
+ *
+ * [[not collective]]
+ */
+ template <class Index, class Assign = Assigner::plus_assign,
+ REQUIRES(Concepts::MultiIndex<Index>)>
+ void insert(Index const& idx, typename Impl::value_type const& value, Assign assign = {})
+ {
+ test_exit_dbg(state_ == VectorStateOf_t<Assign>::value ||
+ state_ == VectorState::unknown ||
+ state_ == VectorState::synchronized,
+ "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
+
+ impl_.insert(idx, value, assign);
+
+ // set the state to insert_values or add_values
+ state_ = VectorStateOf_t<Assign>::value;
+ }
+
+ /// See \ref insert for assignment operation \ref Assigner::assign
+ template <class Index,
+ REQUIRES(Concepts::MultiIndex<Index>)>
+ void set(Index const& idx, typename Impl::value_type const& value)
+ {
+ insert(idx, value, Assigner::assign{});
+ }
+
+ /// See \ref insert for assignment operation \ref Assigner::plus_assign
+ template <class Index,
+ REQUIRES(Concepts::MultiIndex<Index>)>
+ void add(Index const& idx, typename Impl::value_type const& value)
+ {
+ insert(idx, value, Assigner::plus_assign{});
+ }
+
+
+ /// \brief Extract values from the vector referring to the given local indices
+ /// and store it into a buffer
+ /**
+ * Collect value of indices and store them into a buffer. The buffer must be
+ * a vector-like container with `buffer.resize()` and `buffer.begin()`. The
+ * indices must be stored in an iterable container.
+ *
+ * If the vector is not in synchronized state, collects all necessary values
+ * possibly from neighbouring processors.
+ *
+ * [[expects: localView is bound to an element]]
+ * [[expects: node is in localView.tree()]]
+ * [[possibly neighbor-wise collective operation]]
+ */
+ template <class LocalView, class Node, class Buffer,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>()),
+ REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
+ void gather(LocalView const& localView, Node const& node, Buffer& buffer) const
+ {
+ test_exit(state_ == VectorState::unknown ||
+ state_ == VectorState::synchronized,
+ "Vector in invalid state {} for gather operations.", to_string(state_));
+
+ const_cast<Self*>(this)->synchronize();
+
+ buffer.resize(node.size());
+ impl_.gather(nodeIndices(localView, node), buffer.begin());
+ }
+
+ // [[expects: localView is bound to an element]]
+ template <class LocalView, class Buffer,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>())>
+ void gather(LocalView const& localView, Buffer& buffer) const
+ {
+ gather(localView, localView.tree(), buffer);
+ }
+
+ /// Insert a block of values into the vector (add or overwrite to existing values)
+ /**
+ * Inserts or adds values into certain locations of a vector. Insertion indices
+ * are extracted from the given \p localView. The insertion mode is determined
+ * by the \p assign functor. Use \ref Assigner::plus_assign for adding values
+ * (default) or \ref Assigner::assign for overwriting (setting) values. Different
+ * insertion modes can not be mixed! The \p localVector is assumed to be a continuous
+ * memory container with a `data()` method to get a pointer to the beginning.
+ *
+ * The \p mask models a boolean range with at least a `begin()` method. Must
+ * be forward iterable for at least `localVector.size()` elements. Does not
+ * need an `end()` method. See, e.g. \ref FakeContainer.
+ *
+ * Insertion must be closed with a call to \ref finish(). It is not allowed
+ * to switch insertion mode before calling `finish()`.
+ *
+ * [[expects: localView is bound to an element]]
+ * [[expects: node is in localView.tree()]]
+ * [[not collective]]
+ */
+ template <class LocalView, class Node, class NodeVector, class MaskRange, class Assign,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>()),
+ REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
+ void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector,
+ MaskRange const& mask, Assign assign)
+ {
+ test_exit(state_ == VectorStateOf_t<Assign>::value ||
+ state_ == VectorState::unknown ||
+ state_ == VectorState::synchronized,
+ "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
+
+ assert(localVector.size() == node.size());
+
+ // create a range of DOF indices on the node
+ impl_.scatter(nodeIndices(localView, node), localVector, mask, assign);
+
+ // set the state to insert_values or add_values
+ state_ = VectorStateOf_t<Assign>::value;
+ }
+
+ /// Call \ref scatter with `MaskRange` set to \ref FakeContainer.
+ // [[expects: localView is bound to an element]]
+ // [[expects: node is in localView.tree()]]
+ template <class LocalView, class Node, class NodeVector, class Assign,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>()),
+ REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
+ void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector, Assign assign)
+ {
+ scatter(localView, node, localVector, FakeContainer<bool,true>{}, assign);
+ }
+
+ /// Call \ref scatter with `Node` given by the tree of the \p localView.
+ // [[expects: localView is bound to an element]]
+ template <class LocalView, class LocalVector, class Assign,
+ REQUIRES(Dune::models<Dune::Functions::Concept::LocalView<typename LocalView::GlobalBasis>, LocalView>())>
+ void scatter(LocalView const& localView, LocalVector const& localVector, Assign assign)
+ {
+ scatter(localView, localView.tree(), localVector, assign);
+ }
+
+ /// Apply \p func to each value at given indices \p localInd
+ /**
+ * First, synchronizes the values of the vector, then applies the functor to
+ * each local value associated to the local indices \p localInd.
+ *
+ * [[mutable]]
+ **/
+ template <class LocalInd, class Func>
+ void forEach(LocalInd const& localInd, Func&& func)
+ {
+ synchronize();
+ impl_.forEach(localInd, FWD(func));
+ }
+
+ /// Call \ref forEach for all entries in the vector.
+ template <class Func>
+ void forEach(Func&& func)
+ {
+ synchronize();
+ impl_.forEach(FWD(func));
+ }
+
+ /// Apply \p func to each value at given indices \p localInd
+ /**
+ * First, synchronizes the values of the vector, then applies the functor to
+ * each local value associated to the local indices \p localInd.
+ *
+ * [[const]]
+ **/
+ template <class LocalInd, class Func>
+ void forEach(LocalInd const& localInd, Func&& func) const
+ {
+ const_cast<Self*>(this)->synchronize();
+ impl_.forEach(localInd, FWD(func));
+ }
+
+ /// Call \ref forEach for all entries in the vector.
+ template <class Func>
+ void forEach(Func&& func) const
+ {
+ const_cast<Self*>(this)->synchronize();
+ impl_.forEach(FWD(func));
+ }
+
+
+ private:
+ // synchronizes ghost values. Does not touch owned values
+ void synchronize()
+ {
+ if (state_ != VectorState::synchronized)
+ impl_.synchronize();
+
+ state_ = VectorState::synchronized;
+ }
+
+ // print the vector state to string for debugging purposes
+ static std::string to_string(VectorState state)
+ {
+ switch (state) {
+ case VectorState::synchronized: return "synchronized";
+ case VectorState::insert_values: return "insert_values";
+ case VectorState::add_values: return "add_values";
+ default: return "unknown";
+ }
+ }
+
+ private:
+ /// Data backend
+ Impl impl_;
+
+ /// Functor returning the basis dimension
+ SizeInfo sizeInfo_;
+
+ /// The current state of the vector, one of {synchronized, insert_values,
+ /// add_values, unknown}
+ VectorState state_ = VectorState::unknown;
+ };
+
+
+ namespace Impl
+ {
+ template <class T, class Facade>
+ struct VectorTypeImpl;
+
+ template <class T, class S, template <class> class Impl>
+ struct VectorTypeImpl<T,VectorFacade<S,Impl>>
+ {
+ using type = VectorFacade<T,Impl>;
+ };
+
+ } // end namespace Impl
+
+ /// Type transformation changing the value type of the vector
+ template <class T, class Facade>
+ using VectorType_t = typename Impl::VectorTypeImpl<T,Facade>::type;
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/eigen/CMakeLists.txt b/src/amdis/linearalgebra/eigen/CMakeLists.txt
index a7c200cb..908d643b 100644
--- a/src/amdis/linearalgebra/eigen/CMakeLists.txt
+++ b/src/amdis/linearalgebra/eigen/CMakeLists.txt
@@ -3,6 +3,7 @@ install(FILES
DirectRunner.hpp
IterativeRunner.hpp
MatrixBackend.hpp
+ MatrixSize.hpp
PreconConfig.hpp
SolverConfig.hpp
SolverCreator.hpp
diff --git a/src/amdis/linearalgebra/eigen/Constraints.hpp b/src/amdis/linearalgebra/eigen/Constraints.hpp
index f5a6842e..24f7e227 100644
--- a/src/amdis/linearalgebra/eigen/Constraints.hpp
+++ b/src/amdis/linearalgebra/eigen/Constraints.hpp
@@ -11,11 +11,11 @@
namespace AMDiS
{
- template <class Traits>
- struct Constraints<MatrixBackend<Traits>>
+ template <class T, int O>
+ struct Constraints<EigenSparseMatrix<T,O>>
{
- using Matrix = MatrixBackend<Traits>;
- using Vector = VectorBackend<Traits>;
+ using Matrix = EigenSparseMatrix<T,O>;
+ using Vector = EigenVector<T>;
template <class BitVector>
static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
@@ -37,7 +37,7 @@ namespace AMDiS
}
protected:
- template <class T, class BitVector>
+ template <class BitVector>
static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::ColMajor>& mat, BitVector const& nodes, bool setDiagonal)
{
using Mat = Eigen::SparseMatrix<T, Eigen::ColMajor>;
@@ -51,7 +51,7 @@ namespace AMDiS
}
}
- template <class T, class BitVector>
+ template <class BitVector>
static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::RowMajor>& mat, BitVector const& nodes, bool setDiagonal)
{
using Mat = Eigen::SparseMatrix<T, Eigen::RowMajor>;
diff --git a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
index ed184958..b4e45262 100644
--- a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
@@ -5,6 +5,8 @@
#include <string>
#include <vector>
+#include <Eigen/SparseCore>
+
#include <dune/common/timer.hh>
#include <amdis/Output.hpp>
@@ -14,14 +16,12 @@ namespace AMDiS
{
/// \brief The basic container that stores a base matrix and a corresponding
/// row/column feSpace.
- template <class T>
- class MatrixBackend
+ template <class T, int Orientation = Eigen::RowMajor>
+ class EigenSparseMatrix
{
public:
- using Traits = T;
-
/// The matrix type of the underlying base matrix
- using BaseMatrix = typename Traits::Mat;
+ using BaseMatrix = Eigen::SparseMatrix<T, Orientation>;
/// The type of the elements of the DOFMatrix
using value_type = typename BaseMatrix::Scalar;
@@ -32,7 +32,7 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseMatrix.
template <class Basis>
- explicit MatrixBackend(Basis const&, Basis const&) {}
+ explicit EigenSparseMatrix(Basis const&, Basis const&) {}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix()
@@ -74,10 +74,10 @@ namespace AMDiS
/// Create inserter. Assumes that no inserter is currently active on this matrix.
- template <class RowBasis, class ColBasis>
- void init(RowBasis const& rowBasis, ColBasis const& colBasis, SymmetryStructure symmetry)
+ template <class Pattern>
+ void init(Pattern const& pattern, SymmetryStructure symmetry)
{
- matrix_.resize(rowBasis.dimension(), colBasis.dimension());
+ matrix_.resize(pattern.rows(), pattern.cols());
matrix_.setZero();
}
diff --git a/src/amdis/linearalgebra/eigen/MatrixSize.hpp b/src/amdis/linearalgebra/eigen/MatrixSize.hpp
new file mode 100644
index 00000000..18364e93
--- /dev/null
+++ b/src/amdis/linearalgebra/eigen/MatrixSize.hpp
@@ -0,0 +1,60 @@
+#pragma once
+
+#include <functional>
+#include <memory>
+
+#include <amdis/Observer.hpp>
+#include <amdis/common/Index.hpp>
+
+namespace AMDiS
+{
+ class MatrixSize
+ : private ObserverSequence<event::adapt,2>
+ {
+ struct BasisCallback
+ {
+ template <class Basis>
+ explicit BasisCallback(Basis const* basis)
+ : dimension([basis](){ return basis->dimension(); })
+ {}
+
+ std::function<std::size_t()> dimension;
+ };
+
+ public:
+ template <class RowBasis, class ColBasis>
+ MatrixSize(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
+ , rowBasis_(&rowBasis)
+ , colBasis_(&colBasis)
+ {
+ updateImpl(event::adapt{true}, index_t<0>{});
+ updateImpl(event::adapt{true}, index_t<1>{});
+ }
+
+ /// Number of rows in the matrix
+ std::size_t rows() const
+ {
+ return rows_;
+ }
+
+ /// Number of columns in the matrix
+ std::size_t cols() const
+ {
+ return cols_;
+ }
+
+ protected:
+
+ void updateImpl(event::adapt e, index_t<0> i) final { rows_ = rowBasis_.dimension(); }
+ void updateImpl(event::adapt e, index_t<1> i) final { cols_ = colBasis_.dimension(); }
+
+ private:
+ BasisCallback rowBasis_;
+ BasisCallback colBasis_;
+
+ std::size_t rows_;
+ std::size_t cols_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/eigen/Traits.hpp b/src/amdis/linearalgebra/eigen/Traits.hpp
index f655af31..6558e1df 100644
--- a/src/amdis/linearalgebra/eigen/Traits.hpp
+++ b/src/amdis/linearalgebra/eigen/Traits.hpp
@@ -1,9 +1,10 @@
#pragma once
-#include <Eigen/SparseCore>
-
#include <dune/grid/common/partitionset.hh>
#include <amdis/linearalgebra/Communication.hpp>
+#include <amdis/linearalgebra/eigen/MatrixSize.hpp>
+#include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
+#include <amdis/linearalgebra/eigen/VectorBackend.hpp>
namespace AMDiS
{
@@ -13,11 +14,19 @@ namespace AMDiS
template <class Basis, class T = double>
struct BackendTraits
{
- using Mat = Eigen::SparseMatrix<T, Eigen::RowMajor>;
- using Vec = Eigen::Matrix<T, Eigen::Dynamic, 1>;
- using CoefficientType = T;
+ template <class Value>
+ using MatrixImpl = EigenSparseMatrix<Value, Eigen::RowMajor>;
+
+ template <class Value>
+ using VectorImpl = EigenVector<Value>;
+
+ using Mat = typename MatrixImpl<T>::BaseMatrix;
+ using Vec = typename VectorImpl<T>::BaseVector;
using Comm = SequentialCommunication;
+
+ using CoefficientType = T;
using PartitionSet = Dune::Partitions::All;
+ using SparsityPattern = MatrixSize;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/eigen/VectorBackend.hpp b/src/amdis/linearalgebra/eigen/VectorBackend.hpp
index 5efb9506..fe5ed9d7 100644
--- a/src/amdis/linearalgebra/eigen/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/eigen/VectorBackend.hpp
@@ -6,35 +6,37 @@
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
+#include <amdis/linearalgebra/VectorBase.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
namespace AMDiS
{
/// The basic container that stores a base vector and a corresponding basis
template <class T>
- class VectorBackend
+ class EigenVector
+ : public VectorBase<EigenVector<T>>
{
public:
- using Traits = T;
-
/// The type of the base vector
- using BaseVector = typename Traits::Vec;
+ using BaseVector = Eigen::Matrix<T, Eigen::Dynamic, 1>;
/// The type of the elements of the DOFVector
using value_type = typename BaseVector::Scalar;
+ /// The index/size - type
+ using size_type = typename BaseVector::Index;
+
+ private:
/// The type of the elements of the DOFVector
using block_type = value_type;
/// The underlying field type
using field_type = typename Dune::FieldTraits<value_type>::field_type;
- /// The index/size - type
- using size_type = typename BaseVector::Index;
-
public:
/// Constructor. Constructs new BaseVector.
template <class Basis>
- explicit VectorBackend(Basis const&) {}
+ explicit EigenVector(Basis const&) {}
/// Return the data-vector \ref vector_
BaseVector const& vector() const
@@ -49,7 +51,7 @@ namespace AMDiS
}
/// Return the current size of the \ref vector_
- std::size_t size() const
+ size_type size() const
{
return vector_.size();
}
@@ -64,66 +66,23 @@ namespace AMDiS
vector_.setZero();
}
- void finish() { /* do nothing */ }
- void synchronize() { /* do nothing */ }
-
-
- /// Access the entry \p i of the \ref vector with read-access.
- value_type const& at(size_type i) const
- {
- test_exit_dbg(i < vector_.size(),
- "Index {} out of range [0,{})", i, vector_.size());
- return vector_.coeff(i);
- }
/// Access the entry \p i of the \ref vector with write-access.
- template <class Assign>
- void insert(size_type i, value_type value, Assign assign)
+ template <class MultiIndex>
+ value_type& at(MultiIndex const& idx)
{
- test_exit_dbg(i < vector_.size(),
- "Index {} out of range [0,{})", i, vector_.size());
- assign(value, vector_.coeffRef(i));
+ const size_type i = flatMultiIndex(idx);
+ test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
+ return vector_.coeffRef(i);
}
- template <class IndexRange, class OutputIterator>
- void gather(IndexRange const& localInd, OutputIterator buffer) const
- {
- for (size_type i : localInd)
- *(buffer++) = vector_.coeff(i);
- }
-
- template <class IndexRange, class LocalVec, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, FakeContainer<bool,true>, Assign assign)
- {
- auto vec_it = std::begin(vec);
- for (size_type i : localInd)
- assign(*vec_it++, vector_.coeffRef(i));
- }
-
- template <class IndexRange, class LocalVec, class MaskRange, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, MaskRange const& mask, Assign assign)
- {
- auto vec_it = std::begin(vec);
- auto mask_it = std::begin(mask);
- auto ind_it = std::begin(localInd);
- for (; vec_it != std::end(vec); ++vec_it, ++mask_it, ++ind_it) {
- if (*mask_it)
- assign(*vec_it, vector_.coeffRef(*ind_it));
- }
- }
-
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f) const
- {
- for (size_type i : localInd)
- f(i, vector_.coeff(i));
- }
-
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f)
+ /// Access the entry \p i of the \ref vector with read-access.
+ template <class MultiIndex>
+ value_type const& at(MultiIndex const& idx) const
{
- for (size_type i : localInd)
- f(i, vector_.coeffRef(i));
+ const size_type i = flatMultiIndex(idx);
+ test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
+ return vector_.coeff(i);
}
private:
diff --git a/src/amdis/linearalgebra/istl/Constraints.hpp b/src/amdis/linearalgebra/istl/Constraints.hpp
index bbfed667..873a2a71 100644
--- a/src/amdis/linearalgebra/istl/Constraints.hpp
+++ b/src/amdis/linearalgebra/istl/Constraints.hpp
@@ -7,17 +7,15 @@
namespace AMDiS
{
- template <class Traits>
- struct Constraints<MatrixBackend<Traits>>
+ template <class T>
+ struct Constraints<ISTLBCRSMatrix<T>>
{
- using Matrix = MatrixBackend<Traits>;
- using Vector = VectorBackend<Traits>;
+ using Matrix = ISTLBCRSMatrix<T>;
+ using Vector = ISTLBlockVector<T>;
template <class BitVector>
static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
{
- using T = typename Matrix::value_type;
-
// loop over the matrix rows
for (std::size_t i = 0; i < mat.matrix().N(); ++i) {
if (nodes[i]) {
diff --git a/src/amdis/linearalgebra/istl/MatrixBackend.hpp b/src/amdis/linearalgebra/istl/MatrixBackend.hpp
index a62e1135..1ae83e01 100644
--- a/src/amdis/linearalgebra/istl/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/istl/MatrixBackend.hpp
@@ -4,6 +4,7 @@
#include <string>
#include <memory>
+#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <amdis/Output.hpp>
@@ -11,14 +12,24 @@
namespace AMDiS
{
+ template <class T, class = void>
+ struct BlockMatrixType
+ {
+ using type = Dune::FieldMatrix<T,1,1>;
+ };
+
template <class T>
- class MatrixBackend
+ struct BlockMatrixType<T, typename T::field_type>
{
- public:
- using Traits = T;
+ using type = T;
+ };
+ template <class T>
+ class ISTLBCRSMatrix
+ {
+ public:
/// The matrix type of the underlying base matrix
- using BaseMatrix = typename Traits::Mat;
+ using BaseMatrix = Dune::BCRSMatrix<typename BlockMatrixType<T>::type>;
/// The type of the elements of the DOFMatrix
using value_type = typename BaseMatrix::block_type;
@@ -29,7 +40,7 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
template <class Basis>
- explicit MatrixBackend(Basis const&, Basis const&) {}
+ explicit ISTLBCRSMatrix(Basis const&, Basis const&) {}
/// Return the data-vector \ref vector
BaseMatrix const& matrix() const
@@ -44,30 +55,10 @@ namespace AMDiS
}
/// create occupation pattern and apply it to the matrix
- template <class RowBasis, class ColBasis>
- void init(RowBasis const& rowBasis, ColBasis const& colBasis, SymmetryStructure symmetry)
+ template <class Pattern>
+ void init(Pattern const& pattern, SymmetryStructure symmetry)
{
- auto occupationPattern = Dune::MatrixIndexSet{rowBasis.dimension(), colBasis.dimension()};
-
- auto rowLocalView = rowBasis.localView();
- auto colLocalView = colBasis.localView();
- for (const auto& element : elements(rowBasis.gridView())) {
- rowLocalView.bind(element);
- colLocalView.bind(element);
-
- for (std::size_t i = 0; i < rowLocalView.size(); ++i) {
- // The global index of the i-th vertex of the element
- auto row = rowLocalView.index(i);
- 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.exportIdx(matrix_);
-
+ pattern.applyTo(matrix_);
symmetry_ = symmetry;
initialized_ = true;
}
diff --git a/src/amdis/linearalgebra/istl/Traits.hpp b/src/amdis/linearalgebra/istl/Traits.hpp
index fa368680..4372cb5f 100644
--- a/src/amdis/linearalgebra/istl/Traits.hpp
+++ b/src/amdis/linearalgebra/istl/Traits.hpp
@@ -3,54 +3,38 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/partitionset.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioner.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/solver.hh>
+#include <amdis/linearalgebra/SparsityPattern.hpp>
#include <amdis/linearalgebra/istl/Communication.hpp>
#include <amdis/linearalgebra/istl/Creators.hpp>
+#include <amdis/linearalgebra/istl/MatrixBackend.hpp>
+#include <amdis/linearalgebra/istl/VectorBackend.hpp>
namespace AMDiS
{
- template <class T, class = void>
- struct BlockMatrixType
- {
- using type = Dune::FieldMatrix<T,1,1>;
- };
-
- template <class T>
- struct BlockMatrixType<T, typename T::field_type>
- {
- using type = T;
- };
-
- template <class T, class = void>
- struct BlockVectorType
- {
- using type = Dune::FieldVector<T,1>;
- };
-
- template <class T>
- struct BlockVectorType<T, typename T::field_type>
- {
- using type = T;
- };
-
-
/** Traits class for a linear solver for the system AX=B using an FE space described by a dune-functions Basis
* Contains typedefs specific to the ISTL backend.
*/
template <class Basis, class T = double>
struct BackendTraits
{
- using Mat = Dune::BCRSMatrix<typename BlockMatrixType<T>::type>;
- using Vec = Dune::BlockVector<typename BlockVectorType<T>::type>;
- using CoefficientType = T;
+ template <class Value>
+ using MatrixImpl = ISTLBCRSMatrix<Value>;
+
+ template <class Value>
+ using VectorImpl = ISTLBlockVector<Value>;
+
+ using Mat = typename MatrixImpl<T>::BaseMatrix;
+ using Vec = typename VectorImpl<T>::BaseVector;
using Comm = ISTLCommunication_t<Basis>;
+
+ using CoefficientType = T;
using PartitionSet = Dune::Partitions::All;
+ using SparsityPattern = AMDiS::SparsityPattern<Basis, Basis>;
using ScalProd = Dune::ScalarProduct<Vec>;
using LinOp = Dune::AssembledLinearOperator<Mat, Vec, Vec>;
diff --git a/src/amdis/linearalgebra/istl/VectorBackend.hpp b/src/amdis/linearalgebra/istl/VectorBackend.hpp
index 378248a6..ab405cf4 100644
--- a/src/amdis/linearalgebra/istl/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/istl/VectorBackend.hpp
@@ -1,24 +1,40 @@
#pragma once
+#include <dune/istl/bvector.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
+#include <amdis/linearalgebra/VectorBase.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
namespace AMDiS
{
+ template <class T, class = void>
+ struct BlockVectorType
+ {
+ using type = Dune::FieldVector<T,1>;
+ };
+
+ template <class T>
+ struct BlockVectorType<T, typename T::field_type>
+ {
+ using type = T;
+ };
+
template <class T>
- class VectorBackend
+ class ISTLBlockVector
+ : public VectorBase<ISTLBlockVector<T>>
{
public:
- using Traits = T;
-
/// The vector type of the underlying base vector
- using BaseVector = typename Traits::Vec;
+ using BaseVector = Dune::BlockVector<typename BlockVectorType<T>::type>;
/// The type of the elements of the DOFVector
using block_type = typename BaseVector::block_type;
/// The type of the elements of the DOFVector
- using value_type = block_type;
+ using value_type = T;
/// The underlying field type
using field_type = typename block_type::field_type;
@@ -29,7 +45,7 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
template <class Basis>
- explicit VectorBackend(Basis const&) {}
+ explicit ISTLBlockVector(Basis const&) {}
/// Return the data-vector \ref vector
BaseVector const& vector() const
@@ -49,74 +65,28 @@ namespace AMDiS
return vector_.size();
}
- /// Resize the \ref vector_ to the size \p s
+ /// Resize the \ref vector_ to the size given by the \p sizeProvider
+ // NOTE: this resize works recursively on the hierarchy of the vector
template <class SizeInfo>
- void init(SizeInfo const& size, bool clear)
+ void init(SizeInfo const& sizeInfo, bool clear)
{
- vector_.resize(size_type(size));
+ Dune::Functions::istlVectorBackend(vector_).resize(sizeInfo);
if (clear)
vector_ = 0;
}
- void finish() { /* do nothing */ }
- void synchronize() { /* so nothing */ }
-
- /// Access the entry \p i of the \ref vector with read-access.
- block_type const& at(size_type i) const
- {
- test_exit_dbg(i < vector_.size(),
- "Index {} out of range [0,{})", i, vector_.size());
- return vector_[i];
- }
-
- /// Access the entry \p i of the \ref vector with write-access.
- template <class Assign>
- void insert(size_type i, block_type const& value, Assign assign)
- {
- test_exit_dbg(i < vector_.size(),
- "Index {} out of range [0,{})", i, vector_.size());
- assign(value, vector_[i]);
- }
-
- template <class IndexRange, class OutputIterator>
- void gather(IndexRange const& localInd, OutputIterator buffer) const
- {
- for (size_type i : localInd)
- *buffer++ = vector_[i];
- }
-
- template <class IndexRange, class LocalVec, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, FakeContainer<bool,true>, Assign assign)
- {
- auto vec_it = std::begin(vec);
- for (size_type i : localInd)
- assign(*vec_it++, vector_[i]);
- }
-
- template <class IndexRange, class LocalVec, class MaskRange, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, MaskRange const& mask, Assign assign)
- {
- auto vec_it = std::begin(vec);
- auto mask_it = std::begin(mask);
- auto ind_it = std::begin(localInd);
- for (; vec_it != std::end(vec); ++vec_it, ++mask_it, ++ind_it) {
- if (*mask_it)
- assign(*vec_it, vector_[*ind_it]);
- }
- }
-
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f) const
+ /// Access the entry \p idx of the \ref vector_ with write-access.
+ template <class MultiIndex>
+ value_type& at(MultiIndex const& idx)
{
- for (size_type i : localInd)
- f(i, vector_[i]);
+ return Dune::Functions::istlVectorBackend(vector_)[idx];
}
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f)
+ /// Access the entry \p idx of the \ref vector_ with read-access.
+ template <class MultiIndex>
+ value_type const& at(MultiIndex const& idx) const
{
- for (size_type i : localInd)
- f(i, vector_[i]);
+ return Dune::Functions::istlVectorBackend(vector_)[idx];
}
private:
diff --git a/src/amdis/linearalgebra/mtl/CMakeLists.txt b/src/amdis/linearalgebra/mtl/CMakeLists.txt
index 8b79796a..39165a47 100644
--- a/src/amdis/linearalgebra/mtl/CMakeLists.txt
+++ b/src/amdis/linearalgebra/mtl/CMakeLists.txt
@@ -6,6 +6,7 @@ install(FILES
KrylovRunner.hpp
MatrixBackend.hpp
Preconditioner.hpp
+ SlotSize.hpp
SolverPrecon.hpp
Traits.hpp
UmfpackRunner.hpp
diff --git a/src/amdis/linearalgebra/mtl/Constraints.hpp b/src/amdis/linearalgebra/mtl/Constraints.hpp
index 962bc51d..f15487fa 100644
--- a/src/amdis/linearalgebra/mtl/Constraints.hpp
+++ b/src/amdis/linearalgebra/mtl/Constraints.hpp
@@ -15,11 +15,11 @@
namespace AMDiS
{
- template <class Traits>
- struct Constraints<MatrixBackend<Traits>>
+ template <class T>
+ struct Constraints<MTLSparseMatrix<T>>
{
- using Matrix = MatrixBackend<Traits>;
- using Vector = VectorBackend<Traits>;
+ using Matrix = MTLSparseMatrix<T>;
+ using Vector = MTLVector<T>;
template <class BitVector>
static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
@@ -128,11 +128,11 @@ namespace AMDiS
}
- template <class T>
+ template <class Value>
struct Triplet
{
std::size_t row, col;
- T value;
+ Value value;
};
template <class Mat, class Vec, class BitVector, class Associations>
diff --git a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
index 364bf7c9..78ce1d71 100644
--- a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
@@ -5,24 +5,24 @@
#include <string>
#include <vector>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
+#include <amdis/linearalgebra/mtl/SlotSize.hpp>
namespace AMDiS
{
/// \brief The basic container that stores a base matrix
- template <class TraitsType>
- class MatrixBackend
+ template <class T>
+ class MTLSparseMatrix
{
public:
- using Traits = TraitsType;
-
/// The matrix type of the underlying base matrix
- using BaseMatrix = typename Traits::Mat;
+ using BaseMatrix = mtl::compressed2D<T>;
/// The type of the elements of the DOFMatrix
using value_type = typename BaseMatrix::value_type;
@@ -30,13 +30,17 @@ namespace AMDiS
/// The index/size - type
using size_type = typename BaseMatrix::size_type;
+ private:
/// 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>>;
+ /// Type of the matrix non-zero pattern
+ using Pattern = SlotSize;
+
public:
/// Constructor. Constructs new BaseMatrix.
template <class Basis>
- explicit MatrixBackend(Basis const&, Basis const&) {}
+ explicit MTLSparseMatrix(Basis const& /*rowBasis*/, Basis const& /*colBasis*/) {}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix()
@@ -53,16 +57,15 @@ namespace AMDiS
}
/// Create inserter. Assumes that no inserter is currently active on this matrix.
- template <class RowBasis, class ColBasis>
- void init(RowBasis const& rowBasis, ColBasis const& colBasis, SymmetryStructure symmetry)
+ void init(Pattern const& pattern, SymmetryStructure symmetry)
{
test_exit(!inserter_, "Matrix already in insertion mode!");
- calculateSlotSize();
- matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
+ std::size_t slotSize = pattern.avgRowSize(matrix_);
+ matrix_.change_dim(pattern.rows(), pattern.cols());
set_to_zero(matrix_);
- inserter_ = new Inserter(matrix_, slotSize_);
+ inserter_ = new Inserter(matrix_, slotSize);
symmetry_ = symmetry;
}
@@ -109,36 +112,20 @@ namespace AMDiS
return matrix_.nnz();
}
+ /// Symmetry of the matrix entries
SymmetryStructure symmetry() const
{
return symmetry_;
}
- protected:
- // Estimates the slot-size used in the inserter to reserve enough space per row.
- void calculateSlotSize()
- {
- slotSize_ = 0;
-
- 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 data-matrix (can hold a new BaseMatrix or a pointer to external data
BaseMatrix matrix_;
/// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
Inserter* inserter_ = nullptr;
- /// The size of the slots used to insert values per row
- int slotSize_ = MIN_NNZ_PER_ROW;
-
+ /// Symmetry of the matrix entries
SymmetryStructure symmetry_ = SymmetryStructure::unknown;
};
diff --git a/src/amdis/linearalgebra/mtl/SlotSize.hpp b/src/amdis/linearalgebra/mtl/SlotSize.hpp
new file mode 100644
index 00000000..c04d8c31
--- /dev/null
+++ b/src/amdis/linearalgebra/mtl/SlotSize.hpp
@@ -0,0 +1,90 @@
+#pragma once
+
+#include <cassert>
+#include <functional>
+#include <memory>
+#include <dune/common/math.hh>
+
+#include <amdis/Observer.hpp>
+#include <amdis/common/Index.hpp>
+#include <amdis/common/Math.hpp>
+
+namespace AMDiS
+{
+ class SlotSize
+ : private ObserverSequence<event::adapt,2>
+ {
+ struct BasisCallback
+ {
+ template <class Basis>
+ explicit BasisCallback(Basis const* basis)
+ : dimension([basis](){ return basis->dimension(); })
+ , maxSize([basis](){ return basis->localView().maxSize(); })
+ {}
+
+ std::function<std::size_t()> dimension;
+ std::function<std::size_t()> maxSize;
+ };
+
+ public:
+ template <class RowBasis, class ColBasis>
+ SlotSize(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
+ , rowBasis_(&rowBasis)
+ , colBasis_(&colBasis)
+ , surrounding_(Math::pow<RowBasis::GridView::dimension>(2) * Dune::factorial(int(RowBasis::GridView::dimension)))
+ {
+ updateImpl(event::adapt{true}, index_t<0>{});
+ updateImpl(event::adapt{true}, index_t<1>{});
+ }
+
+ /// Number of rows in the matrix
+ std::size_t rows() const
+ {
+ return rows_;
+ }
+
+ /// Number of columns in the matrix
+ std::size_t cols() const
+ {
+ return cols_;
+ }
+
+ /// Average number of non-zeros per row
+ /// In the first time this function is called, use the estRowSize,
+ /// otherwise take `1.5 * nnz/rows`
+ template <class Matrix>
+ std::size_t avgRowSize(Matrix const& matrix) const
+ {
+ assert(rows_ > 0 && cols_ > 0);
+ return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_;
+ }
+
+ protected:
+
+ // update of the row basis, just update the dimension
+ void updateImpl(event::adapt e, index_t<0> i) final
+ {
+ rows_ = rowBasis_.dimension();
+ }
+
+ // estimate the number of columns by multiplying the maximal node size with the
+ // number of element surrounding a vertex. This number is approximated by the
+ // number of simplices surrounding a vertex in a kuhn tringulation
+ void updateImpl(event::adapt e, index_t<1> i) final
+ {
+ cols_ = colBasis_.dimension();
+ estRowSize_ = std::min(cols_, colBasis_.maxSize() * surrounding_);
+ }
+
+ private:
+ BasisCallback rowBasis_;
+ BasisCallback colBasis_;
+
+ std::size_t rows_;
+ std::size_t cols_;
+ std::size_t estRowSize_;
+ std::size_t surrounding_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/mtl/Traits.hpp b/src/amdis/linearalgebra/mtl/Traits.hpp
index f474fb63..e3fa936a 100644
--- a/src/amdis/linearalgebra/mtl/Traits.hpp
+++ b/src/amdis/linearalgebra/mtl/Traits.hpp
@@ -1,10 +1,10 @@
#pragma once
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
#include <dune/grid/common/partitionset.hh>
#include <amdis/linearalgebra/Communication.hpp>
+#include <amdis/linearalgebra/mtl/SlotSize.hpp>
+#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
+#include <amdis/linearalgebra/mtl/VectorBackend.hpp>
namespace AMDiS
{
@@ -14,11 +14,19 @@ namespace AMDiS
template <class Basis, class T = double>
struct BackendTraits
{
- using Mat = mtl::compressed2D<T>;
- using Vec = mtl::dense_vector<T>;
- using CoefficientType = T;
+ template <class Value>
+ using MatrixImpl = MTLSparseMatrix<Value>;
+
+ template <class Value>
+ using VectorImpl = MTLVector<Value>;
+
+ using Mat = typename MatrixImpl<T>::BaseMatrix;
+ using Vec = typename VectorImpl<T>::BaseVector;
using Comm = SequentialCommunication;
+
+ using CoefficientType = T;
using PartitionSet = Dune::Partitions::All;
+ using SparsityPattern = SlotSize;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/mtl/VectorBackend.hpp b/src/amdis/linearalgebra/mtl/VectorBackend.hpp
index c0108522..ada312c6 100644
--- a/src/amdis/linearalgebra/mtl/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/mtl/VectorBackend.hpp
@@ -1,38 +1,42 @@
#pragma once
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
#include <dune/common/ftraits.hh>
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
+#include <amdis/linearalgebra/VectorBase.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
namespace AMDiS
{
/// \brief The basic container that stores a base vector data
- template <class TraitsType>
- class VectorBackend
+ template <class T>
+ class MTLVector
+ : public VectorBase<MTLVector<T>>
{
public:
- using Traits = TraitsType;
-
/// The type of the base vector
- using BaseVector = typename Traits::Vec;
+ using BaseVector = mtl::dense_vector<T>;
/// The type of the elements of the DOFVector
using value_type = typename BaseVector::value_type;
+ /// The index/size - type
+ using size_type = typename BaseVector::size_type;
+
+ private:
/// The type of the elements of the DOFVector
using block_type = value_type;
/// The underlying field type
using field_type = typename Dune::FieldTraits<value_type>::field_type;
- /// The index/size - type
- using size_type = typename BaseVector::size_type;
-
public:
/// Constructor. Constructs new BaseVector.
template <class Basis>
- explicit VectorBackend(Basis const&) {}
+ explicit MTLVector(Basis const&) {}
/// Return the data-vector \ref vector_
BaseVector const& vector() const
@@ -55,71 +59,29 @@ namespace AMDiS
/// Resize the \ref vector_ to the size \p s
template <class SizeInfo>
- void init(SizeInfo const& size, bool clear)
+ void init(SizeInfo const& sizeInfo, bool clear)
{
- vector_.change_dim(size_type(size));
+ vector_.change_dim(sizeInfo({}));
if (clear)
set_to_zero(vector_);
}
- void finish() { /* do nothing */ }
- void synchronize() { /* do nothing */ }
-
/// Access the entry \p i of the \ref vector with read-access.
- value_type const& at(size_type i) const
+ template <class MultiIndex>
+ value_type& at(MultiIndex const& idx)
{
- test_exit_dbg(i < mtl::vec::size(vector_),
- "Index {} out of range [0,{})", i, mtl::vec::size(vector_));
+ const size_type i = flatMultiIndex(idx);
+ test_exit_dbg(i < size(),"Index {} out of range [0,{})", i, size());
return vector_[i];
}
- template <class Assign>
- void insert(size_type i, value_type const& value, Assign assign)
- {
- test_exit_dbg(i < mtl::vec::size(vector_),
- "Index {} out of range [0,{})", i, mtl::vec::size(vector_));
- assign(value, vector_[i]);
- }
-
- template <class IndexRange, class OutputIterator>
- void gather(IndexRange const& localInd, OutputIterator buffer) const
- {
- for (size_type i : localInd)
- *buffer++ = vector_[i];
- }
-
- template <class IndexRange, class LocalVec, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, FakeContainer<bool,true>, Assign assign)
- {
- auto vec_it = std::begin(vec);
- for (size_type i : localInd)
- assign(*vec_it++, vector_[i]);
- }
-
- template <class IndexRange, class LocalVec, class MaskRange, class Assign>
- void scatter(IndexRange const& localInd, LocalVec const& vec, MaskRange const& mask, Assign assign)
- {
- auto vec_it = std::begin(vec);
- auto mask_it = std::begin(mask);
- auto ind_it = std::begin(localInd);
- for (; vec_it != std::end(vec); ++vec_it, ++mask_it, ++ind_it) {
- if (*mask_it)
- assign(*vec_it, vector_[*ind_it]);
- }
- }
-
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f) const
- {
- for (size_type i : localInd)
- f(i, vector_[i]);
- }
-
- template <class IndexRange, class Func>
- void forEach(IndexRange const& localInd, Func&& f)
+ /// Access the entry \p i of the \ref vector with read-access.
+ template <class MultiIndex>
+ value_type const& at(MultiIndex const& idx) const
{
- for (size_type i : localInd)
- f(i, vector_[i]);
+ const size_type i = flatMultiIndex(idx);
+ test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
+ return vector_[i];
}
private:
diff --git a/src/amdis/linearalgebra/petsc/Communication.hpp b/src/amdis/linearalgebra/petsc/Communication.hpp
index 32d21214..8687cc9a 100644
--- a/src/amdis/linearalgebra/petsc/Communication.hpp
+++ b/src/amdis/linearalgebra/petsc/Communication.hpp
@@ -15,7 +15,7 @@
namespace AMDiS
{
- template <class GlobalId, class LocalIndex>
+ template <class GlobalId, class SizeType>
class DistributedCommunication
{
public:
@@ -26,8 +26,8 @@ namespace AMDiS
#else
struct IndexSet
{
- constexpr LocalIndex size() const { return size_; }
- LocalIndex size_ = 0;
+ constexpr SizeType size() const { return size_; }
+ SizeType size_ = 0;
};
struct RemoteIndices {};
@@ -57,7 +57,7 @@ namespace AMDiS
template <class Basis>
void update(Basis const& basis)
{
- indexSet_.emplace();
+ indexSet_ = std::make_unique<IndexSet>();
#if HAVE_MPI
buildParallelIndexSet(basis, *indexSet_);
@@ -71,7 +71,7 @@ namespace AMDiS
}
protected:
- Dune::Std::optional<IndexSet> indexSet_;
+ std::unique_ptr<IndexSet> indexSet_ = nullptr;
RemoteIndices remoteIndices_;
DofMap dofMap_;
};
diff --git a/src/amdis/linearalgebra/petsc/Constraints.hpp b/src/amdis/linearalgebra/petsc/Constraints.hpp
index d00dba32..8dd91646 100644
--- a/src/amdis/linearalgebra/petsc/Constraints.hpp
+++ b/src/amdis/linearalgebra/petsc/Constraints.hpp
@@ -8,11 +8,11 @@
namespace AMDiS
{
- template <class Traits>
- struct Constraints<MatrixBackend<Traits>>
+ template <class DofMap>
+ struct Constraints<PetscMatrix<DofMap>>
{
- using Matrix = MatrixBackend<Traits>;
- using Vector = VectorBackend<Traits>;
+ using Matrix = PetscMatrix<DofMap>;
+ using Vector = PetscVector<DofMap>;
/// Dirichlet boundary conditions.
/**
@@ -27,7 +27,7 @@ namespace AMDiS
{
thread_local std::vector<PetscInt> rows;
rows.clear();
- auto const& dofMap = mat.dofMap();
+ auto const& dofMap = mat.dofMap_;
for (std::size_t i = 0; i < nodes.size(); ++i)
if (nodes[i])
rows.push_back(dofMap.global(i));
@@ -53,7 +53,7 @@ namespace AMDiS
thread_local std::vector<PetscInt> rows, associated_rows;
rows.clear();
associated_rows.clear();
- auto const& dofMap = mat.dofMap();
+ auto const& dofMap = mat.dofMap_;
for (PetscInt i = 0; i < PetscInt(left.size()); ++i) {
if (left[i]) {
// get global row index
diff --git a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
index 2c70d304..e72a33fb 100644
--- a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
@@ -6,23 +6,20 @@
#include <petscmat.h>
+#include <dune/common/timer.hh>
+
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
-#include <amdis/linearalgebra/petsc/MatrixNnzStructure.hpp>
namespace AMDiS
{
/// \brief The basic container that stores a base matrix
- template <class T>
- class MatrixBackend
+ template <class DofMap>
+ class PetscMatrix
{
- friend struct Constraints<MatrixBackend<T>>;
-
- using Communication = typename T::Comm;
+ friend struct Constraints<PetscMatrix<DofMap>>;
public:
- using Traits = T;
-
/// The matrix type of the underlying base matrix
using BaseMatrix = ::Mat;
@@ -35,17 +32,17 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseMatrix.
template <class Basis>
- explicit MatrixBackend(Basis const& rowBasis, Basis const& /*colBasis*/)
- : comm_(rowBasis.communicator())
+ explicit PetscMatrix(Basis const& rowBasis, Basis const& /*colBasis*/)
+ : dofMap_(rowBasis.communicator().dofMap())
{}
// disable copy and move operations
- MatrixBackend(MatrixBackend const&) = delete;
- MatrixBackend(MatrixBackend&&) = delete;
- MatrixBackend& operator=(MatrixBackend const&) = delete;
- MatrixBackend& operator=(MatrixBackend&&) = delete;
+ PetscMatrix(PetscMatrix const&) = delete;
+ PetscMatrix(PetscMatrix&&) = delete;
+ PetscMatrix& operator=(PetscMatrix const&) = delete;
+ PetscMatrix& operator=(PetscMatrix&&) = delete;
- ~MatrixBackend()
+ ~PetscMatrix()
{
destroy();
}
@@ -66,8 +63,8 @@ namespace AMDiS
template <class RowIndex, class ColIndex>
void insert(RowIndex const& r, ColIndex const& c, PetscScalar value)
{
- PetscInt row = dofMap().global(r);
- PetscInt col = dofMap().global(c);
+ PetscInt row = dofMap_.global(r);
+ PetscInt col = dofMap_.global(c);
MatSetValue(matrix_, row, col, value, ADD_VALUES);
}
@@ -80,7 +77,7 @@ namespace AMDiS
// create a vector of global indices from the local indices using the local-to-global map
std::transform(localInd.begin(), localInd.end(), idx.begin(),
- [&](auto const& mi) { return this->dofMap().global(mi); });
+ [this](auto const& mi) { return dofMap_.global(mi); });
MatSetValues(matrix_, idx.size(), idx.data(), idx.size(), idx.data(), localMat.data(), ADD_VALUES);
}
@@ -96,24 +93,19 @@ namespace AMDiS
// create vectors of global indices from the local indices using the local-to-global map
std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
- [&](auto const& mi) { return this->dofMap().global(mi); });
+ [this](auto const& mi) { return dofMap_.global(mi); });
std::transform(colLocalInd.begin(), colLocalInd.end(), ci.begin(),
- [&](auto const& mi) { return this->dofMap().global(mi); });
+ [this](auto const& mi) { return dofMap_.global(mi); });
MatSetValues(matrix_, ri.size(), ri.data(), ci.size(), ci.data(), localMat.data(), ADD_VALUES);
}
/// Create and initialize the matrix
- template <class RowBasis, class ColBasis>
- void init(RowBasis const& rowBasis, ColBasis const& colBasis, SymmetryStructure symmetry)
+ template <class Pattern>
+ void init(Pattern const& pattern, SymmetryStructure symmetry)
{
Dune::Timer t;
- // update the non-zeros structure of the matrix.
- pattern_.update(rowBasis);
- info(3, " update pattern needed {} seconds", t.elapsed());
- t.reset();
-
// destroy an old matrix if created before
destroy();
info(3, " destroy old matrix needed {} seconds", t.elapsed());
@@ -121,14 +113,14 @@ namespace AMDiS
// create a MATAIJ or MATSEQAIJ matrix with provided sparsity pattern
MatCreateAIJ(comm(),
- dofMap().localSize(), dofMap().localSize(),
- dofMap().globalSize(), dofMap().globalSize(),
- 0, pattern_.d_nnz().data(), 0, pattern_.o_nnz().data(), &matrix_);
+ dofMap_.localSize(), dofMap_.localSize(),
+ dofMap_.globalSize(), dofMap_.globalSize(),
+ 0, pattern.d_nnz().data(), 0, pattern.o_nnz().data(), &matrix_);
// keep sparsity pattern even if we delete a row / column with e.g. MatZeroRows
MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
- // set symmetriy properties of the matrix
+ // set symmetry properties of the matrix
switch (symmetry) {
case SymmetryStructure::spd:
MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
@@ -171,12 +163,6 @@ namespace AMDiS
}
private:
- // The local-to-global map
- auto const& dofMap() const
- {
- return comm_.dofMap();
- }
-
// An MPI Communicator or PETSC_COMM_SELF
MPI_Comm comm() const
{
@@ -195,11 +181,11 @@ namespace AMDiS
}
private:
- Communication const& comm_;
+ // The local-to-global map
+ DofMap const& dofMap_;
/// The data-matrix
Mat matrix_;
- MatrixNnzStructure pattern_;
bool initialized_ = false;
};
diff --git a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
index 1a83ec3d..215d8aff 100644
--- a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
@@ -3,6 +3,9 @@
#include <algorithm>
#include <vector>
+#include <amdis/Observer.hpp>
+#include <amdis/common/Index.hpp>
+
#if HAVE_MPI
#include <amdis/common/parallel/Communicator.hpp>
#endif
@@ -10,9 +13,19 @@
namespace AMDiS
{
/// Sparsity pattern used to create PETSc matrices
+ template <class RowBasis, class ColBasis = RowBasis>
class MatrixNnzStructure
+ : private ObserverSequence<event::adapt,2>
{
public:
+ MatrixNnzStructure(RowBasis const& rowBasis, ColBasis const& colBasis)
+ : ObserverSequence<event::adapt,2>(rowBasis, colBasis)
+ , rowBasis_(&rowBasis)
+ , colBasis_(&colBasis)
+ {
+ update();
+ }
+
// Return Number of nonzeros in the diagonal part (owner part)
std::vector<PetscInt> const& d_nnz() const
{
@@ -25,12 +38,34 @@ namespace AMDiS
return onnz_;
}
- /// Calculate the matrix pattern from the local pattern of a basis.
- // NOTE: this assumes a square matrix with row-basis = col-basis.
- template <class Basis>
- void update(Basis const& basis);
+ protected:
+
+ void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
+ void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }
private:
+
+ // Track for each basis wether updated and if both are, call updateImpl3()
+ template <std::size_t I>
+ void updateImpl2(event::adapt, index_t<I>)
+ {
+ assert(!updateCounter_.test(I));
+ updateCounter_.set(I);
+
+ if (updateCounter_.all()) {
+ update();
+ updateCounter_.reset();
+ }
+ }
+
+ // Update pattern when basis is updated
+ void update();
+
+ private:
+ RowBasis const* rowBasis_;
+ ColBasis const* colBasis_;
+ std::bitset<2> updateCounter_ = 0;
+
std::vector<PetscInt> dnnz_; //< number of nonzeros in the diagonal part (owner part)
std::vector<PetscInt> onnz_; //< number of nonzeros in the off-diagonal part (overlap part)
diff --git a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
index cf7a5070..9231ca33 100644
--- a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
@@ -16,19 +16,12 @@
namespace AMDiS {
-template <class Basis>
-void MatrixNnzStructure::update(Basis const& basis)
+template <class RB, class CB>
+void MatrixNnzStructure<RB,CB>::update()
{
- unsigned long newChangeIndex = changeIndex(basis.gridView().grid());
- bool gridChanged = newChangeIndex > changeIndex_;
- changeIndex_ = newChangeIndex;
-
- // update is not necessary
- if (initialized_ && !gridChanged)
- return;
-
Dune::Timer t;
+ auto const& basis = *rowBasis_; // TODO: generalize to row != col basis
auto const& dofMap = basis.communicator().dofMap();
std::size_t localSize = dofMap.localSize();
test_exit(localSize > 0, "Local domain has no owned DOFs!");
@@ -185,8 +178,7 @@ void MatrixNnzStructure::update(Basis const& basis)
#endif // HAVE_MPI
- initialized_ = true;
- msg("update MatrixNnzStructure need {} sec", t.elapsed());
+ info(2,"update MatrixNnzStructure need {} sec", t.elapsed());
}
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/petsc/Traits.hpp b/src/amdis/linearalgebra/petsc/Traits.hpp
index 9272533b..f0c275bb 100644
--- a/src/amdis/linearalgebra/petsc/Traits.hpp
+++ b/src/amdis/linearalgebra/petsc/Traits.hpp
@@ -7,6 +7,9 @@
#include <dune/grid/common/partitionset.hh>
#include <amdis/linearalgebra/petsc/Communication.hpp>
+#include <amdis/linearalgebra/petsc/MatrixNnzStructure.hpp>
+#include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
+#include <amdis/linearalgebra/petsc/VectorBackend.hpp>
namespace AMDiS
{
@@ -20,9 +23,17 @@ namespace AMDiS
using Mat = ::Mat;
using Vec = ::Vec;
- using CoefficientType = PetscScalar;
using Comm = DistributedCommunication_t<Basis>;
+
+ template <class>
+ using MatrixImpl = PetscMatrix<typename Comm::DofMap>;
+
+ template <class>
+ using VectorImpl = PetscVector<typename Comm::DofMap>;
+
+ using CoefficientType = PetscScalar;
using PartitionSet = Dune::Partitions::Interior;
+ using SparsityPattern = MatrixNnzStructure<Basis>;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/petsc/VectorBackend.hpp b/src/amdis/linearalgebra/petsc/VectorBackend.hpp
index 3634c607..a8d67718 100644
--- a/src/amdis/linearalgebra/petsc/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/petsc/VectorBackend.hpp
@@ -4,28 +4,25 @@
#include <dune/common/ftraits.hh>
#include <dune/common/shared_ptr.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
-#include <amdis/functions/Interpolate.hpp>
#include <amdis/operations/Assigner.hpp>
+#include <amdis/typetree/MultiIndex.hpp>
namespace AMDiS
{
/// The basic container that stores a base vector data
- template <class T>
- class VectorBackend
+ template <class DofMap>
+ class PetscVector
{
- using Communication = typename T::Comm;
-
template <class A>
using AssignMode = std::conditional_t<std::is_same<A,Assigner::plus_assign>::value,
std::integral_constant<::InsertMode, ADD_VALUES>,
std::integral_constant<::InsertMode, INSERT_VALUES>>;
public:
- using Traits = T;
-
/// The type of the base vector
using BaseVector = ::Vec;
@@ -38,19 +35,20 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
template <class Basis>
- explicit VectorBackend(Basis const& basis)
- : comm_(Dune::stackobject_to_shared_ptr(basis.communicator()))
+ explicit PetscVector(Basis const& basis)
+ : dofMap_(basis.communicator().dofMap())
{}
/// Move constructor
- VectorBackend(VectorBackend&& other)
+ PetscVector(PetscVector&& other)
+ : dofMap_(other.dofMap_)
{
swap(*this, other);
}
/// Copy constructor
- VectorBackend(VectorBackend const& other)
- : comm_(other.comm_)
+ PetscVector(PetscVector const& other)
+ : dofMap_(other.dofMap_)
, initialized_(other.initialized_)
{
if (other.initialized_) {
@@ -60,22 +58,22 @@ namespace AMDiS
}
/// Destructor
- ~VectorBackend()
+ ~PetscVector()
{
destroy();
}
/// Move assignment operator
- VectorBackend& operator=(VectorBackend&& other)
+ PetscVector& operator=(PetscVector&& other)
{
swap(*this, other);
return *this;
}
/// Copy assignment operator
- VectorBackend& operator=(VectorBackend const& other)
+ PetscVector& operator=(PetscVector const& other)
{
- return *this = VectorBackend(other);
+ return *this = PetscVector(other);
}
/// Return the data-vector \ref vector_
@@ -93,13 +91,13 @@ namespace AMDiS
/// Return the number of entries in the local part of the vector
std::size_t localSize() const
{
- return dofMap().localSize();
+ return dofMap_.localSize();
}
/// Return the number of entries in the global vector
std::size_t globalSize() const
{
- return dofMap().globalSize();
+ return dofMap_.globalSize();
}
/// Resize the \ref vector_ to the size \p s
@@ -120,8 +118,8 @@ namespace AMDiS
VecGetLocalSize(vector_, &localSize);
// re-create the vector
- if (std::size_t(localSize) != dofMap().localSize() ||
- std::size_t(globalSize) != dofMap().globalSize()) {
+ if (std::size_t(localSize) != dofMap_.localSize() ||
+ std::size_t(globalSize) != dofMap_.globalSize()) {
destroy();
create();
info(3, " destroy old and create new vector needed {} seconds", t.elapsed());
@@ -157,20 +155,21 @@ namespace AMDiS
/// Access the entry \p i of the \ref vector with read-only-access.
// [[not collective]]
- PetscScalar at(std::size_t i) const
+ template <class MultiIndex>
+ PetscScalar at(MultiIndex const& idx) const
{
PetscScalar y = 0;
- gather(std::array<std::size_t,1>{i}, &y);
+ gather(std::array<MultiIndex,1>{idx}, &y);
return y;
}
/// Access the entry \p i of the \ref vector with write-access.
// [[not collective]]
- template <class Assign>
- void insert(std::size_t dof, PetscScalar value, Assign)
+ template <class MultiIndex, class Assign>
+ void insert(MultiIndex const& dof, PetscScalar value, Assign)
{
assert(initialized_);
- VecSetValue(vector_, dofMap().global(dof), value, AssignMode<Assign>::value);
+ VecSetValue(vector_, dofMap_.global(flatMultiIndex(dof)), value, AssignMode<Assign>::value);
}
/// Collect all values from this vector to indices given by `localInd` and store it into the output buffer.
@@ -201,7 +200,7 @@ namespace AMDiS
auto ind_it = std::begin(localInd);
auto mask_it = std::begin(mask);
for (; ind_it != std::end(localInd); ++mask_it, ++ind_it)
- globalInd.push_back(*mask_it ? dofMap().global(*ind_it) : -1);
+ globalInd.push_back(*mask_it ? dofMap_.global(flatMultiIndex(*ind_it)) : -1);
if (! std::is_same<MaskRange, FakeContainer<bool,true>>::value)
VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
@@ -226,13 +225,14 @@ namespace AMDiS
for (auto const& dof : localInd)
{
- if (dofMap().owner(dof)) {
+ size_type i = flatMultiIndex(dof);
+ if (dofMap_.owner(i)) {
// the index is a purely local entry
- f(dof, ptr[dofMap().dofToLocal(dof)]);
+ f(dof, ptr[dofMap_.dofToLocal(i)]);
}
else {
// ghost entries are stored at the end of the array
- f(dof, ptr[dofMap().localSize() + dofMap().dofToGhost(dof)]);
+ f(dof, ptr[dofMap_.localSize() + dofMap_.dofToGhost(i)]);
}
}
@@ -240,6 +240,20 @@ namespace AMDiS
VecGhostRestoreLocalForm(vector_, &localGhostedVec);
}
+ template <class Func>
+ void forEach(Func&& f)
+ {
+ forEach(mappedRangeView(Dune::range(localSize()),
+ [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
+ }
+
+ template <class Func>
+ void forEach(Func&& f) const
+ {
+ forEach(mappedRangeView(Dune::range(localSize()),
+ [](auto i) { return Dune::Functions::FlatMultiIndex<decltype(i)>{i}; }), FWD(f));
+ }
+
/// Set all entries to \p value, including the ghost entries
// [[not collective]]
@@ -265,20 +279,14 @@ namespace AMDiS
VecGhostRestoreLocalForm(vector_, &localGhostedVec);
}
- friend void swap(VectorBackend& lhs, VectorBackend& rhs)
+ friend void swap(PetscVector& lhs, PetscVector& rhs)
{
- std::swap(lhs.comm_, rhs.comm_);
+ assert(&lhs.dofMap_ == &rhs.dofMap_);
std::swap(lhs.vector_, rhs.vector_);
std::swap(lhs.initialized_, rhs.initialized_);
}
private:
- // The local-to-global map
- auto const& dofMap() const
- {
- return comm_->dofMap();
- }
-
// An MPI Communicator or PETSC_COMM_SELF
MPI_Comm comm() const
{
@@ -293,7 +301,7 @@ namespace AMDiS
void create()
{
VecCreateGhost(comm(),
- dofMap().localSize(), dofMap().globalSize(), dofMap().ghostSize(), dofMap().ghostIndices().data(),
+ dofMap_.localSize(), dofMap_.globalSize(), dofMap_.ghostSize(), dofMap_.ghostIndices().data(),
&vector_);
}
@@ -307,7 +315,8 @@ namespace AMDiS
}
private:
- std::shared_ptr<Communication const> comm_;
+ // The local-to-global map
+ DofMap const& dofMap_;
/// The data-vector
mutable Vec vector_ = nullptr;
diff --git a/src/amdis/typetree/MultiIndex.hpp b/src/amdis/typetree/MultiIndex.hpp
index 62f9a087..fb997f13 100644
--- a/src/amdis/typetree/MultiIndex.hpp
+++ b/src/amdis/typetree/MultiIndex.hpp
@@ -20,4 +20,13 @@ namespace AMDiS
return idx[0];
}
+ template <class MultiIndex>
+ void multiIndexPushFront(MultiIndex& M, std::size_t M0)
+ {
+ M.resize(M.size() + 1);
+ for (std::size_t i = M.size()-1; i > 0; --i)
+ M[i] = M[i-1];
+ M[0] = M0;
+ }
+
} // end namespace AMDiS
diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp
index e673e38b..ee8feb7b 100644
--- a/test/DOFVectorTest.cpp
+++ b/test/DOFVectorTest.cpp
@@ -22,33 +22,38 @@ using namespace Dune::Functions::BasisFactory;
template <class B, class T>
void test_dofvector(DOFVector<B,T>& vec)
{
- using size_type = typename DOFVector<B,T>::GlobalBasis::LocalView::size_type;
- using value_type = typename DOFVector<B,T>::value_type;
+ using MultiIndex = typename B::MultiIndex;
auto const& basis = *vec.basis();
- AMDIS_TEST( vec.localSize() == basis.dimension() );
- vec.resizeZero();
- vec.init(false);
- vec.insert(0, value_type(1));
+ vec.init(true);
+ AMDIS_TEST_EQ(vec.localSize(), basis.size());
+
+ vec.insert(MultiIndex{0}, T(1));
vec.finish();
- value_type min(10);
- value_type max(-10);
- auto lv = basis.localView();
- auto const& gv = basis.gridView();
- for (auto const& e : elements(gv))
- {
- lv.bind(e);
- for (size_type i = 0; i < lv.size(); ++i)
- {
- auto value = vec.at(lv.index(i));
- min = std::min(min, value);
- max = std::max(max, value);
- }
- }
- AMDIS_TEST( min == T(0) || vec.localSize() == 1 );
- AMDIS_TEST( max == T(1) );
+ AMDIS_TEST_EQ(vec.at(MultiIndex{0}), T(1));
+
+ T min(10);
+ T max(-10);
+ vec.forEach([&min,&max](auto, auto value) {
+ min = std::min(min, value);
+ max = std::max(max, value);
+ });
+
+ if (vec.localSize() > 1)
+ AMDIS_TEST_EQ(min, T(0));
+ AMDIS_TEST_EQ(max, T(1));
+
+ // interpolate constant function on DOFVector
+ vec << T(2);
+ AMDIS_TEST_EQ(vec.at(MultiIndex{0}), T(2));
+
+ vec.interpolate(T(3));
+ AMDIS_TEST_EQ(vec.at(MultiIndex{0}), T(3));
+
+ vec.interpolate_noalias(T(4));
+ AMDIS_TEST_EQ(vec.at(MultiIndex{0}), T(4));
}
int main(int argc, char** argv)
@@ -97,15 +102,15 @@ int main(int argc, char** argv)
// Conversion from Dune::Functions::DefaultGlobalBasis
auto vec2 = makeDOFVector(makeBasis(gridView, preBasis));
- test_dofvector(vec1);
for (auto const& e : elements(gridView))
grid.mark(1, e);
+
grid.preAdapt();
grid.adapt();
grid.postAdapt();
- std::size_t newSize = basis.dimension();
- AMDIS_TEST_EQ(vec1.localSize(), newSize);
- AMDIS_TEST_EQ(vec2.localSize(), newSize);
+
+ test_dofvector(vec1);
+ test_dofvector(vec2);
}
report_errors();
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
index 58d51d00..048803c8 100644
--- a/test/DiscreteFunctionTest.cpp
+++ b/test/DiscreteFunctionTest.cpp
@@ -63,8 +63,8 @@ int main(int argc, char** argv)
auto& U_m = *prob.solutionVector();
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- DiscreteFunction u0_c(U_c);
- DiscreteFunction u0_m(U_m);
+ DiscreteFunction u0_c(U_c.coefficients(), *U_c.basis());
+ DiscreteFunction u0_m(U_m.coefficients(), *U_m.basis());
#endif
auto u1_c = makeDiscreteFunction(U_c);
auto u1_m = makeDiscreteFunction(U_m);
diff --git a/test/ExpressionsTest.cpp b/test/ExpressionsTest.cpp
index d8183517..aa1ecae2 100644
--- a/test/ExpressionsTest.cpp
+++ b/test/ExpressionsTest.cpp
@@ -94,7 +94,7 @@ int main(int argc, char** argv)
// integration of expressions
- auto gv = u.basis()->gridView();
+ auto gv = u.basis().gridView();
DUNE_UNUSED auto int1 = integrate(op1, gv, 5);
DUNE_UNUSED auto int2 = integrate(op2, gv, 5);
diff --git a/test/IntegrateTest.cpp b/test/IntegrateTest.cpp
index ea82e2c9..adb28a08 100644
--- a/test/IntegrateTest.cpp
+++ b/test/IntegrateTest.cpp
@@ -24,7 +24,7 @@ int main(int argc, char** argv)
prob.initialize(INIT_ALL);
auto u = prob.solution(0);
- auto gv = u.basis()->gridView();
+ auto gv = u.basis().gridView();
u << 1.0;
double i1 = integrate(1.0, gv);
diff --git a/test/ObserverTest.cpp b/test/ObserverTest.cpp
index 8f82e888..93e66256 100644
--- a/test/ObserverTest.cpp
+++ b/test/ObserverTest.cpp
@@ -63,8 +63,8 @@ class B
{
public:
B(std::shared_ptr<A> a)
- : Observer<event::fooEvent>(*a)
- , Observer<event::bazEvent>(*a)
+ : Observer<event::fooEvent>(a)
+ , Observer<event::bazEvent>(a)
, a_(a)
{}
@@ -92,8 +92,8 @@ class C
{
public:
C(std::shared_ptr<B> b)
- : Observer<event::barEvent>(*b->a())
- , Observer<event::bazEvent>(*b)
+ : Observer<event::barEvent>(b->a())
+ , Observer<event::bazEvent>(b)
{}
void updateImpl(event::barEvent) override
@@ -128,8 +128,8 @@ class E
{
public:
E(std::shared_ptr<A> a, std::shared_ptr<D> d)
- : Observer<event::fooEvent>(*a)
- , Observer<event::barEvent>(*d)
+ : Observer<event::fooEvent>(a)
+ , Observer<event::barEvent>(d)
{}
void updateImpl(event::fooEvent) override
diff --git a/test/PETScCommTest.cpp b/test/PETScCommTest.cpp
index 8f4b2df1..e712e6d8 100644
--- a/test/PETScCommTest.cpp
+++ b/test/PETScCommTest.cpp
@@ -35,7 +35,7 @@ void test_basis(Basis const& basis, std::string gridName, std::string basisName)
using Comm = typename Traits::Comm;
auto comm = CommunicationCreator<Comm>::create(basis, "communication");
- auto const& dofMap = comm->dofMap();
+ auto const& dofMap = comm.dofMap();
using GlobalIndex = typename Comm::DofMap::GlobalIndex;
std::vector<std::vector<GlobalIndex>> data(world.size());
diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp
index 4d99b534..5f3563fc 100644
--- a/test/ProblemStatTest.cpp
+++ b/test/ProblemStatTest.cpp
@@ -31,6 +31,7 @@ void test()
// use T as coefficient type
using Param = DefaultProblemTraits<Basis, T>;
+ msg("Test<{}>", __PRETTY_FUNCTION__);
ProblemStat<Param> prob("ellipt", grid, basis);
prob.initialize(INIT_ALL);
prob.boundaryManager()->setBoxBoundary({1,1,2,2});
@@ -42,12 +43,12 @@ void test()
AdaptInfo adaptInfo("adapt");
- // test copy constructor of problem
+ msg(" test copy constructor of problem");
auto prob2(prob);
AdaptStationary adaptStat2("adapt", prob2, adaptInfo);
adaptStat2.adapt();
- // test move constructor of problem
+ msg(" test move constructor of problem");
auto prob3(std::move(prob2));
AdaptStationary adaptStat3("adapt", prob3, adaptInfo);
adaptStat3.adapt();
@@ -57,19 +58,19 @@ int main(int argc, char** argv)
{
Environment env(argc, argv);
-#if !HAVE_PETSC || PETSC_USE_REAL_SINGLE
+#if !AMDIS_HAS_PETSC || PETSC_USE_REAL_SINGLE
test<float>();
#endif
-#if !HAVE_PETSC || PETSC_USE_REAL_DOUBLE
+#if !AMDIS_HAS_PETSC || PETSC_USE_REAL_DOUBLE
test<double>();
#endif
-#if !HAVE_PETSC
+#if !AMDIS_HAS_PETSC
test<long double>();
#endif
-#if HAVE_QUADMATH && (!HAVE_PETSC || PETSC_USE_REAL___FLOAT128)
+#if HAVE_QUADMATH && (!AMDIS_HAS_PETSC || PETSC_USE_REAL___FLOAT128) && (!AMDIS_HAS_EIGEN)
test<Dune::Float128>();
#endif
--
GitLab