diff --git a/CMakeLists.txt b/CMakeLists.txt
index eeeae340262a5440c90e67568502f6b1c307a647..0d198f4bed15a57a500683d7e693f9e5987035bf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,8 @@ if(NOT (dune-common_DIR OR dune-common_ROOT OR
       ${PROJECT_BINARY_DIR})
 endif()
 
+set(CXX_MAX_STANDARD 14)
+
 #find dune-common and set the module path
 find_package(dune-common REQUIRED)
 list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
diff --git a/dune/amdis/Assembler.hpp b/dune/amdis/Assembler.hpp
index 684c053c36646be488e8e7d82bf499d475b1e2b1..5cd91279552bb2e32576571b9c57bc99173bd3a9 100644
--- a/dune/amdis/Assembler.hpp
+++ b/dune/amdis/Assembler.hpp
@@ -6,139 +6,89 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
-#include <dune/amdis/FiniteElementSpaces.hpp>
 #include <dune/amdis/LinearAlgebra.hpp>
 #include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/common/TypeDefs.hpp>
 
 namespace AMDiS
 {
-  template <class FeSpaces>
+  template <class GlobalBasis>
   class Assembler
   {
-    template <std::size_t I>
-    using FeSpace = std::tuple_element_t<I, FeSpaces>;
-
-    /// Number of problem components
-    static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
-
     /// The grid view the global FE basis lives on
-    using GridView = typename FeSpace<0>::GridView;
-
-    using SystemMatrixType = SystemMatrix<FeSpaces>;
-    using SystemVectorType = SystemVector<FeSpaces>;
-
-    template <class T>
-    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
+    using GridView = typename GlobalBasis::GridView;
 
-    template <class T>
-    using VectorEntries = Dune::FieldVector<T, nComponents>;
+    using SystemMatrixType = SystemMatrix<GlobalBasis>;
+    using SystemVectorType = SystemVector<GlobalBasis>;
 
-    using ElementVector = Impl::ElementVector;
-    using ElementMatrix = Impl::ElementMatrix;
+    using TypeTree = typename GlobalBasis::LocalView::Tree;
+    using
 
   public:
     /// Constructor, stores a shared-pointer to the feSpaces
-    Assembler(std::shared_ptr<FeSpaces> const& feSpaces, bool asmMatrix, bool asmVector)
-      : globalBases_(feSpaces)
+    Assembler(std::shared_ptr<GlobalBasis> const& globalBasis, bool asmMatrix, bool asmVector)
+      : globalBasis_(globalBasis)
       , asmMatrix_(asmMatrix)
       , asmVector_(asmVector)
     {}
 
     /// Assemble the linear system
-    template <class Operators>
+    template <class MatrixOperators, class VectorOperators>
     void assemble(
         GridView const& gv,
         SystemMatrixType& matrix,
         SystemVectorType& solution,
         SystemVectorType& rhs,
-        MatrixEntries<Operators>& matrix_operators,
-        VectorEntries<Operators>& rhs_operators);
+        MatrixOperators& matrixOperators,
+        VectorOperators& rhsOperators);
 
   private:
     /// Sets the system to zero and initializes all operators and boundary conditions
-    template <class Operators>
+    template <class MatrixOperators, class VectorOperators>
     void initMatrixVector(
         SystemMatrixType& matrix,
         SystemVectorType& solution,
         SystemVectorType& rhs,
-        MatrixEntries<Operators>& matrix_operators,
-        VectorEntries<Operators>& rhs_operators) const;
+        MatrixOperators& matrixOperators,
+        VectorOperators& rhsOperators) const;
 
 
-    /// Assemble a block-matrix of element-matrices and return a matrix of flags, whether
-    /// a given block has received some entries.
-    template <class Operators>
-    MatrixEntries<bool> assembleElementMatrices(
-        MatrixEntries<Operators>& operators,
-        MatrixEntries<ElementMatrix>& elementMatrix) const;
-
-    /// Assemble a block-vector of element-vectors and return a vector of flags, whether
-    /// a given block has received some entries.
-    template <class Operators>
-    VectorEntries<bool> assembleElementVectors(
-        VectorEntries<Operators>& operators,
-        VectorEntries<ElementVector>& elementVector) const;
+    template <class MatrixOperators, class VectorOperators>
+    void assembleElement(
+        SystemMatrixType& matrix,
+        SystemVectorType& rhs,
+        MatrixOperators& matrixOperators,
+        VectorOperators& rhsOperators) const;
 
 
-    /// Assemble one block of the block-element-matrix
-    // The MatrixData argument stores all matrix-operators
-    template <class Operators, class RowView, class ColView>
-    bool assembleElementMatrix(
+    template <class ElementContainer, class Container, class Operators, class... Bases>
+    void assembleElementOperators(
+        ElementContainer& elementContainer,
+        Container& container,
         Operators& operators,
-        ElementMatrix& elementMatrix,
-        RowView const& rowLocalView,
-        ColView const& colLocalView) const;
-
-    /// Assemble one block of the block-element-vector
-    // The VectorData argument stores all vector-operators
-    template <class Operators, class RowView>
-    bool assembleElementVector(
-        Operators& operators,
-        ElementVector& elementVector,
-        RowView const& rowLocalView) const;
-
-
-    /// Add the block-element-matrix to the system-matrix
-    void addElementMatrices(
-        SystemMatrixType& dofmatrix,
-        MatrixEntries<bool> const& addMat,
-        MatrixEntries<ElementMatrix> const& elementMatrix) const;
+        Bases const&... subBases) const;
 
-    /// Add the block-element-vector to the system-vector
-    void addElementVectors(
-        SystemVectorType& dofvector,
-        VectorEntries<bool> const& addVec,
-        VectorEntries<ElementVector> const& elementVector) const;
 
     /// Finish insertion into the matrix and assembles boundary conditions
     /// Return the number of nonzeros assembled into the matrix
-    template <class Operators>
+    template <class MatrixOperators, class VectorOperators>
     std::size_t finishMatrixVector(
         SystemMatrixType& matrix,
         SystemVectorType& solution,
         SystemVectorType& rhs,
-        MatrixEntries<Operators>& matrix_operators,
-        VectorEntries<Operators>& rhs_operators) const;
+        MatrixOperators& matrixOperators,
+        VectorOperators& rhsOperators) const;
 
 
     /// Return whether the matrix-block needs to be assembled
-    template <std::size_t R, std::size_t C, class Operators>
-    bool assembleMatrix(const index_t<R>, const index_t<C>,
-                        MatrixEntries<Operators> const& matrix_operators) const
-    {
-      return asmMatrix_ && (!matrix_operators[R][C].assembled || matrix_operators[R][C].changing);
-    }
-
-    /// Return whether the vector-block needs to be assembled
-    template <std::size_t R, class Operators>
-    bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators) const
+    template <class Basis0, class... Bases>
+    auto getElement(Basis0 const& basis0, Bases const&... bases) const
     {
-      return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
+      return basis0.localView().element();
     }
 
   private:
-    FiniteElementSpaces<FeSpaces> globalBases_;
+    std::shared_ptr<GlobalBasis> globalBasis_;
     bool asmMatrix_;
     bool asmVector_;
   };
