diff --git a/dune.module b/dune.module
index 123a1d68389b2070a8bc1791685b84a2da1045e7..0a19748dfd64b9b1c23ea8dff5c83d155df460b2 100644
--- a/dune.module
+++ b/dune.module
@@ -7,5 +7,5 @@ Module: amdis
 Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on
-Depends: dune-common dune-geometry dune-localfunctions dune-typetree dune-grid dune-functions
+Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-localfunctions (>= 2.6) dune-typetree (>= 2.6) dune-grid (>= 2.6) dune-functions (>= 2.6)
 #Suggests: dune-uggrid dune-alugrid dune-foamgrid
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 5b620582a0ced2f7488db986b575f613e2b7b78a..67b30e6f1d594455e1cdce1329d4ae9f1a95fbf7 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -36,15 +36,11 @@ int main(int argc, char** argv)
     return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
   };
 
-  using Grid = Dune::YaspGrid<AMDIS_DIM>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
-  auto basis = makeBasis(gridPtr->leafGridView(), lagrange<1>());
-
-  auto prob = makeProblemStat("ellipt", *gridPtr, basis);
+  ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
-  prob.addMatrixOperator(opL);
+  prob.addMatrixOperator(opL, _0, _0);
 
   auto opForce = makeOperator(tag::test{}, f, 6);
   prob.addVectorOperator(opForce, _0);
@@ -60,7 +56,7 @@ int main(int argc, char** argv)
   std::vector<double> widths; widths.reserve(numLevels);
   for (int i = 0; i < numLevels; ++i) {
     prob.getGrid()->globalRefine(1);
-    auto gridView = prob.getGlobalBasis()->gridView();
+    auto gridView = prob.gridView();
 
     double h = 0;
     for (auto const& e : edges(gridView))
@@ -76,7 +72,7 @@ int main(int argc, char** argv)
     errH1.push_back(std::sqrt(errorH1));
 
 #if WRITE_FILES
-    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
+    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
     vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
     vtkWriter.write("u_" + std::to_string(i));
 #endif
diff --git a/examples/heat.cc b/examples/heat.cc
index f09f7505c41c20ad7a50720ae0a753bd8779ed05..d232264159f3ca35290fbaf364d9b3c06b032400 100644
--- a/examples/heat.cc
+++ b/examples/heat.cc
@@ -13,44 +13,44 @@
 
 using namespace AMDiS;
 
+// 1 component with polynomial degree 1
+//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
+using HeatParam   = YaspGridBasis<AMDIS_DIM, 2>;
+using HeatProblem = ProblemStat<HeatParam>;
+using HeatProblemInstat = ProblemInstat<HeatParam>;
+
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-  using Grid = Dune::YaspGrid<AMDIS_DIM>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
-  auto basis = makeBasis(gridPtr->leafGridView(), lagrange<2>());
-
-  auto prob = makeProblemStat("heat", *gridPtr, basis);
+  HeatProblem prob("heat");
   prob.initialize(INIT_ALL);
 
-  auto probInstat = makeProblemInstat("heat", prob);
+  HeatProblemInstat probInstat("heat", prob);
   probInstat.initialize(INIT_UH_OLD);
 
   AdaptInfo adaptInfo("adapt");
 
   double* invTau = probInstat.getInvTau();
 
-  auto _ = Dune::TypeTree::hybridTreePath();
-
   auto opTimeLhs = makeOperator(tag::test_trial{}, std::ref(*invTau));
-  prob.addMatrixOperator(opTimeLhs, _, _);
+  prob.addMatrixOperator(opTimeLhs, 0, 0);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
-  prob.addMatrixOperator(opL, _, _);
+  prob.addMatrixOperator(opL, 0, 0);
 
   auto opTimeRhs = makeOperator(tag::test{},
-    invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution()), 2);
-  prob.addVectorOperator(opTimeRhs, _);
+    invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution(0)), 2);
+  prob.addVectorOperator(opTimeRhs, 0);
 
   auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
-  prob.addVectorOperator(opForce, _);
+  prob.addVectorOperator(opForce, 0);
 
 
   // set boundary condition
   auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
   auto dbcValues = [](auto const& p){ return 0.0; };
-  prob.addDirichletBC(predicate, _, _, dbcValues);
+  prob.addDirichletBC(predicate, 0, 0, dbcValues);
 
   AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
   adapt.adapt();
diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 4ac468a07641690d09fa85066773107758b99316..4fcf941ba26aa0d16f10dc0ef26832dc01a881d7 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -12,20 +12,23 @@
 
 using namespace AMDiS;
 
