diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 4735e5f3fcbb2b495341be394463c56677db9ad6..86fea59d461bd00da837149013b15258a36b794f 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -7,17 +7,21 @@
using namespace AMDiS;
-using Grid = Dune::YaspGrid<2>;
-using Basis = TaylorHoodBasis<Grid>;
+struct NavierStokesBasis
+{
+ using Grid = Dune::YaspGrid<GRIDDIM>;
+ using GlobalBasis = typename TaylorHoodBasis<Grid>::GlobalBasis;
+ using CoefficientType = double;
+};
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
- ProblemStat<Basis> prob("stokes");
+ ProblemStat<NavierStokesBasis> prob("stokes");
prob.initialize(INIT_ALL);
- ProblemInstat<Basis> probInstat("stokes", prob);
+ ProblemInstat<NavierStokesBasis> probInstat("stokes", prob);
probInstat.initialize(INIT_UH_OLD);
double viscosity = 1.0;
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 59d521695cbce90acacd93b6b22ae2d2dea681ae..dad37e90430b380cdca1867ba025156edec714b9 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -16,19 +16,18 @@ namespace AMDiS
{
/// Implementation of interface \ref AssemblerBase
/**
- * \tparam LC LocalContext
+ * \tparam Traits See \ref DefaultAssemblerTraits
**/
- template <class LC, class Operator, class... Nodes>
+ template <class Traits, class Operator, class... Nodes>
class Assembler
- : public AssemblerInterface<LC, Nodes...>
+ : public AssemblerInterface<Traits, Nodes...>
{
- using Super = AssemblerInterface<LC, Nodes...>;
-
- private:
+ using Super = AssemblerInterface<Traits, Nodes...>;
+ using LocalContext = typename Traits::LocalContext;
using Element = typename Super::Element;
using Geometry = typename Super::Geometry;
- using ElementMatrixVector = typename Super::ElementMatrixVector;
+ using ElementContainer = typename Traits::ElementContainer;
public:
@@ -77,11 +76,11 @@ namespace AMDiS
* \ref calculateElementVector or \ref calculateElementMatrix on the
* vector or matrix operator, respectively.
**/
- void assemble(LC const& localContext, Nodes const&... nodes,
- ElementMatrixVector& elementMatrixVector) final
+ void assemble(LocalContext const& localContext, Nodes const&... nodes,
+ ElementContainer& ElementContainer) final
{
- ContextGeometry<LC> contextGeo{localContext, element(), geometry()};
- assembleImpl(contextGeo, nodes..., elementMatrixVector);
+ ContextGeometry<LocalContext> contextGeo{localContext, element(), geometry()};
+ assembleImpl(contextGeo, nodes..., ElementContainer);
}
@@ -133,10 +132,10 @@ namespace AMDiS
/// Generate a \ref Assembler on a given LocalContext `LC` (element or intersection)
- template <class LC, class Operator, class... Nodes>
+ template <class Traits, class Operator, class... Nodes>
auto makeAssembler(Operator&& op, Nodes const&...)
{
- return Assembler<LC, Underlying_t<Operator>, Nodes...>{FWD(op)};
+ return Assembler<Traits, Underlying_t<Operator>, Nodes...>{FWD(op)};
}
} // end namespace AMDiS
diff --git a/src/amdis/AssemblerInterface.hpp b/src/amdis/AssemblerInterface.hpp
index 68e9b5d289c1732cbaeaaceaa3b41fb77fd3eb5d..5cfee9affb986feddc078d5ccc4d84edd221dd4f 100644
--- a/src/amdis/AssemblerInterface.hpp
+++ b/src/amdis/AssemblerInterface.hpp
@@ -9,10 +9,18 @@
namespace AMDiS
{
+ template <class LC, class C>
+ struct DefaultAssemblerTraits
+ {
+ using LocalContext = LC;
+ using ElementContainer = C;
+ };
+
/// Abstract base-class of a \ref Assembler
- template <class LocalContext, class... Nodes>
+ template <class Traits, class... Nodes>
class AssemblerInterface
{
+ using LocalContext = typename Traits::LocalContext;
using ContextType = Impl::ContextTypes<LocalContext>;
public:
@@ -25,14 +33,6 @@ namespace AMDiS
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
- using ElementMatrix = Dune::DynamicMatrix<double>; // TODO: choose correct value_type
- using ElementVector = Dune::DynamicVector<double>;
-
- /// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
- using ElementMatrixVector = std::conditional_t<
- (sizeof...(Nodes)==1), ElementVector, std::conditional_t<
- (sizeof...(Nodes)==2), ElementMatrix, void>>;
-
public:
/// Virtual destructor
virtual ~AssemblerInterface() = default;
@@ -44,9 +44,9 @@ namespace AMDiS
virtual void unbind() = 0;
/// Assemble an element matrix or element vector on the test- (and trial-) function node(s)
- virtual void assemble(LocalContext const& localContext,
+ virtual void assemble(typename Traits::LocalContext const& localContext,
Nodes const&... nodes,
- ElementMatrixVector& elementMatrixVector) = 0;
+ typename Traits::ElementContainer& elementMatrixVector) = 0;
};
} // end namespace AMDiS
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index b6f01e9377ea03d15b7f68806fd13ff409db6e21..bf2efd7834a493243a3665923ff59bc072c814d8 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -156,7 +156,8 @@ namespace AMDiS
// Interpolate from possibly vanishing elements
if (mightCoarsen) {
auto maxLevel = gv.grid().maxLevel();
- ctype const checkInsideTolerance = std::sqrt(std::numeric_limits<ctype>::epsilon());
+ using std::sqrt;
+ ctype const checkInsideTolerance = sqrt(std::numeric_limits<ctype>::epsilon());
for (auto const& e : elements(gv))
{
auto father = e;
diff --git a/src/amdis/Mesh.hpp b/src/amdis/Mesh.hpp
index 1d1d0948b1e7397d733682770188af35741a7bea..be361b7392f7b0e49fce52da89bf93473e1aad39 100644
--- a/src/amdis/Mesh.hpp
+++ b/src/amdis/Mesh.hpp
@@ -141,7 +141,7 @@ namespace AMDiS
static std::unique_ptr<Grid> create(std::string const& meshName)
{
- Dune::FieldVector<double, dim> L; L = 1.0; // extension of the domain
+ Dune::FieldVector<T, dim> L; L = 1.0; // extension of the domain
Parameters::get(meshName + "->dimension", L);
auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
@@ -161,8 +161,8 @@ namespace AMDiS
static std::unique_ptr<Grid> create(std::string const& meshName)
{
- Dune::FieldVector<double, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain
- Dune::FieldVector<double, dim> upperright; upperright = 1.0; // Upper right corner of the domain
+ Dune::FieldVector<T, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain
+ Dune::FieldVector<T, dim> upperright; upperright = 1.0; // Upper right corner of the domain
Parameters::get(meshName + "->min corner", lowerleft);
Parameters::get(meshName + "->max corner", upperright);
diff --git a/src/amdis/OperatorList.hpp b/src/amdis/OperatorList.hpp
index bf780920a3f2a977174ec7d84725a187fda0d532..929634127dcc9c496d13e7508f5351aae4b187ba 100644
--- a/src/amdis/OperatorList.hpp
+++ b/src/amdis/OperatorList.hpp
@@ -20,7 +20,7 @@ namespace AMDiS
};
}
- template <class GridView>
+ template <class GridView, class Container>
class OperatorLists
{
using Element = typename GridView::template Codim<0>::Entity;
@@ -102,11 +102,14 @@ namespace AMDiS
auto const& onBoundary() const { return boundary_; }
private:
+ using ElementTraits = DefaultAssemblerTraits<Element,Container>;
+ using IntersectionTraits = DefaultAssemblerTraits<Intersection,Container>;
+
/// The type of local operators associated with grid elements
- using ElementAssembler = AssemblerInterface<Element, Nodes...>;
+ using ElementAssembler = AssemblerInterface<ElementTraits, Nodes...>;
/// The type of local operators associated with grid intersections
- using IntersectionAssembler = AssemblerInterface<Intersection, Nodes...>;
+ using IntersectionAssembler = AssemblerInterface<IntersectionTraits, Nodes...>;
/// List of operators to be assembled on grid elements
std::vector<std::shared_ptr<ElementAssembler>> element_;
@@ -136,12 +139,12 @@ namespace AMDiS
};
- template <class RowBasis, class ColBasis>
+ template <class RowBasis, class ColBasis, class ElementMatrix>
using MatrixOperators
- = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView>::template MatData>;
+ = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView,ElementMatrix>::template MatData>;
- template <class GlobalBasis>
+ template <class GlobalBasis, class ElementVector>
using VectorOperators
- = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
+ = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView,ElementVector>::template VecData>;
} // end namespace AMDiS
diff --git a/src/amdis/PeriodicBC.inc.hpp b/src/amdis/PeriodicBC.inc.hpp
index 94aa74693cb826dc2c1a8f8cedc85d5c64bad8cd..0c92fecd8ba3fbbe67237bb39c7f1bd59de33d1c 100644
--- a/src/amdis/PeriodicBC.inc.hpp
+++ b/src/amdis/PeriodicBC.inc.hpp
@@ -43,7 +43,8 @@ template <class D, class MI>
void PeriodicBC<D,MI>::
initAssociations(Basis const& basis)
{
- static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ using std::sqrt;
+ static const auto tol = sqrt(std::numeric_limits<typename D::field_type>::epsilon());
periodicNodes_.clear();
periodicNodes_.resize(basis.dimension(), false);
@@ -107,8 +108,9 @@ namespace Impl
bool operator()(D const& lhs, D const& rhs) const
{
+ using std::abs;
for (int i = 0; i < D::dimension; i++) {
- if (std::abs(lhs[i] - rhs[i]) < tol_)
+ if (abs(lhs[i] - rhs[i]) < tol_)
continue;
return lhs[i] < rhs[i];
}
@@ -127,7 +129,8 @@ template <class D, class MI>
void PeriodicBC<D,MI>::
initAssociations2(Basis const& basis)
{
- static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ using std::sqrt;
+ static const auto tol = sqrt(std::numeric_limits<typename D::field_type>::epsilon());
periodicNodes_.clear();
periodicNodes_.resize(basis.dimension(), false);
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 97046397ba315ccdddc8a42294718ec9ef0e5c4c..35cea992b5fd4f2ab7f801e3402f835409ee219c 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -70,8 +70,8 @@ namespace AMDiS
/// Dimension of the world
static constexpr int dow = Grid::dimensionworld;
- using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
- using SystemVector = DOFVector<GlobalBasis, double>;
+ using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, typename Traits::CoefficientType>;
+ using SystemVector = DOFVector<GlobalBasis, typename Traits::CoefficientType>;
using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index 0be129b51b7e447dfcc5b3733ed6860783f60ec4..6bd0c49b0cc0d38d372f1c9471103e055790fd8d 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -64,7 +64,7 @@ namespace AMDiS
};
- template <class Grid, class PreBasisCreator>
+ template <class Grid, class PreBasisCreator, class T = double>
struct DefaultBasisCreator
{
using GridView = typename Grid::LeafGridView;
@@ -75,6 +75,7 @@ namespace AMDiS
}
using GlobalBasis = decltype(create(std::declval<GridView>()));
+ using CoefficientType = T;
};
/// \brief ProblemStatTraits for a (composite) basis composed of
diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
index 3a367a3a94d95b925a7de0db9c9d721753c43fac..9323112d25ea27f1dc52351974dba3a373d14c1c 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -32,7 +32,7 @@ namespace AMDiS
template <class GB, class VT, class TP>
class DiscreteFunction
{
- static_assert(std::is_arithmetic<VT>::value, "");
+ // static_assert(Dune::IsNumber<VT>::value, "");
private:
using GlobalBasis = GB;
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.hpp b/src/amdis/linearalgebra/DOFMatrixBase.hpp
index 213759edfdd00d6a8c02c47a5f7d3dda79e87581..a197930c5cc0b91ca437e4e24bb4cc46f76317c7 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.hpp
@@ -37,12 +37,13 @@ namespace AMDiS
/// The index/size - type
using size_type = typename RowBasis::size_type;
+ using value_type = typename Backend::value_type;
/// The type of the data matrix used in the backend
using BaseMatrix = typename Backend::BaseMatrix;
/// The type of the matrix filled on an element with local contributions
- using ElementMatrix = Dune::DynamicMatrix<double>;
+ using ElementMatrix = Dune::DynamicMatrix<value_type>;
public:
/// Constructor. Stores the shared_ptr to the bases.
@@ -166,7 +167,7 @@ namespace AMDiS
ElementMatrix elementMatrix_;
/// List of operators associated to row/col node
- MatrixOperators<RowBasis,ColBasis> operators_;
+ MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
index f06bbde64f2d68fc52042ece6d62ebc0fe93e57e..e37a52c2ab05d68aa35c0eaa7df22121c1895e0f 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
@@ -32,10 +32,11 @@ void DOFMatrixBase<RB,CB,B>::
insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
ElementMatrix const& elementMatrix)
{
+ using std::abs;
for (size_type i = 0; i < rowLocalView.size(); ++i) {
size_type const row = flatMultiIndex(rowLocalView.index(i));
for (size_type j = 0; j < colLocalView.size(); ++j) {
- if (std::abs(elementMatrix[i][j]) > threshold<double>) {
+ if (abs(elementMatrix[i][j]) > threshold<double>) {
size_type const col = flatMultiIndex(colLocalView.index(j));
backend_.insert(row, col, elementMatrix[i][j]);
}
@@ -59,8 +60,9 @@ addOperator(ContextTag contextTag, Expr const& expr,
auto j = child(colBasis_->localView().tree(), makeTreePath(col));
using LocalContext = typename ContextTag::type;
+ using Traits = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
auto op = makeLocalOperator<LocalContext>(expr, rowBasis_->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i, j));
+ auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i, j));
operators_[i][j].push(contextTag, std::move(localAssembler));
}
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index 62603d6617af9fa050ce1ef79d7965c9bd265b45..9f17f9aa9f3462b6183656ced9dd57aa81e23d44 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -73,7 +73,7 @@ namespace AMDiS
using BaseVector = typename Backend::BaseVector;
/// The type of the vector filled on an element with local contributions
- using ElementVector = Dune::DynamicVector<double>;
+ using ElementVector = Dune::DynamicVector<value_type>;
/// Defines an interface to transfer the data during grid adaption
using DataTransfer = DataTransferInterface<Self>;
@@ -290,7 +290,7 @@ namespace AMDiS
ElementVector elementVector_;
/// List of operators associated to nodes, filled in \ref addOperator().
- VectorOperators<GlobalBasis> operators_;
+ VectorOperators<GlobalBasis,ElementVector> operators_;
/// Data interpolation when the grid changes, set by default
/// to \ref DataTransferOperation::INTERPOLATE.
diff --git a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
index 4795c768f9f58a9579096d41a0dc40514dbcd3eb..783eb8998bced4271c82a18f3414245b6481a838 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
@@ -30,8 +30,9 @@ template <class GB, class B>
void DOFVectorBase<GB,B>::
insert(LocalView const& localView, ElementVector const& elementVector)
{
+ using std::abs;
for (size_type i = 0; i < localView.size(); ++i) {
- if (std::abs(elementVector[i]) > threshold<double>) {
+ if (abs(elementVector[i]) > threshold<double>) {
size_type const idx = flatMultiIndex(localView.index(i));
backend_[idx] += elementVector[i];
}
@@ -50,8 +51,9 @@ addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
auto i = child(basis_->localView().tree(), makeTreePath(path));
using LocalContext = typename ContextTag::type;
+ using Traits = DefaultAssemblerTraits<LocalContext, ElementVector>;
auto op = makeLocalOperator<LocalContext>(expr, basis_->gridView());
- auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i));
+ auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i));
operators_[i].push(contextTag, std::move(localAssembler));
}
diff --git a/src/amdis/linearalgebra/eigen/DirectRunner.hpp b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
index afcfa23453e06223bc329d669ff32dc5ed01c6f4..f37208b9bfb65f29c9c0c586fd5e1e555242e680 100644
--- a/src/amdis/linearalgebra/eigen/DirectRunner.hpp
+++ b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
@@ -10,23 +10,24 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::DirectRunner
- * \brief \implements RunnerInterface for the (external) direct solvers
- */
- template <class LU, class Matrix, class VectorX, class VectorB>
+ * \ingroup Solver
+ * \class AMDiS::DirectRunner
+ * \brief \implements RunnerInterface for the (external) direct solvers
+ */
+ template <template <class> class Solver, class Matrix, class VectorX, class VectorB>
class DirectRunner
: public RunnerInterface<Matrix, VectorX, VectorB>
{
protected:
using Super = RunnerInterface<Matrix, VectorX, VectorB>;
+ using EigenSolver = Solver<Matrix>;
public:
/// Constructor.
DirectRunner(std::string const& prefix)
: solver_{}
{
- SolverConfig<LU>::init(prefix, solver_);
+ SolverConfig<Solver>::init(prefix, solver_);
Parameters::get(prefix + "->reuse pattern", reusePattern_);
}
@@ -64,7 +65,7 @@ namespace AMDiS
}
private:
- LU solver_;
+ EigenSolver solver_;
bool reusePattern_ = false;
bool initialized_ = false;
};
diff --git a/src/amdis/linearalgebra/eigen/SolverCreator.hpp b/src/amdis/linearalgebra/eigen/SolverCreator.hpp
index 4110f9efd959db502d73224fe6ae3b8292af8135..8d28383777227e808336470f5db592caf48f5bb8 100644
--- a/src/amdis/linearalgebra/eigen/SolverCreator.hpp
+++ b/src/amdis/linearalgebra/eigen/SolverCreator.hpp
@@ -100,7 +100,7 @@ namespace AMDiS
using SolverCreator
= IterativeSolverCreator<Matrix, VectorX, VectorB, IterativeSolver>;
- template <class DirectSolver>
+ template <template <class> class DirectSolver>
using DirectSolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;
@@ -148,21 +148,27 @@ namespace AMDiS
auto dgmres = new SolverCreator<DGMRES>;
Map::addCreator("dgmres", dgmres);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#if HAVE_SUITESPARSE_UMFPACK
// sparse LU factorization and solver based on UmfPack
- auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU<Matrix>>;
+ auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU>;
Map::addCreator("umfpack", umfpack);
#endif
#if HAVE_SUPERLU
// sparse direct LU factorization and solver based on the SuperLU library
- auto superlu = new DirectSolverCreator<Eigen::SuperLU<Matrix>>;
+ auto superlu = new DirectSolverCreator<Eigen::SuperLU>;
Map::addCreator("superlu", superlu);
#endif
- // default iterative solver
- Map::addCreator("default", gmres);
-
// default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
Map::addCreator("direct", umfpack);
diff --git a/src/amdis/linearalgebra/istl/DirectRunner.hpp b/src/amdis/linearalgebra/istl/DirectRunner.hpp
index f2785b644291c421e8d50eae16e0014a79c12a8d..6a3f7852fa72a1cdf1e8cf094ebedca81a9011e1 100644
--- a/src/amdis/linearalgebra/istl/DirectRunner.hpp
+++ b/src/amdis/linearalgebra/istl/DirectRunner.hpp
@@ -3,6 +3,8 @@
#include <algorithm>
#include <string>
+#include <dune/common/ftraits.hh>
+
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
@@ -10,17 +12,18 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::DirectRunner
- * \brief \implements RunnerInterface for the (external) direct solvers
- */
- template <class Solver, class Matrix, class VectorX, class VectorB>
+ * \ingroup Solver
+ * \class AMDiS::DirectRunner
+ * \brief \implements RunnerInterface for the (external) direct solvers
+ */
+ template <template <class> class Solver, class Mat, class Sol, class Rhs>
class DirectRunner
- : public RunnerInterface<Matrix, VectorX, VectorB>
+ : public RunnerInterface<Mat, Sol, Rhs>
{
protected:
- using Super = RunnerInterface<Matrix, VectorX, VectorB>;
- using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
+ using Super = RunnerInterface<Mat, Sol, Rhs>;
+ using BaseSolver = Dune::InverseOperator<Sol, Rhs>;
+ using ISTLSolver = Solver<Mat>;
public:
/// Constructor.
@@ -29,7 +32,7 @@ namespace AMDiS
{}
/// Implementes \ref RunnerInterface::init()
- void init(Matrix const& A) override
+ void init(Mat const& A) override
{
solver_ = solverCreator_.create(A);
}
@@ -41,14 +44,14 @@ namespace AMDiS
}
/// Implementes \ref RunnerInterface::solve()
- int solve(Matrix const& A, VectorX& x, VectorB const& b,
- SolverInfo& solverInfo) override
+ int solve(Mat const& A, Sol& x, Rhs const& b,
+ SolverInfo& solverInfo) override
{
// storing some statistics
Dune::InverseOperatorResult statistics;
// solve the linear system
- VectorB _b = b;
+ Rhs _b = b;
solver_->apply(x, _b, statistics);
solverInfo.setRelResidual(statistics.reduction);
@@ -57,7 +60,7 @@ namespace AMDiS
}
private:
- ISTLSolverCreator<Solver> solverCreator_;
+ ISTLSolverCreator<Solver<Mat>> solverCreator_;
std::shared_ptr<BaseSolver> solver_;
};
}
diff --git a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
index 923e1d006de12e9819af378e41d04fc00087ead9..ca8c16e42731ec3d7aec42a2555044f98f8f10e1 100644
--- a/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
+++ b/src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp
@@ -102,19 +102,31 @@ namespace AMDiS
auto ssor = new PreconCreator<Dune::SeqSSOR<Mat, Sol, Rhs>>;
Map::addCreator("ssor", ssor);
+ init_ilu(std::is_arithmetic<typename Dune::FieldTraits<Mat>::field_type>{});
+ init_amg(std::is_same<typename Dune::FieldTraits<Mat>::real_type, double>{});
+
+ auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
+ Map::addCreator("richardson", richardson);
+ Map::addCreator("default", richardson);
+ }
+
+ static void init_ilu(std::false_type) {}
+ static void init_ilu(std::true_type)
+ {
auto ilu = new PreconCreator<Dune::SeqILU<Mat, Sol, Rhs>>;
- Map::addCreator("ilu", ilu);
- Map::addCreator("ilu0", ilu);
+ Map::addCreator("ilu", id(ilu));
+ Map::addCreator("ilu0", id(ilu));
+ }
+ static void init_amg(std::false_type) {}
+ static void init_amg(std::true_type)
+ {
auto amg = new AMGPreconCreator<Dune::Amg::AMG, Mat, Sol, Rhs>;
Map::addCreator("amg", amg);
auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Mat, Sol, Rhs>;
Map::addCreator("fastamg", fastamg);
-
- auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
- Map::addCreator("richardson", richardson);
- Map::addCreator("default", richardson);
}
+
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
index 35bcd6bf8986461e3dee1213deeb0c4f053d9334..700cb1fa31978230f069cbd37d1446d882c56e4a 100644
--- a/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
+++ b/src/amdis/linearalgebra/istl/ISTL_Solver.hpp
@@ -143,15 +143,15 @@ namespace AMDiS
using Matrix = Dune::BCRSMatrix<Block,Alloc>;
using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
- template <class ISTLSolver>
+ template <class Solver>
using SolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
- ISTLRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+ ISTLRunner<Solver, Matrix, VectorX, VectorB>>::Creator;
- template <class ISTLSolver>
+ template <template <class M> class Solver>
using DirectSolverCreator
= typename LinearSolver<Matrix, VectorX, VectorB,
- DirectRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+ DirectRunner<Solver, Matrix, VectorX, VectorB>>::Creator;
using Map = CreatorMap<SolverBase>;
@@ -175,25 +175,31 @@ namespace AMDiS
auto fcg = new SolverCreator<Dune::GeneralizedPCGSolver<VectorX>>;
Map::addCreator("fcg", fcg);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<Matrix>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#if HAVE_SUITESPARSE_UMFPACK
- auto umfpack = new DirectSolverCreator<Dune::UMFPack<Matrix>>;
+ auto umfpack = new DirectSolverCreator<Dune::UMFPack>;
Map::addCreator("umfpack", umfpack);
#endif
#if HAVE_SUPERLU
- auto superlu = new DirectSolverCreator<Dune::SuperLU<Matrix>>;
+ auto superlu = new DirectSolverCreator<Dune::SuperLU>;
Map::addCreator("superlu", superlu);
#endif
- // default direct solvers
+ // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
Map::addCreator("direct", umfpack);
#elif HAVE_SUPERLU
Map::addCreator("direct", superlu);
#endif
-
- // default iterative solver
- Map::addCreator("default", gmres);
}
};
diff --git a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
index 2cd2df2b53c1fc6bbc6b51b7affa0f7233786c9a..5f2946d27a28cfc2b993819cc36179b643d0c1fe 100644
--- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
@@ -453,6 +453,15 @@ namespace AMDiS
auto preonly = new SolverCreator<preonly_type>;
Map::addCreator("preonly", preonly);
+ // default iterative solver
+ Map::addCreator("default", gmres);
+
+ init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{});
+ }
+
+ static void init_direct(std::false_type) {}
+ static void init_direct(std::true_type)
+ {
#ifdef HAVE_UMFPACK
auto umfpack = new UmfpackSolverCreator;
Map::addCreator("umfpack", umfpack);
@@ -460,9 +469,6 @@ namespace AMDiS
// default direct solvers
Map::addCreator("direct", umfpack);
#endif
-
- // default iterative solver
- Map::addCreator("default", gmres);
}
};
diff --git a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
index e51a1a746a3fa1c6616fab2933c3d40cedf9983d..2e5c62c18611ac9a40ab090d4d35026aabd3d947 100644
--- a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
@@ -89,9 +89,9 @@ namespace AMDiS
r -= A * x;
// print information about the solution process
- itl::cyclic_iteration<double> iter(r, maxIter_, solverInfo.relTolerance(),
- solverInfo.absTolerance(),
- printCycle_);
+ using T = typename Dune::FieldTraits<typename Vector::value_type>::real_type;
+ itl::cyclic_iteration<T> iter(two_norm(r),
+ maxIter_, solverInfo.relTolerance(), solverInfo.absTolerance(), printCycle_);
iter.set_quite(solverInfo.info() == 0);
iter.suppress_resume(solverInfo.info() == 0);
diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
index f1f382db56fa0ede0967fd78352fb9bcea18a298..01a45ef99f6cd06d6c0c81978ccdb0c1923737ae 100644
--- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
@@ -15,37 +15,49 @@
namespace AMDiS
{
/**
- * \ingroup Solver
- * \class AMDiS::UmfpackRunner
- * \brief \implements RunnerInterface
- * Wrapper for the external
- * [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
- *
- * This is a direct solver for large sparse matrices.
- */
+ * \ingroup Solver
+ * \class AMDiS::UmfpackRunner
+ * \brief \implements RunnerInterface
+ * Wrapper for the external
+ * [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
+ *
+ * This is a direct solver for large sparse matrices.
+ */
template <class Matrix, class Vector>
class UmfpackRunner;
- template <class Matrix, class Vector>
- class UmfpackRunnerBase
- : public RunnerInterface<Matrix, Vector, Vector>
+ template <class T, class Param, class Vector>
+ class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
+ : public RunnerInterface<mtl::compressed2D<T, Param>, Vector, Vector>
{
- protected:
+ using Matrix = mtl::compressed2D<T, Param>;
using Super = RunnerInterface<Matrix, Vector>;
using PreconBase = typename Super::PreconBase;
- using FullMatrix = Matrix;
- using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
+ using SolverType = mtl::mat::umfpack::solver<Matrix>;
public:
/// Constructor. Reads UMFPACK parameters from initfile
- UmfpackRunnerBase(std::string const& prefix)
+ UmfpackRunner(std::string const& prefix)
{
Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ?
Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
Parameters::get(prefix + "->alloc init", allocInit_);
}
+ /// Implementation of \ref RunnerInterface::init()
+ void init(Matrix const& matrix) override
+ {
+ try {
+ solver_.reset(new SolverType(matrix, symmetricStrategy_, allocInit_));
+ } catch (mtl::mat::umfpack::error const& e) {
+ error_exit("UMFPACK_ERROR(factorize, {}) = {}", e.code, e.what());
+ }
+ }
+
+ /// Implementation of \ref RunnerInterface::exit()
+ void exit() override {}
+
/// Implementation of \ref RunnerBase::solve()
int solve(Matrix const& A, Vector& x, Vector const& b,
SolverInfo& solverInfo) override
@@ -55,9 +67,9 @@ namespace AMDiS
int code = 0;
try {
- code = (*solver_)(x, b);
+ code = (*solver_)(x, b);
} catch (mtl::mat::umfpack::error& e) {
- error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what());
+ error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what());
}
auto r = Vector(b);
@@ -70,9 +82,6 @@ namespace AMDiS
return code;
}
- /// Implementation of \ref RunnerInterface::exit()
- void exit() override {}
-
protected:
std::shared_ptr<SolverType> solver_;
@@ -80,33 +89,6 @@ namespace AMDiS
int symmetricStrategy_ = 0;
double allocInit_ = 0.7;
};
-
-
- // specialization for full-matrices
- template <class T, class Param, class Vector>
- class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
- : public UmfpackRunnerBase<mtl::compressed2D<T, Param>, Vector>
- {
- using FullMatrix = mtl::compressed2D<T, Param>;
- using Super = UmfpackRunnerBase<FullMatrix, Vector>;
-
- using SolverType = typename Super::SolverType;
-
- public:
- UmfpackRunner(std::string const& prefix)
- : Super(prefix)
- {}
-
- /// Implementation of \ref RunnerInterface::init()
- void init(FullMatrix const& fullMatrix) override
- {
- try {
- Super::solver_.reset(new SolverType(fullMatrix, Super::symmetricStrategy_, Super::allocInit_));
- } catch (mtl::mat::umfpack::error const& e) {
- error_exit("UMFPACK_ERROR(factorize, {}) = {}", e.code, e.what());
- }
- }
- };
}
#endif // HAVE_UMFPACK
diff --git a/src/amdis/linearalgebra/mtl/itl/details.hpp b/src/amdis/linearalgebra/mtl/itl/details.hpp
index 1d199bc35ef10206b8523c5e1ce80433fd8782d9..310b19fcdf518aa7dedcbe5f28e8185fe13cfb16 100644
--- a/src/amdis/linearalgebra/mtl/itl/details.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/details.hpp
@@ -49,7 +49,8 @@ namespace itl
// }
// }
- inline void rotmat(const double& a, const double& b , double& c, double& s)
+ template <class T>
+ inline void rotmat(const T& a, const T& b , T& c, T& s)
{
using std::abs;
using std::sqrt;
@@ -60,13 +61,13 @@ namespace itl
}
else if ( abs(b) > abs(a) )
{
- double temp = a / b;
+ T temp = a / b;
s = 1.0 / sqrt( 1.0 + temp*temp );
c = temp * s;
}
else
{
- double temp = b / a;
+ T temp = b / a;
c = 1.0 / sqrt( 1.0 + temp*temp );
s = temp * c;
}
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
index 101ff799539c5d03c6514fc049e9913866ac61e2..21e4216b4b43060d870edc37107dd68587756a3c 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
@@ -139,7 +139,8 @@ namespace itl
H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
H[k+1][k] = 0.0;
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
index 2609cda13905a0cde83ee26e82adbaf3ee94b5b1..2f066ee7472f46b3f0e32b010cee9be92dd9a6a6 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
@@ -152,7 +152,8 @@ namespace itl
irange range(num_rows(H));
H[iall][k] = w[range];
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
index e41806020dd2b051e3bf61991698d9eb2d7597ce..c65c7afe97c7dabd7393d48ad365b4dba2ea09bd 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
@@ -140,7 +140,8 @@ namespace itl
H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
H[k+1][k] = 0.0;
- rho = std::abs(s[k+1]);
+ using std::abs;
+ rho = abs(s[k+1]);
}
// reduce k, to get regular matrix
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
index c8bdb83643e6302bee33d7397269c3b30074559e..eea759692e35db89045145db006892faef69998a 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
@@ -148,7 +148,8 @@ namespace itl
irange range(num_rows(H));
H[iall][k] = w[range];
- rho = std::abs(s[k+1]) / bnrm2;
+ using std::abs;
+ rho = abs(s[k+1]) / bnrm2;
}
// reduce k, to get regular matrix
diff --git a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
index ef7360dd115653cebf3d8a80bebb4f00f26e1b5f..95e7761fcebd726f2727b95748c11aea49e6e4df 100644
--- a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
@@ -42,7 +42,7 @@ namespace AMDiS
"RN must be Leaf-Node and CN must be a Power-Node.");
static const std::size_t CHILDREN = CN::CHILDREN;
- static_assert( std::is_same<typename GridFct::Range, Dune::FieldVector<double, int(CHILDREN)>>::value, "" );
+ static_assert( Size_v<typename GridFct::Range> == CHILDREN, "" );
auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
diff --git a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
index dab0476f1527d4eedb87831db68664a0fe7cf85e..2d04c6f3b9b854880f8d696cf569f47400b1ace5 100644
--- a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
@@ -154,7 +154,7 @@ namespace AMDiS
Mat& elementMatrix, tag::matrix)
{
static const std::size_t CHILDREN = RN::CHILDREN;
- static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)>>::value, "" );
+ static_assert(Rows_v<expr_value_type> == CHILDREN && Cols_v<expr_value_type> == CHILDREN, "" );
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
@@ -170,7 +170,7 @@ namespace AMDiS
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
- const Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)> exprValue = Super::coefficient(local);
+ const auto exprValue = Super::coefficient(local);
auto const& rowShapeValues = rowShapeValuesCache[iq];
auto const& colShapeValues = colShapeValuesCache[iq];
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 51040f8f464d1f90f8e52102f184cbb32e14fc00..51f2e08cf2cf30108b0f09f1f7da182be0b23317 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -49,6 +49,9 @@ dune_add_test(SOURCES OperationsTest.cpp
dune_add_test(SOURCES OperatorsTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES ProblemStatTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES RangeTypeTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ad40cb96c064ebd418b80b5aafe3f47102763ce
--- /dev/null
+++ b/test/ProblemStatTest.cpp
@@ -0,0 +1,36 @@
+#include <amdis/AMDiS.hpp>
+#include <amdis/LocalOperators.hpp>
+#include <amdis/ProblemStat.hpp>
+
+using namespace AMDiS;
+
+using Grid = Dune::YaspGrid<2>;
+
+template <class T>
+struct Param
+ : DefaultBasisCreator<Grid, Impl::LagrangePreBasisCreator<1,1>, T>
+{};
+
+template <class T>
+void test()
+{
+ ProblemStat<Param<T>> prob("ellipt");
+ prob.initialize(INIT_ALL);
+ prob.boundaryManager()->setBoxBoundary({1,1,2,2});
+
+ prob.addMatrixOperator(sot(T(1)), 0, 0);
+ prob.addVectorOperator(zot(T(1)), 0);
+ prob.addDirichletBC(BoundaryType{1}, 0,0, T(0));
+}
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ test<float>();
+ test<double>();
+ test<long double>();
+
+ AMDiS::finalize();
+ return 0;
+}