diff --git a/dune/amdis/Assembler.inc.hpp b/dune/amdis/Assembler.inc.hpp
index 5b4112a806de1188d120bc7e1e55c2299bd0944a..5ba790d2bf0a35ba2ecc26f712411c0f0fd03fc9 100644
--- a/dune/amdis/Assembler.inc.hpp
+++ b/dune/amdis/Assembler.inc.hpp
@@ -2,224 +2,131 @@
 
 namespace AMDiS {
 
-template <class FeSpaces>
-  template <class Operators>
-void Assembler<FeSpaces>::assemble(
+template <class GlobalBasis>
+  template <class MatrixOperators, class VectorOperators>
+void Assembler<GlobalBasis>::assemble(
     GridView const& gv,
     SystemMatrixType& matrix,
     SystemVectorType& solution,
     SystemVectorType& rhs,
-    MatrixEntries<Operators>& matrix_operators,
-    VectorEntries<Operators>& rhs_operators)
+    MatrixOperators& matrixOperators,
+    VectorOperators& rhsOperators)
 {
   // 1. Update global bases
-  globalBases_.update(gv);
+  globalBasis_.update(gv);
 
   // 2. init matrix and rhs vector and initialize dirichlet boundary conditions
-  initMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
+  initMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
 
   // 3. traverse grid and assemble operators on the elements
-  for (auto const& element : elements(globalBases_.gridView()))
+  for (auto const& element : elements(globalBasis_.gridView()))
   {
-    globalBases_.bind(element);
-
-    MatrixEntries<ElementMatrix> elementMatrix;
-    VectorEntries<ElementVector> elementVector;
-
-    MatrixEntries<bool> addMat = assembleElementMatrices(matrix_operators, elementMatrix);
-    VectorEntries<bool> addVec = assembleElementVectors(rhs_operators, elementVector);
-
-    addElementMatrices(matrix, addMat, elementMatrix);
-    addElementVectors(rhs, addVec, elementVector);
-
-    globalBases_.unbind();
+    globalBasis_.bind(element);
+    assembleElement(matrix, rhs, matrixOperators, rhsOperators);
+    globalBasis_.unbind();
   }
 
   // 4. finish matrix insertion and apply dirichlet boundary conditions
-  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
+  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
 
   msg("fillin of assembled matrix: ", nnz);
 }
 
 
-template <class FeSpaces>
-  template <class Operators>
-void Assembler<FeSpaces>::initMatrixVector(
+template <class GlobalBasis>
+  template <class MatrixOperators, class VectorOperators>
+void Assembler<GlobalBasis>::initMatrixVector(
     SystemMatrixType& matrix,
     SystemVectorType& solution,
     SystemVectorType& rhs,
-    MatrixEntries<Operators>& matrix_operators,
-    VectorEntries<Operators>& rhs_operators) const
+    MatrixOperators& matrixOperators,
+    VectorOperators& rhsOperators) const
 {
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
-  {
-    static const int R = decltype(_r)::value;
-    auto const& rowFeSpace = globalBases_.basis(_r);
-
-    msg(rowFeSpace.size(), " DOFs for FeSpace[", R, "]");
-
-    if (this->assembleVector(_r, rhs_operators)) {
-      rhs.compress(_r);
-      rhs.getVector(_r) = 0.0;
-
-      // init vector operators
-      for (auto& scaled : rhs_operators[R].element)
-        scaled.op->init(rowFeSpace);
-      for (auto& scaled : rhs_operators[R].boundary)
-        scaled.op->init(rowFeSpace);
-      for (auto& scaled : rhs_operators[R].intersection)
-        scaled.op->init(rowFeSpace);
-    }
+  matrix.init(globalBasis_);
+  solution.init(globalBasis_);
+  rhs.init(globalBasis_);
 
-    forEach(range_<0, nComponents>, [&,this](auto const _c)
-    {
-      static const int C = decltype(_c)::value;
-      auto const& colFeSpace = globalBases_.basis(_c);
-
-      bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
-      matrix(_r, _c).init(asmMatrix);
-
-      if (asmMatrix) {
-        // init matrix operators
-        for (auto& scaled : matrix_operators[R][C].element)
-          scaled.op->init(rowFeSpace, colFeSpace);
-        for (auto& scaled : matrix_operators[R][C].boundary)
-          scaled.op->init(rowFeSpace, colFeSpace);
-        for (auto& scaled : matrix_operators[R][C].intersection)
-          scaled.op->init(rowFeSpace, colFeSpace);
-
-        // init boundary condition
-        for (int c = 0; c < nComponents; ++c)
-          for (auto bc : matrix_operators[R][c].dirichlet)
-            bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
-      }
-    });
-  });
-}
+  auto localView = globalBasis_.localView();
+  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
+  {
+    if (rowNode.isLeaf)
+      msg(globalBasis_.size(rowTreePath), " DOFs for Basis[", Dune::TypeTree::treePathIndex(rowTreePath), "]");
 
+    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
+    if (rhsOperators[rowNode].assemble(asmVector_))
+      rhsOperators[rowNode].init(rowBasis);
 
-template <class FeSpaces>
-  template <class Operators>
-typename Assembler<FeSpaces>::template MatrixEntries<bool>
-Assembler<FeSpaces>::assembleElementMatrices(
-    MatrixEntries<Operators>& operators,
-    MatrixEntries<ElementMatrix>& elementMatrix) const
-{
-  MatrixEntries<bool> addMat;
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
-  {
-    static const std::size_t R = decltype(_r)::value;
-    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
     {
-      static const std::size_t C = decltype(_c)::value;
-
-      // assemble block of element matrix
-      addMat[R][C] = this->assembleMatrix(_r, _c, operators)
-        ? this->assembleElementMatrix(operators[R][C], elementMatrix[R][C],
-                                      globalBases_.localView(_r), globalBases_.localView(_c))
-        : false;
+      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
+      if (matrixOperators[rowNode][colNode].assemble(asmMatrix_))
+        matrixOperators[rowNode][colNode].init(rowBasis, colBasis);
+
+      // init boundary condition
+      // for (int c = 0; c < nComponents; ++c)
+      //   for (auto bc : matrixOperators[R][c].dirichlet)
+      //     bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
     });
   });
-
-  return addMat;
 }
 
 
-template <class FeSpaces>
-  template <class Operators>
-typename Assembler<FeSpaces>::template VectorEntries<bool>
-  Assembler<FeSpaces>::assembleElementVectors(
-    VectorEntries<Operators>& operators,
-    VectorEntries<ElementVector>& elementVector) const
+template <class GlobalBasis>
+  template <class MatrixOperators, class VectorOperators>
+void Assembler<GlobalBasis>::assembleElement(
+    SystemMatrixType& matrix,
+    SystemVectorType& rhs,
+    MatrixOperators& matrixOperators,
+    VectorOperators& rhsOperators) const
 {
-  VectorEntries<bool> addVec;
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  auto localView = globalBasis_.localView();
+  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
   {
-    static const std::size_t R = decltype(_r)::value;
-
-    // assemble block of element vector
-    addVec[R] = this->assembleVector(_r, operators)
-      ? this->assembleElementVector(operators[R], elementVector[R], globalBases_.localView(_r))
-      : false;
-  });
-
-  return addVec;
-}
-
-
-template <class FeSpaces>
-  template <class Operators, class RowView, class ColView>
-bool Assembler<FeSpaces>::assembleElementMatrix(
-    Operators& operators,
-    ElementMatrix& elementMatrix,
-    RowView const& rowLocalView,
-    ColView const& colLocalView) const
-{
-  if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
-    return false; // nothing to do
-
-  auto const nRows = rowLocalView.tree().finiteElement().size();
-  auto const nCols = colLocalView.tree().finiteElement().size();
+    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
 
-  auto const& element = rowLocalView.element();
-  auto const& gridView = rowLocalView.globalBasis().gridView();
+    auto& rhsOp = rhsOperators[rowNode];
+    if (rhsOp.assemble(asmVector_) && !rhsOp.empty()) {
+      auto elementVector = makeElementVector(rowNode);
+      set_to_zero(elementVector);
 
-  // fills the entire matrix with zeroes
-  elementMatrix.change_dim(nRows, nCols);
-  set_to_zero(elementMatrix);
-
-  bool add = false;
-
-  auto assemble_operators = [&](auto& e, auto& operator_list) {
-    for (auto scaled : operator_list) {
-      bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
-      add = add || add_op;
+      assembleElementOperators(elementVector, rhs, rhsOp, rowBasis);
     }
-  };
 
-  // assemble element operators
-  assemble_operators(element, operators.element);
+    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
+    {
+      auto& matOp = matrixOperators[rowNode][colNode];
+      if (matOp.assemble(asmMatrix_) && !matOp.empty()) {
+        auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
 
-  // assemble intersection operators
-  if (!operators.intersection.empty()
-      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
-  {
-    for (auto const& intersection : intersections(gridView, element)) {
-      if (intersection.boundary())
-        assemble_operators(intersection, operators.boundary);
-      else
-        assemble_operators(intersection, operators.intersection);
-    }
-  }
+        auto elementMatrix = makeElementMatrix(rowNode, colNode);
+        set_to_zero(elementMatrix);
 
-  return add;
+        assembleElementOperators(elementMatrix, matrix, matOp, rowBasis, colBasis);
+      }
+    });
+  });
 }
 
 
-template <class FeSpaces>
-  template <class Operators, class RowView>
-bool Assembler<FeSpaces>::assembleElementVector(
+template <class GlobalBasis>
+  template <class ElementContainer, class Container, class Operators, class... Bases>
+void Assembler<GlobalBasis>::assembleElementOperators(
+    ElementContainer& elementContainer,
+    Container& container,
     Operators& operators,
-    ElementVector& elementVector,
-    RowView const& rowLocalView) const
+    Bases const&... subBases) const
 {
-  if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
+  if (operators.empty())
     return false;
 
-  auto const nRows = rowLocalView.tree().finiteElement().size();
-
-  auto const& element = rowLocalView.element();
-  auto const& gridView = rowLocalView.globalBasis().gridView();
-
-  // Set all entries to zero
-  elementVector.change_dim(nRows);
-  set_to_zero(elementVector);
+  auto const& element = getElement(subBases...);
+  auto const& gridView = subBasis.gridView();
 
   bool add = false;
 
-  auto assemble_operators = [&](auto& e, auto& operator_list) {
+  auto assemble_operators = [&](auto const& context, auto& operator_list) {
     for (auto scaled : operator_list) {
-      bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
+      bool add_op = scaled.op->assemble(gridView, context, subBases.localView()..., elementContainer, scaled.factor);
       add = add || add_op;
     }
   };
@@ -239,109 +146,50 @@ bool Assembler<FeSpaces>::assembleElementVector(
     }
   }
 
-  return add;
-}
+  if (!add)
+    return;
 
-
-template <class FeSpaces>
-void Assembler<FeSpaces>::addElementMatrices(
-    SystemMatrixType& dofmatrix,
-    MatrixEntries<bool> const& addMat,
-    MatrixEntries<ElementMatrix> const& elementMatrix) const
-{
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
-  {
-    static const std::size_t R = decltype(_r)::value;
-    forEach(range_<0, nComponents>, [&,this](auto const _c)
-    {
-      static const std::size_t C = decltype(_c)::value;
-      if (!addMat[R][C])
-        return;
-
-      auto const& rowIndexSet = globalBases_.localIndexSet(_r);
-      auto const& colIndexSet = globalBases_.localIndexSet(_c);
-
-      // NOTE: current implementation does not utilize the multi-indices that we get from localIndexSet.
-
-      for (std::size_t i = 0; i < num_rows(elementMatrix[R][C]); ++i) {
-        // The global index of the i−th vertex of the element
-        auto const row = rowIndexSet.index(i);
-        for (std::size_t j = 0; j < num_cols(elementMatrix[R][C]); ++j) {
-          // The global index of the j−th vertex of the element
-          auto const col = colIndexSet.index(j);
-          dofmatrix(_r,_c)(row,col) += elementMatrix[R][C](i,j);
-        }
-      }
-    });
-  });
-}
-
-
-template <class FeSpaces>
-void Assembler<FeSpaces>::addElementVectors(
-    SystemVectorType& dofvector,
-    VectorEntries<bool> const& addVec,
-    VectorEntries<ElementVector> const& elementVector) const
-{
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
-  {
-    static const std::size_t R = decltype(_r)::value;
-    if (!addVec[R])
-      return;
-
-    auto const& localIndexSet = globalBases_.localIndexSet(_r);
-
-    for (std::size_t i = 0; i < size(elementVector[R]); ++i) {
-      // The global index of the i-th vertex
-      auto const row = localIndexSet.index(i);
-      dofvector[_r][row] += elementVector[R][i];
-    }
-  });
+  elementContainer.apply(subBases.localIndexSet()..., container);
 }
 
 