-int main(int argc, char** argv)
+struct NavierStokesBasis
 {
-  AMDiS::init(argc, argv);
-
   using Grid = Dune::YaspGrid<AMDIS_DIM>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
+  using GlobalBasis = typename TaylorHoodBasis<Grid::LeafGridView>::GlobalBasis;
+};
+
+using StokesProblem = ProblemStat<NavierStokesBasis>;
+using StokesProblemInstat = ProblemInstat<NavierStokesBasis>;
 
-  auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
-  auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
 
-  auto prob = makeProblemStat("stokes", *gridPtr, basis);
+  StokesProblem prob("stokes");
   prob.initialize(INIT_ALL);
 
-  auto probInstat = makeProblemInstat("stokes", prob);
+  StokesProblemInstat probInstat("stokes", prob);
   probInstat.initialize(INIT_UH_OLD);
 
   double viscosity = 1.0;
diff --git a/examples/stokes0.cc b/examples/stokes0.cc
index 37588d94ca729cf1f39664d204dc73cb9c7101d2..fbcc5fb80bd5d7bb703046cd30a48bbc8c933d34 100644
--- a/examples/stokes0.cc
+++ b/examples/stokes0.cc
@@ -13,16 +13,17 @@
 
 using namespace AMDiS;
 
+// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
+
+using Grid = Dune::YaspGrid<AMDIS_DIM>;
+using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
+using StokesProblem = ProblemStat<StokesParam>;
+
 int main(int argc, char** argv)
 {
     AMDiS::init(argc, argv);
 
-    using Grid = Dune::YaspGrid<AMDIS_DIM>;
-    auto gridPtr = MeshCreator<Grid>::create("mesh");
-    auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
-    auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
-
-    auto prob = makeProblemStat("stokes", *gridPtr, basis);
+    StokesProblem prob("stokes");
     prob.initialize(INIT_ALL);
 
     double viscosity = 1.0;
diff --git a/examples/stokes1.cc b/examples/stokes1.cc
index cee81273568942ad1f2b18ff7734d08ee7f2c87d..cd62b4cdb60d9c1598bfd74c2ba8484f630a1cb7 100644
--- a/examples/stokes1.cc
+++ b/examples/stokes1.cc
@@ -13,17 +13,17 @@
 
 using namespace AMDiS;
 
+// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
+
+using Grid = Dune::YaspGrid<AMDIS_DIM>;
+using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
+using StokesProblem = ProblemStat<StokesParam>;
+
 int main(int argc, char** argv)
 {
     AMDiS::init(argc, argv);
 
-    using Grid = Dune::YaspGrid<AMDIS_DIM>;
-    auto gridPtr = MeshCreator<Grid>::create("mesh");
-
-    auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
-    auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
-
-    auto prob = makeProblemStat("stokes", *gridPtr, basis);
+    StokesProblem prob("stokes");
     prob.initialize(INIT_ALL);
 
     double viscosity = 1.0;
diff --git a/examples/stokes3.cc b/examples/stokes3.cc
index f7b29e875722f022e0a7e2c655b9724d3c8b48e3..54c75d0da556897759920a782f3cf0b029875902 100644
--- a/examples/stokes3.cc
+++ b/examples/stokes3.cc
@@ -9,17 +9,15 @@
 
 using namespace AMDiS;
 
+using Grid = Dune::YaspGrid<AMDIS_DIM>;
+using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
+using StokesProblem = ProblemStat<StokesParam>;
+
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-  using Grid = Dune::YaspGrid<AMDIS_DIM>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
-
-  auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
-  auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
-
-  auto prob = makeProblemStat("stokes", *gridPtr, basis);
+  StokesProblem prob("stokes");
   prob.initialize(INIT_ALL);
 
   double viscosity = 1.0;
@@ -40,12 +38,12 @@ int main(int argc, char** argv)
   auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
 
   // define boundary values
-  auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,2>
+  auto parabolic_y = [](auto const& x) -> FieldVector<double,2>
   {
     return {0.0, x[1]*(1.0 - x[1])};
   };
 
-  auto zero = [](auto const& x) -> Dune::FieldVector<double,2>
+  auto zero = [](auto const& x) -> FieldVector<double,2>
   {
     return {0.0, 0.0};
   };
diff --git a/examples/vecellipt.cc b/examples/vecellipt.cc
index 78aaa0e383dfbbcbc0d07d0864e50632cbe114bb..8ea8636049a3f5c7a91be55b85f310d0996bb25c 100644
--- a/examples/vecellipt.cc
+++ b/examples/vecellipt.cc
@@ -6,17 +6,16 @@
 
 using namespace AMDiS;
 
+// 1 component with polynomial degree 1
+//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
+using ElliptParam   = YaspGridBasis<AMDIS_DIM, 2,2>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-  using Grid = Dune::YaspGrid<AMDIS_DIM>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
-
-  auto preBasis = repeat<2>(lagrange<2>());
-  auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
-
-  auto prob = makeProblemStat("stokes", *gridPtr, basis);
+  ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
   AdaptInfo adaptInfo("adapt");
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index 7afe0dfba60e8276bc73d413cc5ee9ad9799fa36..5ab4ed185154626e52901f22613b2e4409fa8144 100644
--- a/src/amdis/Assembler.inc.hpp
+++ b/src/amdis/Assembler.inc.hpp
@@ -37,7 +37,7 @@ void Assembler<Traits>::assemble(
     auto geometry = element.geometry();
 
     // traverse type-tree of global-basis
-    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
+    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
     {
       auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
       auto rowLocalView = rowBasis.localView();
@@ -55,7 +55,7 @@ void Assembler<Traits>::assemble(
         this->assembleElementOperators(element, rhsOp, vecAssembler);
       }
 
-      AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
+      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
       {
         auto& matOp = matrixOperators_[rowNode][colNode];
         if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
@@ -95,9 +95,9 @@ void Assembler<Traits>::assemble(
     }
 
     // unbind all operators
-    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
+    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
       rhsOperators_[rowNode].unbind();
-      AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
+      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
         matrixOperators_[rowNode][colNode].unbind();
       });
     });
@@ -151,7 +151,7 @@ void Assembler<Traits>::initMatrixVector(
     rhs = 0;
 
   auto localView = globalBasis_.localView();
-  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
+  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
   {
 #ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
     if (rowNode.isLeaf)
@@ -160,7 +160,7 @@ void Assembler<Traits>::initMatrixVector(
 
     auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
 
-    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
+    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
     {
       auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
 
@@ -183,14 +183,14 @@ std::size_t Assembler<Traits>::finishMatrixVector(
   matrix.finish();
 
   auto localView = globalBasis_.localView();
-  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
+  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
   {
     auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
     auto& rhsOp = rhsOperators_[rowNode];
     if (rhsOp.doAssemble(asmVector))
       rhsOp.assembled = true;
 
-    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
+    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
     {
       auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
       auto& matOp = matrixOperators_[rowNode][colNode];
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index d456189e62747eff3bad6265f70af886c41c4e73..34a281e1c992a40192133583403bac3484a1fa91 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -5,9 +5,13 @@
 #include <type_traits>
 #include <vector>
 
+#include <dune/common/hybridutilities.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+
 #include <amdis/Output.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/ValueCategory.hpp>
+#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
 #include <amdis/utility/RangeType.hpp>
 #include <amdis/utility/TreeData.hpp>
 
@@ -75,7 +79,18 @@ namespace AMDiS
       test_exit_dbg(initialized_, "Boundary condition not initialized!");
       auto columns = matrix.applyDirichletBC(dirichletNodes_);
 
-      finishImpl(matrix, rhs, solution, rowBasis, colBasis, ValueCategory_t<Range>{});
+      using Dune::Functions::interpolate;
+      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
+        [&](auto id) {
+          auto rhsWrapper = wrapper(rhs.getVector());
+          interpolate(id(rowBasis), rhsWrapper, values_, dirichletNodes_);
+        });
+
+      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
+        [&](auto id) {
+          auto solutionWrapper = wrapper(solution.getVector());
+          interpolate(id(colBasis), solutionWrapper, values_, dirichletNodes_);
+        });
 
       // subtract columns of dirichlet nodes from rhs
       // for (auto const& triplet : columns)
@@ -96,16 +111,6 @@ namespace AMDiS
     template <class RB, class RowNodeTag>
     void initImpl(RB const&, RowNodeTag) {}
 
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
-    void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
-                    RowBasis const& rowBasis, ColBasis const& colBasis,
-                    ValueCat);
-
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
-    void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
-                    RowBasis const& rowBasis, ColBasis const& colBasis,
-                    tag::unknown) {}
-
   private:
     std::function<bool(Domain)> predicate_;
     std::function<Range(Domain)> values_;
diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp
index 1c50ea3abdec93202731875aa72f2757b213778f..4fcaf888fac8e546355d5a6b59bcb21fe1643cc6 100644
--- a/src/amdis/DirichletBC.inc.hpp
+++ b/src/amdis/DirichletBC.inc.hpp
@@ -4,6 +4,7 @@
 #include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <amdis/LinearAlgebra.hpp>
 #include <amdis/linear_algebra/HierarchicWrapper.hpp>
+#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
 
 namespace AMDiS
 {
@@ -13,7 +14,7 @@ namespace AMDiS
       RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
   {
     using Dune::Functions::interpolate;
-    interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
+    interpolate(rowBasis, dirichletNodes_, predicate_);
   }
 
   template <class WorldVector, class Range>
@@ -46,19 +47,4 @@ namespace AMDiS
     });
   }
 
-
-  template <class WorldVector, class Range>
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
-  void DirichletBC<WorldVector, Range>::finishImpl(
-      Matrix& matrix, VectorX& solution, VectorB& rhs,
-      RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
-  {
-    using Dune::Functions::interpolate;
-    interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values_,
-      hierarchicVectorWrapper(dirichletNodes_));
-    interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values_,
-      hierarchicVectorWrapper(dirichletNodes_));
-  }
-
-
 } // end namespace AMDiS
