diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index a086f54aa1a410d497e72a6dcd289094c19dbdbd..2e8021df7777c87867511b3045ed85665b30421a 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -4,6 +4,7 @@
 #include <string>
 #include <utility>
 
+#include <dune/common/hybridutilities.hh>
 #include <dune/common/timer.hh>
 #include <dune/typetree/childextraction.hh>
 
@@ -12,7 +13,6 @@
 #include <amdis/Assembler.hpp>
 #include <amdis/GridFunctionOperator.hpp>
 #include <amdis/GridTransferManager.hpp>
-#include <amdis/common/Loops.hpp>
 
 namespace AMDiS {
 
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index 5d6f05e109fb715658987f4ea0d1f856e3074560..27a4610dafce018344467bb502a27228f5194b7f 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -7,7 +7,7 @@
 #include <dune/functions/functionspacebases/pqknodalbasis.hh>
 #include <dune/grid/yaspgrid.hh>
 
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Logical.hpp>
 
 namespace AMDiS
 {
diff --git a/src/amdis/common/Apply.hpp b/src/amdis/common/Apply.hpp
index 7cade13417b9106709376ac7d98b496a76aadb11..21e2fef0993dc8a7c5c31095c8eb89a8704afc8a 100644
--- a/src/amdis/common/Apply.hpp
+++ b/src/amdis/common/Apply.hpp
@@ -4,7 +4,8 @@
 #include <type_traits>
 #include <utility>
 
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Index.hpp>
+#include <amdis/common/StaticSize.hpp>
 #include <amdis/common/TypeTraits.hpp>
 
 namespace AMDiS
@@ -13,14 +14,15 @@ namespace AMDiS
   {
     namespace Impl_
     {
-      template <class F, class Tuple, std::size_t... I>
-      constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, std::index_sequence<I...>)
+      template <class Functor, class Tuple, std::size_t... I>
+      constexpr decltype(auto) apply_impl(Functor&& f, Tuple&& t, std::index_sequence<I...>)
       {
-        return f(std::get<I>(FWD(t))...);
+        using std::get;
+        return f(get<I>(FWD(t))...);
       }
 
-      template <class F, std::size_t I0, std::size_t... I>
-      constexpr decltype(auto) apply_indices_impl(F&& f, index_t<I0>, std::index_sequence<I...>)
+      template <class Functor, std::size_t I0, std::size_t... I>
+      constexpr decltype(auto) apply_indices_impl(Functor&& f, index_t<I0>, std::index_sequence<I...>)
       {
         return f(index_t<I0+I>{}...);
       }
@@ -29,25 +31,25 @@ namespace AMDiS
     template <class F, class Tuple>
     constexpr decltype(auto) apply(F&& f, Tuple&& t)
     {
-        return Impl_::apply_impl(FWD(f), FWD(t),
-            std::make_index_sequence<std::tuple_size<std::remove_reference_t<Tuple>>::value>{});
+      return Impl_::apply_impl(FWD(f), FWD(t),
+          std::make_index_sequence<Size_v<std::remove_reference_t<Tuple>>>{});
     }
 
-    template <class F, class... Args>
-    constexpr decltype(auto) apply_variadic(F&& f, Args&&... args)
+    template <class Functor, class... Args>
+    constexpr decltype(auto) apply_variadic(Functor&& f, Args&&... args)
     {
-        return apply(FWD(f), std::forward_as_tuple(args...));
+      return apply(FWD(f), std::forward_as_tuple(args...));
     }
 
-    template <class F, std::size_t N>
-    constexpr decltype(auto) apply_indices(F&& f, index_t<N>)
+    template <class Functor, std::size_t N>
+    constexpr decltype(auto) apply_indices(Functor&& f, index_t<N>)
     {
       return Impl_::apply_indices_impl(FWD(f), index_t<0>{},
           std::make_index_sequence<N>{});
     }
 
-    template <class F, std::size_t I0, std::size_t I1>
-    constexpr decltype(auto) apply_indices(F&& f, index_t<I0>, index_t<I1>)
+    template <class Functor, std::size_t I0, std::size_t I1>
+    constexpr decltype(auto) apply_indices(Functor&& f, index_t<I0>, index_t<I1>)
     {
       return Impl_::apply_indices_impl(FWD(f), index_t<I0>{},
           std::make_index_sequence<I1-I0>{});
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index c3b5a867cf54695effa45dc949904984920d9c89..1eec42e568a6e7b63a9135be22071466aae8f489 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -12,12 +12,15 @@ install(FILES
     FieldMatVec.hpp
     FieldMatVec.inc.hpp
     Filesystem.hpp
+    ForEach.hpp
+    Index.hpp
     Literals.hpp
+    Logical.hpp
     Loops.hpp
     Math.hpp
-    Mpl.hpp
     MultiTypeMatrix.hpp
     MultiTypeVector.hpp
+    Range.hpp
     StaticSize.hpp
     String.hpp
     Tags.hpp
diff --git a/src/amdis/common/Concepts.hpp b/src/amdis/common/Concepts.hpp
index d17060df9a384a82fbee31f18e9f1b5bc0327853..fc4aba7ad538be05eb359d8f4c6a73029e0bb8df 100644
--- a/src/amdis/common/Concepts.hpp
+++ b/src/amdis/common/Concepts.hpp
@@ -5,7 +5,7 @@
 #include <dune/functions/common/functionconcepts.hh>
 
 #include <amdis/common/ConceptsBase.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Logical.hpp>
 #include <amdis/common/TypeTraits.hpp>
 
 namespace AMDiS
@@ -33,6 +33,18 @@ namespace AMDiS
     struct IsSimilar<Types<A0,As...>, Types<B0,Bs...>>
       : and_t<IsSimilar<A0, B0>::value, IsSimilar<Types<As...>, Types<Bs...>>::value> {};
 
+
+    template <class... Ts>
+    struct IsSame;
+
+    template <class T0, class... Ts>
+    struct IsSame<T0, Ts...>
+      : std::is_same<Types<T0,Ts...>, Types<Ts...,T0>> {};
+
+    template <>
+    struct IsSame<> : std::true_type {};
+
+
     template <class T>
     struct IsReferenceWrapper
       : std::false_type {};
@@ -71,6 +83,9 @@ namespace AMDiS
     } // end namespace Definition
 #endif // DOXYGEN
 
+    /// Types are the same
+    template <class... Ts>
+    constexpr bool Same = Traits::IsSame<Ts...>::value;
 
     /// Types are the same, up to decay of qualifiers
     template <class A, class B>
diff --git a/src/amdis/common/ForEach.hpp b/src/amdis/common/ForEach.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dbe64d7d535132cd9833d2b81f42bc7d595885c2
--- /dev/null
+++ b/src/amdis/common/ForEach.hpp
@@ -0,0 +1,50 @@
+#pragma once
+
+#include <initializer_list>
+
+#include <amdis/common/Apply.hpp>
+#include <amdis/common/Index.hpp>
+#include <amdis/common/Range.hpp>
+
+namespace AMDiS
+{
+  namespace Tools
+  {
+    namespace Impl_
+    {
+      template <class T>
+      void ignored_evaluation(std::initializer_list<T>&&) { /* do nothing */ }
+    }
+
+    template <class Tuple, class Functor>
+    constexpr void for_each(Tuple&& tuple, Functor&& f)
+    {
+  #if AMDIS_HAS_CXX_FOLD_EXPRESSIONS
+      Tools::apply([f=std::move(f)](auto&&... t) { (f(FWD(t)),...); }, tuple);
+  #else
+      Tools::apply([f=std::move(f)](auto&&... t) {
+        Impl_::ignored_evaluation<int>({0, (f(FWD(t)), 0)...});
+      }, tuple);
+  #endif
+    }
+
+    template <std::size_t I0, std::size_t I1, class Functor>
+    constexpr void for_range(index_t<I0> i0, index_t<I1> i1, Functor&& f)
+    {
+      Tools::for_each(range_t<I0,I1>{}, FWD(f));
+    }
+
+    template <std::size_t N, class Functor>
+    constexpr void for_range(index_t<N>, Functor&& f)
+    {
+      Tools::for_each(range_t<0,N>{}, FWD(f));
+    }
+
+    template <std::size_t I0, std::size_t I1, class Functor>
+    constexpr void for_range(Functor&& f)
+    {
+      Tools::for_each(range_t<I0,I1>{}, FWD(f));
+    }
+
+  } // end namespace Tools
+} // end namespace AMDiS
diff --git a/src/amdis/common/Index.hpp b/src/amdis/common/Index.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ab44576f89fe15d5fa44dc2f5e6e431cf18e75b4
--- /dev/null
+++ b/src/amdis/common/Index.hpp
@@ -0,0 +1,45 @@
+#pragma once
+
+// std c++ headers
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace AMDiS
+{
+  // introduce some shortcuts for integral constants
+  // ---------------------------------------------------------------------------
+
+  /// A wrapper for int type
+  template <int I>
+  using int_t = std::integral_constant<int, I>;
+
+  /// Variable template to generate int-type
+  template <int I>
+  constexpr int_t<I> int_ = {};
+
+  /// class that represents a sequence of integers
+  template <int... I>
+  using Ints = std::integer_sequence<int, I...>;
+
+  template <std::size_t I, int... J>
+  auto get(Ints<J...>) { return std::get<I>(std::make_tuple(int_<J>...)); }
+
+
+  /// A wrapper for std::size_t type
+  template <std::size_t I>
+  using index_t = std::integral_constant<std::size_t, I>;
+
+  /// Variable template to generate std::size_t-type
+  template <std::size_t I>
+  constexpr index_t<I> index_ = {};
+
+
+  /// class that represents a sequence of indices
+  template <std::size_t... I>
+  using Indices = std::index_sequence<I...>;
+
+  template <std::size_t I, std::size_t... J>
+  auto get(Indices<J...>) { return std::get<I>(std::make_tuple(index_<J>...)); }
+
+} // end namespace AMDiS
diff --git a/src/amdis/common/Logical.hpp b/src/amdis/common/Logical.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..98e8fac4ee8f32cdc857a53475adea9f5631b366
--- /dev/null
+++ b/src/amdis/common/Logical.hpp
@@ -0,0 +1,97 @@
+#pragma once
+
+// std c++ headers
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace AMDiS
+{
+  /// A wrapper for bool types
+  template <bool B>
+  using bool_t  = std::integral_constant<bool, B>;
+
+  /// Variable template to generate bool type
+  template <bool B>
+  constexpr bool_t<B> bool_ = {};
+
+  // some boolean operations
+  // ---------------------------------------------------------------------------
+
+  namespace Impl
+  {
+    template <bool...> struct all_helper {};
+
+  } // end namespace Impl
+
+  template <bool... Bs>
+  using all_of_t = std::is_same<Impl::all_helper<true, Bs...>, Impl::all_helper<Bs..., true>>;
+
+  template <bool... Bs>
+  constexpr bool all_of_v = all_of_t<Bs...>::value;
+
+  template <bool... Bs>
+  using and_t = all_of_t<Bs...>;
+
+  template <bool... Bs>
+  using none_of_t = std::is_same<Impl::all_helper<false, Bs...>, Impl::all_helper<Bs..., false>>;
+
+  template <bool... Bs>
+  constexpr bool none_of_v = none_of_t<Bs...>::value;
+
+  template <bool... Bs>
+  using any_of_t = bool_t<not none_of_t<Bs...>::value>;
+
+  template <bool... Bs>
+  constexpr bool any_of_v = any_of_t<Bs...>::value;
+
+  template <bool... Bs>
+  using or_t = any_of_t<Bs...>;
+
+
+  template <bool... Bs>
+  constexpr bool_t<and_t<Bs...>::value> and_ = {};
+
+  template <bool B0, bool B1>
+  constexpr bool_t<B0 && B1> operator&&(bool_t<B0>, bool_t<B1>) { return {}; }
+
+
+  template <bool... Bs>
+  constexpr bool_t<or_t<Bs...>::value> or_ = {};
+
+  template <bool B0, bool B1>
+  constexpr bool_t<B0 || B1> operator||(bool_t<B0>, bool_t<B1>) { return {}; }
+
+
+  template <bool B>
+  using not_t = bool_t<!B>;
+
+  template <bool B>
+  constexpr bool_t<not_t<B>::value> not_ = {};
+
+  template <bool B>
+  constexpr bool_t<!B> operator!(bool_t<B>) { return {}; }
+
+
+  namespace Impl
+  {
+    template <class T, T... Is>
+    struct IsEqualImpl;
+
+    template <class T, T I0, T... Is>
+    struct IsEqualImpl<T, I0, Is...>
+        : public std::is_same<std::integer_sequence<T,I0,Is...>,
+                              std::integer_sequence<T,Is...,I0>> {};
+
+    template <class T>
+    struct IsEqualImpl<T> { enum { value = true }; };
+
+  } // end namespace Impl
+
+  template <class T, T... values>
+  using IsEqual = Impl::IsEqualImpl<T,values...>;
+
+  template <class T, class... Ts>
+  using is_one_of = or_t<std::is_same<T, Ts>::value...>;
+
+} // end namespace AMDiS
diff --git a/src/amdis/common/Loops.hpp b/src/amdis/common/Loops.hpp
deleted file mode 100644
index ca004320659762fc02aa97d4053ea1e07e9ce9da..0000000000000000000000000000000000000000
--- a/src/amdis/common/Loops.hpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#pragma once
-
-#include <dune/common/hybridutilities.hh>
-#include <dune/common/rangeutilities.hh>
-
-namespace AMDiS
-{
-  using Dune::Hybrid::forEach;
-
-} // end namespace AMDiS
diff --git a/src/amdis/common/Mpl.hpp b/src/amdis/common/Mpl.hpp
deleted file mode 100644
index 2e49f5a56fbed2962a82976307eccc46271567eb..0000000000000000000000000000000000000000
--- a/src/amdis/common/Mpl.hpp
+++ /dev/null
@@ -1,188 +0,0 @@
-#pragma once
-
-// std c++ headers
-#include <tuple>
-#include <type_traits>
-#include <utility>
-
-namespace AMDiS
-{
-  // introduce some shortcuts for integral constants
-  // ---------------------------------------------------------------------------
-
-  /// A wrapper for int type
-  template <int I>
-  using int_t = std::integral_constant<int, I>;
-
-  /// Variable template to generate int-type
-  template <int I>
-  constexpr int_t<I> int_ = {};
-
-  /// class that represents a sequence of integers
-  template <int... I>
-  using Ints = std::integer_sequence<int, I...>;
-
-  template <std::size_t I, int... J>
-  auto get(Ints<J...>) { return std::get<I>(std::make_tuple(int_<J>...)); }
-
-
-  /// A wrapper for std::size_t type
-  template <std::size_t I>
-  using index_t = std::integral_constant<std::size_t, I>;
-
-  /// Variable template to generate std::size_t-type
-  template <std::size_t I>
-  constexpr index_t<I> index_ = {};
-
-
-  /// class that represents a sequence of indices
-  template <std::size_t... I>
-  using Indices = std::index_sequence<I...>;
-
-  template <std::size_t I, std::size_t... J>
-  auto get(Indices<J...>) { return std::get<I>(std::make_tuple(index_<J>...)); }
-
-
-  /// A wrapper for bool types
-  template <bool B>
-  using bool_t  = std::integral_constant<bool, B>;
-
-  /// Variable template to generate bool type
-  template <bool B>
-  constexpr bool_t<B> bool_ = {};
-
-
-  namespace Impl
-  {
-    /// A range of indices [I,J)
-    template <class Int, Int I, Int J>
-    struct range_impl
-    {
-      using type = range_impl;
-
-      /// Return the first element in the range
-      static constexpr auto begin() { return std::integral_constant<Int, I>{}; }
-
-      /// Returns the element after the last element in the range
-      static constexpr auto end() { return std::integral_constant<Int, J>{}; }
-
-      /// Returns the ith index in the range as integral constant
-      template <std::size_t i>
-      constexpr auto operator[](index_t<i>) const
-      {
-        return std::integral_constant<Int, I+Int(i)>{};
-      }
-
-      /// Return whether the range is empty
-      static constexpr bool empty() { return I >= J; }
-
-      /// Returns the size of the range
-      static constexpr auto size() { return std::integral_constant<Int, J-I>{}; }
-    };
-
-    /// Extracts the Ith element element from the range
-    template <std::size_t I, class Int, Int begin, Int end>
-    constexpr auto get(range_impl<Int, begin, end> const& r) { return r[index_<I>]; }
-
-  } // end namespace Impl
-
-  template <std::size_t I, std::size_t J>
-  using range_t = Impl::range_impl<std::size_t, I, J>;
-
-  template <std::size_t I, std::size_t J>
-  constexpr range_t<I,J> range_ = {};
-
-  // some boolean operations
-  // ---------------------------------------------------------------------------
-
-  namespace Impl
-  {
-    template <bool...> struct all_helper {};
-
-  } // end namespace Impl
-
-  template <bool... Bs>
-  using all_of_t = std::is_same<Impl::all_helper<true, Bs...>, Impl::all_helper<Bs..., true>>;
-
-  template <bool... Bs>
-  constexpr bool all_of_v = all_of_t<Bs...>::value;
-
-  template <bool... Bs>
-  using and_t = all_of_t<Bs...>;
-
-  template <bool... Bs>
-  using none_of_t = std::is_same<Impl::all_helper<false, Bs...>, Impl::all_helper<Bs..., false>>;
-
-  template <bool... Bs>
-  constexpr bool none_of_v = none_of_t<Bs...>::value;
-
-  template <bool... Bs>
-  using any_of_t = bool_t<not none_of_t<Bs...>::value>;
-
-  template <bool... Bs>
-  constexpr bool any_of_v = any_of_t<Bs...>::value;
-
-  template <bool... Bs>
-  using or_t = any_of_t<Bs...>;
-
-
-  template <bool... Bs>
-  constexpr bool_t<and_t<Bs...>::value> and_ = {};
-
-  template <bool B0, bool B1>
-  constexpr bool_t<B0 && B1> operator&&(bool_t<B0>, bool_t<B1>) { return {}; }
-
-
-  template <bool... Bs>
-  constexpr bool_t<or_t<Bs...>::value> or_ = {};
-
-  template <bool B0, bool B1>
-  constexpr bool_t<B0 || B1> operator||(bool_t<B0>, bool_t<B1>) { return {}; }
-
-
-  template <bool B>
-  using not_t = bool_t<!B>;
-
-  template <bool B>
-  constexpr bool_t<not_t<B>::value> not_ = {};
-
-  template <bool B>
-  constexpr bool_t<!B> operator!(bool_t<B>) { return {}; }
-
-
-  namespace Impl
-  {
-    template <class T, T... Is>
-    struct IsEqualImpl;
-
-    template <class T, T I0, T... Is>
-    struct IsEqualImpl<T, I0, Is...>
-        : public std::is_same<std::integer_sequence<T,I0,Is...>,
-                              std::integer_sequence<T,Is...,I0>> {};
-
-    template <class T>
-    struct IsEqualImpl<T> { enum { value = true }; };
-
-    template <class... Ts>
-    struct IsSameImpl;
-
-    template <class T0, class... Ts>
-    struct IsSameImpl<T0, Ts...>
-      : public std::is_same<std::tuple<T0,Ts...>,
-                            std::tuple<Ts...,T0>> {};
-
-    template <>
-    struct IsSameImpl<> { enum { value = true }; };
-
-  } // end namespace Impl
-
-  template <class T, T... values>
-  using IsEqual = Impl::IsEqualImpl<T,values...>;
-
-  template <class T, class... Ts>
-  using is_one_of = or_t<std::is_same<T, Ts>::value...>;
-
-  template <class... Ts>
-  using IsSame = Impl::IsSameImpl<Ts...>;
-
-} // end namespace AMDiS
diff --git a/src/amdis/common/MultiTypeMatrix.hpp b/src/amdis/common/MultiTypeMatrix.hpp
index 392005aefa4c8f896b57b28a97ace714cf4c6f85..e8a46ec8eac568358f53b50f1b53fdd100a36087 100644
--- a/src/amdis/common/MultiTypeMatrix.hpp
+++ b/src/amdis/common/MultiTypeMatrix.hpp
@@ -6,10 +6,11 @@
 #include <dune/common/ftraits.hh>
 
 #include <amdis/common/Concepts.hpp>
-#include <amdis/common/Loops.hpp>
+#include <amdis/common/ForEach.hpp>
+#include <amdis/common/Index.hpp>
 #include <amdis/common/Math.hpp>
-#include <amdis/common/Mpl.hpp>
 #include <amdis/common/MultiTypeVector.hpp>
+#include <amdis/common/Range.hpp>
 #include <amdis/common/StaticSize.hpp>
 #include <amdis/common/TypeTraits.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
@@ -71,35 +72,35 @@ namespace AMDiS
     /// Assignment of real number to all tuple elements
     MultiTypeMatrix& operator=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv = value; });
+      Tools::for_each(*this, [value](auto& fv) { fv = value; });
       return *this;
     }
 
     // Compound assignment operator +=
     MultiTypeMatrix& operator+=(MultiTypeMatrix const& that)
     {
-      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
+      Tools::for_range<0,rows>([&that,this](auto const _i) { (*this)[_i] += that[_i]; });
       return *this;
     }
 
     // Compound assignment operator -=
     MultiTypeMatrix& operator-=(MultiTypeMatrix const& that)
     {
-      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
+      Tools::for_range<0,rows>([&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
       return *this;
     }
 
     // Scaling of all tuple elements by a constant value
     MultiTypeMatrix& operator*=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv *= value; });
+      Tools::for_each(*this, [value](auto& fv) { fv *= value; });
       return *this;
     }
 
     // Scaling of all tuple elements by the inverse of a constant value
     MultiTypeMatrix& operator/=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv /= value; });