-template <class FeSpaces>
-  template <class Operators>
-std::size_t Assembler<FeSpaces>::finishMatrixVector(
+template <class GlobalBasis>
+  template <class MatrixOperators, class VectorOperators>
+std::size_t Assembler<GlobalBasis>::finishMatrixVector(
     SystemMatrixType& matrix,
     SystemVectorType& solution,
     SystemVectorType& rhs,
-    MatrixEntries<Operators>& matrix_operators,
-    VectorEntries<Operators>& rhs_operators) const
+    MatrixOperators& matrixOperators,
+    VectorOperators& rhsOperators) const
 {
-  std::size_t nnz = 0;
-  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  matrix.finish();
+
+  auto localView = globalBasis_.localView();
+  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
   {
-    static const int R = decltype(_r)::value;
-    if (this->assembleVector(_r, rhs_operators))
-      rhs_operators[R].assembled = true;
+    auto& rhsOp = rhsOperators[rowNode];
+    if (rhsOp.assemble(asmVector_))
+      rhsOp.assembled = true;
 
-    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
     {
-      static const int C = decltype(_c)::value;
-      bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
-
-      matrix(_r, _c).finish();
-
-      if (asmMatrix) {
-        matrix_operators[R][C].assembled = true;
-
-        // finish boundary condition
-        for (int c = 0; c < nComponents; ++c) {
-          for (int r = 0; r < nComponents; ++r) {
-            if (r != R && c != C)
-              continue;
-            for (auto bc : matrix_operators[r][c].dirichlet)
-              bc->finish(r == R, c == C, matrix(_r, _c), solution[_c], rhs[_r]);
-          }
-        }
-
-        nnz += matrix(_r, _c).getNnz();
-      }
+      auto& matOp = matrixOperators[rowNode][colNode];
+      if (matOp.assemble(asmMatrix_))
+        matOp.assembled = true;
+
+      // finish boundary condition
+      // for (int c = 0; c < nComponents; ++c) {
+      //   for (int r = 0; r < nComponents; ++r) {
+      //     if (r != R && c != C)
+      //       continue;
+      //     for (auto bc : matrixOperators[r][c].dirichlet)
+      //       bc->finish(r == R, c == C, matrix(_r, _c), solution[_c], rhs[_r]);
+      //   }
+      // }
     });
   });
 
-  return nnz;
+  return matrix.getNnz();
 }
 
 } // end namespace AMDiS
diff --git a/dune/amdis/CMakeLists.txt b/dune/amdis/CMakeLists.txt
index a2f3a368c9b7b3e4f81fd7427c66b9389fdb2f2c..a05425672c06048d5675db94e7c54d7329eb6939 100644
--- a/dune/amdis/CMakeLists.txt
+++ b/dune/amdis/CMakeLists.txt
@@ -19,7 +19,7 @@ dune_add_library("duneamdis" NO_EXPORT
     )
 add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
 target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)
-target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0)
+target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0 -Wall -pedantic -Wno-unused-parameter)
 
 find_package(MTL REQUIRED
              PATHS /usr/local/lib/mtl4)
diff --git a/dune/amdis/CreatorMap.hpp b/dune/amdis/CreatorMap.hpp
index 375ecd488f27583d68849ef72982c5fed56b8557..24b40c7e1d05bda457deb4f8379628e16031c60b 100644
--- a/dune/amdis/CreatorMap.hpp
+++ b/dune/amdis/CreatorMap.hpp
@@ -12,7 +12,7 @@ namespace AMDiS
   // All classes that need creators must specialize this class and implement
   // a static void init() method.
   template <class BaseClass>