diff --git a/src/amdis/ProblemInstat.inc.hpp b/src/amdis/ProblemInstat.inc.hpp
index 4c268f5a52181be316539b4959d25ced87acd3c5..9005f11fb067bbb2c1c40a233f2adc47e99aab1e 100644
--- a/src/amdis/ProblemInstat.inc.hpp
+++ b/src/amdis/ProblemInstat.inc.hpp
@@ -45,7 +45,7 @@ void ProblemInstat<Traits>::createUhOld()
   if (oldSolution)
     warning("oldSolution already created\n");
   else // create oldSolution
-    oldSolution.reset(new SystemVector(*problemStat.getGlobalBasis(), name + "_uOld"));
+    oldSolution.reset(new SystemVector(*problemStat.globalBasis(), name + "_uOld"));
 }
 
 
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index fc5f60480ec54c35ca6de6b2593d88035c3bf7c2..9fd6901acebb699597d612c0c235117b5346cc3a 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -67,7 +67,6 @@ namespace AMDiS
     using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
 
   public:
-#if 0
     /**
      * \brief Constructor. Takes the name of the problem that is used to
      * access values corresponding to this problem in the parameter file.
@@ -94,7 +93,6 @@ namespace AMDiS
       globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
       initGlobalBasis(*globalBasis_);
     }
-#endif
 
     /// \brief Constructor taking a grid reference and a basis reference.
     /// Stores pointers to both.
@@ -211,7 +209,6 @@ namespace AMDiS
     /// Return a pointer to the grid, \ref grid
     std::shared_ptr<Grid> getGrid() { return grid_; }
 
-#if 0
     /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
     /// matrices and vectors, as well as the file-writer.
     void setGrid(Grid& grid)
@@ -222,13 +219,12 @@ namespace AMDiS
       createMatricesAndVectors();
       createFileWriter();
     }
-#endif
 
     /// Return the gridView of the leaf-level
     auto const& gridView() { return globalBasis_->gridView(); }
 
     /// Return the \ref feSpaces
-    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
+    std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
 
 
     /// Implementation of \ref ProblemStatBase::getName
@@ -239,6 +235,27 @@ namespace AMDiS
 
   protected: // initialization methods
 
+    template <class T, class GV>
+    using HasCreate = decltype(T::create(std::declval<GV>()));
+
+    void createGlobalBasis()
+    {
+      createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
+      initGlobalBasis(*globalBasis_);
+    }
+
+    void createGlobalBasis(std::true_type)
+    {
+      assert( bool(grid) );
+      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
+      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid->leafGridView()));
+    }
+
+    void createGlobalBasis(std::false_type)
+    {
+      error_exit("Can not create GlobalBasis from type. Pass a BasisCreator instead!");
+    }
+
     void createGrid()
     {
       gridName_ = "";
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 1db46d2b60ae4506318e3fe79e7d50a0d66eb706..92a7b2155a095f204eea5f93351a73925bec6fc7 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -57,12 +57,10 @@ void ProblemStat<Traits>::initialize(
     warning("globalBasis already created");
   }
   else {
-#if 0
     if (initFlag.isSet(INIT_FE_SPACE) ||
         (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
       createGlobalBasis();
     }
-#endif
 
     if (adoptProblem &&
         (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index 6dcb9e170483293871c98084c5f43641c43487b0..5d6f05e109fb715658987f4ea0d1f856e3074560 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -1,114 +1,87 @@
 #pragma once
 
-#include <tuple>
-
 #include <dune/functions/functionspacebases/basistags.hh>
 #include <dune/functions/functionspacebases/compositebasis.hh>
-#include <dune/functions/functionspacebases/flatmultiindex.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/powerbasis.hh>
 #include <dune/functions/functionspacebases/pqknodalbasis.hh>
 #include <dune/grid/yaspgrid.hh>
 
-namespace Dune {
-namespace Functions {
-namespace BasisFactory {
-
-  // define the basis-build (factory-tag) for a given pre-basis
-  template <class PreBasis>
-  struct FactoryTag;
-
-  template <class MultiIndex, class IndexMergingStrategy, class... ChildPreBasisFactory>
-  struct FactoryTag<CompositePreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasisFactory...>>
-  {
-    using type = Imp::CompositePreBasisFactory<IndexMergingStrategy, typename FactoryTag<ChildPreBasisFactory>::type...>;
-  };
-
-  template <class MultiIndex, class IndexMergingStrategy, class ChildPreBasisFactory, std::size_t C>
-  struct FactoryTag<PowerPreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasisFactory, C>>
-  {
-    using type = Imp::PowerPreBasisFactory<C, IndexMergingStrategy, typename FactoryTag<ChildPreBasisFactory>::type>;
-  };
-
-  template <class GridView, int k, class MultiIndex>
-  struct FactoryTag<PQkPreBasis<GridView, k, MultiIndex>>
-  {
-    using type = Imp::PQkPreBasisFactory<k>;
-  };
-
-}}} // end namespace Dune::Functions::BasisFactory
+#include <amdis/common/Mpl.hpp>
 
 namespace AMDiS
 {
   namespace Impl
   {
-    template <class GridView, bool same, int... deg>
-    struct LagrangeBasisCreatorImpl;
+    template <bool same, int... degs>
+    struct LagrangePreBasisCreatorImpl;
 
     // specialization for composite basis
-    template <class GridView, int... deg>
-    struct LagrangeBasisCreatorImpl<GridView, false, deg...>
+    template <int... degs>
+    struct LagrangePreBasisCreatorImpl<false, degs...>
     {
-      using MultiIndex = Dune::ReservedVector<std::size_t,1>;
-      using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
-
-      using type = Dune::Functions::CompositePreBasis<MultiIndex, IndexTag,
-        Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>...>;
+      static auto create()
+      {
+        using namespace Dune::Functions::BasisFactory;
+        return composite(lagrange<degs>()..., flatLexicographic());
+      }
     };
 
     // specialization for power basis
-    template <class GridView, int deg, int... degs>
-    struct LagrangeBasisCreatorImpl<GridView, true, deg, degs...>
+    template <int deg, int... degs>
+    struct LagrangePreBasisCreatorImpl<true, deg, degs...>
     {
-      using MultiIndex = Dune::ReservedVector<std::size_t,1>;
-      using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
-
-      using type = Dune::Functions::PowerPreBasis<MultiIndex, IndexTag,
-        Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>, 1 + sizeof...(degs)>;
+      static auto create()
+      {
+        using namespace Dune::Functions::BasisFactory;
+        return power<1+sizeof...(degs)>(lagrange<deg>(), flatLexicographic());
+      }
     };
 
     // factory to construct a global basis of several lagrange bases, with flat indexing.
-    template <class GridView, int deg, int... degs>
-    struct LagrangeBasisCreator
-    {
-      using type = typename LagrangeBasisCreatorImpl<GridView, and_t<(deg == degs)...>::value, deg, degs...>::type;
-    };
+    template <int deg, int... degs>
+    struct LagrangePreBasisCreator
+        : public LagrangePreBasisCreatorImpl<all_of_v<(deg == degs)...>, deg, degs...>
+    {};
 
 
-    template <class GridView, int k>
-    struct TaylorHoodBasisCreator
+    template <int dow, int k = 1>
+    struct TaylorHoodPreBasisCreator
     {
-      using MultiIndex = Dune::ReservedVector<std::size_t,1>;
-      using IndexTagV = Dune::Functions::BasisFactory::FlatInterleaved;
-
-      using VelocityPreBasis = Dune::Functions::PowerPreBasis<MultiIndex, IndexTagV,
-        Dune::Functions::PQkPreBasis<GridView,k+1,MultiIndex>, GridView::dimensionworld>;
-
-      using PressurePreBasis = Dune::Functions::PQkPreBasis<GridView,k,MultiIndex>;
-
-      using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
-      using type = Dune::Functions::CompositePreBasis<MultiIndex, IndexTag,
-        VelocityPreBasis, PressurePreBasis>;
+      static auto create()
+      {
+        using namespace Dune::Functions::BasisFactory;
+        return composite(power<dow>(lagrange<k+1>(), flatInterleaved()), lagrange<k>(), flatLexicographic());
+      }
     };
 
   } // end namespace Impl
 
   /// An Exemplar for ProblemStatTraits
   template <class GlobalBasisType>
-  class DefaultProblemTraits
+  struct DefaultProblemTraits
   {
-  public:
-    /// A global-basis of type dune-functions DefaultGlobalBasis<NodeFactory>
     using GlobalBasis = GlobalBasisType;
   };
 
+
+  template <class GridView, class PreBasisCreator>
+  struct DefaultBasisCreator
+  {
+    static auto create(GridView const& gridView)
+    {
+      using namespace Dune::Functions::BasisFactory;
+      return makeBasis(gridView, PreBasisCreator::create());
+    }
+
+    using GlobalBasis = decltype(create(std::declval<GridView>()));
+  };
+
   /// \brief ProblemStatTraits for a (composite) basis composed of
   /// lagrange bases of different degree.
   template <class GridView, int... degrees>
   struct LagrangeBasis
-      : public DefaultProblemTraits<
-          Dune::Functions::DefaultGlobalBasis<
-            typename Impl::LagrangeBasisCreator<GridView, degrees...>::type> >
+      : public DefaultBasisCreator<GridView, Impl::LagrangePreBasisCreator<degrees...>>
   {};
 
   /// \brief Specialization of \ref LagrangeBasis for Grid type
@@ -122,29 +95,7 @@ namespace AMDiS
   /// with pressure degree k
   template <class GridView, int k = 1>
   struct TaylorHoodBasis
-      : public DefaultProblemTraits<
-          Dune::Functions::DefaultGlobalBasis<
-            typename Impl::TaylorHoodBasisCreator<GridView,k>::type> >
+      : public DefaultBasisCreator<GridView, Impl::TaylorHoodPreBasisCreator<GridView::dimensionworld,k>>
   {};
 
-  using Dune::Functions::BasisFactory::makeBasis;
-  using Dune::Functions::BasisFactory::lagrange;
-  using Dune::Functions::BasisFactory::pq;
-
-  template <std::size_t k, class ChildPBF>
-  auto repeat(ChildPBF&& childPreBasisFactory)
-  {
-    return Dune::Functions::BasisFactory::power<k>(
-      std::forward<ChildPBF>(childPreBasisFactory),
-      Dune::Functions::BasisFactory::flatInterleaved());
-  }
-
-  template <class... Args>
-  auto compose(Args&&... args)
-  {
-    return Dune::Functions::BasisFactory::composite(
-      std::forward<Args>(args)...,
-      Dune::Functions::BasisFactory::flatLexicographic());
-  }
-
 } // end namespace AMDiS
diff --git a/src/amdis/gridfunctions/DOFVectorView.inc.hpp b/src/amdis/gridfunctions/DOFVectorView.inc.hpp
index 5d1222632a6e5c3eb45905850996543d02aaf0aa..a97dc1a47b1922b0db421f80b19468858a4f2076 100644
--- a/src/amdis/gridfunctions/DOFVectorView.inc.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.inc.hpp
@@ -13,7 +13,7 @@ LocalFunction::operator()(LocalDomain const& x) const
 
   auto&& coefficients = *globalFunction_->dofVector_;
   auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
-  AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto)
+  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
   {
     auto&& fe = node.finiteElement();
     auto&& localBasis = fe.localBasis();
@@ -22,7 +22,7 @@ LocalFunction::operator()(LocalDomain const& x) const
     localBasis.evaluateFunction(x, shapeFunctionValues);
 
     // Get range entry associated to this node
-    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, y));
+    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
 
     for (std::size_t i = 0; i < localBasis.size(); ++i) {
       auto&& multiIndex = localView_.index(node.localIndex(i));
@@ -59,7 +59,7 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
 
   auto&& coefficients = *globalFunction_->dofVector_;
   auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
-  AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto)
+  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
   {
     // TODO: may DOFVectorView::Range to FieldVector type if necessary
     using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
@@ -80,7 +80,7 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
       multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
 
     // Get range entry associated to this node
-    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, dy));
+    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
 
     for (std::size_t i = 0; i < localBasis.size(); ++i) {
       auto&& multiIndex = localView_.index(node.localIndex(i));
diff --git a/src/amdis/linear_algebra/HierarchicWrapper.hpp b/src/amdis/linear_algebra/HierarchicWrapper.hpp
index 9da189724f5c80806897008e2a8d233e7935c155..94617e6de355bc2fbd9b7bb8ea3953fad3abb998 100644
--- a/src/amdis/linear_algebra/HierarchicWrapper.hpp
+++ b/src/amdis/linear_algebra/HierarchicWrapper.hpp
@@ -64,6 +64,11 @@ namespace AMDiS
       resizeImpl(vector_, s);
     }
 
+    std::size_t size() const
+    {
+      return sizeImpl(vector_);
+    }
+
   protected:
 
     template <class Vec, class Basis>
@@ -80,6 +85,22 @@ namespace AMDiS
       vec.change_dim(std::size_t(s));
     }
 
+    template <class... Params>
+    std::size_t sizeImpl(mtl::dense_vector<Params...> const& v) const
+    {
+      return mtl::size(v);
+    }
+
+    template <class V>
+    using HasSizeMember = decltype(std::declval<V>().size());
+
+    template <class V,
+      std::enable_if_t<Dune::Std::is_detected<HasSizeMember, V>::value, int> = 0>
+    std::size_t sizeImpl(V const& v) const
+    {
+      return v.size();
+    }
+
   private:
     Vector& vector_;
   };
diff --git a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp b/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
index 87dd50a933ca6c79625b875498304e5738bf044e..4e08e4698dca62edaafda5c74647266decfc4911 100644
--- a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
+++ b/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
@@ -23,7 +23,7 @@ namespace AMDiS
       using size_type = typename Super::size_type;
 
       template <class Vector_>
-      explicit MTLWrapper(Vector_&& vec)
+      MTLWrapper(Vector_&& vec)
         : Super(std::forward<Vector_>(vec))
       {}
 
diff --git a/src/amdis/utility/Traversal.hpp b/src/amdis/utility/Traversal.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3a02e7422776912b5cc025741a5375111ecad5af
--- /dev/null
+++ b/src/amdis/utility/Traversal.hpp
@@ -0,0 +1,137 @@
+#pragma once
+
+#include <dune/typetree/nodetags.hh>
+#include <dune/typetree/treepath.hh>
+#include <dune/typetree/visitor.hh>
+
+#include <amdis/common/Loops.hpp>
+
+namespace AMDiS
+{
+  // forward declaration of main engine struct
+  template <typename NodeTag, bool visit = true>
+  struct TraverseTree;
+
+
+  // Do not visit nodes the visitor is not interested in
+  template <typename NodeTag>
+  struct TraverseTree<NodeTag, false>
+  {
+    template <typename Node, typename Visitor, typename TreePath>
+    static void apply(const Node& node, const Visitor& visitor, TreePath const& tp)
+    {}
+  };
+
+#ifndef DOXYGEN
+
+  // some implementation details
+
+  template <class Node, class Index>
+  struct HybridChildType
+      : HybridChildType<std::remove_const_t<Node>, std::remove_const_t<Index>> {};
+
+  template <class Node>
+  struct HybridChildType<Node, std::size_t>
+  {
+    using type = typename Node::template Child<0>::Type;
+  };
+
+  template <class Node, std::size_t K>
+  struct HybridChildType<Node, Dune::index_constant<K>>
+  {
+    using type = typename Node::template Child<K>::Type;
+  };
+
+  template <class NodeTag, class Node>
+  constexpr std::size_t hybridDegree(NodeTag, Node const& node)
+  {
+    return Dune::TypeTree::degree(node);
+  }
+
+  template <class Node>
+  constexpr auto hybridDegree(Dune::TypeTree::CompositeNodeTag, Node const& node)
+  {
+    return Dune::index_constant<Node::CHILDREN>{};
+  }
+
+
+  template <std::size_t k, std::size_t n>
+  constexpr bool notLastChild(Dune::index_constant<k> const&, Dune::index_constant<n> const&)
+  {
+    return k < n-1;
+  }
+
+  constexpr bool notLastChild(std::size_t k, std::size_t n)
+  {
+    return k < n-1;
+  }
+
+#endif
+
+
+  template <class NodeTag>
+  struct TraverseTree<NodeTag, true>
+  {
+    template <typename N, typename V, typename TreePath>
+    static void apply(N&& n, V&& v, TreePath const& tp)
+    {
+      using Node = std::remove_reference_t<N>;
+      using Visitor = std::remove_reference_t<V>;
+
+      v.pre(std::forward<N>(n),tp);
+
+      auto const deg = hybridDegree(NodeTag{}, n);
+      forEach(Dune::range(deg), [&](auto const _k)
+      {
+        // always call beforeChild(), regardless of the value of visit
+        v.beforeChild(std::forward<N>(n),n.child(_k),tp,_k);
+
+        // descend to child
+        using C = typename HybridChildType<Node, decltype(_k)>::type;
+        const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value;
+        TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(_k),std::forward<V>(v),push_back(tp, _k));
+
+        // always call afterChild(), regardless of the value of visit
+        v.afterChild(std::forward<N>(n),n.child(_k),tp,_k);
+
+        // if this is not the last child, call infix callback
+        if (notLastChild(_k, deg))
+          v.in(std::forward<N>(n),tp);
+      });
+
+      v.post(std::forward<N>(n),tp);
+    }
+  };
+
+  // LeafNode - just call the leaf() callback
+  template <>
+  struct TraverseTree<Dune::TypeTree::LeafNodeTag, true>
+  {
+    template <typename N, typename V, typename TreePath>
+    static void apply(N&& n, V&& v, TreePath const& tp)
+    {
+      v.leaf(std::forward<N>(n),tp);
+    }
+  };
+
+  //! Apply visitor to TypeTree.
+  /**
+    * This function applies the given visitor to the given tree. Both visitor and tree may be const
+    * or non-const (if the compiler supports rvalue references, they may even be a non-const temporary).
+    *
+    * \note The visitor must implement the interface laid out by DefaultVisitor (most easily achieved by
+    *       inheriting from it).
+    *
+    * \param tree    The tree the visitor will be applied to.
+    * \param visitor The visitor to apply to the tree.
+    */
+  template <typename Tree, typename Visitor>
+  void traverseTree(Tree&& tree, Visitor&& visitor)
+  {
+    using Node = std::remove_reference_t<Tree>;
+    using NodeTag = Dune::TypeTree::NodeTag<Node>;
+    using TreePath = Dune::TypeTree::HybridTreePath<>;
+    TraverseTree<NodeTag>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor), TreePath{});
+  }
+
+} // end namespace AMDiS
diff --git a/src/amdis/utility/TreeData.hpp b/src/amdis/utility/TreeData.hpp
index b913ad464a80ffb60b4181ea27fa27ba35a7d399..9325cf1398d50431a3ebb8301cc4a8bd704f4963 100644
--- a/src/amdis/utility/TreeData.hpp
+++ b/src/amdis/utility/TreeData.hpp
@@ -5,8 +5,9 @@
 #include <utility>
 #include <vector>
 
