diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 6c1c83620c32b9cb1e1e295e9f038eada426a96f..59d521695cbce90acacd93b6b22ae2d2dea681ae 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -15,11 +15,14 @@
 namespace AMDiS
 {
   /// Implementation of interface \ref AssemblerBase
-  template <class LocalContext, class Operator, class... Nodes>
+  /**
+   * \tparam LC  LocalContext
+   **/
+  template <class LC, class Operator, class... Nodes>
   class Assembler
-      : public AssemblerInterface<LocalContext, Nodes...>
+      : public AssemblerInterface<LC, Nodes...>
   {
-    using Super = AssemblerInterface<LocalContext, Nodes...>;
+    using Super = AssemblerInterface<LC, Nodes...>;
 
   private:
 
@@ -74,12 +77,11 @@ namespace AMDiS
      * \ref calculateElementVector or \ref calculateElementMatrix on the
      * vector or matrix operator, respectively.
      **/
-    void assemble(LocalContext const& localContext,
-                  Nodes const&... nodes,
+    void assemble(LC const& localContext, Nodes const&... nodes,
                   ElementMatrixVector& elementMatrixVector) final
     {
-      ContextGeometry<LocalContext> context{localContext, element(), geometry()};
-      assembleImpl(context, nodes..., elementMatrixVector);
+      ContextGeometry<LC> contextGeo{localContext, element(), geometry()};
+      assembleImpl(contextGeo, nodes..., elementMatrixVector);
     }
 
 
@@ -88,21 +90,17 @@ namespace AMDiS
   protected: // implementation detail
 
     // matrix assembling
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void assembleImpl(Context const& context,
-                      RowNode const& rowNode, ColNode const& colNode,
-                      ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void assembleImpl(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      op_->calculateElementMatrix(context, rowNode, colNode, elementMatrix);
+      op_->calculateElementMatrix(contextGeo, rowNode, colNode, elementMatrix);
     }
 
     // vector assembling
-    template <class Context, class Node, class ElementVector>
-    void assembleImpl(Context const& context,
-                      Node const& node,
-                      ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void assembleImpl(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
-      op_->calculateElementVector(context, node, elementVector);
+      op_->calculateElementVector(contextGeo, node, elementVector);
     }
 
 #endif // DOXYGEN
@@ -134,11 +132,11 @@ namespace AMDiS
   };
 
 
-  /// Generate a \ref Assembler on a given `LocalContext` (element or intersection)
-  template <class LocalContext, class Operator, class... Nodes>
+  /// Generate a \ref Assembler on a given LocalContext `LC` (element or intersection)
+  template <class LC, class Operator, class... Nodes>
   auto makeAssembler(Operator&& op, Nodes const&...)
   {
-    return Assembler<LocalContext, Underlying_t<Operator>, Nodes...>{FWD(op)};
+    return Assembler<LC, Underlying_t<Operator>, Nodes...>{FWD(op)};
   }
 
 } // end namespace AMDiS
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index 3946f17f9c552ecf19d348ee507020b7a46e9b8c..a181ffe1d722b2c146d2e128967c86aa0ed074bc 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -128,8 +128,8 @@ namespace AMDiS
     NodeDataTransferContainer nodeDataTransfer_;
   };
 
-  template <class Container, class Basis>
-  void DataTransfer<Container, Basis>::preAdapt(Container const& coeff, bool mightCoarsen)
+  template <class C, class B>
+  void DataTransfer<C,B>::preAdapt(C const& coeff, bool mightCoarsen)
   {
     auto gv = basis_->gridView();
     auto lv = basis_->localView();
@@ -217,8 +217,8 @@ namespace AMDiS
     }
   }
 
-  template <class Container, class Basis>
-  void DataTransfer<Container, Basis>::postAdapt(Container& coeff, bool refined)
+  template <class C, class B>
+  void DataTransfer<C,B>::postAdapt(C& coeff, bool refined)
   {
     coeff.resize(*basis_);
     auto gv = basis_->gridView();
@@ -370,9 +370,9 @@ namespace AMDiS
     NodeElementData fatherDOFsTemp_;
   };
 
-  template <class Node, class Container, class Basis>
+  template <class Node, class C, class B>
     template <class Trafo>
-  bool NodeDataTransfer<Node, Container, Basis>::
+  bool NodeDataTransfer<Node,C,B>::
     restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
                   NodeElementData const& childDOFs, bool init)
   {
@@ -426,9 +426,9 @@ namespace AMDiS
                            std::logical_and<bool>());
   }
 
-  template <class Node, class Container, class Basis>
+  template <class Node, class C, class B>
     template <class Trafo>
-  void NodeDataTransfer<Node, Container, Basis>::
+  void NodeDataTransfer<Node,C,B>::
     prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
                  Trafo const& trafo, bool init)
   {
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index 10e932b947ec299bdbc3f6687ef26696ad5597dc..a74b0eabd797d0c4778822dd78757734fc89ba9c 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -53,22 +53,28 @@ namespace AMDiS
       , values_(FWD(values))
     {}
 
-    template <class Intersection>
-    bool onBoundary(Intersection const& intersection) const
+    template <class IS>
+    bool onBoundary(IS const& intersection) const
     {
       return predicate_ ? (*predicate_)(intersection.geometry().center())
                         : Super::onBoundary(intersection);
     }
 
     // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
-    template <class RowBasis, class ColBasis>
-    void init(RowBasis const& rowBasis, ColBasis const& colBasis);
+    template <class RB, class CB>
+    void init(RB const& rowBasis, CB const& colBasis);
 
     /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
     /// to the matrix and optionally delete the corresponding matrix-column.
-    template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
-    void fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
-                               RowNode const& rowNode, ColNode const& colNode);
+    /**
+     * \tparam Mat  Matrix
+     * \tparam Sol  Vector of solution
+     * \tparam Rhs  Vector of rhs
+     * \tparam RN   RowNode
+     * \tparam CN   ColNode
+     **/
+    template <class Mat, class Sol, class Rhs, class RN, class CN>
+    void fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, CN const& colNode);
 
   private:
     Dune::Std::optional<std::function<bool(Domain)>> predicate_;
@@ -82,8 +88,8 @@ namespace AMDiS
   {
     using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
 
-    template <class RowNode, class ColNode>
-    using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RowNode>>> >;
+    template <class RN, class CN>
+    using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RN>>> >;
   };
 
   template <class RowBasis, class ColBasis>
diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp
index cdbed4cccfe0a5343593dab7125c039667f94269..69f873ccb63a5a74514371ea0dd7a14295189f00 100644
--- a/src/amdis/DirichletBC.inc.hpp
+++ b/src/amdis/DirichletBC.inc.hpp
@@ -13,29 +13,28 @@
 namespace AMDiS {
 
 template <class D, class R>
-  template <class RowBasis, class ColBasis>
+  template <class RB, class CB>
 void DirichletBC<D,R>::
-init(RowBasis const& rowBasis, ColBasis const& colBasis)
+init(RB const& rowBasis, CB const& colBasis)
 {
   using Dune::Functions::forEachBoundaryDOF;
-  using LV = typename ColBasis::LocalView;
-  using I = typename ColBasis::GridView::Intersection;
+  using LV = typename CB::LocalView;
+  using IS = typename CB::GridView::Intersection;
   dirichletNodes_.resize(colBasis.dimension());
-  forEachBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, I const& intersection) {
+  forEachBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, IS const& intersection) {
     dirichletNodes_[localView.index(localIndex)] = onBoundary(intersection);
   });
 }
 
 
 template <class D, class R>
-  template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+  template <class Mat, class Sol, class Rhs, class RN, class CN>
 void DirichletBC<D,R>::
-fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
-                      RowNode const& rowNode, ColNode const& colNode)
+fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, CN const& colNode)
 {
-  auto columns = Constraints<Matrix>::dirichletBC(matrix, dirichletNodes_);
+  auto columns = Constraints<Mat>::dirichletBC(matrix, dirichletNodes_);
 
-  Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, R>{},
+  Dune::Hybrid::ifElse(std::is_same<RangeType_t<RN>, R>{},
     [&](auto id) {
       auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
       Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
@@ -45,7 +44,7 @@ fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
   for (auto const& triplet : columns)
     rhs[triplet.row] -= triplet.value * solution[triplet.col];
 
-  Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, R>{},
+  Dune::Hybrid::ifElse(std::is_same<RangeType_t<CN>, R>{},
     [&](auto id) {
       auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
       Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
diff --git a/src/amdis/GridFunctionOperator.hpp b/src/amdis/GridFunctionOperator.hpp
index 13fcc88641c7196637a093c7151a0778799256de..d427cd0e96826286f1af58dbafa2d98ec522499d 100644
--- a/src/amdis/GridFunctionOperator.hpp
+++ b/src/amdis/GridFunctionOperator.hpp
@@ -26,28 +26,30 @@ namespace AMDiS
    *
    * The class is specialized, by deriving from it, in \ref GridFunctionOperator.
    *
-   * \tparam Derived      The class derived from GridFunctionOperatorBase
-   * \tparam LocalContext The Element or Intersection type
-   * \tparam GridFunction The GridFunction, a LocalFunction is created from, and
-   *                      that is evaluated at quadrature points.
+   * \tparam Derived  The class derived from GridFunctionOperatorBase
+   * \tparam LC       The Element or Intersection type
+   * \tparam GF       The GridFunction, a LocalFunction is created from, and
+   *                  that is evaluated at quadrature points.
    *
    * **Requirements:**
-   * - `LocalContext` models either Entity (of codim 0) or Intersection
-   * - `GridFunction` models the \ref Concepts::GridFunction
+   * - `LC` models either Entity (of codim 0) or Intersection
+   * - `GF` models the \ref Concepts::GridFunction
    **/
-  template <class Derived, class LocalContext, class GridFunction>
+  template <class Derived, class LC, class GF>
   class GridFunctionOperatorBase
-      : public LocalOperator<Derived, LocalContext>
+      : public LocalOperator<Derived, LC>
   {
-    template <class D, class LC>
+    template <class, class>
     friend class LocalOperator;
 
-    using ContextType = Impl::ContextTypes<LocalContext>;
-    using Super = LocalOperator<Derived, LocalContext>;
+    using ContextType = Impl::ContextTypes<LC>;
+    using Super = LocalOperator<Derived, LC>;
 
   private:
+    using GridFunction = GF;
+
     /// The type of the localFunction associated with the GridFunction
-    using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
+    using LocalFunction = decltype(localFunction(std::declval<GF>()));
 
     /// The Codim=0 entity of the grid, the localFunction can be bound to
     using Element = typename ContextType::Entity;
@@ -56,10 +58,10 @@ namespace AMDiS
     using Geometry = typename Element::Geometry;
 
     /// The type of the local coordinates in the \ref Element
-    using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate;
+    using LocalCoordinate = typename GF::EntitySet::LocalCoordinate;
 
     /// A factory for QuadratureRules that incooperate the order of the LocalFunction
-    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>;
+    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LC::mydimension, LocalFunction>;
 
   public:
     /// \brief Constructor. Stores a copy of `gridFct`.
@@ -68,19 +70,20 @@ namespace AMDiS
      * differentiation order of the operator, to calculate the
      * quadrature degree in \ref getDegree.
      **/
-    template <class GF>
-    GridFunctionOperatorBase(GF&& gridFct, int termOrder)
+    template <class GridFct>
+    GridFunctionOperatorBase(GridFct&& gridFct, int termOrder)
       : gridFct_(FWD(gridFct))
       , termOrder_(termOrder)
     {}
 
     /// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
-    template <class PreQuadFactory>
-    void setQuadFactory(PreQuadFactory&& pre)
+    /// \tparam  PQF  A PreQuadratureFactory
+    template <class PQF>
+    void setQuadFactory(PQF&& pre)
     {
       using ctype = typename Geometry::ctype;
       quadFactory_.emplace(
-        makeQuadratureFactory<ctype, LocalContext::mydimension, LocalFunction>(FWD(pre)));
+        makeQuadratureFactory<ctype, LC::mydimension, LocalFunction>(FWD(pre)));
     }
 
   protected:
@@ -147,7 +150,7 @@ namespace AMDiS
   class GridFunctionOperatorTransposed
       : public LocalOperator<Derived, typename Transposed::LocalContext>
   {
-    template <class D, class LC>
+    template <class, class>
     friend class LocalOperator;
 
     template <class T, class... Args>
@@ -161,8 +164,9 @@ namespace AMDiS
     {}
 
     /// Redirects the setQuadFactory call top the transposed operator
-    template <class PreQuadFactory>
-    void setQuadFactory(PreQuadFactory&& pre)
+    /// \tparam  PQF  A PreQuadratureFactory
+    template <class PQF>
+    void setQuadFactory(PQF&& pre)
     {
       transposedOp_.setQuadFactory(FWD(pre));
     }
@@ -182,15 +186,18 @@ namespace AMDiS
     }
 
     /// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    /**
+     * \tparam CG  ContextGeometry
+     * \tparam RN  RowNode
+     * \tparam CN  ColNode
+     * \tparam Mat ElementMatrix
+     **/
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeometry, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
       auto elementMatrixTransposed = transposed(elementMatrix);
       transposedOp_.getElementMatrix(
-        context, colNode, rowNode, elementMatrixTransposed);
+        contextGeometry, colNode, rowNode, elementMatrixTransposed);
     }
 
   private:
@@ -199,12 +206,12 @@ namespace AMDiS
 
 
 #ifndef DOXYGEN
-  template <class Tag, class PreGridFct, class PreQuadFactory>
+  template <class Tag, class PreGridFct, class PQF>
   struct PreGridFunctionOperator
   {
     Tag tag;
     PreGridFct expr;
-    PreQuadFactory quadFactory;
+    PQF quadFactory;
   };
 #endif
 
@@ -228,12 +235,12 @@ namespace AMDiS
    * other effect.
    *
    * \tparam Tag  An Identifier for this GridFunctionOperator
-   * \tparam LocalContext  An Element or Intersection the operator is evaluated on
+   * \tparam LC  An Element or Intersection the operator is evaluated on
    * \tparam GridFct  A GridFunction evaluated in local coordinates on the bound element
    **/
-  template <class Tag, class LocalContext, class GridFct>
+  template <class Tag, class LC, class GridFct>
   class GridFunctionOperator
-      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LocalContext, GridFct>, LocalContext, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LC, GridFct>, LC, GridFct>
   {};
 
   template <class Context, class Tag, class GF, class QF>
diff --git a/src/amdis/GridTransferManager.hpp b/src/amdis/GridTransferManager.hpp
index 00f2f64dada9374218ea1442570a1eee31cde518..5721bf3ba9e176a7fad427c51390ae0cae3ff1ff 100644
--- a/src/amdis/GridTransferManager.hpp
+++ b/src/amdis/GridTransferManager.hpp
@@ -17,7 +17,7 @@ namespace AMDiS
   /// Static administration class for automatic handling of DOFVectors during grid adaption
   class GridTransferManager
   {
-    template <class T1, class T2>
+    template <class, class>
     friend class DOFVectorBase;
 
     using Self = GridTransferManager;
diff --git a/src/amdis/Integrate.hpp b/src/amdis/Integrate.hpp
index 2a9522689838189068d18e4d57af1a0ce39b706c..a94e9781250b739ab7eba6f596b3c3408f2139ab 100644
--- a/src/amdis/Integrate.hpp
+++ b/src/amdis/Integrate.hpp
@@ -8,8 +8,8 @@ namespace AMDiS
 {
   namespace Impl
   {
-    template <class GF, class GridView, class QuadProvider>
-    auto integrateImpl(GF&& gf, GridView const& gridView, QuadProvider makeQuad)
+    template <class GF, class GV, class QP>
+    auto integrateImpl(GF&& gf, GV const& gridView, QP makeQuad)
     {
       auto localFct = localFunction(FWD(gf));
 
@@ -31,8 +31,8 @@ namespace AMDiS
       return result;
     }
 
-    template <class GF, class GridView, class QuadProvider>
-    auto integrateImpl(GF&& gf, GridView const& gv, QuadProvider makeQuad, std::true_type)
+    template <class GF, class GV, class QP>
+    auto integrateImpl(GF&& gf, GV const& gv, QP makeQuad, std::true_type)
     {
       return integrateImpl(FWD(gf), gv, makeQuad);
     }
diff --git a/src/amdis/LocalOperator.hpp b/src/amdis/LocalOperator.hpp
index ad3772674336649ca838b51961db931a69d26d9e..5cf504149e979ffda8523465e7d92be4dfe1d59f 100644
--- a/src/amdis/LocalOperator.hpp
+++ b/src/amdis/LocalOperator.hpp
@@ -73,34 +73,29 @@ namespace AMDiS
 
     /// \brief Assemble a local element matrix on the element that is bound.
     /**
-     * \param context Container for geometry related data, \ref ContextGeometry
+     * \param contextGeo Container for geometry related data, \ref ContextGeometry
      * \param rowNode The node of the test-basis
      * \param colNode The node of the trial-basis
      * \param elementMatrix The output element matrix
      **/
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void calculateElementMatrix(Context const& context,
-                                RowNode const& rowNode,
-                                ColNode const& colNode,
-                                ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void calculateElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
       assert( bound_ );
-      derived().getElementMatrix(context, rowNode, colNode, elementMatrix);
+      derived().getElementMatrix(contextGeo, rowNode, colNode, elementMatrix);
     }
 
     /// \brief Assemble a local element vector on the element that is bound.
     /**
-     * \param context Container for geometry related data \ref ContextGeometry
+     * \param contextGeo Container for geometry related data \ref ContextGeometry
      * \param node The node of the test-basis
      * \param elementVector The output element vector
      **/
-    template <class Context, class Node, class ElementVector>
-    void calculateElementVector(Context const& context,
-                                Node const& node,
-                                ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void calculateElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       assert( bound_ );
-      derived().getElementVector(context, node, elementVector);
+      derived().getElementVector(contextGeo, node, elementVector);
     }
 
 
@@ -121,33 +116,17 @@ namespace AMDiS
     // default implementation. Can be overridden in the derived classes
     void unbind_impl() {}
 
-    // default implementation called by \ref calculateElementMatrix
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& /*context*/,
-                          RowNode const& /*rowNode*/,
-                          ColNode const& /*colNode*/,
-                          ElementMatrix& /*elementMatrix*/)
-    {
-      error_exit("Needs to be implemented by derived class!");
-    }
-
-    // default implementation called by \ref calculateElementVector
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& /*context*/,
-                          Node const& /*node*/,
-                          ElementVector& /*elementVector*/)
-    {
-      error_exit("Needs to be implemented by derived class!");
-    }
-
     /// \brief Return the quadrature degree for a matrix operator.
     /**
      * The quadrature degree that is necessary, to integrate the expression
      * multiplied with (possibly derivatives of) shape functions. Thus it needs
      * the order of derivatives, this operator implements.
+     *
+     * \tparam RN  RowNode
+     * \tparam CN  ColNode
      **/
-    template <class RowNode, class ColNode>
-    int getDegree(int order, int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
+    template <class RN, class CN>
+    int getDegree(int order, int coeffDegree, RN const& rowNode, CN const& colNode) const
     {
       assert( bound_ );
       test_warning(coeffDegree >= 0,
@@ -198,8 +177,13 @@ namespace AMDiS
 
 
   /// Generate a \ref LocalOperator from a PreOperator.
-  template <class Derived, class Context, class GridView>
-  auto makeLocalOperator(LocalOperator<Derived, Context> const& localOp, GridView const& /*gridView*/)
+  /**
+   * \tparam Derived  Implementation of LocalOperator
+   * \tparam LC       LocalContext
+   * \tparam GV       GridView
+   **/
+  template <class Derived, class LC, class GV>
+  auto makeLocalOperator(LocalOperator<Derived, LC> const& localOp, GV const& /*gridView*/)
   {
     return localOp.derived();
   }
diff --git a/src/amdis/OperatorList.hpp b/src/amdis/OperatorList.hpp
index 35adb2c99e0960f9000f46b44b2835f9714c45b2..bf780920a3f2a977174ec7d84725a187fda0d532 100644
--- a/src/amdis/OperatorList.hpp
+++ b/src/amdis/OperatorList.hpp
@@ -26,10 +26,10 @@ namespace AMDiS
     using Element = typename GridView::template Codim<0>::Entity;
     using Intersection = typename GridView::Intersection;
 
-    template <class OperatorType>
+    template <class Op>
     struct DataElement
     {
-      std::shared_ptr<OperatorType> op;
+      std::shared_ptr<Op> op;
       BoundaryCondition bc;
     };
 
diff --git a/src/amdis/PeriodicBC.hpp b/src/amdis/PeriodicBC.hpp
index 65541856963d02c5d426a4c0560f8e4d91f66bbd..4b8229e52a7ea384fcc6cb005caed50c2320550e 100644
--- a/src/amdis/PeriodicBC.hpp
+++ b/src/amdis/PeriodicBC.hpp
@@ -40,6 +40,10 @@ namespace AMDiS
 
 
   /// Implements a periodic boundary condition
+  /**
+   * \tparam D  Domain
+   * \tparam MI Type of the MultiIndex
+   **/
   template <class D, class MI>
   class PeriodicBC
       : public BoundaryCondition
@@ -58,11 +62,11 @@ namespace AMDiS
       , faceTrafo_(std::move(faceTrafo))
     {}
 
-    template <class RowBasis, class ColBasis>
-    void init(RowBasis const& rowBasis, ColBasis const& colBasis);
+    template <class RB, class CB>
+    void init(RB const& rowBasis, CB const& colBasis);
 
-    template <class Matrix, class X, class B, class RN, class CN>
-    void fillBoundaryCondition(Matrix& A, X& x, B& b, RN const& rowNode, CN const& colNode);
+    template <class Mat, class Sol, class Rhs, class RN, class CN>
+    void fillBoundaryCondition(Mat& A, Sol& x, Rhs& b, RN const& rowNode, CN const& colNode);
 
     auto const& associations() const
     {
diff --git a/src/amdis/PeriodicBC.inc.hpp b/src/amdis/PeriodicBC.inc.hpp
index a0c03627349d0686ccc042fd3daddfd5745400f1..2b449f3bd4a108ddbf91ce4f259624831802e891 100644
--- a/src/amdis/PeriodicBC.inc.hpp
+++ b/src/amdis/PeriodicBC.inc.hpp
@@ -11,9 +11,9 @@
 namespace AMDiS {
 
 template <class D, class MI>
-  template <class Basis, class ColBasis>
+  template <class RB, class CB>
 void PeriodicBC<D,MI>::
-init(Basis const& basis, ColBasis const& /*colBasis*/)
+init(RB const& basis, CB const& /*colBasis*/)
 {
   if (!bool(hasConnectedIntersections_)) {
     hasConnectedIntersections_ = true;
@@ -231,12 +231,11 @@ coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
 
 
 template <class D, class MI>
-  template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+  template <class Mat, class Sol, class Rhs, class RN, class CN>
 void PeriodicBC<D,MI>::
-fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
-                      RowNode const& rowNode, ColNode const& colNode)
+fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, CN const& colNode)
 {
-  Constraints<Matrix>::periodicBC(matrix, periodicNodes_, associations_);
+  Constraints<Mat>::periodicBC(matrix, periodicNodes_, associations_);
 
   for (auto const& a : associations_) {
     rhs[a.second] += rhs[a.first];
diff --git a/src/amdis/gridfunctions/AnalyticGridFunction.hpp b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
index c0994d4f947187c9bba1fd3d8224908172c8466b..1c48728c46705f66910062fb8a6582b4c9bbfd09 100644
--- a/src/amdis/gridfunctions/AnalyticGridFunction.hpp
+++ b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
@@ -16,8 +16,8 @@ namespace AMDiS
   class AnalyticLocalFunction;
 #endif
 
-  template <class R, class D, class LocalContext, class Function>
-  class AnalyticLocalFunction<R(D), LocalContext, Function>
+  template <class R, class D, class LC, class Function>
+  class AnalyticLocalFunction<R(D), LC, Function>
   {
   public:
     /// The LocalDomain this LocalFunction can be evaluated in
@@ -30,7 +30,7 @@ namespace AMDiS
     enum { hasDerivative = true };
 
   private:
-    using Geometry = typename LocalContext::Geometry;
+    using Geometry = typename LC::Geometry;
 
   public:
     /// \brief Constructor. stores the function `fct`.
@@ -39,7 +39,7 @@ namespace AMDiS
     {}
 
     /// \brief Create a geometry object from the element
-    void bind(LocalContext const& element)
+    void bind(LC const& element)
     {
       geometry_.emplace(element.geometry());
     }
@@ -79,9 +79,9 @@ namespace AMDiS
    * **Requirements:**
    * - The functor `F` must fulfill the concept \ref Concepts::HasFunctorOrder
    **/
-  template <class Sig, class LocalContext, class F,
+  template <class Sig, class LC, class F,
     REQUIRES(Concepts::HasFunctorOrder<F,1>)>
-  int order(AnalyticLocalFunction<Sig,LocalContext,F> const& lf)
+  int order(AnalyticLocalFunction<Sig,LC,F> const& lf)
   {
     return order(lf.fct(),1);
   }
@@ -97,8 +97,8 @@ namespace AMDiS
    * **Requirements:**
    * - The functor `F` must fulfill the concept \ref Concepts::HasPartial
    **/
-  template <class R, class D, class LocalContext, class F>
-  auto derivative(AnalyticLocalFunction<R(D),LocalContext,F> const& lf)
+  template <class R, class D, class LC, class F>
+  auto derivative(AnalyticLocalFunction<R(D),LC,F> const& lf)
   {
     static_assert(Concepts::HasPartial<F>,
       "No partial(...,_0) defined for Functor F of AnalyticLocalFunction.");
@@ -107,7 +107,7 @@ namespace AMDiS
 
     using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
     using DerivativeSignature = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range(D);
-    return AnalyticLocalFunction<DerivativeSignature,LocalContext,decltype(df)>{df};
+    return AnalyticLocalFunction<DerivativeSignature,LC,decltype(df)>{df};
   }
 
 
diff --git a/src/amdis/gridfunctions/ConstantGridFunction.hpp b/src/amdis/gridfunctions/ConstantGridFunction.hpp
index 936f02edc6109a3157d09c3c67259b3990db8353..ed1ed7973ad895b4535084d81bb3e275ed86c1f7 100644
--- a/src/amdis/gridfunctions/ConstantGridFunction.hpp
+++ b/src/amdis/gridfunctions/ConstantGridFunction.hpp
@@ -20,8 +20,8 @@ namespace AMDiS
 #endif
 
   /// \brief LocalFunction of a Gridfunction returning a constant value.
-  template <class R, class D, class LocalContext, class T>
-  class ConstantLocalFunction<R(D), LocalContext, T>
+  template <class R, class D, class LC, class T>
+  class ConstantLocalFunction<R(D), LC, T>
   {
   public:
     /// The LocalDomain this LocalFunction can be evaluated in
@@ -34,7 +34,7 @@ namespace AMDiS
     enum { hasDerivative = true };
 
   private:
-    using Geometry = typename LocalContext::Geometry;
+    using Geometry = typename LC::Geometry;
 
   public:
     /// \brief Constructor. Stores the constant value.
@@ -42,7 +42,7 @@ namespace AMDiS
       : value_(value)
     {}
 
-    void bind(LocalContext const& /*element*/) { /* do nothing */ }
+    void bind(LC const& /*element*/) { /* do nothing */ }
     void unbind() { /* do nothing */ }
 
     /// Return the constant `value_`.
@@ -58,7 +58,7 @@ namespace AMDiS
       using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
       using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
       DerivativeRange diff(0);
-      return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
+      return ConstantLocalFunction<DerivativeRange(D),LC,DerivativeRange>{diff};
     }
 
     /// Return the constant polynomial order 0.
@@ -72,15 +72,15 @@ namespace AMDiS
   };
 
   /// \relates ConstantLocalFunction
-  template <class R, class D, class LocalContext, class T>
-  auto derivative(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
+  template <class R, class D, class LC, class T>
+  auto derivative(ConstantLocalFunction<R(D), LC, T> const& lf)
   {
     return lf.derivative();
   }
 
   /// \relates ConstantLocalFunction
-  template <class R, class D, class LocalContext, class T>
-  int order(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
+  template <class R, class D, class LC, class T>
+  int order(ConstantLocalFunction<R(D), LC, T> const& lf)
   {
     return 0;
   }
diff --git a/src/amdis/linearalgebra/Constraints.hpp b/src/amdis/linearalgebra/Constraints.hpp
index f539113b9b32a3aaa4ede4e788351bfdfec5857d..7dc6a0617805e4ef1eef5bf7c875b3d1e9f801db 100644
--- a/src/amdis/linearalgebra/Constraints.hpp
+++ b/src/amdis/linearalgebra/Constraints.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <amdis/linearalgebra/Common.hpp>
 #include <amdis/linearalgebra/DOFMatrixBase.hpp>
 
 namespace AMDiS
@@ -7,15 +8,34 @@ namespace AMDiS
   template <class Matrix>
   struct Constraints
   {
-    template <class R, class C, class B, class BitVector>
-    static auto dirichletBC(DOFMatrixBase<R,C,B>& matrix, BitVector const& nodes, bool setDiagonal = true)
+    template <class BitVec>
+    static auto dirichletBC(Matrix& matrix, BitVec const& nodes, bool setDiagonal = true)
+    {
+      /* do nothing */
+      return std::array<Triplet<typename Matrix::value_type>,0>{};
+    }
+
+    template <class BitVec, class Assoc>
+    static auto periodicBC(Matrix& matrix, BitVec const& left, Assoc const& association, bool setDiagonal = true)
+    {
+      /* do nothing */
+      return std::array<Triplet<typename Matrix::value_type>,0>{};
+    }
+  };
+
+  template <class RB, class CB, class Backend>
+  struct Constraints<DOFMatrixBase<RB,CB,Backend>>
+  {
+    using Matrix = DOFMatrixBase<RB,CB,Backend>;
+
+    template <class BitVec>
+    static auto dirichletBC(Matrix& matrix, BitVec const& nodes, bool setDiagonal = true)
     {
       return Constraints<typename Matrix::BaseMatrix>::dirichletBC(matrix.matrix(), nodes, setDiagonal);
     }
 
-    template <class R, class C, class B, class BitVector, class Association>
-    static auto periodicBC(DOFMatrixBase<R,C,B>& matrix, BitVector const& left, Association const& association,
-      bool setDiagonal = true)
+    template <class BitVec, class Assoc>
+    static auto periodicBC(Matrix& matrix, BitVec const& left, Assoc const& association, bool setDiagonal = true)
     {
       return Constraints<typename Matrix::BaseMatrix>::periodicBC(matrix.matrix(), left, association, setDiagonal);
     }
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.hpp b/src/amdis/linearalgebra/DOFMatrixBase.hpp
index 116032d57dee58a52ecf1256936a16065919b15c..620d92d033ae298508c8ce64e2e7e3c20d2be104 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.hpp
@@ -20,16 +20,16 @@ namespace AMDiS
    * \tparam CB  Basis of matrix columns
    * \tparam Backend  A linear-algebra backend for the matrix storage
    **/
-  template <class RowBasisType, class ColBasisType, class Backend>
+  template <class RB, class CB, class Backend>
   class DOFMatrixBase
   {
   public:
     /// The type of the finite element space / basis of the row
-    using RowBasis = RowBasisType;
+    using RowBasis = RB;
     using RowLocalView = typename RowBasis::LocalView;
 
     /// The type of the finite element space / basis of the column
-    using ColBasis = ColBasisType;
+    using ColBasis = CB;
     using ColLocalView = typename ColBasis::LocalView;
 
     using Element = typename RowLocalView::Element;
@@ -118,7 +118,7 @@ namespace AMDiS
      *
      * \tparam ContextTag  One of \ref tag::element_operator, \ref tag::intersection_operator
      *                     or \ref tag::boundary_operator indicating where to assemble this operator.
-     * \tparam PreOperator  An pre-operator that can be bound to a gridView, or a valid
+     * \tparam Expr        An pre-operator that can be bound to a gridView, or a valid
      *                      GridOperator.
      * \tparam row  A tree-path for the RowBasis
      * \tparam col  A tree-path for the ColBasis
@@ -126,9 +126,9 @@ namespace AMDiS
      * [[expects: row is valid tree-path in RowBasis]]
      * [[expects: col is valid tree-path in ColBasis]]
      **/
-    template <class ContextTag, class Operator,
+    template <class ContextTag, class Expr,
               class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
-    void addOperator(ContextTag contextTag, Operator const& preOp,
+    void addOperator(ContextTag contextTag, Expr const& expr,
                      RowTreePath row = {}, ColTreePath col = {});
 
     /// Assemble the matrix operators on the bound element.
diff --git a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
index 647d2d7ff3ef67eeccd8203112151cc8c42025cb..7088f001067996bb4ad8622267fa9f6e7eafe53c 100644
--- a/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFMatrixBase.inc.hpp
@@ -45,9 +45,9 @@ insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
 
 
 template <class RB, class CB, class B>
-  template <class ContextTag, class PreOperator, class RowTreePath, class ColTreePath>
+  template <class ContextTag, class Expr, class RowTreePath, class ColTreePath>
 void DOFMatrixBase<RB,CB,B>::
-addOperator(ContextTag contextTag, PreOperator const& preOp,
+addOperator(ContextTag contextTag, Expr const& expr,
             RowTreePath row, ColTreePath col)
 {
   static_assert( Concepts::PreTreePath<RowTreePath>,
@@ -58,9 +58,9 @@ addOperator(ContextTag contextTag, PreOperator const& preOp,
   auto i = child(rowBasis_->localView().tree(), makeTreePath(row));
   auto j = child(colBasis_->localView().tree(), makeTreePath(col));
 
-  using Context = typename ContextTag::type;
-  auto op = makeLocalOperator<Context>(preOp, rowBasis_->gridView());
-  auto localAssembler = makeUniquePtr(makeAssembler<Context>(std::move(op), i, j));
+  using LocalContext = typename ContextTag::type;
+  auto op = makeLocalOperator<LocalContext>(expr, rowBasis_->gridView());
+  auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i, j));
 
   operators_[i][j].push(contextTag, std::move(localAssembler));
 }
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index 3dc723eddd86722aceee872f11af0e495f50a3fb..1b438b0fd88c10511c21b1e7f1b3702db3e2b7dc 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -22,10 +22,10 @@ namespace AMDiS
    * assembled Operators indexed with DOF indices. The vector data is associated
    * to a global basis.
    *
-   * \tparam B  Basis of the vector
+   * \tparam GB  Basis of the vector
    * \tparam Backend  A linear algebra backend implementing the storage and operations.
    **/
-  template <class BasisType, class Backend>
+  template <class GB, class Backend>
   class DOFVectorBase
       : public DOFVectorInterface
   {
@@ -33,14 +33,14 @@ namespace AMDiS
 
   public:
     /// The type of the functionspace basis associated to this vector
-    using Basis = BasisType;
-    using LocalView = typename Basis::LocalView;
+    using GlobalBasis = GB;
+    using LocalView = typename GB::LocalView;
 
     using Element = typename LocalView::Element;
     using Geometry = typename Element::Geometry;
 
     /// The index/size - type
-    using size_type  = typename BasisType::size_type;
+    using size_type  = typename GB::size_type;
 
     /// The type of the elements of the DOFVector
     using value_type = typename Backend::value_type;
@@ -57,7 +57,7 @@ namespace AMDiS
 
   public:
     /// Constructor. Constructs new BaseVector.
-    DOFVectorBase(BasisType const& basis, DataTransferOperation op = NO_OPERATION)
+    DOFVectorBase(GlobalBasis const& basis, DataTransferOperation op = NO_OPERATION)
       : basis_(&basis)
       , dataTransfer_(DataTransferFactory::create(op, basis))
     {
@@ -119,7 +119,7 @@ namespace AMDiS
     }
 
     /// Return the basis \ref basis_ associated to the vector
-    Basis const& basis() const
+    GlobalBasis const& basis() const
     {
       return *basis_;
     }
@@ -142,7 +142,7 @@ namespace AMDiS
       return basis_->dimension();
     }
 
-    void resize(Dune::Functions::SizeInfo<Basis> const& s)
+    void resize(Dune::Functions::SizeInfo<GB> const& s)
     {
       backend_.resize(size_type(s));
     }
@@ -180,8 +180,8 @@ namespace AMDiS
     void insert(LocalView const& localView, ElementVector const& elementVector);
 
     /// Associate a local operator with this DOFVector
-    template <class ContextTag, class Operator, class TreePath = RootTreePath>
-    void addOperator(ContextTag contextTag, Operator const& preOp, TreePath path = {});
+    template <class ContextTag, class Expr, class TreePath = RootTreePath>
+    void addOperator(ContextTag contextTag, Expr const& expr, TreePath path = {});
 
     /// Assemble the vector operators on the bound element.
     void assemble(LocalView const& localView);
@@ -227,7 +227,7 @@ namespace AMDiS
 
   private:
     /// The finite element space / basis associated with the data vector
-    Basis const* basis_;
+    GlobalBasis const* basis_;
 
     /// Data backend
     Backend backend_;
@@ -236,7 +236,7 @@ namespace AMDiS
     ElementVector elementVector_;
 
     /// List of operators associated to nodes
-    VectorOperators<Basis> operators_;
+    VectorOperators<GlobalBasis> operators_;
 
     /// Data interpolation when the grid changes
     std::shared_ptr<DataTransfer> dataTransfer_;
diff --git a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
index 47264d540934cc02500dd37fb3b170396e20c836..3be8021a685366a95376096fad69044fff64029b 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
@@ -7,8 +7,8 @@
 
 namespace AMDiS {
 
-template <class BS, class B>
-void DOFVectorBase<BS,B>::
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
 init(bool asmVector)
 {
   backend_.resize(size());
@@ -20,14 +20,14 @@ init(bool asmVector)
 }
 
 
-template <class BS, class B>
-void DOFVectorBase<BS,B>::
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
 finish(bool /*asmVector*/)
 {}
 
 
-template <class BS, class B>
-void DOFVectorBase<BS,B>::
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
 insert(LocalView const& localView, ElementVector const& elementVector)
 {
   for (size_type i = 0; i < localView.size(); ++i) {
@@ -39,26 +39,26 @@ insert(LocalView const& localView, ElementVector const& elementVector)
 }
 
 
-template <class BS, class B>
-  template <class ContextTag, class PreOperator, class TreePath>
-void DOFVectorBase<BS,B>::
-addOperator(ContextTag contextTag, PreOperator const& preOp, TreePath path)
+template <class GB, class B>
+  template <class ContextTag, class Expr, class TreePath>
+void DOFVectorBase<GB,B>::
+addOperator(ContextTag contextTag, Expr const& expr, TreePath path)
 {
   static_assert( Concepts::PreTreePath<TreePath>,
       "path must be a valid treepath, or an integer/index-constant");
 
   auto i = child(basis_->localView().tree(), makeTreePath(path));
 
-  using Context = typename ContextTag::type;
-  auto op = makeLocalOperator<Context>(preOp, basis_->gridView());
-  auto localAssembler = makeUniquePtr(makeAssembler<Context>(std::move(op), i));
+  using LocalContext = typename ContextTag::type;
+  auto op = makeLocalOperator<LocalContext>(expr, basis_->gridView());
+  auto localAssembler = makeUniquePtr(makeAssembler<LocalContext>(std::move(op), i));
 
   operators_[i].push(contextTag, std::move(localAssembler));
 }
 
 
-template <class BS, class B>
-void DOFVectorBase<BS,B>::
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
 assemble(LocalView const& localView)
 {
   elementVector_ = 0;
@@ -79,8 +79,8 @@ assemble(LocalView const& localView)
 }
 
 
-template <class BS, class B>
-void DOFVectorBase<BS,B>::
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
 assemble()
 {
   auto localView = basis_->localView();
diff --git a/src/amdis/linearalgebra/LinearSolver.hpp b/src/amdis/linearalgebra/LinearSolver.hpp
index 7f338f7ab6a7211fdae486c6db92d48224ea6b05..fc62ce2f1515a19a9b2c8925ed4e63c71d0d5b60 100644
--- a/src/amdis/linearalgebra/LinearSolver.hpp
+++ b/src/amdis/linearalgebra/LinearSolver.hpp
@@ -20,12 +20,12 @@ namespace AMDiS
    * solvers where the backend provides an interface, can be assigned
    * by different Runner objects.
    **/
-  template <class Matrix, class VectorX, class VectorB, class Runner>
+  template <class Mat, class Sol, class Rhs, class Runner>
   class LinearSolver
-      : public LinearSolverInterface<Matrix, VectorX, VectorB>
+      : public LinearSolverInterface<Mat, Sol, Rhs>
   {
     using Self = LinearSolver;
-    using Super = LinearSolverInterface<Matrix, VectorX, VectorB>;
+    using Super = LinearSolverInterface<Mat, Sol, Rhs>;
     using RunnerBase = typename Super::RunnerBase;
 
   public:
@@ -52,8 +52,7 @@ namespace AMDiS
 
   private:
     /// Implements \ref LinearSolverInterface::solveSystemImpl()
-    void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
-                   SolverInfo& solverInfo) override
+    void solveImpl(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) override
     {
       Dune::Timer t;
       if (solverInfo.doCreateMatrixData()) {
diff --git a/src/amdis/linearalgebra/LinearSolverInterface.hpp b/src/amdis/linearalgebra/LinearSolverInterface.hpp
index 285d7dfa8f6f87273d8bcef3cfa0fad994e77874..4e0474bf6bae64d164b95f0b883966083ca09d84 100644
--- a/src/amdis/linearalgebra/LinearSolverInterface.hpp
+++ b/src/amdis/linearalgebra/LinearSolverInterface.hpp
@@ -20,11 +20,11 @@
 namespace AMDiS
 {
   /// Abstract base class for linear solvers
-  template <class Matrix, class VectorX, class VectorB = VectorX>
+  template <class Mat, class Sol, class Rhs = Sol>
   class LinearSolverInterface
   {
   protected:
-    using RunnerBase = RunnerInterface<Matrix, VectorX, VectorB>;
+    using RunnerBase = RunnerInterface<Mat, Sol, Rhs>;
 
   public:
     /// Destructor.
@@ -40,8 +40,7 @@ namespace AMDiS
      *  \p x     A [block-]vector for the unknown components.
      *  \p b     A [block-]vector for the right-hand side of the linear system.
      **/
-    void solve(Matrix const& A, VectorX& x, VectorB const& b,
-               SolverInfo& solverInfo)
+    void solve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo)
     {
       solveImpl(A, x, b, solverInfo);
     }
@@ -51,8 +50,7 @@ namespace AMDiS
 
   private:
     /// main methods that all solvers must implement
-    virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
-                           SolverInfo& solverInfo) = 0;
+    virtual void solveImpl(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) = 0;
   };
 
 
diff --git a/src/amdis/linearalgebra/PreconditionerInterface.hpp b/src/amdis/linearalgebra/PreconditionerInterface.hpp
index 409e8270584c6a1c005197997f2d01c1f64d2301..512b924c71afa04aee437e7327b8836581621579 100644
--- a/src/amdis/linearalgebra/PreconditionerInterface.hpp
+++ b/src/amdis/linearalgebra/PreconditionerInterface.hpp
@@ -5,23 +5,23 @@
 namespace AMDiS
 {
   /// Interface for Preconditioner types
-  template <class Matrix, class VectorX, class VectorB = VectorX>
+  template <class Mat, class Sol, class Rhs = Sol>
   struct PreconditionerInterface
   {
     /// Virtual destructor.
     virtual ~PreconditionerInterface() = default;
 
     /// Is called a the beginning of a solution procedure
-    virtual void init(Matrix const& fullMatrix) = 0;
+    virtual void init(Mat const& fullMatrix) = 0;
 
     /// Is called at the end of a solution procedure
     virtual void exit() = 0;
 
     /// Apply the preconditioner to a vector \p b and store the result in \p x
-    virtual void solve(VectorB const& b, VectorX& x) const = 0;
+    virtual void solve(Rhs const& b, Sol& x) const = 0;
 
     /// Apply the transposed preconditioner to a vector \p b and store the result in \p x
-    virtual void adjoint_solve(VectorB const& b, VectorX& x) const
+    virtual void adjoint_solve(Rhs const& b, Sol& x) const
     {
       error_exit("Must be implemented by derived class.");
     }
diff --git a/src/amdis/linearalgebra/RunnerInterface.hpp b/src/amdis/linearalgebra/RunnerInterface.hpp
index a9e733d4da61d0c1d44ef603675e5ad539fe5253..4aebd88460e50b68cce440a8cd5479971bd906ba 100644
--- a/src/amdis/linearalgebra/RunnerInterface.hpp
+++ b/src/amdis/linearalgebra/RunnerInterface.hpp
@@ -8,29 +8,27 @@ namespace AMDiS
   class SolverInfo;
 
   /// Interface for Runner / Worker types used in solver classes
-  template <class Matrix, class VectorX, class VectorB = VectorX>
+  template <class Mat, class Sol, class Rhs = Sol>
   class RunnerInterface
   {
   protected:
-    using PreconBase = PreconditionerInterface<Matrix, VectorX, VectorB>;
+    using PreconBase = PreconditionerInterface<Mat, Sol, Rhs>;
 
   public:
     /// virtual destructor
     virtual ~RunnerInterface() = default;
 
     /// Is called at the beginning of a solution procedure
-    virtual void init(Matrix const& A) = 0;
+    virtual void init(Mat const& A) = 0;
 
     /// Is called at the end of a solution procedure
     virtual void exit() = 0;
 
     /// Solve the system A*x = b
-    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
-                      SolverInfo& solverInfo) = 0;
+    virtual int solve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) = 0;
 
     /// Solve the adjoint system A^T*x = b
-    virtual int adjointSolve(Matrix const& A, VectorX& x, VectorB const& b,
-                             SolverInfo& solverInfo)
+    virtual int adjointSolve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo)
     {
       error_exit("Must be implemented by derived class.");
       return 0;
diff --git a/src/amdis/localoperators/ConvectionDiffusionOperator.hpp b/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
index 3022a5b16dcedbc613e0485f876bb120f940ddaa..a03c70a3857d71c55af65b9f232ebef47d57cb3c 100644
--- a/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
+++ b/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
@@ -20,12 +20,11 @@ namespace AMDiS
   /// convection-diffusion operator, see \ref convectionDiffusion
   /// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or
   /// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving)
-  template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
+  template <class LC, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
   class ConvectionDiffusionOperator
-      : public LocalOperator<ConvectionDiffusionOperator<LocalContext, GridFctA, GridFctB, GridFctC, GridFctF, conserving>,
-                             LocalContext>
+      : public LocalOperator<ConvectionDiffusionOperator<LC, GridFctA, GridFctB, GridFctC, GridFctF, conserving>, LC>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
 
     using A_range = typename GridFctA::Range;
     static_assert( Size_v<A_range> == 1 || (Rows_v<A_range> == dow && Cols_v<A_range> == dow),
@@ -49,22 +48,20 @@ namespace AMDiS
       , gridFctF_(gridFctF)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode, ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only." );
 
-      static_assert(std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{},
+      static_assert(std::is_same<FiniteElementType_t<RN>, FiniteElementType_t<CN>>{},
         "Galerkin operator requires equal finite elements for test and trial space." );
 
-      using RangeFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
 
-      auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element());
-      auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element());
-      auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element());
+      auto localFctA = localFunction(gridFctA_); localFctA.bind(contextGeo.element());
+      auto localFctB = localFunction(gridFctB_); localFctB.bind(contextGeo.element());
+      auto localFctC = localFunction(gridFctC_); localFctC.bind(contextGeo.element());
 
       auto const& localFE = rowNode.finiteElement();
       std::size_t numLocalFe = localFE.size();
@@ -73,22 +70,22 @@ namespace AMDiS
                               this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
                               this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)});
 
-      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>;
-      auto const& quad = QuadratureRules::rule(context.type(), quadDeg);
+      using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::LocalContext::mydimension>;
+      auto const& quad = QuadratureRules::rule(contextGeo.type(), quadDeg);
 
-      NodeQuadCache<RowNode> cache(localFE.localBasis());
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context, quad);
+      NodeQuadCache<RN> cache(localFE.localBasis());
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo, quad);
 
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // the values of the shape functions on the reference element at the quadrature point
         auto const& shapeValues = shapeValuesCache[iq];
@@ -97,7 +94,7 @@ namespace AMDiS
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
 
@@ -139,36 +136,34 @@ namespace AMDiS
     }
 
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isLeaf,
         "Operator can be applied to Leaf-Nodes only." );
 
-      auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
+      auto localFctF = localFunction(gridFctF_); localFctF.bind(contextGeo.element());
 
       auto const& localFE = node.finiteElement();
       std::size_t numLocalFe = localFE.size();
 
       int quad_order = this->getDegree(0,coeffOrder(localFctF),node);
 
-      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>;
-      auto const& quad = QuadratureRules::rule(context.type(), quad_order);
+      using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::LocalContext::dimension>;
+      auto const& quad = QuadratureRules::rule(contextGeo.type(), quad_order);
 
       NodeQuadCache<Node> cache(localFE.localBasis());
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context, quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo, quad);
 
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // the values of the shape functions on the reference element at the quadrature point
         auto const& shapeValues = shapeValuesCache[iq];
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = localFctF(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = localFctF(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         for (std::size_t i = 0; i < numLocalFe; ++i) {
           const auto local_i = node.localIndex(i);
diff --git a/src/amdis/localoperators/FirstOrderDivTestvec.hpp b/src/amdis/localoperators/FirstOrderDivTestvec.hpp
index d48f3c4f0563fbd1787da6fb180732d792fea713..fd749e470b776ec4682e0e5f1fd1811cba425bd2 100644
--- a/src/amdis/localoperators/FirstOrderDivTestvec.hpp
+++ b/src/amdis/localoperators/FirstOrderDivTestvec.hpp
@@ -21,13 +21,12 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\nabla\cdot\Psi, c\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::divtestvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -36,15 +35,13 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& context, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isPower, "Node must be a Power-Node.");
 
       static const std::size_t CHILDREN = Node::CHILDREN;
-      static_assert( CHILDREN == Context::dow, "" );
+      static_assert( CHILDREN == CG::dow, "" );
 
       auto const& quad = this->getQuadratureRule(context.type(), node);
       auto const& localFE = node.child(0).finiteElement();
@@ -68,7 +65,7 @@ namespace AMDiS
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/FirstOrderDivTestvecTrial.hpp b/src/amdis/localoperators/FirstOrderDivTestvecTrial.hpp
index 655be6c0f6388f48a128be94802fb00e6bf55229..45f85c223069422ef95a35a9039c6502e21b92aa 100644
--- a/src/amdis/localoperators/FirstOrderDivTestvecTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderDivTestvecTrial.hpp
@@ -18,13 +18,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\nabla\cdot\Psi, c\,\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::divtestvec_trial, LocalContext, GridFct>
-      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::divtestvec_trial, LocalContext, GridFct>,
-                                              GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::divtestvec_trial, LC, GridFct>
+      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::divtestvec_trial, LC, GridFct>,
+                                              GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>>
   {
     using Self = GridFunctionOperator;
-    using Transposed = GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>;
+    using Transposed = GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>;
     using Super = GridFunctionOperatorTransposed<Self, Transposed>;
 
   public:
diff --git a/src/amdis/localoperators/FirstOrderGradTest.hpp b/src/amdis/localoperators/FirstOrderGradTest.hpp
index c7060f5109b76040ee7dbb04c542357538700861..150cbc0dc1a2f743ada361437ef93ec4595c806e 100644
--- a/src/amdis/localoperators/FirstOrderGradTest.hpp
+++ b/src/amdis/localoperators/FirstOrderGradTest.hpp
@@ -21,14 +21,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\nabla\psi, b\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::gradtest, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::gradtest, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
 
@@ -37,37 +36,35 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isLeaf, "Node must be Leaf-Node.");
 
-      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
       auto const& localFE = node.finiteElement();
       std::size_t feSize = localFE.size();
 
       using RangeFieldType = typename NodeQuadCache<Node>::Traits::RangeFieldType;
       NodeQuadCache<Node> cache(localFE.localBasis());
 
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
         const auto factor = Super::coefficient(local);
-        const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto dx = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // The gradients of the shape functions on the reference element
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/FirstOrderGradTestTrial.hpp b/src/amdis/localoperators/FirstOrderGradTestTrial.hpp
index c742f28fe60e715169fc7096b361929aa01852ce..41ec01635510397804a7ce9e7a0a822a2f41ec0f 100644
--- a/src/amdis/localoperators/FirstOrderGradTestTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderGradTestTrial.hpp
@@ -19,13 +19,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\nabla\psi, \mathbf{b}\,\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::gradtest_trial, LocalContext, GridFct>
-      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trial, LocalContext, GridFct>,
-                                              GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::gradtest_trial, LC, GridFct>
+      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trial, LC, GridFct>,
+                                              GridFunctionOperator<tag::test_gradtrial, LC, GridFct>>
   {
     using Self = GridFunctionOperator;
-    using Transposed = GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>;
+    using Transposed = GridFunctionOperator<tag::test_gradtrial, LC, GridFct>;
     using Super = GridFunctionOperatorTransposed<Self, Transposed>;
 
   public:
diff --git a/src/amdis/localoperators/FirstOrderGradTestTrialvec.hpp b/src/amdis/localoperators/FirstOrderGradTestTrialvec.hpp
index 4649db6fb51a5fe37a0b73b92b5157151d682d78..ab26534faa1e579b10727a5449b580ffb549e5cc 100644
--- a/src/amdis/localoperators/FirstOrderGradTestTrialvec.hpp
+++ b/src/amdis/localoperators/FirstOrderGradTestTrialvec.hpp
@@ -18,13 +18,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\nabla\psi, c\,\Phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::gradtest_trialvec, LocalContext, GridFct>
-      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trialvec, LocalContext, GridFct>,
-                                              GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::gradtest_trialvec, LC, GridFct>
+      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trialvec, LC, GridFct>,
+                                              GridFunctionOperator<tag::testvec_gradtrial, LC, GridFct>>
   {
     using Self = GridFunctionOperator;
-    using Transposed = GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>;
+    using Transposed = GridFunctionOperator<tag::testvec_gradtrial, LC, GridFct>;
     using Super = GridFunctionOperatorTransposed<Self, Transposed>;
 
   public:
diff --git a/src/amdis/localoperators/FirstOrderPartialTest.hpp b/src/amdis/localoperators/FirstOrderPartialTest.hpp
index ca72328b17c73b4baf8c80d5846d2b4b0d08decd..46fcb9f34d61965ed601cbca84d740e25d765370 100644
--- a/src/amdis/localoperators/FirstOrderPartialTest.hpp
+++ b/src/amdis/localoperators/FirstOrderPartialTest.hpp
@@ -22,13 +22,12 @@ namespace AMDiS
 
   /// first-order operator \f$ \langle\partial_i\psi, c\rangle \f$
 
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::partialtest, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::partialtest, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -38,31 +37,29 @@ namespace AMDiS
       , comp_(tag.comp)
     {}
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");
 
-      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
       auto const& localFE = node.finiteElement();
       std::size_t feSize = localFE.size();
 
       using RangeFieldType = typename NodeQuadCache<Node>::Traits::RangeFieldType;
       NodeQuadCache<Node> cache(localFE.localBasis());
 
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
-        assert(jacobian.M() == Context::dim);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
+        assert(jacobian.M() == CG::dim);
 
         // The multiplicative factor in the integral transformation formula
-        const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto dx = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto c = Super::coefficient(local);
 
         // The gradients of the shape functions on the reference element
diff --git a/src/amdis/localoperators/FirstOrderPartialTestTrial.hpp b/src/amdis/localoperators/FirstOrderPartialTestTrial.hpp
index 246994760204e741f172d9083a2f67e65514e4fc..ca62e1947f9542bfc2fb04cc2de2f4483d7ffa0b 100644
--- a/src/amdis/localoperators/FirstOrderPartialTestTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderPartialTestTrial.hpp
@@ -26,13 +26,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\partial_i\psi, c\,\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::partialtest_trial, LocalContext, GridFct>
-      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::partialtest_trial, LocalContext, GridFct>,
-                                              GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::partialtest_trial, LC, GridFct>
+      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::partialtest_trial, LC, GridFct>,
+                                              GridFunctionOperator<tag::test_partialtrial, LC, GridFct>>
   {
     using Self = GridFunctionOperator;
-    using Transposed = GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>;
+    using Transposed = GridFunctionOperator<tag::test_partialtrial, LC, GridFct>;
     using Super = GridFunctionOperatorTransposed<Self, Transposed>;
 
   public:
diff --git a/src/amdis/localoperators/FirstOrderTestDivTrialvec.hpp b/src/amdis/localoperators/FirstOrderTestDivTrialvec.hpp
index 07e31f1d40948abab4fb6d65154b35c856d19042..94ca60f94f66b78626f7e4f3e869886e76cab6b8 100644
--- a/src/amdis/localoperators/FirstOrderTestDivTrialvec.hpp
+++ b/src/amdis/localoperators/FirstOrderTestDivTrialvec.hpp
@@ -21,13 +21,12 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\psi, c\,\nabla\cdot\Phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -36,40 +35,37 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isPower,
-        "RowNode must be Leaf-Node and ColNode must be a Power-Node.");
+      static_assert(RN::isLeaf && CN::isPower,
+        "row-node must be Leaf-Node and col-node must be a Power-Node.");
 
-      static const std::size_t CHILDREN = ColNode::CHILDREN;
-      static_assert( CHILDREN == Context::dow, "" );
+      static const std::size_t CHILDREN = CN::CHILDREN;
+      static_assert( CHILDREN == CG::dow, "" );
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<typename ColNode::ChildType>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<typename CN::ChildType>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<typename ColNode::ChildType> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<typename CN::ChildType> colCache(colLocalFE.localBasis());
 
-      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(context, quad);
-      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo, quad);
+      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // the values of the shape functions on the reference element at the quadrature point
         auto const& shapeValues = shapeValuesCache[iq];
@@ -78,7 +74,7 @@ namespace AMDiS
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> colGradients;
         colGradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/FirstOrderTestGradTrial.hpp b/src/amdis/localoperators/FirstOrderTestGradTrial.hpp
index fa161142757843a3cba123936d6f3a47e71becfc..f59f7e39b14a0a5fa5f8692f6d22e5b8d670c75f 100644
--- a/src/amdis/localoperators/FirstOrderTestGradTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderTestGradTrial.hpp
@@ -22,14 +22,13 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\psi, \mathbf{b}\cdot\nabla\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test_gradtrial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_gradtrial, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
 
@@ -38,37 +37,34 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.finiteElement();
       auto const& colLocalFE = colNode.finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<ColNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<CN> colCache(colLocalFE.localBasis());
 
-      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(context, quad);
-      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo, quad);
+      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto b = Super::coefficient(local);
 
         // the values of the shape functions on the reference element at the quadrature point
@@ -78,7 +74,7 @@ namespace AMDiS
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> colGradients;
         colGradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp b/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
index 33628e2fe054a361d5bc594aec53eecbb684e521..97f647a111a29ec10de22cf21df54c96dbd4c8c1 100644
--- a/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
@@ -29,13 +29,12 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\psi, c\,\partial_i\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test_partialtrial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_partialtrial, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -45,38 +44,35 @@ namespace AMDiS
       , comp_(tag.comp)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.finiteElement();
       auto const& colLocalFE = colNode.finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<ColNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<CN> colCache(colLocalFE.localBasis());
 
-      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(context, quad);
-      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo, quad);
+      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
-        assert(jacobian.M() == Context::dim);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
+        assert(jacobian.M() == CG::dim);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto c = Super::coefficient(local);
 
         // the values of the shape functions on the reference element at the quadrature point
diff --git a/src/amdis/localoperators/FirstOrderTestvecGradTrial.hpp b/src/amdis/localoperators/FirstOrderTestvecGradTrial.hpp
index 35949d1aca82f5458ad4bbd6f83aef4ce27ea277..5386586062c907921d96953285287995f888f9db 100644
--- a/src/amdis/localoperators/FirstOrderTestvecGradTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderTestvecGradTrial.hpp
@@ -21,13 +21,12 @@ namespace AMDiS
 
 
   /// first-order operator \f$ \langle\Psi, c\,\nabla\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::testvec_gradtrial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_gradtrial, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -36,40 +35,37 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isPower && ColNode::isLeaf,
-        "ColNode must be Leaf-Node and RowNode must be a Power-Node.");
+      static_assert(RN::isPower && CN::isLeaf,
+        "col-node must be Leaf-Node and row-node must be a Power-Node.");
 
-      static const std::size_t CHILDREN = RowNode::CHILDREN;
-      static_assert( CHILDREN == Context::dow, "" );
+      static const std::size_t CHILDREN = RN::CHILDREN;
+      static_assert( CHILDREN == CG::dow, "" );
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<ColNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
 
-      NodeQuadCache<typename RowNode::ChildType> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
+      NodeQuadCache<typename RN::ChildType> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<CN> colCache(colLocalFE.localBasis());
 
-      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(context, quad);
-      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo, quad);
+      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // the values of the shape functions on the reference element at the quadrature point
         auto const& shapeValues = shapeValuesCache[iq];
@@ -78,7 +74,7 @@ namespace AMDiS
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> colGradients;
         colGradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/SecondOrderDivTestvecDivTrialvec.hpp b/src/amdis/localoperators/SecondOrderDivTestvecDivTrialvec.hpp
index f0cfa6cf90197f83bd79d417a0622a2d1ede4a8d..25692c7f6b4d00254e1c18ebed6e64de295ab8de 100644
--- a/src/amdis/localoperators/SecondOrderDivTestvecDivTrialvec.hpp
+++ b/src/amdis/localoperators/SecondOrderDivTestvecDivTrialvec.hpp
@@ -21,13 +21,12 @@ namespace AMDiS
 
 
   /// second-order operator \f$ \langle\nabla\cdot\Psi, c\,\nabla\cdot\Phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::divtestvec_divtrialvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec_divtrialvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::divtestvec_divtrialvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec_divtrialvec, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -36,72 +35,67 @@ namespace AMDiS
       : Super(expr, 2)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isPower && ColNode::isPower,
+      static_assert(RN::isPower && CN::isPower,
         "Operator can be applied to Power-Nodes only.");
 
-      const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+      const bool sameFE = std::is_same<FiniteElementType_t<RN>, FiniteElementType_t<CN>>::value;
       const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       if (sameFE && sameNode)
-        getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix);
+        getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix);
       else
-        getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+        getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
     }
 
   protected:
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixStandard(Context const& context,
-                                  QuadratureRule const& quad,
-                                  RowNode const& rowNode,
-                                  ColNode const& colNode,
-                                  ElementMatrix& elementMatrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
+                                  RN const& rowNode, CN const& colNode,
+                                  Mat& elementMatrix)
     {
-      static const std::size_t CHILDREN = RowNode::CHILDREN;
-      static_assert( RowNode::CHILDREN == ColNode::CHILDREN, "" );
+      static const std::size_t CHILDREN = RN::CHILDREN;
+      static_assert( RN::CHILDREN == CN::CHILDREN, "" );
 
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RowFieldType = typename NodeQuadCache<typename RowNode::ChildType>::Traits::RangeFieldType;
-      using ColFieldType = typename NodeQuadCache<typename ColNode::ChildType>::Traits::RangeFieldType;
+      using RowFieldType = typename NodeQuadCache<typename RN::ChildType>::Traits::RangeFieldType;
+      using ColFieldType = typename NodeQuadCache<typename CN::ChildType>::Traits::RangeFieldType;
 
-      NodeQuadCache<typename RowNode::ChildType> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<typename ColNode::ChildType> colCache(colLocalFE.localBasis());
+      NodeQuadCache<typename RN::ChildType> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<typename CN::ChildType> colCache(colLocalFE.localBasis());
 
-      auto const& rowGradientsCache = rowCache.evaluateJacobianAtQP(context, quad);
-      auto const& colGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& rowGradientsCache = rowCache.evaluateJacobianAtQP(contextGeo, quad);
+      auto const& colGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // The gradients of the shape functions on the reference element
         auto const& rowShapeGradients = rowGradientsCache[iq];
         auto const& colShapeGradients = colGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using RowWorldVector = FieldVector<RowFieldType,Context::dow>;
+        using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
         thread_local std::vector<RowWorldVector> rowGradients;
         rowGradients.resize(rowShapeGradients.size());
 
         for (std::size_t i = 0; i < rowGradients.size(); ++i)
           jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
 
-        using ColWorldVector = FieldVector<ColFieldType,Context::dow>;
+        using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
         thread_local std::vector<ColWorldVector> colGradients;
         colGradients.resize(colShapeGradients.size());
 
@@ -121,12 +115,10 @@ namespace AMDiS
     }
 
 
-    template <class Context, class QuadratureRule, class Node, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   Node const& node,
-                                   Node const& /*colNode*/,
-                                   ElementMatrix& elementMatrix)
+    template <class CG, class QR, class Node, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   Node const& node, Node const& /*colNode*/,
+                                   Mat& elementMatrix)
     {
       static const std::size_t CHILDREN = Node::CHILDREN;
 
@@ -138,22 +130,22 @@ namespace AMDiS
 
       NodeQuadCache<LeafNode> cache(localFE.localBasis());
 
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         // The gradients of the shape functions on the reference element
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
 
diff --git a/src/amdis/localoperators/SecondOrderGradTestGradTrial.hpp b/src/amdis/localoperators/SecondOrderGradTestGradTrial.hpp
index 3e9416fb19bbffcb9ec9191e0df6ce893df2a0c0..29b71ec06b55857adce7eeb5f77be941a29696c5 100644
--- a/src/amdis/localoperators/SecondOrderGradTestGradTrial.hpp
+++ b/src/amdis/localoperators/SecondOrderGradTestGradTrial.hpp
@@ -22,14 +22,13 @@ namespace AMDiS
 
 
   /// second-order operator \f$ \langle\nabla\psi, c\,\nabla\phi\rangle \f$, or \f$ \langle\nabla\psi, A\,\nabla\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::gradtest_gradtrial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest_gradtrial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::gradtest_gradtrial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest_gradtrial, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     using expr_value_type = typename GridFct::Range;
     static_assert( Size_v<expr_value_type> == 1 || (Rows_v<expr_value_type> == dow && Cols_v<expr_value_type> == dow),
@@ -40,61 +39,56 @@ namespace AMDiS
       : Super(expr, 2)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
-      const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+      const bool sameFE = std::is_same<FiniteElementType_t<RN>, FiniteElementType_t<CN>>::value;
       const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
       using Category = ValueCategory_t<typename GridFct::Range>;
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       if (sameFE && sameNode)
-        getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix, Category{});
+        getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
       else if (sameFE)
-        getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+        getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
       else
         error_exit("Not implemented: currently only the implementation for equal fespaces available");
     }
 
   protected:
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixStandard(Context const& context,
-                                  QuadratureRule const& quad,
-                                  RowNode const& rowNode,
-                                  ColNode const& colNode,
-                                  ElementMatrix& elementMatrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
+                                  RN const& rowNode, CN const& colNode,
+                                  Mat& elementMatrix)
     {
       auto const& localFE = rowNode.finiteElement();
       std::size_t size = localFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> cache(localFE.localBasis());
+      NodeQuadCache<RN> cache(localFE.localBasis());
 
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto exprValue = Super::coefficient(local);
 
         // The gradients of the shape functions on the reference element
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
 
@@ -111,37 +105,34 @@ namespace AMDiS
       }
     }
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   RowNode const& node,
-                                   ColNode const& /*colNode*/,
-                                   ElementMatrix& elementMatrix,
-                                   tag::scalar)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   RN const& node, CN const& /*colNode*/,
+                                   Mat& elementMatrix, tag::scalar)
     {
       auto const& localFE = node.finiteElement();
       std::size_t size = localFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> cache(localFE.localBasis());
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      NodeQuadCache<RN> cache(localFE.localBasis());
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position())
                                                       * quad[iq].weight();
 
         // The gradients of the shape functions on the reference element
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
         for (std::size_t i = 0; i < gradients.size(); ++i)
@@ -163,37 +154,34 @@ namespace AMDiS
       }
     }
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   RowNode const& node,
-                                   ColNode const& /*colNode*/,
-                                   ElementMatrix& elementMatrix,
-                                   tag::matrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   RN const& node, CN const& /*colNode*/,
+                                   Mat& elementMatrix, tag::matrix)
     {
       auto const& localFE = node.finiteElement();
       std::size_t size = localFE.size();
 
-      using RangeFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
+      using RangeFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> cache(localFE.localBasis());
-      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
+      NodeQuadCache<RN> cache(localFE.localBasis());
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto exprValue = Super::coefficient(local);
 
         // The gradients of the shape functions on the reference element
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
-        using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
         for (std::size_t i = 0; i < gradients.size(); ++i)
diff --git a/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp b/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
index b0c78dd487676c840fbc0f535e1e4564e069c088..ae9d2b42dc170793dae6dbf11483b40810ac8e4c 100644
--- a/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
+++ b/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
@@ -25,13 +25,12 @@ namespace AMDiS
 
 
   /// second-order operator \f$ \langle\partial_i\psi, c\,\partial_j\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::partialtest_partialtrial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest_partialtrial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::partialtest_partialtrial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest_partialtrial, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -42,38 +41,35 @@ namespace AMDiS
       , compTrial_(tag.comp_trial)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.finiteElement();
       auto const& colLocalFE = colNode.finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      using RowFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
-      using ColFieldType = typename NodeQuadCache<ColNode>::Traits::RangeFieldType;
+      using RowFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
+      using ColFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<CN> colCache(colLocalFE.localBasis());
 
-      auto const& rowGradientsCache = rowCache.evaluateJacobianAtQP(context, quad);
-      auto const& colGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
+      auto const& rowGradientsCache = rowCache.evaluateJacobianAtQP(contextGeo, quad);
+      auto const& colGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto c = Super::coefficient(local);
 
         // The gradients of the shape functions on the reference element
diff --git a/src/amdis/localoperators/StokesOperator.hpp b/src/amdis/localoperators/StokesOperator.hpp
index d4a8e94a72e0f17a114e2475254d616a6583f173..a8f03f03709b550a6b0533f6b9189904f661cc79 100644
--- a/src/amdis/localoperators/StokesOperator.hpp
+++ b/src/amdis/localoperators/StokesOperator.hpp
@@ -21,13 +21,12 @@ namespace AMDiS
 
 
   /// Stokes operator \f$ \langle\nabla\mathbf{v}, \nu \nabla\mathbf{u}\rangle + \langle\nabla\cdot\mathbf{v}, p\rangle + \langle q, \langle\nabla\cdot\mathbf{u}\rangle \f$
-  template <class LocalContext, class ViscosityExpr>
-  class GridFunctionOperator<tag::stokes, LocalContext, ViscosityExpr>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::stokes, LocalContext, ViscosityExpr>,
-                                        LocalContext, ViscosityExpr>
+  template <class LC, class ViscosityExpr>
+  class GridFunctionOperator<tag::stokes, LC, ViscosityExpr>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::stokes, LC, ViscosityExpr>, LC, ViscosityExpr>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, ViscosityExpr>;
+    using Super = GridFunctionOperatorBase<Self, LC, ViscosityExpr>;
 
     static_assert( Size_v<typename ViscosityExpr::Range> == 1, "Viscosity must be of scalar type." );
 
@@ -36,11 +35,8 @@ namespace AMDiS
       : Super(expr, 1)
     {}
 
-    template <class Context, class Node, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          Node const& tree,
-                          Node const& /*colNode*/,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class Node, class Mat>
+    void getElementMatrix(CG const& contextGeo, Node const& tree, Node const& /*colNode*/, Mat& elementMatrix)
     {
       using VelocityNode = typename Node::template Child<0>::type;
       using VelocityCompNode = typename VelocityNode::template Child<0>::type;
@@ -49,7 +45,7 @@ namespace AMDiS
 
       using namespace Dune::Indices;
 
-      auto const& quad = this->getQuadratureRule(context.type(), tree, tree);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), tree, tree);
 
       auto const& velocityLocalFE = tree.child(_0,0).finiteElement();
       std::size_t numVelocityLocalFE = velocityLocalFE.size();
@@ -60,24 +56,24 @@ namespace AMDiS
       NodeQuadCache<VelocityCompNode> velocityCache(velocityLocalFE.localBasis());
       NodeQuadCache<PressureNode> pressureCache(pressureLocalFE.localBasis());
 
-      auto const& shapeGradientsCache = velocityCache.evaluateJacobianAtQP(context, quad);
-      auto const& pressureValuesCache = pressureCache.evaluateFunctionAtQP(context, quad);
+      auto const& shapeGradientsCache = velocityCache.evaluateJacobianAtQP(contextGeo, quad);
+      auto const& pressureValuesCache = pressureCache.evaluateFunctionAtQP(contextGeo, quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
-        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto vel_factor = Super::coefficient(local) * factor;
 
         auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
         using RangeFieldType = typename NodeQuadCache<VelocityCompNode>::Traits::RangeFieldType;
-        using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
+        using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
         thread_local std::vector<WorldVector> gradients;
         gradients.resize(shapeGradients.size());
         for (std::size_t i = 0; i < gradients.size(); ++i)
@@ -88,14 +84,14 @@ namespace AMDiS
         // <viscosity * grad(u), grad(v)>
         for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
           const auto value_ii = vel_factor * (gradients[i] * gradients[i]);
-          for (std::size_t k = 0; k < Context::dow; ++k) {
+          for (std::size_t k = 0; k < CG::dow; ++k) {
             const auto local_ki = tree.child(_0,k).localIndex(i);
             elementMatrix[local_ki][local_ki] += value_ii;
           }
 
           for (std::size_t j = i+1; j < numVelocityLocalFE; ++j) {
             const auto value = vel_factor * (gradients[i] * gradients[j]);
-            for (std::size_t k = 0; k < Context::dow; ++k) {
+            for (std::size_t k = 0; k < CG::dow; ++k) {
               const auto local_ki = tree.child(_0,k).localIndex(i);
               const auto local_kj = tree.child(_0,k).localIndex(j);
               elementMatrix[local_ki][local_kj] += value;
@@ -108,7 +104,7 @@ namespace AMDiS
         for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
           for (std::size_t j = 0; j < numPressureLocalFE; ++j) {
             const auto value = factor * pressureValues[j];
-            for (std::size_t k = 0; k < Context::dow; ++k) {
+            for (std::size_t k = 0; k < CG::dow; ++k) {
               const auto local_v = tree.child(_0,k).localIndex(i);
               const auto local_p = tree.child(_1).localIndex(j);
 
diff --git a/src/amdis/localoperators/ZeroOrderTest.hpp b/src/amdis/localoperators/ZeroOrderTest.hpp
index 9aad9f277def095d13420f9f186c864c5d24e988..a6bf83bf301d96fc2d4569967ce96a0871d13112 100644
--- a/src/amdis/localoperators/ZeroOrderTest.hpp
+++ b/src/amdis/localoperators/ZeroOrderTest.hpp
@@ -20,13 +20,12 @@ namespace AMDiS
 
 
   /// zero-order vector-operator \f$ (c\, \psi) \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -35,26 +34,24 @@ namespace AMDiS
       : Super(expr, 0)
     {}
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only");
 
-      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
       auto const& localFE = node.finiteElement();
       std::size_t size = localFE.size();
 
       NodeQuadCache<Node> cache(localFE.localBasis());
 
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context,quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position())
                                                       * quad[iq].weight();
 
         auto const& shapeValues = shapeValuesCache[iq];
diff --git a/src/amdis/localoperators/ZeroOrderTestTrial.hpp b/src/amdis/localoperators/ZeroOrderTestTrial.hpp
index 127ef16d564f0b7c239e1e042ce5173ea1de49d4..a9f4687b9fc70a2fea9e04fbe7612e766f66e965 100644
--- a/src/amdis/localoperators/ZeroOrderTestTrial.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestTrial.hpp
@@ -20,13 +20,12 @@ namespace AMDiS
 
 
   /// zero-order operator \f$ \langle\psi, c\,\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test_trial, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trial, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test_trial, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trial, LC, GridFct>, LC, GridFct>
   {
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
 
@@ -35,42 +34,37 @@ namespace AMDiS
       : Super(expr, 0)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+      const bool sameFE = std::is_same<FiniteElementType_t<RN>, FiniteElementType_t<CN>>::value;
       const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       if (sameFE && sameNode)
-        getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix);
+        getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix);
       else
-        getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+        getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
     }
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only");
 
-      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
       auto const& localFE = node.finiteElement();
       std::size_t size = localFE.size();
 
       NodeQuadCache<Node> cache(localFE.localBasis());
 
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context,quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position())
                                                       * quad[iq].weight();
 
         auto const& shapeValues = shapeValuesCache[iq];
@@ -83,15 +77,12 @@ namespace AMDiS
 
 
   protected:
-
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixStandard(Context const& context,
-                                  QuadratureRule const& quad,
-                                  RowNode const& rowNode,
-                                  ColNode const& colNode,
-                                  ElementMatrix& elementMatrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
+                                  RN const& rowNode, CN const& colNode,
+                                  Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
       auto const& rowLocalFE = rowNode.finiteElement();
@@ -99,17 +90,17 @@ namespace AMDiS
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<CN> colCache(colLocalFE.localBasis());
 
-      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(context,quad);
-      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(context,quad);
+      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo,quad);
+      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         auto const& rowShapeValues = rowShapeValuesCache[iq];
         auto const& colShapeValues = colShapeValuesCache[iq];
@@ -127,29 +118,26 @@ namespace AMDiS
     }
 
 
-    template <class Context, class QuadratureRule,
-              class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   RowNode const& node,
-                                   ColNode const& /*colNode*/,
-                                   ElementMatrix& elementMatrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   RN const& node, CN const& /*colNode*/,
+                                   Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isLeaf,
+      static_assert(RN::isLeaf && CN::isLeaf,
         "Operator can be applied to Leaf-Nodes only.");
 
       auto const& localFE = node.finiteElement();
       std::size_t size = localFE.size();
 
-      NodeQuadCache<RowNode> cache(localFE.localBasis());
+      NodeQuadCache<RN> cache(localFE.localBasis());
 
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context,quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         auto const& shapeValues = shapeValuesCache[iq];
 
diff --git a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
index 670b6d95a4c0e17bf34bad3b3279e9137e1aed39..ef7360dd115653cebf3d8a80bebb4f00f26e1b5f 100644
--- a/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestTrialvec.hpp
@@ -20,14 +20,13 @@ namespace AMDiS
 
 
   /// zero-order operator \f$ \langle\psi, \mathbf{b}\cdot\Phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::test_trialvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trialvec, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
 
@@ -36,35 +35,32 @@ namespace AMDiS
       : Super(expr, 0)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isLeaf && ColNode::isPower,
-        "RowNode must be Leaf-Node and ColNode must be a Power-Node.");
+      static_assert(RN::isLeaf && CN::isPower,
+        "RN must be Leaf-Node and CN must be a Power-Node.");
 
-      static const std::size_t CHILDREN = ColNode::CHILDREN;
+      static const std::size_t CHILDREN = CN::CHILDREN;
       static_assert( std::is_same<typename GridFct::Range, Dune::FieldVector<double, int(CHILDREN)>>::value, "" );
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       auto const& rowLocalFE = rowNode.finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
       std::size_t rowSize = rowLocalFE.size();
       std::size_t colSize = colLocalFE.size();
 
-      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<typename ColNode::ChildType> colCache(colLocalFE.localBasis());
+      NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<typename CN::ChildType> colCache(colLocalFE.localBasis());
 
-      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(context,quad);
-      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(context,quad);
+      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo,quad);
+      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto b = Super::coefficient(local);
 
         auto const& rowShapeValues = rowShapeValuesCache[iq];
diff --git a/src/amdis/localoperators/ZeroOrderTestvec.hpp b/src/amdis/localoperators/ZeroOrderTestvec.hpp
index e5ff5403d2efd4719068abfd25d71beb7da164bd..4b09b688fd250ee6f69ce1c383d42635794c76f1 100644
--- a/src/amdis/localoperators/ZeroOrderTestvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestvec.hpp
@@ -20,14 +20,13 @@ namespace AMDiS
 
 
   /// zero-order vector-operator \f$ (\mathbf{b}\cdot\Psi) \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::testvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::testvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
 
@@ -36,29 +35,27 @@ namespace AMDiS
       : Super(expr, 0)
     {}
 
-    template <class Context, class Node, class ElementVector>
-    void getElementVector(Context const& context,
-                          Node const& node,
-                          ElementVector& elementVector)
+    template <class CG, class Node, class Vec>
+    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
     {
       static_assert(Node::isPower,
         "Operator can be applied to Power-Nodes only.");
 
       static const std::size_t CHILDREN = Node::CHILDREN;
 
-      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
       auto const& localFE = node.child(0).finiteElement();
       std::size_t size = localFE.size();
 
       NodeQuadCache<typename Node::ChildType> cache(localFE.localBasis());
 
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context,quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto exprValue = Super::coefficient(local);
 
         auto const& shapeValues = shapeValuesCache[iq];
diff --git a/src/amdis/localoperators/ZeroOrderTestvecTrial.hpp b/src/amdis/localoperators/ZeroOrderTestvecTrial.hpp
index 4b3a6da4a1352e8e705f15734300e3f4373c4083..96938b52edc1ddbb00c2ca5d4679a9125c7f57dd 100644
--- a/src/amdis/localoperators/ZeroOrderTestvecTrial.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestvecTrial.hpp
@@ -18,13 +18,13 @@ namespace AMDiS
 
 
   /// zero-order operator \f$ \langle\Psi, \mathbf{b}\,\phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::testvec_trial, LocalContext, GridFct>
-      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::testvec_trial, LocalContext, GridFct>,
-                                              GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::testvec_trial, LC, GridFct>
+      : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::testvec_trial, LC, GridFct>,
+                                              GridFunctionOperator<tag::test_trialvec, LC, GridFct>>
   {
     using Self = GridFunctionOperator;
-    using Transposed = GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>;
+    using Transposed = GridFunctionOperator<tag::test_trialvec, LC, GridFct>;
     using Super = GridFunctionOperatorTransposed<Self, Transposed>;
 
   public:
diff --git a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
index 1a34286c232164f791d158e262c76bc64655de07..dab0476f1527d4eedb87831db68664a0fe7cf85e 100644
--- a/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
+++ b/src/amdis/localoperators/ZeroOrderTestvecTrialvec.hpp
@@ -21,14 +21,13 @@ namespace AMDiS
 
 
   /// zero-order operator \f$ \langle\Psi, c\,\Phi\rangle \f$, or \f$ \langle\Psi, A\,\Phi\rangle \f$
-  template <class LocalContext, class GridFct>
-  class GridFunctionOperator<tag::testvec_trialvec, LocalContext, GridFct>
-      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_trialvec, LocalContext, GridFct>,
-                                        LocalContext, GridFct>
+  template <class LC, class GridFct>
+  class GridFunctionOperator<tag::testvec_trialvec, LC, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_trialvec, LC, GridFct>, LC, GridFct>
   {
-    static const int dow = ContextGeometry<LocalContext>::dow;
+    static const int dow = ContextGeometry<LC>::dow;
     using Self = GridFunctionOperator;
-    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;
 
     using expr_value_type = typename GridFct::Range;
     static_assert( Size_v<expr_value_type> == 1 || (Rows_v<expr_value_type> == dow && Cols_v<expr_value_type> == dow),
@@ -39,56 +38,50 @@ namespace AMDiS
       : Super(expr, 0)
     {}
 
-    template <class Context, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrix(Context const& context,
-                          RowNode const& rowNode,
-                          ColNode const& colNode,
-                          ElementMatrix& elementMatrix)
+    template <class CG, class RN, class CN, class Mat>
+    void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
     {
-      static_assert(RowNode::isPower && ColNode::isPower,
+      static_assert(RN::isPower && CN::isPower,
         "Operator can be applied to Power-Nodes only.");
 
-      static_assert( (RowNode::CHILDREN == ColNode::CHILDREN), "" );
+      static_assert( (RN::CHILDREN == CN::CHILDREN), "" );
       // theoretically the restriction of the same nr of childs would not be necessary
 
-      const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+      const bool sameFE = std::is_same<FiniteElementType_t<RN>, FiniteElementType_t<CN>>::value;
       const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
 
       using Category = ValueCategory_t<expr_value_type>;
 
-      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+      auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
       if (sameFE && sameNode && std::is_same<Category,tag::scalar>::value)
-        getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix, Category{});
+        getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
       else
-        getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix, Category{});
+        getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
     }
 
   protected: // specialization for different value-categories
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixStandard(Context const& context,
-                                  QuadratureRule const& quad,
-                                  RowNode const& rowNode,
-                                  ColNode const& colNode,
-                                  ElementMatrix& elementMatrix,
-                                  tag::scalar)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
+                                  RN const& rowNode, CN const& colNode,
+                                  Mat& elementMatrix, tag::scalar)
     {
-      static const std::size_t CHILDREN = RowNode::CHILDREN;
+      static const std::size_t CHILDREN = RN::CHILDREN;
 
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
 
-      NodeQuadCache<typename RowNode::ChildType> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<typename ColNode::ChildType> colCache(colLocalFE.localBasis());
+      NodeQuadCache<typename RN::ChildType> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<typename CN::ChildType> colCache(colLocalFE.localBasis());
 
-      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(context,quad);
-      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(context,quad);
+      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo,quad);
+      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         auto const& rowShapeValues = rowShapeValuesCache[iq];
         auto const& colShapeValues = colShapeValuesCache[iq];
@@ -111,27 +104,24 @@ namespace AMDiS
     }
 
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   RowNode const& node,
-                                   ColNode const& /*colNode*/,
-                                   ElementMatrix& elementMatrix,
-                                   tag::scalar)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   RN const& node, CN const& /*colNode*/,
+                                   Mat& elementMatrix, tag::scalar)
     {
-      static const std::size_t CHILDREN = RowNode::CHILDREN;
+      static const std::size_t CHILDREN = RN::CHILDREN;
 
       auto const& localFE = node.child(0).finiteElement();
 
-      NodeQuadCache<typename RowNode::ChildType> cache(localFE.localBasis());
+      NodeQuadCache<typename RN::ChildType> cache(localFE.localBasis());
 
-      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context,quad);
+      auto const& shapeValuesCache = cache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
 
         auto const& shapeValues = shapeValuesCache[iq];
 
@@ -158,31 +148,28 @@ namespace AMDiS
       }
     }
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixStandard(Context const& context,
-                                  QuadratureRule const& quad,
-                                  RowNode const& rowNode,
-                                  ColNode const& colNode,
-                                  ElementMatrix& elementMatrix,
-                                  tag::matrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
+                                  RN const& rowNode, CN const& colNode,
+                                  Mat& elementMatrix, tag::matrix)
     {
-      static const std::size_t CHILDREN = RowNode::CHILDREN;
+      static const std::size_t CHILDREN = RN::CHILDREN;
       static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)>>::value, "" );
 
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
 