-  struct DefaultCreators;
+  class DefaultCreators;
 
 
   /** \ingroup Common
diff --git a/dune/amdis/FileWriter.hpp b/dune/amdis/FileWriter.hpp
index b03d735a4648aa4b7b62448fd7b3cd38736f9648..756eee7d9e7d69e351106bb477228fd50ccc919c 100644
--- a/dune/amdis/FileWriter.hpp
+++ b/dune/amdis/FileWriter.hpp
@@ -28,7 +28,7 @@ namespace AMDiS
     static constexpr int dim = Traits::dim;
 
     /// Dimension of the world
-    static constexpr int dow = Traits::dimworld;
+    static constexpr int dow = Traits::dow;
 
     /// Number of problem components
     static constexpr int nComponents = Traits::nComponents;
diff --git a/dune/amdis/FiniteElementSpaces.hpp b/dune/amdis/FiniteElementSpaces.hpp
index c1ac22b9057113bba82819ee70e42ec1f3fda12f..482f6dfeafbcd42a48bec0785801d7dee010d262 100644
--- a/dune/amdis/FiniteElementSpaces.hpp
+++ b/dune/amdis/FiniteElementSpaces.hpp
@@ -2,6 +2,7 @@
 
 #include <dune/amdis/common/TupleUtility.hpp>
 #include <dune/amdis/common/IndexSeq.hpp>
+#include <dune/amdis/common/Loops.hpp>
 
 namespace AMDiS
 {
@@ -61,10 +62,15 @@ namespace AMDiS
       forEach(range_<0, nComponents>, [&,this](auto const _i)
       {
         static const int I = decltype(_i)::value;
-        std::get<I>(localViews_).bind(element);
-        std::get<I>(localIndexSets_).bind(std::get<I>(localViews_));
+
+        auto& localView = std::get<I>(localViews_);
+        localView.bind(element);
+
+        auto& localIndexSet = std::get<I>(localIndexSets_);
+        localIndexSet.bind(localView);
       });
       // NOTE: maybe create element-geometry here
+      bound_ = true;
     }
 
     /// Unbind from the current element
@@ -76,10 +82,11 @@ namespace AMDiS
         std::get<I>(localIndexSets_).unbind();
         std::get<I>(localViews_).unbind();
       });
+      bound_ = false;
     }
 
     template <std::size_t I>
-    auto const& basis(const index_t<I> _i = {}) const
+    auto const& feSpace(const index_t<I> _i = {}) const
     {
       return std::get<I>(*feSpaces_);
     }
@@ -87,28 +94,40 @@ namespace AMDiS
     template <std::size_t I>
     auto const& localView(const index_t<I> _i = {}) const
     {
+      assert( bound_ && "localViews must be bound to an element." );
       return std::get<I>(localViews_);
     }
 
     template <std::size_t I>
     auto const& localIndexSet(const index_t<I> _i = {}) const
     {
+      assert( bound_ && "localIndexSets must be bound to a localView." );
       return std::get<I>(localIndexSets_);
     }
 
+    auto const& element() const
+    {
+      assert( bound_ && "localViews must be bound to an element." );
+      return std::get<0>(localViews_).element();
+    }
+
     auto const& gridView() const
     {
       return std::get<0>(*feSpaces_).gridView();
     }
 
   private:
+    /// Tuple of global functionspace bases
     std::shared_ptr<FeSpaces> feSpaces_;
 
-    /// tuple of localView objects, obtained from the tuple of global bases
+    /// Tuple of localView objects, obtained from the tuple of global bases
     LocalViews localViews_;
 
-    // tuple of localIndexSet objects, obtained from the tuple of global bases
+    /// Tuple of localIndexSet objects, obtained from the tuple of global bases
     LocalIndexSets localIndexSets_;
+
+    /// True, if localViews and localIndexSets are bound to an element
+    bool bound_ = false;
   };
 
 } // end namespace AMDiS
diff --git a/dune/amdis/LocalAssembler.hpp b/dune/amdis/LocalAssembler.hpp
index 61efb848f84e2cd2d6eab5fd5c05e38cb442790e..27b079e23ea945fecd408ef64a6c235c65ec787f 100644
--- a/dune/amdis/LocalAssembler.hpp
+++ b/dune/amdis/LocalAssembler.hpp
@@ -70,8 +70,6 @@ namespace AMDiS
 
     /// geometry() for entities or geometryInInside() for intersections
     Geometry geometry;
-
-    // TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
   };
 
 } // end namespace AMDiS
diff --git a/dune/amdis/Operator.hpp b/dune/amdis/Operator.hpp
index 34e675a77781f0c29b2419e5ecad8067be3ac1ea..7d52274dd203ff85527925fd2de2dfe0850d0ae0 100644
--- a/dune/amdis/Operator.hpp
+++ b/dune/amdis/Operator.hpp
@@ -18,20 +18,15 @@ namespace AMDiS
 {
 
   // the LocalContext is either an Codim=0-EntityType or an IntersectionType
-  template <class MeshView, class LocalContext>
+  template <class GridView, class LocalContext>
   class Operator
   {
     using Self = Operator;
     using OperatorTermType = OperatorTerm<LocalContext>;
 
-    using ElementVector = Impl::ElementVector;
-    using ElementMatrix = Impl::ElementMatrix;
-
-    using IdType = typename MeshView::Grid::LocalIdSet::IdType;
-
-    template <class, class> friend class ZeroOrderAssembler;
+    template <class, class>                 friend class ZeroOrderAssembler;
     template <class, class, FirstOrderType> friend class FirstOrderAssembler;
-    template <class, class> friend class SecondOrderAssembler;
+    template <class, class>                 friend class SecondOrderAssembler;
 
   public:
     /// \brief Add coefficients for zero-order operator < C(phi), psi >.
@@ -114,13 +109,13 @@ namespace AMDiS
 
   public: // initialize and assemble operator on element
 
-    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
-    template <class RowFeSpace, class ColFeSpace>
-    void init(RowFeSpace const& rowFeSpace,
-              ColFeSpace const& colFeSpace);
+    /// Extract the polynomial degree from \p rowBasis and \p colBasis.
+    template <class RowBasis, class ColBasis>
+    void init(RowBasis const& rowBasis,
+              ColBasis const& colBasis);
 
-    template <class RowFeSpace>
-    void init(RowFeSpace const& rowFeSpace);
+    template <class RowBasis>
+    void init(RowBasis const& rowBasis);
 
 
     /// Calculates the needed quadrature degree for the given term-order \p order.
@@ -131,18 +126,24 @@ namespace AMDiS
                             FirstOrderType firstOrderType = GRD_PHI);
 
 
-    template <class RowView, class ColView>
-    bool getElementMatrix(LocalContext const& element,
-                          RowView const& rowView,
-                          ColView const& colView,
-                          ElementMatrix& elementMatrix,
-                          double* factor = NULL);
+    // assemble matrix operator
+    template <class RowView, class ColView, class ElementMatrix>
+    bool assemble(
+        GridView const& gv,
+        LocalContext const& element,
+        RowView const& rowView,
+        ColView const& colView,
+        ElementMatrix& elementMatrix,
+        double* factor = NULL);
 
-    template <class RowView>
-    bool getElementVector(LocalContext const& element,
-                          RowView const& rowView,
-                          ElementVector& elementVector,
-                          double* factor = NULL);
+    // assemble vector operator
+    template <class LocalView, class ElementVector>
+    bool assemble(
+        GridView const& gv,
+        LocalContext const& element,
+        LocalView const& localView,
+        ElementVector& elementVector,
+        double* factor = NULL);
 
 
   private: // implementation details
@@ -228,7 +229,10 @@ namespace AMDiS
   private:
 
     template <class LocalView>
-    IdType getElementId(LocalView const& localView);
+    std::size_t getElementIndex(GridView const& gridView, LocalView const& localView) const
+    {
+      return gridView.indexSet().index(localView.element());
+    }
 
 
   private:
@@ -248,13 +252,154 @@ namespace AMDiS
     int psiDegree = 1;
     int phiDegree = 1;
 
-    IdType lastMatrixId = 0;
-    IdType lastVectorId = 0;
+    std::size_t lastMatrixIndex = 0;
+    std::size_t lastVectorIndex = 0;
 
     ElementMatrix cachedElementMatrix;
     ElementVector cachedElementVector;
   };
 
+
+  template <class GlobalBasis>
+  class TreeOperators
+  {
+  protected:
+    using GridView = typename GlobalBasis::GridView;
+    using Element = GridView::template Codim<0>::Entity;
+
+    using ElementOperator = Operator<GridView, Element>;
+    using IntersectionOperator = Operator<GridView, typename GridView::Intersection>;
+
+    using Tree = typename Basis::LocalView::Tree;
+
+    template <class OperatorType>
+    struct Scaled
+    {
+      std::shared_ptr<OperatorType> op;
+      double* factor = nullptr;
+      double* estFactor = nullptr;
+      BoundaryType b = {0};
+    };
+
+    template <class Node>
+    struct NodeData
+    {
+      std::list<Scaled<ElementOperator>> element;
+      std::list<Scaled<IntersectionOperator>> boundary;
+      std::list<Scaled<IntersectionOperator>> intersection;
+
+      bool assembled = false; // if false, do reassemble
+      bool changing = false; // if true, or assembled false, do reassemble
+
+      template <class.. Bases>
+      void init(Bases&&... bases)
+      {
+        for (auto scaled : element)
+          scaled.op->init(std::forward<Bases>(bases)...);
+        for (auto scaled : boundary)
+          scaled.op->init(std::forward<Bases>(bases)...);
+        for (auto scaled : intersection)
+          scaled.op->init(std::forward<Bases>(bases)...);
+      }
+
+      bool empty() const
+      {
+        return element.empty() && boundary.empty() && intersection.empty();
+      }
+
+      bool assemble(bool flag) const
+      {
+        return flag && (!assembled || changing);
+      }
+    };
+  };
+
+  template <class GlobalBasis>
+  class MatrixOperators
+      : public TreeOperators<GlobalBasis>
+  {
+    using Tree = typename TreeOperators<GlobalBasis>::Tree;
+    using NodeData = typename TreeOperators<GlobalBasis>::NodeData;
+
+    template <class Node>
+    using RowNodeData = Dune::Functions::TreeData<Tree, NodeData, false>;
+
+    using Data = Dune::Functions::TreeData<Tree, RowNodeData, false>;
+
+  public:
+    MatrixOperators(GlobalBasis const& globalBasis)
+    {
+      auto localView = globalBasis.localView();
+      data_.init(localView.tree());
+      forEachNode(localView.tree(), [&,this](auto const& node, auto const treePath)
+      {
+        data_[node].init(localView.tree());
+      });
+    }
+
+    template <class Node>
+    auto& operator[](Node const& node)
+    {
+      return data_[node];
+    }
+
+    template <class Node>
+    auto const& operator[](Node const& node) const
+    {
+      return data_[node];
+    }
+
+    template <class RowNode, class ColNode>
+    auto& operator()(RowNode const& rowNode, ColNode const& colNode)
+    {
+      return data_[rowNode][colNode];
+    }
+
+    template <class RowNode, class ColNode>
+    auto const& operator()(RowNode const& rowNode, ColNode const& colNode) const
+    {
+      return data_[rowNode][colNode];
+    }
+
+  private:
+    Data data_;
+  };
+
+  template <class GlobalBasis>
+  class VectorOperators
+      : public TreeOperators<GlobalBasis>
+  {
+    using Tree = typename TreeOperators<GlobalBasis>::Tree;
+    using NodeData = typename TreeOperators<GlobalBasis>::NodeData;
+    using Data = Dune::Functions::TreeData<Tree, NodeData, false>;
+
+  public:
+    VectorOperators(GlobalBasis const& globalBasis)
+    {
+      auto localView = globalBasis.localView();
+      data_.init(localView.tree());
+    }
+
+    template <class Node>
+    auto& operator[](Node const& node)
+    {
+      return data_[node];
+    }
+
+    template <class Node>
+    auto const& operator[](Node const& node) const
+    {
+      return data_[node];
+    }
+
+  private:
+    Data data_;
+  };
+
+
+
+
+
 } // end namespace AMDiS
 
 #include "Operator.inc.hpp"
diff --git a/dune/amdis/Operator.inc.hpp b/dune/amdis/Operator.inc.hpp
index 8a784852ca4d6035de7d4f7d3585590de57df98e..7e06852a37a09c4dcbcd53143a7564de1488f20b 100644
--- a/dune/amdis/Operator.inc.hpp
+++ b/dune/amdis/Operator.inc.hpp
@@ -1,15 +1,14 @@
 #pragma once
 
-namespace AMDiS
+namespace AMDiS {
+
+template <class GridView, class Element>
+  template <class RowBasis, class ColBasis>
+void Operator<GridView, Element>::
+init(RowBasis const& rowBasis, ColBasis const& colBasis)
 {
-  template <class MeshView, class Element>
-    template <class RowFeSpace, class ColFeSpace>
-  void Operator<MeshView, Element>::
-  init(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace)
-  {
-    using IdType = typename Operator<MeshView, Element>::IdType;
-    lastMatrixId = std::numeric_limits<IdType>::max();
-    lastVectorId = std::numeric_limits<IdType>::max();
+  lastMatrixIndex = std::numeric_limits<std::size_t>::max();
+  lastVectorIndex = std::numeric_limits<std::size_t>::max();
 
 //     auto const& rowFE = rowView.tree().finiteElement();
 //     auto const& colFE = colView.tree().finiteElement();
@@ -17,238 +16,224 @@ namespace AMDiS
 //     psiDegree = rowFE.localBasis().order();
 //     phiDegree = colFE.localBasis().order();
 
-    psiDegree = getPolynomialDegree<RowFeSpace>;
-    phiDegree = getPolynomialDegree<ColFeSpace>;
+  psiDegree = getPolynomialDegree<RowBasis>;
+  phiDegree = getPolynomialDegree<ColBasis>;
 
-    // TODO: calc quadrature degree here.
-  }
+  // TODO: calc quadrature degree here.
+}
 
 
-  template <class MeshView, class Element>
-    template <class RowFeSpace>
-  void Operator<MeshView, Element>::
-  init(RowFeSpace const& rowFeSpace)
-  {
-    using IdType = typename Operator<MeshView, Element>::IdType;
-    lastVectorId = std::numeric_limits<IdType>::max();
+template <class GridView, class Element>
+  template <class RowBasis>
+void Operator<GridView, Element>::
+init(RowBasis const& rowBasis)
+{
+  lastVectorIndex = std::numeric_limits<std::size_t>::max();
+  psiDegree = getPolynomialDegree<RowBasis>;
+}
 
-    psiDegree = getPolynomialDegree<RowFeSpace>;
-  }
 
+template <class GridView, class Element>
+  template <class Geometry>
+int Operator<GridView, Element>::
+getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType type)
+{
+  std::list<OperatorTermType*>* terms = NULL;
 
-  template <class MeshView, class Element>
-    template <class Geometry>
-  int Operator<MeshView, Element>::
-  getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType type)
+  switch(order)
   {
-    std::list<OperatorTermType*>* terms = NULL;
-
-    switch(order)
-    {
-    case 0:
-      terms = &zeroOrder;
-      break;
-    case 1:
-      if (type == GRD_PHI)
-        terms = &firstOrderGrdPhi;
-      else
-        terms = &firstOrderGrdPsi;
-      break;
-    case 2:
-      terms = &secondOrder;
-      break;
-    }
-    int maxTermDegree = 0;
+  case 0:
+    terms = &zeroOrder;
+    break;
+  case 1:
+    if (type == GRD_PHI)
+      terms = &firstOrderGrdPhi;
+    else
+      terms = &firstOrderGrdPsi;
+    break;
+  case 2:
+    terms = &secondOrder;
+    break;
+  }
+  int maxTermDegree = 0;
 
-    for (OperatorTermType* term : *terms)
-      maxTermDegree = std::max(maxTermDegree, term->getDegree(t));
+  for (OperatorTermType* term : *terms)
+    maxTermDegree = std::max(maxTermDegree, term->getDegree(t));
 
-    int degree = psiDegree + phiDegree + maxTermDegree;
-    if (t.isSimplex())
-      degree -= order;
+  int degree = psiDegree + phiDegree + maxTermDegree;
+  if (t.isSimplex())
+    degree -= order;
 
-    if (!geometry.affine())
-      degree += 1; // oder += (order+1)
+  if (!geometry.affine())
+    degree += 1; // oder += (order+1)
 
-    return degree;
-  }
+  return degree;
+}
 
 
-  template <class MeshView, class Element>
-    template <class LocalView>
-  typename Operator<MeshView, Element>::IdType
-  Operator<MeshView, Element>::getElementId(LocalView const& localView)
-  {
-    auto const& element = localView.element();
-    auto const& grid = localView.globalBasis().gridView().grid();
+template <class GridView, class Element>
+  template <class RowView, class ColView, class ElementMatrix>
+bool Operator<GridView, Element>::assemble(
+    GridView const& gridView,
+    Element const& element,
+    RowView const& rowView,
+    ColView const& colView,
+    ElementMatrix& elementMatrix,
+    double* factor)
+{
+  double fac = factor ? *factor : 1.0;
 
-    auto const& idSet = grid.localIdSet();
-    return idSet.id(element);
-  }
+  if (fac == 0.0)
+    return false;
 
+  const std::size_t nRows = rowView.tree().finiteElement().size();
+  const std::size_t nCols = colView.tree().finiteElement().size();
 
-  template <class MeshView, class Element>
-    template <class RowView, class ColView>
-  bool Operator<MeshView, Element>::
-  getElementMatrix(Element const& element,
-                   RowView const& rowView,
-                   ColView const& colView,
-                   ElementMatrix& elementMatrix,
-                   double* factor)
-  {
-    double fac = factor ? *factor : 1.0;
-
-    if (fac == 0.0)
-      return false;
-
-    const auto nRows = rowView.tree().finiteElement().size();
-    const auto nCols = colView.tree().finiteElement().size();
-
-    auto currentId = getElementId(rowView);
-    bool useCachedMatrix = (lastMatrixId == currentId
-                            && num_rows(cachedElementMatrix) == nRows
-                            && num_cols(cachedElementMatrix) == nCols);
-
-    bool assign = true;
-    if (!useCachedMatrix) {
-
-      cachedElementMatrix.change_dim(nRows, nCols);
-      set_to_zero(cachedElementMatrix);
-
-      if (!zeroOrder.empty()) {
-        ZeroOrderAssembler<MeshView, Element> assembler(*this, element);
-        assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
-      }
-      if (!firstOrderGrdPhi.empty()) {
-        FirstOrderAssembler<MeshView, Element, GRD_PHI> assembler(*this, element);
-        assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
-      }
-      if (!firstOrderGrdPsi.empty()) {
-        FirstOrderAssembler<MeshView, Element, GRD_PSI> assembler(*this, element);
-        assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
-      }
-      if (!secondOrder.empty()) {
-        SecondOrderAssembler<MeshView, Element> assembler(*this, element);
-        assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
-      }
-
-      assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
-               !firstOrderGrdPsi.empty() || !secondOrder.empty();
-    }
+  auto currentIndex = getElementIndex(gridView, rowView.element());
+  bool useCachedMatrix = (lastMatrixIndex == currentIndex
+                          && num_rows(cachedElementMatrix) == nRows
+                          && num_cols(cachedElementMatrix) == nCols);
+
+  bool assign = true;
+  if (!useCachedMatrix) {
 
-    if (assign)
-      elementMatrix += fac * cachedElementMatrix;
+    cachedElementMatrix.change_dim(nRows, nCols);
+    set_to_zero(cachedElementMatrix);
 
-    lastMatrixId = currentId;
+    if (!zeroOrder.empty()) {
+      ZeroOrderAssembler<GridView, Element> assembler(*this, element);
+      assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
+    }
+    if (!firstOrderGrdPhi.empty()) {
+      FirstOrderAssembler<GridView, Element, GRD_PHI> assembler(*this, element);
+      assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
+    }
+    if (!firstOrderGrdPsi.empty()) {
+      FirstOrderAssembler<GridView, Element, GRD_PSI> assembler(*this, element);
+      assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
+    }
+    if (!secondOrder.empty()) {
+      SecondOrderAssembler<GridView, Element> assembler(*this, element);
+      assembler.calculateElementMatrix(*this, rowView, colView, cachedElementMatrix);
+    }
 
-    return true;
+    assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
+              !firstOrderGrdPsi.empty() || !secondOrder.empty();
   }
 
+  if (assign)
+    elementMatrix += fac * cachedElementMatrix;
 
-  template <class MeshView, class Element>
-    template <class RowView>
-  bool Operator<MeshView, Element>::
-  getElementVector(Element const& element,
-                   RowView const& rowView,
-                   ElementVector& elementVector,
-                   double* factor)
-  {
-    double fac = factor ? *factor : 1.0;
+  lastMatrixIndex = currentIndex;
 
-    if (fac == 0.0)
-      return false;
+  return true;
+}
 
-    test_exit( firstOrderGrdPhi.empty() && secondOrder.empty(),
-               "Derivatives on ansatz-functions not allowed on the vector-side!");
 
-    const auto nRows = rowView.tree().finiteElement().size();
+template <class GridView, class Element>
+  template <class LocalView, class ElementVector>
+bool Operator<GridView, Element>::assemble(
+    Element const& element,
+    LocalView const& localView,
+    ElementVector& elementVector,
+    double* factor)
+{
+  double fac = factor ? *factor : 1.0;
 
-    auto currentId = getElementId(rowView);
-    bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);
+  if (fac == 0.0)
+    return false;
 
-    bool assign = true;
-    if (!useCachedVector) {
-      cachedElementVector.change_dim(nRows);
-      set_to_zero(cachedElementVector);
+  test_exit( firstOrderGrdPhi.empty() && secondOrder.empty(),
+              "Derivatives on ansatz-functions not allowed on the vector-side!");
 
-      if (!zeroOrder.empty()) {
-        ZeroOrderAssembler<MeshView, Element> assembler(*this, element);
-        assembler.calculateElementVector(*this, rowView, cachedElementVector);
-      }
-      if (!firstOrderGrdPsi.empty()) {
-        FirstOrderAssembler<MeshView, Element, GRD_PSI> assembler(*this, element);
-        assembler.calculateElementVector(*this, rowView, cachedElementVector);
-      }
+  const std::size_t nRows = rowView.tree().finiteElement().size();
 
-      assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
-    }
+  auto currentIndex = getElementIndex(rowView);
+  bool useCachedVector = (lastVectorIndex == currentIndex && size(cachedElementVector) == nRows);
 
-    if (assign)
-      elementVector += fac * cachedElementVector;
+  bool assign = true;
+  if (!useCachedVector) {
+    cachedElementVector.change_dim(nRows);
+    set_to_zero(cachedElementVector);
 
-    lastVectorId = currentId;
+    if (!zeroOrder.empty()) {
+      ZeroOrderAssembler<GridView, Element> assembler(*this, element);
+      assembler.calculateElementVector(*this, rowView, cachedElementVector);
+    }
+    if (!firstOrderGrdPsi.empty()) {
+      FirstOrderAssembler<GridView, Element, GRD_PSI> assembler(*this, element);
+      assembler.calculateElementVector(*this, rowView, cachedElementVector);
+    }
 
-    return true;
+    assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
   }
 
+  if (assign)
+    elementVector += fac * cachedElementVector;
 
-  template <class MeshView, class Element>
-    template <class Term>
-  Operator<MeshView, Element>& Operator<MeshView, Element>::
-  addZOTImpl(Term const& term)
-  {
-    zeroOrder.push_back(new GenericOperatorTerm<Element, Term>(term));
-    return *this;
-  }
+  lastVectorIndex = currentIndex;
 
+  return true;
+}
 
-  template <class MeshView, class Element>
-    template <class Term>
-  Operator<MeshView, Element>& Operator<MeshView, Element>::
-  addFOTImpl(Term const& term, FirstOrderType type)
-  {
-    using OpTerm = GenericOperatorTerm<Element, Term>;
-    if (type == GRD_PHI)
-      firstOrderGrdPhi.push_back(new OpTerm(term));
-    else
-      firstOrderGrdPsi.push_back(new OpTerm(term));
-    return *this;
-  }
 
+template <class GridView, class Element>
+  template <class Term>
+Operator<GridView, Element>& Operator<GridView, Element>::
+addZOTImpl(Term const& term)
+{
+  zeroOrder.push_back(new GenericOperatorTerm<Element, Term>(term));
+  return *this;
+}
 
-  template <class MeshView, class Element>
-    template <class Term>
-  Operator<MeshView, Element>& Operator<MeshView, Element>::
-  addFOTImpl(Term const& term, std::size_t i, FirstOrderType type)
-  {
-    using OpTerm = GenericOperatorTerm<Element, Term, VectorComponent>;
 
-    if (type == GRD_PHI)
-      firstOrderGrdPhi.push_back(new OpTerm(term, {i}));
-    else
-      firstOrderGrdPsi.push_back(new OpTerm(term, {i}));
-    return *this;
-  }
+template <class GridView, class Element>
+  template <class Term>
+Operator<GridView, Element>& Operator<GridView, Element>::
+addFOTImpl(Term const& term, FirstOrderType type)
+{
+  using OpTerm = GenericOperatorTerm<Element, Term>;
+  if (type == GRD_PHI)
+    firstOrderGrdPhi.push_back(new OpTerm(term));
+  else
+    firstOrderGrdPsi.push_back(new OpTerm(term));
+  return *this;
+}
+
+
+template <class GridView, class Element>
+  template <class Term>
+Operator<GridView, Element>& Operator<GridView, Element>::
+addFOTImpl(Term const& term, std::size_t i, FirstOrderType type)
+{
+  using OpTerm = GenericOperatorTerm<Element, Term, VectorComponent>;
 
+  if (type == GRD_PHI)
+    firstOrderGrdPhi.push_back(new OpTerm(term, {i}));
+  else
+    firstOrderGrdPsi.push_back(new OpTerm(term, {i}));
+  return *this;
+}
 
-  template <class MeshView, class Element>
-    template <class Term>
-  Operator<MeshView, Element>& Operator<MeshView, Element>::
-  addSOTImpl(Term const& term)
-  {
-    secondOrder.push_back(new GenericOperatorTerm<Element, Term>(term));
-    return *this;
-  }
 
+template <class GridView, class Element>
+  template <class Term>
+Operator<GridView, Element>& Operator<GridView, Element>::
+addSOTImpl(Term const& term)
+{
+  secondOrder.push_back(new GenericOperatorTerm<Element, Term>(term));
+  return *this;
+}
 
-  template <class MeshView, class Element>
-    template <class Term>
-  Operator<MeshView, Element>& Operator<MeshView, Element>::
-  addSOTImpl(Term const& term, std::size_t i, std::size_t j)
-  {
-    using OpTerm = GenericOperatorTerm<Element, Term, MatrixComponent>;
-    secondOrder.push_back(new OpTerm(term, {i,j}));
-    return *this;
-  }
+
+template <class GridView, class Element>
+  template <class Term>
+Operator<GridView, Element>& Operator<GridView, Element>::
+addSOTImpl(Term const& term, std::size_t i, std::size_t j)
+{
+  using OpTerm = GenericOperatorTerm<Element, Term, MatrixComponent>;
+  secondOrder.push_back(new OpTerm(term, {i,j}));
+  return *this;
+}
 
 } // end namespace AMDiS
diff --git a/dune/amdis/OperatorTermBase.hpp b/dune/amdis/OperatorTermBase.hpp
index 3989dc5845d0562bb7a5354e1e8a25cd90b0991c..974c303bf7277ce769ed764a38370406a5e015ee 100644
--- a/dune/amdis/OperatorTermBase.hpp
+++ b/dune/amdis/OperatorTermBase.hpp
@@ -20,7 +20,7 @@ namespace AMDiS
   class OperatorTerm
   {
   protected:
-    static constexpr int dim = LocalContext::dimension;
+    static constexpr int dim = Impl::Get<LocalContext>::Entity::dimension;
     static constexpr int dow = LocalContext::Geometry::coorddimension;
 
     using QuadratureRule = QuadratureRuleFactory_t<LocalContext,double,dim>;
diff --git a/dune/amdis/ProblemStat.hpp b/dune/amdis/ProblemStat.hpp
index aae9c49ee063f341b08b1f764cd96608d895aabf..4954cc03a0d27e381b11e7491f3417a74a39c271 100644
--- a/dune/amdis/ProblemStat.hpp
+++ b/dune/amdis/ProblemStat.hpp
@@ -37,7 +37,7 @@ namespace AMDiS
 {
   struct BoundaryType { int b; };
 
-  template <class Traits>
+  template <class GlobalBasis>
   class ProblemStat
       : public ProblemStatBase
       , public StandardProblemIteration
@@ -46,34 +46,29 @@ namespace AMDiS
 
   public: // typedefs and static constants
 
-    using FeSpaces = typename Traits::FeSpaces;
-    using Mesh     = typename Traits::Mesh;
-    using MeshView = typename Mesh::LeafGridView;
+    using GridView = typename GlobalBasis::GridView;
+    using Grid     = typename GridView::Grid;
 
-    using Codim0   = typename MeshView::template Codim<0>;
-    using Element  = typename Codim0::Entity;
+    using Element  = typename GridView::template Codim<0>::Entity;
 
 
     /// Dimension of the mesh
-    static constexpr int dim = Traits::dim;
+    static constexpr int dim = Grid::dimension;
 
     /// Dimension of the world
-    static constexpr int dow = Traits::dow;
+    static constexpr int dow = Grid::dimensionworld;
 
     /// Number of problem components
-    static constexpr int nComponents = Traits::nComponents;
+    static constexpr int nComponents = Traits::nComponents; /* number of leaf nodes in the basis-tree */
 
 
-    template <std::size_t I>
-    using FeSpace = std::tuple_element_t<I, FeSpaces>;
+    using WorldVector = typename Element::Geometry::GlobalCoordinate;
 