+      Tools::for_each(*this, [value](auto& fv) { fv /= value; });
       return *this;
     }
 
@@ -150,4 +151,15 @@ namespace AMDiS
     }
   };
 
+  namespace Impl
+  {
+    template <class... Rows>
+    struct RowsImpl<MultiTypeMatrix<Rows...>>
+        : std::integral_constant<std::size_t, sizeof...(Rows) > {};
+
+    template <class Row0, class... Rows>
+    struct ColsImpl<MultiTypeMatrix<Row0, Rows...>>
+        : SizeImpl<Row0> {};
+  }
+
 } // end namespace AMDiS
diff --git a/src/amdis/common/MultiTypeVector.hpp b/src/amdis/common/MultiTypeVector.hpp
index 1f3e38f57a54795f0cd3deb1ccf6a85846ff8e3f..7bec7543913af53d4bd26c48fc2f267e87869d8d 100644
--- a/src/amdis/common/MultiTypeVector.hpp
+++ b/src/amdis/common/MultiTypeVector.hpp
@@ -7,8 +7,9 @@
 #include <dune/functions/common/indexaccess.hh>
 
 #include <amdis/common/Concepts.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/ForEach.hpp>
+#include <amdis/common/Index.hpp>
+#include <amdis/common/Range.hpp>
 #include <amdis/common/StaticSize.hpp>
 #include <amdis/common/TypeTraits.hpp>
 
