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;
+}