-    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
+    using SystemVectorType = SystemVector<GlobalBasis>;
+    using SystemMatrixType = SystemMatrix<GlobalBasis>;
 
-    using SystemVectorType = SystemVector<FeSpaces>;
-    using SystemMatrixType = SystemMatrix<FeSpaces>;
-
-    using ElementOperator = Operator<MeshView, Element>;
-    using IntersectionOperator = Operator<MeshView, typename MeshView::Intersection>;
+    using ElementOperator = Operator<GridView, Element>;
+    using IntersectionOperator = Operator<GridView, typename GridView::Intersection>;
 
     using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                    typename SystemVectorType::MultiVector>;
@@ -97,11 +92,11 @@ namespace AMDiS
 
     /// Constructor taking additionally a reference to a mesh that is used
     /// instead of the default created mesh, \ref ProblemStat
-    ProblemStat(std::string name, Mesh& mesh)
+    ProblemStat(std::string name, Grid& grid)
       : ProblemStat(name)
     {
-      this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
-      componentMeshes.resize(nComponents, this->mesh.get());
+      this->grid = std::shared_ptr<Grid>(&grid, optional_delete(false));
+      componentMeshes.resize(nComponents, this->grid.get());
 
       meshName = "";
       Parameters::get(name + "->mesh", meshName);
@@ -257,15 +252,15 @@ namespace AMDiS
       linearSolver = solver;
     }
 