@@ -64,35 +65,35 @@ namespace AMDiS
     /// Assignment of real number to all tuple elements
     MultiTypeVector& operator=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv = value; });
+      Tools::for_each(*this, [value](auto& fv) { fv = value; });
       return *this;
     }
 
     // Compound assignment operator +=
     MultiTypeVector& operator+=(MultiTypeVector const& that)
     {
-      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
+      Tools::for_range<0,dimension>([&that,this](auto const _i) { (*this)[_i] += that[_i]; });
       return *this;
     }
 
     // Compound assignment operator -=
     MultiTypeVector& operator-=(MultiTypeVector const& that)
     {
-      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
+      Tools::for_range<0,dimension>([&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
       return *this;
     }
 
     // Scaling of all tuple elements by a constant value
     MultiTypeVector& operator*=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv *= value; });
+      Tools::for_each(*this, [value](auto& fv) { fv *= value; });
       return *this;
     }
 
     // Scaling of all tuple elements by the inverse of a constant value
     MultiTypeVector& operator/=(real_type value)
     {
-      forEach(*this, [value](auto& fv) { fv /= value; });
+      Tools::for_each(*this, [value](auto& fv) { fv /= value; });
       return *this;
     }
 
@@ -133,4 +134,11 @@ namespace AMDiS
     }
   };
 