-      NodeQuadCache<typename RowNode::ChildType> rowCache(rowLocalFE.localBasis());
-      NodeQuadCache<typename ColNode::ChildType> colCache(colLocalFE.localBasis());
+      NodeQuadCache<typename RN::ChildType> rowCache(rowLocalFE.localBasis());
+      NodeQuadCache<typename CN::ChildType> colCache(colLocalFE.localBasis());
 
-      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(context,quad);
-      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(context,quad);
+      auto const& rowShapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo,quad);
+      auto const& colShapeValuesCache = colCache.evaluateFunctionAtQP(contextGeo,quad);
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(quad[iq].position());
+        decltype(auto) local = contextGeo.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
         const Dune::FieldMatrix<double, int(CHILDREN), int(CHILDREN)> exprValue = Super::coefficient(local);
 
         auto const& rowShapeValues = rowShapeValuesCache[iq];
@@ -208,13 +195,10 @@ namespace AMDiS
       }
     }
 
-    template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
-    void getElementMatrixOptimized(Context const& context,
-                                   QuadratureRule const& quad,
-                                   RowNode const& node,
-                                   ColNode const& /*colNode*/,
-                                   ElementMatrix& elementMatrix,
-                                   tag::matrix)
+    template <class CG, class QR, class RN, class CN, class Mat>
+    void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
+                                   RN const& node, CN const& /*colNode*/,
+                                   Mat& elementMatrix, tag::matrix)
     {
       assert(false && "Not implemented.");
     }