-#include <dune/typetree/traversal.hh>
+//#include <dune/typetree/traversal.hh>
 #include <dune/typetree/typetree.hh>
+#include <amdis/utility/Visitor.hpp>
 
 namespace AMDiS
 {
@@ -207,10 +208,10 @@ namespace AMDiS
     }
 
     template <class Func>
-    void applyImpl(Func&& func, std::true_type) { AMDiS::forEachLeafNode_(*tree_, func); }
+    void applyImpl(Func&& func, std::true_type) { forEachLeafNode_(*tree_, func); }
 
     template <class Func>
-    void applyImpl(Func&& func, std::false_type) { AMDiS::forEachNode_(*tree_, func); }
+    void applyImpl(Func&& func, std::false_type) { forEachNode_(*tree_, func); }
 
   protected:
     const Tree* tree_ = nullptr;
@@ -251,7 +252,11 @@ namespace AMDiS
     void init(Tree const& tree, tag::store)
     {
       Super::init(tree, tag::store{});
+<<<<<<< c786927261d1c07223ce3f9feabd89e439f2f4bb
       AMDiS::forEachNode_(tree, [&,this](auto const& node, auto&&)
+=======
+      forEachNode_(tree, [&,this](auto const& node, auto&&)
+>>>>>>> compatibility with dune 2.6.0 of all dependend modules
       {
         (*this)[node].init(tree, tag::store{});
       });
diff --git a/test/ExpressionsTest.cpp b/test/ExpressionsTest.cpp
index 8c01c202972f1d04ec429c03f9067526c7f5be46..bc544975c96726fd944e82b33202119f31dcf1d9 100644
--- a/test/ExpressionsTest.cpp
+++ b/test/ExpressionsTest.cpp
@@ -13,20 +13,19 @@
 
 using namespace AMDiS;
 
+using ElliptParam   = YaspGridBasis<2, 2>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
   using namespace Dune::Indices;
 
-  using Grid = Dune::YaspGrid<2>;
-  auto gridPtr = MeshCreator<Grid>::create("mesh");
-  auto basis = makeBasis(gridPtr->leafGridView(), lagrange<2>());
-
-  auto prob = makeProblemStat("ellipt", *gridPtr, basis);
+  ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
-  auto u = prob.getSolution();
+  auto u = prob.getSolution(_0);
 
   // eval a functor at global coordinates (at quadrature points)
   auto expr1 = evalAtQP([](Dune::FieldVector<double, 2> const& x) { return x[0] + x[1]; });
@@ -43,7 +42,7 @@ int main(int argc, char** argv)
   auto expr8 = X(0);
 
   // Evaluation of the DOFVector (component) at quadrature points
-  auto expr9 = prob.getSolution();
+  auto expr9 = prob.getSolution(_0);
 
   // ---------------------------------