+  namespace Impl
+  {
+    template <class... Ts>
+    struct SizeImpl<MultiTypeVector<Ts...>>
+        : std::integral_constant<std::size_t, sizeof...(Ts) > {};
+  }
+
 } // end namespace AMDiS
diff --git a/src/amdis/common/Range.hpp b/src/amdis/common/Range.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..69f3d79de60d404f0931a52b9162f768fc8d1fc8
--- /dev/null
+++ b/src/amdis/common/Range.hpp
@@ -0,0 +1,69 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/common/Index.hpp>
+#include <amdis/common/StaticSize.hpp>
+
+namespace AMDiS
+{
+  namespace Impl
+  {
+    /// A range of indices [I,J)
+    template <class Int, Int I, Int J>
+    struct range_impl
+    {
+      using type = range_impl;
+
+      /// Return the first element in the range
+      static constexpr auto begin() { return std::integral_constant<Int, I>{}; }
+
+      /// Returns the element after the last element in the range
+      static constexpr auto end() { return std::integral_constant<Int, J>{}; }
+
+      /// Returns the ith index in the range as integral constant
+      template <std::size_t i>
+      constexpr auto operator[](index_t<i>) const
+      {
+        return std::integral_constant<Int, I+Int(i)>{};
+      }
+
+      /// Return whether the range is empty
+      static constexpr bool empty() { return I >= J; }
+
+      /// Returns the size of the range
+      static constexpr auto size() { return std::integral_constant<Int, J-I>{}; }
+    };
+
+    /// Extracts the Ith element element from the range
+    template <std::size_t I, class Int, Int begin, Int end>
+    constexpr auto get(range_impl<Int, begin, end> const& r) { return r[index_<I>]; }
+
+    /// Return the size of the range
+    template <class Int, Int I, Int J>
+    struct SizeImpl<range_impl<Int,I,J>>
+        : std::integral_constant<std::size_t, std::size_t(J-I)> {};
+
+  } // end namespace Impl
+
+  template <std::size_t I, std::size_t J>
+  using range_t = Impl::range_impl<std::size_t, I, J>;
+
+  template <std::size_t I, std::size_t J>
+  constexpr range_t<I,J> range_ = {};
+
+} // end namespace AMDiS
+
+namespace std
+{
+  template <class Int, Int I0, Int I1>
+  struct tuple_size<AMDiS::Impl::range_impl<Int,I0,I1>>
+      : std::integral_constant<std::size_t, std::size_t(I1-I0)> {};
+
+  template <std::size_t I, class Int, Int I0, Int I1>
+  struct tuple_element<I,AMDiS::Impl::range_impl<Int,I0,I1>>
+  {
+    using type = Int;
+  };
+
+} // end namespace std
diff --git a/src/amdis/common/StaticSize.hpp b/src/amdis/common/StaticSize.hpp
index d1b399e15f4c750b4e2e8c7421d74bc397fdb5ec..ec346a0ec753bceebb3f4c2fef67f239ff74a0b0 100644
--- a/src/amdis/common/StaticSize.hpp
+++ b/src/amdis/common/StaticSize.hpp
@@ -4,6 +4,8 @@
 #include <tuple>
 #include <type_traits>
 