diff --git a/src/amdis/utility/AssembleOperators.hpp b/src/amdis/utility/AssembleOperators.hpp
index 955b85d489029bb995abdd8ee49d5c7241623e90..f252bdbe9a4ba08089bb272295a9d80067a6c12b 100644
--- a/src/amdis/utility/AssembleOperators.hpp
+++ b/src/amdis/utility/AssembleOperators.hpp
@@ -17,12 +17,10 @@ namespace AMDiS
   }
 
 
-  template <class GridView, class Element, class Operators, class EA, class IA, class BA>
-  void assembleOperators(
-      GridView const& gridView,
-      Element const& element,
-      Operators& operators,
-      AssemblerTriple<EA,IA,BA> const& assemblerTriple)
+  template <class GV, class Element, class Operators, class EA, class IA, class BA>
+  void assembleOperators(GV const& gridView, Element const& element,
+                         Operators& operators,
+                         AssemblerTriple<EA,IA,BA> const& assemblerTriple)
   {
     // assemble element operators
     assemblerTriple.elementAssembler(element, operators.onElement());
@@ -41,8 +39,8 @@ namespace AMDiS
   }
 
 
-  template <class Node, class ElementVector>
-  auto makeVectorAssembler(Node const& node, ElementVector& elementVector)
+  template <class Node, class Vec>
+  auto makeVectorAssembler(Node const& node, Vec& elementVector)
   {
     return makeAssemblerTriple(
       [&](auto const& element, auto& operators) {
@@ -62,8 +60,8 @@ namespace AMDiS
   }
 
 
-  template <class RowNode, class ColNode, class ElementMatrix>
-  auto makeMatrixAssembler(RowNode const& rowNode, ColNode const& colNode, ElementMatrix& elementMatrix)
+  template <class RN, class CN, class Mat>
+  auto makeMatrixAssembler(RN const& rowNode, CN const& colNode, Mat& elementMatrix)
   {
     return makeAssemblerTriple(
       [&](auto const& element, auto& operators) {
diff --git a/src/amdis/utility/LocalBasisCache.hpp b/src/amdis/utility/LocalBasisCache.hpp
index 98ae9313ae24b38fe3daa50d364246f6dcb9a177..1ee176560a5e66a696de62a4c0e82c5893832e14 100644
--- a/src/amdis/utility/LocalBasisCache.hpp
+++ b/src/amdis/utility/LocalBasisCache.hpp
@@ -185,29 +185,29 @@ namespace AMDiS
     }
 
     /// Evaluate local basis functions at all quadrature points of a quadrature rule.
-    template <class LocalContext, class QuadratureRule>
-    ShapeValuesContainer const& evaluateFunctionAtQP(LocalContext const& context, QuadratureRule const& quad) const
+    template <class CG, class QR>
+    ShapeValuesContainer const& evaluateFunctionAtQP(CG const& contextGeo, QR const& quad) const
     {
-      Key key{context.type().id(), quad.order(), quad.size()};
+      Key key{contextGeo.type().id(), quad.order(), quad.size()};
       return shapeValuesContainer_.get(key, [&](Key const&)
       {
         ShapeValuesContainer data(quad.size());
         for (std::size_t iq = 0; iq < quad.size(); ++iq)
-          localBasis_->evaluateFunction(context.local(quad[iq].position()), data[iq]);
+          localBasis_->evaluateFunction(contextGeo.local(quad[iq].position()), data[iq]);
         return data;
       });
     }
 
     /// Evaluate local basis jacobians at all quadrature points of a quadrature rule.
-    template <class LocalContext, class QuadratureRule>
-    ShapeGradientsContainer const& evaluateJacobianAtQP(LocalContext const& context, QuadratureRule const& quad) const
+    template <class CG, class QR>
+    ShapeGradientsContainer const& evaluateJacobianAtQP(CG const& contextGeo, QR const& quad) const
     {
-      Key key{context.type().id(), quad.order(), quad.size()};
+      Key key{contextGeo.type().id(), quad.order(), quad.size()};
       return shapeGradientsContainer_.get(key, [&](Key const&)
       {
         ShapeGradientsContainer data(quad.size());
         for (std::size_t iq = 0; iq < quad.size(); ++iq)
-          localBasis_->evaluateJacobian(context.local(quad[iq].position()), data[iq]);
+          localBasis_->evaluateJacobian(contextGeo.local(quad[iq].position()), data[iq]);
         return data;
       });
     }
diff --git a/src/amdis/utility/QuadratureFactory.hpp b/src/amdis/utility/QuadratureFactory.hpp
index 74d8c21df223f0813d7a8a49c3b8ca0492e07a7e..a243367f4825ba67b7ce78b04b488240f3d5c10b 100644
--- a/src/amdis/utility/QuadratureFactory.hpp
+++ b/src/amdis/utility/QuadratureFactory.hpp
@@ -10,11 +10,11 @@
 namespace AMDiS
 {
   /// Base class for quadrature factories for localFunctions
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   class QuadratureFactory
   {
   protected:
-    using QuadratureRule = Dune::QuadratureRule<ctype, dimension>;
+    using QuadratureRule = Dune::QuadratureRule<ctype, dim>;
 
   public:
     virtual ~QuadratureFactory() = default;
@@ -33,7 +33,7 @@ namespace AMDiS
     /// polynomial degree of the integrand function.
     virtual QuadratureRule const& rule(Dune::GeometryType const& type, int degree) const
     {
-      using QuadratureRules = Dune::QuadratureRules<ctype, dimension>;
+      using QuadratureRules = Dune::QuadratureRules<ctype, dim>;
       return QuadratureRules::rule(type, degree);
     }
   };
@@ -45,9 +45,9 @@ namespace AMDiS
 
   /// \brief Factory for quadrature rule, that calculates the coefficient order from
   /// a localFunction passed to the bind method.
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   class QuadFactoryFromLocalFunction
-      : public QuadratureFactory<ctype, dimension, LocalFunction>
+      : public QuadratureFactory<ctype, dim, LocalFunction>
   {
     using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
     static_assert(Concept::value,
@@ -75,16 +75,16 @@ namespace AMDiS
   }
 
   /// Generator for \ref QuadFactoryFromLocalFunction. \relates QuadFactoryFromLocalFunction
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   auto makeQuadratureFactory(PreQuadFactoryFromLocalFunction const& /*pre*/)
   {
-    return QuadFactoryFromLocalFunction<ctype, dimension, LocalFunction>{};
+    return QuadFactoryFromLocalFunction<ctype, dim, LocalFunction>{};
   }
 
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   auto makeQuadratureFactoryPtr(PreQuadFactoryFromLocalFunction const& /*pre*/)
   {
-    using Factory = QuadFactoryFromLocalFunction<ctype, dimension, LocalFunction>;
+    using Factory = QuadFactoryFromLocalFunction<ctype, dim, LocalFunction>;
     return std::make_unique<Factory>();
   }
 
@@ -92,11 +92,11 @@ namespace AMDiS
 
   /// \brief Factory for quadrature rule, that takes to order of the localFunction
   /// and a quadrature-type
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   class QuadFactoryFromOrder
-      : public QuadratureFactory<ctype, dimension, LocalFunction>
+      : public QuadratureFactory<ctype, dim, LocalFunction>
   {
-    using Super = QuadratureFactory<ctype, dimension, LocalFunction>;
+    using Super = QuadratureFactory<ctype, dim, LocalFunction>;
     using QuadratureRule = typename Super::QuadratureRule;
 
   public:
@@ -110,7 +110,7 @@ namespace AMDiS
 
     QuadratureRule const& rule(Dune::GeometryType const& type, int degree) const final
     {
-      using QuadratureRules = Dune::QuadratureRules<ctype, dimension>;
+      using QuadratureRules = Dune::QuadratureRules<ctype, dim>;
       return QuadratureRules::rule(type, degree, qt_);
     }
 
@@ -131,27 +131,27 @@ namespace AMDiS
   }
 
   /// Generator for \ref QuadFactoryFromOrder. \relates QuadFactoryFromOrder
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   auto makeQuadratureFactory(PreQuadFactoryFromOrder const& pre)
   {
-    return QuadFactoryFromOrder<ctype, dimension, LocalFunction>{pre.order, pre.qt};
+    return QuadFactoryFromOrder<ctype, dim, LocalFunction>{pre.order, pre.qt};
   }
 
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   auto makeQuadratureFactoryPtr(PreQuadFactoryFromOrder const& pre)
   {
-    using Factory = QuadFactoryFromOrder<ctype, dimension, LocalFunction>;
+    using Factory = QuadFactoryFromOrder<ctype, dim, LocalFunction>;
     return std::make_unique<Factory>(pre.order, pre.qt);
   }
 
 
 
   /// \brief Factory for quadrature rule, that is based on an existing rule
-  template <class ctype, int dimension, class LocalFunction>
+  template <class ctype, int dim, class LocalFunction>
   class QuadFactoryFromRule
-      : public QuadratureFactory<ctype, dimension, LocalFunction>
+      : public QuadratureFactory<ctype, dim, LocalFunction>
   {
-    using Super = QuadratureFactory<ctype, dimension, LocalFunction>;
+    using Super = QuadratureFactory<ctype, dim, LocalFunction>;
     using QuadratureRule = typename Super::QuadratureRule;
 
   public:
@@ -174,23 +174,23 @@ namespace AMDiS
     QuadRule const& rule;
   };
 
-  template <class ctype, int dimension>
-  auto makePreQuadratureFactory(Dune::QuadratureRule<ctype, dimension> const& rule)
+  template <class ctype, int dim>
+  auto makePreQuadratureFactory(Dune::QuadratureRule<ctype, dim> const& rule)
   {
-    return PreQuadFactoryFromRule<Dune::QuadratureRule<ctype, dimension>>{rule};
+    return PreQuadFactoryFromRule<Dune::QuadratureRule<ctype, dim>>{rule};
   }
 
   /// Generator for \ref QuadFactoryFromRule. \relates QuadFactoryFromRule
-  template <class ctype, int dimension, class LocalFunction, class QuadRule>
+  template <class ctype, int dim, class LocalFunction, class QuadRule>
   auto makeQuadratureFactory(PreQuadFactoryFromRule<QuadRule> const& pre)
   {
-    return QuadFactoryFromRule<ctype, dimension, LocalFunction>{pre.rule};
+    return QuadFactoryFromRule<ctype, dim, LocalFunction>{pre.rule};
   }
 
-  template <class ctype, int dimension, class LocalFunction, class QuadRule>
+  template <class ctype, int dim, class LocalFunction, class QuadRule>
   auto makeQuadratureFactoryPtr(PreQuadFactoryFromRule<QuadRule> const& pre)
   {
-    using Factory = QuadFactoryFromRule<ctype, dimension, LocalFunction>;
+    using Factory = QuadFactoryFromRule<ctype, dim, LocalFunction>;
     return std::make_unique<Factory>(pre.rule);
   }