-    /// Return a pointer to the mesh, \ref mesh
-    std::shared_ptr<Mesh> getMesh() { return mesh; }
+    /// Return a pointer to the grid, \ref grid
+    std::shared_ptr<Grid> getGrid() { return grid; }
 
     /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
     /// matrices and vectors, as well as the file-writer.
-    void setMesh(Mesh& mesh_)
+    void setGrid(Grid& grid_)
     {
-      mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
-      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());
+      grid = std::shared_ptr<Grid>(&grid_, optional_delete(false));
+      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
 
       createFeSpaces();
       createMatricesAndVectors();
@@ -274,26 +269,26 @@ namespace AMDiS
 
 
     /// Return the gridView of the leaf-level
-    auto leafGridView() { return mesh->leafGridView(); }
+    auto leafGridView() { return grid->leafGridView(); }
 
     /// Return the gridView of levle `level`
-    auto levelGridView(int level) { return mesh->levelGridView(level); }
+    auto levelGridView(int level) { return grid->levelGridView(level); }
 
     /// Return the \ref feSpaces
-    auto getFeSpaces() { return feSpaces; }
+    auto getGlobalBasis() { return globalBasis; }
 
 
     /// Return the I'th \ref feSpaces component
     template <std::size_t I = 0>
-    FeSpace<I> const& getFeSpace(const index_t<I> = {}) const
+    FeSpace<I> const& getComponentBasis(const index_t<I> = {}) const
     {
-      return std::get<I>(*feSpaces);
+      return std::get<I>(*globalBasis);
     }
 
     template <std::size_t I = 0>
-    FeSpace<I>& getFeSpace(const index_t<I> = {})
+    FeSpace<I>& getComponentBasis(const index_t<I> = {})
     {
-      return std::get<I>(*feSpaces);
+      return std::get<I>(*globalBasis);
     }
 
     /// Implementation of \ref ProblemStatBase::getName
@@ -310,33 +305,33 @@ namespace AMDiS
 
   protected: // initialization methods
 
-    void createMesh()
+    void createGrid()
     {
-      meshName = "";
-      Parameters::get(name + "->mesh", meshName);
-      test_exit(!meshName.empty(), "No mesh name specified for '", name, "->mesh'!");
+      gridName = "";
+      Parameters::get(name + "->mesh", gridName);
+      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
 
-      mesh = MeshCreator<Mesh>::create(meshName);
+      grid = MeshCreator<Grid>::create(gridName);
 
-      msg("Create mesh:");
-      msg("#elements = "   , mesh->size(0));
-      msg("#faces/edges = ", mesh->size(1));
-      msg("#vertices = "   , mesh->size(dim));
+      msg("Create grid:");
+      msg("#elements = "   , grid->size(0));
+      msg("#faces/edges = ", grid->size(1));
+      msg("#vertices = "   , grid->size(dim));
       msg("");
     }
 
-    void createFeSpaces()
+    void createGlobalBasis()
     {
-      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
+      globalBasis = std::make_shared<GlobalBasis>(new GlobalBasis(leafGridView()));
     }
 
     void createMatricesAndVectors()
     {
-      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
-      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
+      systemMatrix = std::make_shared<SystemMatrixType>(*globalBasis);
+      solution = std::make_shared<SystemVectorType>(*globalBasis, componentNames);
 
       auto rhsNames = std::vector<std::string>(nComponents, "rhs");
-      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
+      rhs = std::make_shared<SystemVectorType>(*globalBasis, rhsNames);
     }
 
     void createSolver()
@@ -364,17 +359,17 @@ namespace AMDiS
     /// vectors, \ref solution.
     std::vector<std::string> componentNames;
 
-    /// Mesh of this problem.
-    std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
+    /// Grid of this problem.
+    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
 
     /// Name of the mesh
-    std::string meshName = "none";
+    std::string gridName = "none";
 
     /// Pointer to the meshes for the different problem components
-    std::vector<Mesh*> componentMeshes;
+    std::vector<Grid*> componentGrids;
 
     /// FE spaces of this problem.
-    std::shared_ptr<FeSpaces> feSpaces; // eventuell const
+    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
 
     /// A FileWriter object
     std::shared_ptr<FileWriter<Traits>> filewriter;
diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp
index 8d90f6c078f6f55ccd3d792e649b321edaf1f2c9..cbcb30ef70f0138cc0340d04b271ccd4d948ec20 100644
--- a/dune/amdis/ProblemStat.inc.hpp
+++ b/dune/amdis/ProblemStat.inc.hpp
@@ -6,15 +6,15 @@
 
 namespace AMDiS {
 
-template <class Traits>
-void ProblemStat<Traits>::initialize(
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::initialize(
     Flag initFlag,
     Self* adoptProblem,
     Flag adoptFlag)
 {
-  // create meshes
-  if (mesh) {
-    warning("mesh already created");
+  // create grides
+  if (grid) {
+    warning("grid already created");
   }
   else {
     if (initFlag.isSet(CREATE_MESH) ||
@@ -27,40 +27,40 @@ void ProblemStat<Traits>::initialize(
         (adoptFlag.isSet(INIT_MESH) ||
         adoptFlag.isSet(INIT_SYSTEM) ||
         adoptFlag.isSet(INIT_FE_SPACE))) {
-      mesh = adoptProblem->getMesh();
+      grid = adoptProblem->getMesh();
     }
 
-    componentMeshes.resize(nComponents, mesh.get());
+    componentMeshes.resize(nComponents, grid.get());
   }
 
-  if (!mesh)
-    warning("no mesh created");
+  if (!grid)
+    warning("no grid created");
 
   int globalRefinements = 0;
-  Parameters::get(meshName + "->global refinements", globalRefinements);
+  Parameters::get(gridName + "->global refinements", globalRefinements);
   if (globalRefinements > 0) {
-    mesh->globalRefine(globalRefinements);
+    grid->globalRefine(globalRefinements);
   }
 
 
   // create fespace
-  if (feSpaces) {
-    warning("feSpaces already created");
+  if (globalBasis) {
+    warning("globalBasis already created");
   }
   else {
     if (initFlag.isSet(INIT_FE_SPACE) ||
         (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
-      createFeSpaces();
+      createGlobalBasis();
     }
 
     if (adoptProblem &&
         (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
-      feSpaces = adoptProblem->getFeSpaces();
+      globalBasis = adoptProblem->getGlobalBasis();
     }
   }
 
-  if (!feSpaces)
-    warning("no feSpaces created\n");
+  if (!globalBasis)
+    warning("no globalBasis created\n");
 
 
   // create system
@@ -103,32 +103,32 @@ void ProblemStat<Traits>::initialize(
 
 // add matrix/vector operator terms
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addMatrixOperator(ElementOperator const& op, int i, int j, double* factor, double* estFactor)
 {
   matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
   matrixOperators[i][j].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addMatrixOperator(IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
 {
   matrixOperators[i][j].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
   matrixOperators[i][j].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
 {
   matrixOperators[i][j].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
   matrixOperators[i][j].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
 {
   for (auto op : ops){
@@ -139,8 +139,8 @@ addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
   }
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
 {
   for (auto op : ops) {
@@ -152,32 +152,32 @@ addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool>
 }
 
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addVectorOperator(ElementOperator const& op, int i, double* factor, double* estFactor)
 {
   rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
   rhsOperators[i].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addVectorOperator(IntersectionOperator const& op, int i, double* factor, double* estFactor)
 {
   rhsOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
   rhsOperators[i].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
 {
   rhsOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
   rhsOperators[i].changing = true;
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addVectorOperator(std::map< int, ElementOperator> ops)
 {
   for (auto op : ops) {
@@ -186,8 +186,8 @@ addVectorOperator(std::map< int, ElementOperator> ops)
   }
 }
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
 {
   for (auto op : ops) {
@@ -198,9 +198,9 @@ addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
 
 
 // Adds a Dirichlet boundary condition
-template <class Traits>
+template <class GlobalBasis>
   template <class Predicate, class Values>
-void ProblemStat<Traits>::
+void ProblemStat<GlobalBasis>::
 addDirichletBC(Predicate const& predicate, int row, int col, Values const& values)
 {
   static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
@@ -219,8 +219,8 @@ addDirichletBC(Predicate const& predicate, int row, int col, Values const& value
 }
 
 
-template <class Traits>
-void ProblemStat<Traits>::
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
 solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
 {
   Timer t;
@@ -259,13 +259,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
 }
 
 
-template <class Traits>
-void ProblemStat<Traits>::
-buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix, bool asmVector)
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
+buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
 {
   Timer t;
 
-  Assembler<FeSpaces> assembler(getFeSpaces(), asmMatrix, asmVector);
+  Assembler<GlobalBasis> assembler(getGlobalBasis(), asmMatrix, asmVector);
 
   auto gridView = leafGridView();
   assembler.assemble(gridView,
@@ -276,9 +276,9 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix, bool asmV
 }
 
 
-template <class Traits>
-void ProblemStat<Traits>::
-writeFiles(AdaptInfo& adaptInfo, bool force)
+template <class GlobalBasis>
+void ProblemStat<GlobalBasis>::
+writeFiles(AdaptInfo& adaptInfo, bool /*force*/)
 {
   Timer t;
   filewriter->write(adaptInfo.getTime(), *solution);
diff --git a/dune/amdis/SecondOrderAssembler.hpp b/dune/amdis/SecondOrderAssembler.hpp
index 8760df11debee7299320281e74ae1876c0289c72..46149865a1bc2dd5c50621ceb8d02671f94f4596 100644
--- a/dune/amdis/SecondOrderAssembler.hpp
+++ b/dune/amdis/SecondOrderAssembler.hpp
@@ -48,7 +48,7 @@ namespace AMDiS
       auto const& quad = Super::getQuadrature().getRule();
 
       // TODO: currently only the implementation for equal fespaces
-      assert( psiDegree == phiDegree );
+      test_exit( psiDegree == phiDegree, "" );
 
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
diff --git a/dune/amdis/linear_algebra/ElementMatrix.hpp b/dune/amdis/linear_algebra/ElementMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..434a7a0d34e5487018f691e65b5510f184b8b357
--- /dev/null
+++ b/dune/amdis/linear_algebra/ElementMatrix.hpp
@@ -0,0 +1,181 @@
+#pragma once
+
+namespace AMDiS
+{
+  namespace Impl
+  {
+    using BaseElementMatrix = mtl::dense2D<double>;
+    using PowerIndex = std::array<std::size_t, 2>;
+
+    template <class RowNode, class ColNode, class RowTag, class ColTag>
+    class ElementMatrix;
+
+    template <class RowNode, class ColNode>
+    class ElementMatrix<RowNode, ColNode, Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag>
+    {
+    public:
+      ElementMatrix() = default;
+
+      ElementMatrix(RowNode const& rowNode, ColNode const& colNode)
+      {
+        resize(rowNode, colNode);
+      }
+
+      void resize(RowNode const& rowNode, ColNode const& colNode)
+      {
+        matrix_.change_dim(rowNode.finiteElement().size(),
+                           colNode.finiteElement().size());
+      }
+
+      double const& operator()(std::size_t i, std::size_t j) const
+      {
+        return matrix_[i][j];
+      }
+
+      double& operator()(std::size_t i, std::size_t j)
+      {
+        return matrix_[i][j];
+      }
+
+      template <class RowIndexSet, class ColIndexSet, class Matrix>
+      void apply(RowIndexSet const& rowIndexSet, ColIndexSet const& colIndexSet, Matrix& mat) const
+      {
+        for (std::size_t i = 0; i < num_rows(matrix_); ++i) {
+          // The global index of the i−th vertex of the element
+          auto const row = rowIndexSet.index(i);
+          for (std::size_t j = 0; j < num_cols(matrix_); ++j) {
+            // The global index of the j−th vertex of the element
+            auto const col = colIndexSet.index(j);
+            mat(row,col) += matrix_(i,j);
+          }
+        }
+      }
+
+    private:
+      BaseElementMatrix matrix_;
+    };
+
+
+    template <class RowNode, class ColNode>
+    class ElementMatrix<RowNode, ColNode, Dune::TypeTree::PowerNodeTag, Dune::TypeTree::LeafNodeTag>
+    {
+      static_assert( RowNode::Child<0>::type::isLeaf,
+        "Requires near-leaf power nodes." );
+
+    public:
+      ElementMatrix() = default;
+
+      ElementMatrix(RowNode const& rowNode, ColNode const& colNode)
+      {
+        resize(rowNode, colNode);
+      }
+
+      void resize(RowNode const& rowNode, ColNode const& colNode)
+      {
+        for (std::size_t i = 0; i < RowNode::CHILDREN; ++i)
+          matrix_[i].change_dim(rowNode.child(0).finiteElement().size(),
+                                colNode.finiteElement().size());
+      }
+
+      double const& operator()(PowerIndex ii, std::size_t j) const
+      {
+        return matrix_[ii[0]][ii[1]][j];
+      }
+
+      double& operator()(PowerIndex ii, std::size_t j)
+      {
+        return matrix_[ii[0]][ii[1]][j];
+      }
+
+    private:
+      std::array<BaseElementMatrix, RowNode::CHILDREN> matrix_;
+    };
+
+
+    template <class RowNode, class ColNode>
+    class ElementMatrix<RowNode, ColNode, Dune::TypeTree::LeafNodeTag, Dune::TypeTree::PowerNodeTag>
+    {
+      static_assert( ColNode::Child<0>::type::isLeaf,
+        "Requires near-leaf power nodes." );
+
+    public:
+      ElementMatrix() = default;
+
+      ElementMatrix(RowNode const& rowNode, ColNode const& colNode)
+      {
+        resize(rowNode, colNode);
+      }
+
+      void resize(RowNode const& rowNode, ColNode const& colNode)
+      {
+        for (std::size_t i = 0; i < ColNode::CHILDREN; ++i)
+          matrix_[i].change_dim(rowNode.finiteElement().size(),
+                                colNode.child(0).finiteElement().size());
+      }
+
+      double const& operator()(std::size_t i, PowerIndex jj) const
+      {
+        return matrix_[jj[0]][i][jj[1]];
+      }
+
+      double& operator()(std::size_t i, PowerIndex jj)
+      {
+        return matrix_[jj[0]][i][jj[1]];
+      }
+
+    private:
+      std::array<BaseElementMatrix, ColNode::CHILDREN> matrix_;
+    };
+
+
+    template <class RowNode, class ColNode>
+    class ElementMatrix<RowNode, ColNode, Dune::TypeTree::PowerNodeTag, Dune::TypeTree::PowerNodeTag>
+    {
+      static_assert( RowNode::Child<0>::type::isLeaf &&
+                     ColNode::Child<0>::type::isLeaf,
+        "Requires near-leaf power nodes." );
+
+    public:
+      ElementMatrix() = default;
+
+      ElementMatrix(RowNode const& rowNode, ColNode const& colNode)
+      {
+        resize(rowNode, colNode);
+      }
+
+      void resize(RowNode const& rowNode, ColNode const& colNode)
+      {
+        for (std::size_t i = 0; i < RowNode::CHILDREN; ++i)
+          for (std::size_t j = 0; j < ColNode::CHILDREN; ++j)
+            matrix_[i][j].change_dim(rowNode.child(0).finiteElement().size(),
+                                     colNode.child(0).finiteElement().size());
+      }
+
+      double const& operator()(PowerIndex ii, PowerIndex jj) const
+      {
+        return matrix_[ii[0]][jj[0]][ii[1]][jj[1]];
+      }
+
+      double& operator()(PowerIndex ii, PowerIndex jj)
+      {
+        return matrix_[ii[0]][jj[0]][ii[1]][jj[1]];
+      }
+
+    private:
+      std::array<std::array<BaseElementMatrix, ColNode::CHILDREN>, RowNode::CHILDREN> matrix_;
+    };
+
+  } // end namespace Impl
+
+
+  template <class RowNode, class ColNode>
+  using ElementMatrix
+    = Impl::ElementMatrix<RowNode, ColNode, typename RowNode::NodeTag, typename ColNode::NodeTag>;
+
+  template <class RowNode, class ColNode>
+  ElementMatrix<RowNode, ColNode> makeElementMatrix(RowNode const& rowNode, ColNode const& colNode)
+  {
+    return {rowNode, colNode};
+  }
+
+} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/ElementVector.hpp b/dune/amdis/linear_algebra/ElementVector.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef63cd14e423b97c882b184dc00ae8d83691763c
--- /dev/null
+++ b/dune/amdis/linear_algebra/ElementVector.hpp
@@ -0,0 +1,94 @@
+#pragma once
+
+namespace AMDiS
+{
+  namespace Impl
+  {
+    using BaseElementVector = mtl::dense_vector<double>;
+    using PowerIndex = std::array<std::size_t, 2>;
+
+    template <class Node, class NodeTag>
+    class ElementVector;
+
+    template <class Node>
+    class ElementVector<Node, Dune::TypeTree::LeafNodeTag>
+    {
+    public:
+      ElementVector() = default;
+
+      ElementVector(Node const& node)
+      {
+        resize(node);
+      }
+
+      void resize(Node const& node)
+      {
+        vector_.change_dim(node.finiteElement().size());
+      }
+
+      double const& operator[](std::size_t i) const { return vector_[i]; }
+      double& operator[](std::size_t i) { return vector_[i]; }
+
+
+      template <class LocalIndexSet, class Vector>
+      void apply(LocalIndexSet const& localIndexSet, Vector& vec) const
+      {
+        for (std::size_t i = 0; i < num_rows(vector_); ++i) {
+          // The global index of the i−th vertex of the element
+          auto const idx = localIndexSet.index(i);
+          vec[idx] += vector_[i];
+        }
+      }
+
+    private:
+      BaseElementVector vector_;
+    };
+
+
+    template <class Node>
+    class ElementVector<Node, Dune::TypeTree::PowerNodeTag>
+    {
+      static_assert( Node::Child<0>::type::isLeaf,
+        "Requires near-leaf power nodes." );
+
+    public:
+      ElementVector() = default;
+
+      ElementVector(Node const& node)
+      {
+        resize(node);
+      }
+
+      void resize(Node const& node)
+      {
+        for (std::size_t i = 0; i < Node::CHILDREN; ++i)
+          vector_[i].change_dim(node.child(0).finiteElement().size());
+      }
+
+      double const& operator()(PowerIndex ii) const
+      {
+        return vector_[ii[0]][ii[1]];
+      }
+
+      double& operator()(PowerIndex ii)
+      {
+        return vector_[ii[0]][ii[1]];
+      }
+
+    private:
+      std::array<BaseElementVector, Node::CHILDREN> vector_;
+    };
+
+  } // end namespace Impl
+
+  template <class Node>
+  using ElementVector = Impl::ElementVector<Node, typename Node::NodeTag>;
+
+  template <class Node>
+  ElementVector<Node> makeElementVector(Node const& node)
+  {
+    return {node};
+  }
+
+
+} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp b/dune/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
index 901e67e663baa24f4da62eab00dbd8ee2677a157..f042e53c52d262d047abb7a05cbe3849db883351 100644
--- a/dune/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
+++ b/dune/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
@@ -9,7 +9,7 @@ namespace AMDiS
 {
 
   template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
-  struct DefaultCreators<PreconditionerInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>>
+  class DefaultCreators<PreconditionerInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>>
   {
     using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
     using PreconBase = PreconditionerInterface<Matrix, Vector>;
@@ -20,6 +20,7 @@ namespace AMDiS
 
     using Map = CreatorMap<PreconBase>;
 
+  public:
     static void init()
     {
       auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
diff --git a/dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp b/dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
index 878caeb3db7efeb96cdd3326890838ec9cb2b4b5..c5e88cbb4796200a2392e731b20fcec9dd59e177 100644
--- a/dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
+++ b/dune/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
@@ -74,46 +74,47 @@ namespace AMDiS
    */
   template <class Matrix>
   using ICPreconditioner = itl::pc::ic_0<Matrix>;
-               
-  
-  
+
+
+
   template <class T, class Param, class Vector>
-  struct DefaultCreators<PreconditionerInterface<mtl::compressed2D<T, Param>, Vector>>
+  class DefaultCreators<PreconditionerInterface<mtl::compressed2D<T, Param>, Vector>>
   {
     using Matrix = mtl::compressed2D<T, Param>;
     using PreconBase = PreconditionerInterface<Matrix, Vector>;
-    
+
     template <template<class> class ITLPrecon>
-    using PreconCreator 
+    using PreconCreator
       = typename Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>::Creator;
-    
+
     using Map = CreatorMap<PreconBase>;
-  
+
+  public:
     static void init()
     {
       auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
       Map::addCreator("diag", pc_diag);
       Map::addCreator("jacobi", pc_diag);
-      
+
       auto pc_mass = new PreconCreator<MassLumpingPreconditioner>;
       Map::addCreator("masslumping", pc_mass);
-      
+
       auto pc_ilu = new PreconCreator<ILUPreconditioner>;
       Map::addCreator("ilu", pc_ilu);
       Map::addCreator("ilu0", pc_ilu);
-      
+
       auto pc_ic = new PreconCreator<ICPreconditioner>;
       Map::addCreator("ic", pc_ic);
       Map::addCreator("ic0", pc_ic);
-      
+
       auto pc_id = new PreconCreator<IdentityPreconditioner>;
       Map::addCreator("identity", pc_id);
       Map::addCreator("no", pc_id);
-      
+
       Map::addCreator("default", pc_id);
     }
   };
-  
+
 
   template <class Matrix, class Vector>
   itl::pc::solver<PreconditionerInterface<Matrix, Vector>, Vector, false>
@@ -128,5 +129,5 @@ namespace AMDiS
   {
     return {P, vin};
   }
-  
+
 } // namespace AMDiS
diff --git a/dune/amdis/utility/GetEntity.hpp b/dune/amdis/utility/GetEntity.hpp
index 9f345a2df2816ee827c036664e85de6776f5f402..879ab752a841d94d8fa5e8903f1883521258200b 100644
--- a/dune/amdis/utility/GetEntity.hpp
+++ b/dune/amdis/utility/GetEntity.hpp
@@ -25,7 +25,7 @@ namespace AMDiS
     template <class I>
     struct Get<I, Void_t<decltype(std::declval<I>().inside())>>
     {
-      static_assert( int(I::mydimension) < int(I::dimension), "" );
+      static_assert( int(I::mydimension) < int(I::Entity::dimension), "" );
 
       using Entity = typename I::Entity;
       using Geometry = typename I::LocalGeometry;
diff --git a/dune/amdis/utility/Visitor.hpp b/dune/amdis/utility/Visitor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..010b4c06d04c3a065e2bbdabc120519a35e967e5
--- /dev/null
+++ b/dune/amdis/utility/Visitor.hpp
@@ -0,0 +1,92 @@
+#pragma once
+
+#include <dune/typetree/visitor.hh>
+
+namespace AMDiS
+{
+  namespace filter
+  {
+    // all interior nodes
+    struct All
+    {
+      template <class Fct, class... Args>
+      void apply(Fct&&, Args&&...) const {};
+    };
+
+    // restrict interior nodes to PowerNodes with direct leaf childs
+    struct PowerNodes
+    {
+      template <class Fct, class... Args>
+      void apply(Fct&&, Args&&...) const {};
+
+      template <class Fct, class T, class TreePath>
+      void apply(Fct&& fct, T&& t, TreePath treePath) const
+      {
+        using Node = std::decay_t<T>;
+        apply(std::forward<Fct>(fct), std::forward<T>(t), treePath, Node::NodeTag{});
+      };
+
+      template <class Fct, class T, class TreePath>
+      void apply(Fct&& fct, T&& t, TreePath treePath, Dune::TypeTree::PowerNodeTag tag) const
+      {
+        using Node = std::decay_t<T>;
+        apply(std::forward<Fct>(fct), std::forward<T>(t), treePath, tag, bool_<Node::Child<0>::isLeaf>);
+      }
+
+      template <class Fct, class T, class TreePath, bool ChildIsLeaf>
+      void apply(Fct&& fct, T&& t, TreePath treePath, Dune::TypeTree::PowerNodeTag, bool_t<true> /*childIsLeaf*/) const
+      {
+        fct(std::forward<T>(t), treePath);
+      }
+    };
+
+  } // end namespace filter
+
+
+  template <class Filter, class Fct>
+  class FilteredVisitor
+      : public Dune::TypeTree::TreeVisitor
+  {
+  public:
+    FilteredVisitor(Filter const& filter, Fct const& fct)
+      : filter_(filter)
+      , fct_(fct)
+    {}
+
+    /// Method for infix tree traversal.
+    template<typename T, typename TreePath>
+    void in(T&& t, TreePath treePath) const
+    {
+      filter_.apply(fct_, std::forward<T>(t), treePath);
+    }
+
+    /// Method for leaf traversal.
+    template<typename T, typename TreePath>
+    void leaf(T&& t, TreePath treePath) const
+    {
+      fct_(std::forward<T>(t), treePath);
+    }
+
+  private:
+    Filter filter_;
+    Fct fct_;
+  };
+
+
+  template <class Tree, class Fct>
+  void forEachNode(Tree&& tree, Fct&& fct)
+  {
+    using Visitor = FilteredVisitor<filter::PowerNodes, std::decay_t<Fct>>;
+    Dune::TypeTree::applyToTree(std::forward<Tree>(tree),
+      Visitor{filter::PowerNodes{}, std::forward<Fct>(fct)});
+  }
+
+  template <class Tree, class Filter, class Fct>
+  void forEachNode(Tree&& tree, Filter&& filter, Fct&& fct)
+  {
+    using Visitor = FilteredVisitor<std::decay_t<Filter>, std::decay_t<Fct>>;
+    Dune::TypeTree::applyToTree(std::forward<Tree>(tree),
+      Visitor{std::forward<Filter>(filter), std::forward<Fct>(fct)});
+  }
+
+} // end namespace AMDiS