+#include <amdis/common/TypeTraits.hpp>
+
 namespace Dune
 {
   // forward declarations
@@ -25,14 +27,6 @@ namespace Dune
 
 namespace AMDiS
 {
-  // forward declarations
-  template <class... Rows>
-  class MultiTypeMatrix;
-
-  template <class... Ts>
-  class MultiTypeVector;
-
-
   namespace Impl
   {
     template <class Tuple, class = void>
@@ -58,10 +52,6 @@ namespace AMDiS
     struct SizeImpl<Dune::FieldMatrix<T,N,M>>
         : std::integral_constant<std::size_t, N*M> {};
 
-    template <class... Ts>
-    struct SizeImpl<MultiTypeVector<Ts...>>
-        : std::integral_constant<std::size_t, sizeof...(Ts) > {};
-
     template <class... Ts>
     struct SizeImpl<Dune::TupleVector<Ts...>>
         : std::integral_constant<std::size_t, sizeof...(Ts)> {};
@@ -94,10 +84,6 @@ namespace AMDiS
     struct RowsImpl<Dune::FieldMatrix<T,N,M>>
         : std::integral_constant<std::size_t, N> {};
 
-    template <class... Rows>
-    struct RowsImpl<MultiTypeMatrix<Rows...>>
-        : std::integral_constant<std::size_t, sizeof...(Rows) > {};
-
     template <class... Rows>
     struct RowsImpl<Dune::MultiTypeBlockMatrix<Rows...>>
         : std::integral_constant<std::size_t, sizeof...(Rows)> {};
@@ -126,10 +112,6 @@ namespace AMDiS
     struct ColsImpl<Dune::FieldMatrix<T,N,M>>
         : std::integral_constant<std::size_t, M> {};
 
-    template <class Row0, class... Rows>
-    struct ColsImpl<MultiTypeMatrix<Row0, Rows...>>
-        : SizeImpl<Row0> {};
-
     template <class Row0, class... Rows>
     struct ColsImpl<Dune::MultiTypeBlockMatrix<Row0, Rows...>>
         : SizeImpl<Row0> {};
diff --git a/src/amdis/common/TypeTraits.hpp b/src/amdis/common/TypeTraits.hpp
index 721ec6f035b50dc43446919ea1d9b4d3ab538924..5d63c0775d8869ea38258351ac6e57bd3a88ed17 100644
--- a/src/amdis/common/TypeTraits.hpp
+++ b/src/amdis/common/TypeTraits.hpp
@@ -2,13 +2,6 @@
 
 #include <memory>
 #include <type_traits>
-#include <vector>
-
-#if AMDIS_HAS_CXX_CONSTEXPR_IF
-#define IF_CONSTEXPR if constexpr
-#else
-#define IF_CONSTEXPR if
-#endif
 
 namespace AMDiS
 {
diff --git a/src/amdis/gridfunctions/FunctorGridFunction.hpp b/src/amdis/gridfunctions/FunctorGridFunction.hpp
index 628cc050116118ec556e159838049c13414e0c39..2054b342e6be260c380743b86f086675d1f34bba 100644
--- a/src/amdis/gridfunctions/FunctorGridFunction.hpp
+++ b/src/amdis/gridfunctions/FunctorGridFunction.hpp
@@ -6,8 +6,8 @@
 #include <dune/common/std/apply.hh>
 
 #include <amdis/Operations.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/ForEach.hpp>
+#include <amdis/common/Logical.hpp>
 #include <amdis/gridfunctions/GridFunctionConcepts.hpp>
 
 namespace AMDiS
@@ -74,7 +74,7 @@ namespace AMDiS
     template <class Element>
     void bind(Element const& element)
     {
-      forEach(localFcts_, [&](auto& localFct) {
+      Tools::for_each(localFcts_, [&](auto& localFct) {
         (*localFct).bind(element);
       });
     }
@@ -82,7 +82,7 @@ namespace AMDiS
     /// Calls \ref unbind for all localFunctions
     void unbind()
     {
-      forEach(localFcts_, [&](auto& localFct) {
+      Tools::for_each(localFcts_, [&](auto& localFct) {
         (*localFct).unbind();
       });
     }
@@ -90,8 +90,7 @@ namespace AMDiS
     /// Applies the functor \ref fct_ to the evaluated localFunctions
     Range operator()(Domain const& x) const
     {
-      using Dune::Std::apply;
-      return apply([&](auto&&... localFct) { return fct_((*localFct)(x)...); }, localFcts_);
+      return Tools::apply([&](auto&&... localFct) { return fct_((*localFct)(x)...); }, localFcts_);
     }
 
   public:
@@ -130,8 +129,7 @@ namespace AMDiS
     // d_i(f)[lgfs...] * lgfs_i
     auto term_i = [&](auto const _i)
     {
-      using Dune::Std::apply;
-      auto di_f = apply([&](auto const&... lgfs) {
+      auto di_f = Tools::apply([&](auto const&... lgfs) {
         return makeFunctorGridFunction(partial(lf.fct(), _i), (*lgfs)...);
       }, lf.localFcts());
 
@@ -141,8 +139,7 @@ namespace AMDiS
     };
 
     // sum_i [ d_i(f)[lgfs...] * derivative(lgfs_i)
-    using Dune::Std::apply;
-    auto gridFct = apply([&](auto const... _i)
+    auto gridFct = Tools::apply([&](auto const... _i)
     {
       return makeFunctorGridFunction(Operation::Plus{}, term_i(_i)...);
     }, index_seq);
@@ -164,8 +161,7 @@ namespace AMDiS
              && all_of_v<Concepts::HasOrder<LFs>...>)>
   int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
   {
-    using Dune::Std::apply;
-    return apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(*lgfs)...); },
+    return Tools::apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(*lgfs)...); },
       lf.localFcts());
   }
 
@@ -221,21 +217,18 @@ namespace AMDiS
     /// Applies the functor to the evaluated gridfunctions
     Range operator()(Domain const& x) const
     {
-      using Dune::Std::apply;
-      return apply([&](auto&&... gridFct) { return fct_(gridFct(x)...); }, gridFcts_);
+      return Tools::apply([&](auto&&... gridFct) { return fct_(gridFct(x)...); }, gridFcts_);
     }
 
     /// Return the stored \ref EntitySet of the first GridFunction
     EntitySet const& entitySet() const
     {
-      using std::get;
-      return get<0>(gridFcts_).entitySet();
+      return std::get<0>(gridFcts_).entitySet();
     }
 
     LocalFunction localFunction() const
     {
-      using Dune::Std::apply;
-      return apply([&](auto const&... gridFcts) { return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
+      return Tools::apply([&](auto const&... gridFcts) { return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
     }
 
   private:
@@ -278,7 +271,7 @@ namespace AMDiS
       template <class GridView>
       static auto create(Self const& self, GridView const& gridView)
       {
-        return Dune::Std::apply([&](auto const&... pgf) {
+        return Tools::apply([&](auto const&... pgf) {
           return makeFunctorGridFunction(self.fct_,
             makeGridFunction(pgf, gridView)...);
         }, self.preGridFcts_);
diff --git a/src/amdis/gridfunctions/GridFunctionConcepts.hpp b/src/amdis/gridfunctions/GridFunctionConcepts.hpp
index 59d7e2346302a8cb8b1adfe3c71adb5503f9d470..772716f8d5bba1f6cfd722ab490bcbca5c319ace 100644
--- a/src/amdis/gridfunctions/GridFunctionConcepts.hpp
+++ b/src/amdis/gridfunctions/GridFunctionConcepts.hpp
@@ -3,7 +3,7 @@
 #include <dune/common/typeutilities.hh>
 
 #include <amdis/common/Concepts.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Logical.hpp>
 
 namespace AMDiS
 {
diff --git a/src/amdis/linearalgebra/eigen/Constraints.hpp b/src/amdis/linearalgebra/eigen/Constraints.hpp
index f183586f9edf3ecf758fe7979e3d6403fea02102..46b99b92cb1700a42f5c9e12102b8c48587ea1b4 100644
--- a/src/amdis/linearalgebra/eigen/Constraints.hpp
+++ b/src/amdis/linearalgebra/eigen/Constraints.hpp
@@ -4,7 +4,7 @@
 
 #include <Eigen/SparseCore>
 
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Index.hpp>
 #include <amdis/linearalgebra/Common.hpp>
 #include <amdis/linearalgebra/Constraints.hpp>
 
diff --git a/src/amdis/localoperators/ConvectionDiffusionOperator.hpp b/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
index e70558d220b5919fd6c8ef499bf89d983f835f37..3022a5b16dcedbc613e0485f876bb120f940ddaa 100644
--- a/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
+++ b/src/amdis/localoperators/ConvectionDiffusionOperator.hpp
@@ -2,6 +2,8 @@
 
 #include <type_traits>
 
+#include <dune/common/hybridutilities.hh>
+
 #include <amdis/utility/LocalBasisCache.hpp>
 #include <amdis/LocalOperator.hpp>
 #include <amdis/Output.hpp>
@@ -106,7 +108,7 @@ namespace AMDiS
         WorldVector b = makeB(localFctB(local));
         const auto c = localFctC(local);
 
-        IF_CONSTEXPR(conserving) {
+        if (conserving) {
           WorldVector gradAi;
           for (std::size_t i = 0; i < numLocalFe; ++i) {
             const auto local_i = rowNode.localIndex(i);
diff --git a/src/amdis/operations/Basic.hpp b/src/amdis/operations/Basic.hpp
index 304023d1af9fef68e46f171088aa51b813204b0e..9e214b89f6e4c790ab7ae27b85c3952f81052800 100644
--- a/src/amdis/operations/Basic.hpp
+++ b/src/amdis/operations/Basic.hpp
@@ -3,7 +3,7 @@
 #include <algorithm>
 #include <utility>
 
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Index.hpp>
 #include <amdis/common/ConceptsBase.hpp>
 
 namespace AMDiS
diff --git a/src/amdis/operations/Composer.hpp b/src/amdis/operations/Composer.hpp
index 586cc45b0fbae9a8bdc88aea904aeff6b1a5d574..0cdc49cdb98d23e9b363323bd9c91b0419cd53d4 100644
--- a/src/amdis/operations/Composer.hpp
+++ b/src/amdis/operations/Composer.hpp
@@ -5,8 +5,9 @@
 
 #include <dune/common/std/apply.hh>
 
+#include <amdis/common/Apply.hpp>
 #include <amdis/common/Concepts.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/Logical.hpp>
 #include <amdis/operations/Basic.hpp>
 
 namespace AMDiS
@@ -40,7 +41,7 @@ namespace AMDiS
       constexpr auto operator()(Ts const&... args) const
       {
         auto eval = [&](auto const& g) { return g(args...); };
-        return Dune::Std::apply([&,this](auto const&... gs) { return f_(eval(gs)...); }, gs_);
+        return Tools::apply([&,this](auto const&... gs) { return f_(eval(gs)...); }, gs_);
       }
 
       F f_;
@@ -75,7 +76,7 @@ namespace AMDiS
     int order(Composer<F,Gs...> const& c, Int... degrees)
     {
       auto deg = [&](auto const& g) { return order(g, int(degrees)...); };
-      return Dune::Std::apply([&](auto const&... gs) { return order(c.f_, deg(gs)...); }, c.gs_);
+      return Tools::apply([&](auto const&... gs) { return order(c.f_, deg(gs)...); }, c.gs_);
     }
 
     /// Partial derivative of composed function:
diff --git a/src/amdis/typetree/FiniteElementType.hpp b/src/amdis/typetree/FiniteElementType.hpp
index a2a79defb90052b085c082a96f06d4d848456dae..682a74472f2912cbf7f3d169ce40c1fb8fa45474 100644
--- a/src/amdis/typetree/FiniteElementType.hpp
+++ b/src/amdis/typetree/FiniteElementType.hpp
@@ -5,8 +5,8 @@
 #include <dune/common/typetraits.hh>
 #include <dune/typetree/nodetags.hh>
 
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/Mpl.hpp>
+#include <amdis/common/ForEach.hpp>
+#include <amdis/common/Index.hpp>
 #include <amdis/common/Tags.hpp>
 
 namespace AMDiS
@@ -73,7 +73,7 @@ namespace AMDiS
       static int order(Node const& node)
       {
         int degree = 0;
-        forEach(range_<0,Node::CHILDREN>, [&](auto const _i) {
+        Tools::for_range<0,Node::CHILDREN>([&](auto const _i) {
           degree = std::max(degree, polynomialDegree(node.child(_i)));
         });
         return degree;
diff --git a/src/amdis/typetree/Traversal.hpp b/src/amdis/typetree/Traversal.hpp
index a0255eb3112c2f1c6176f68f09c40db090be7320..45751bcdd23495dc9efcf53dcf4e8ba823e543ef 100644
--- a/src/amdis/typetree/Traversal.hpp
+++ b/src/amdis/typetree/Traversal.hpp
@@ -1,10 +1,13 @@
 #pragma once
 
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/rangeutilities.hh>
+
 #include <dune/typetree/nodetags.hh>
 #include <dune/typetree/treepath.hh>
 #include <dune/typetree/visitor.hh>
 
-#include <amdis/common/Loops.hpp>
+#include <amdis/common/ForEach.hpp>
 #include <amdis/common/TypeTraits.hpp>
 
 namespace AMDiS
@@ -82,7 +85,7 @@ namespace AMDiS
       v.pre(FWD(n),tp);
 
       auto const deg = hybridDegree(NodeTag{}, n);
-      forEach(Dune::range(deg), [&](auto const _k)
+      Dune::Hybrid::forEach(Dune::range(deg), [&](auto const _k)
       {
         // always call beforeChild(), regardless of the value of visit
         v.beforeChild(FWD(n),n.child(_k),tp,_k);