diff --git a/dune/solvers/common/algorithm.hh b/dune/solvers/common/algorithm.hh
index d40eab96cbb340fc7e432c7239a2518cc8db9a93..99d238ca811809a49aac9e11e1e3c3120ad942e2 100644
--- a/dune/solvers/common/algorithm.hh
+++ b/dune/solvers/common/algorithm.hh
@@ -140,68 +140,66 @@ constexpr auto elementAt(Container&& c, Index&& i)
 
 namespace Imp {
 
-template<class Begin, class End>
-class StaticIntegralRange
-{
-public:
-
-  template<std::size_t i>
-  constexpr auto operator[](Dune::index_constant<i>) const
+  template<class Begin, class End>
+  class StaticIntegralRange
   {
-    return Dune::index_constant<Begin::value+i>();
-  }
+  public:
+
+    template<std::size_t i>
+    constexpr auto operator[](Dune::index_constant<i>) const
+    {
+      return Dune::index_constant<Begin::value+i>();
+    }
+
+    static constexpr auto size()
+    {
+      return std::integral_constant<typename Begin::value_type, End::value - Begin::value>();
+    }
+  };
 
-  static constexpr auto size()
+  template<class T>
+  class DynamicIntegralRange
   {
-    return std::integral_constant<typename Begin::value_type, End::value - Begin::value>();
-  }
-};
+  public:
+    constexpr DynamicIntegralRange(const T& begin, const T& end):
+      begin_(begin),
+      end_(end)
+    {}
 
-template<class T>
-class DynamicIntegralRange
-{
-public:
-  constexpr DynamicIntegralRange(const T& begin, const T& end):
-    begin_(begin),
-    end_(end)
-  {}
+    const T& begin() const
+    { return begin_; }
 
-  const T& begin() const
-  { return begin_; }
+    const T& end() const
+    { return end_; }
 
-  const T& end() const
-  { return end_; }
+    constexpr auto size() const
+    {
+      return end() - begin();
+    }
 
-  constexpr auto size() const
+    constexpr T operator[](const T&i) const
+    { return begin()+i; }
+
+  private:
+    T begin_;
+    T end_;
+  };
+
+  template<class Begin, class End,
+    std::enable_if_t<IsIntegralConstant<Begin>::value and IsIntegralConstant<End>::value, int> = 0>
+  constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<1>&)
   {
-    return end() - begin();
+    static_assert(Begin::value <= End::value, "You cannot create an integralRange where end<begin");
+    return Imp::StaticIntegralRange<Begin,End>();
   }
 
-  constexpr T operator[](const T&i) const
+  template<class Begin, class End>
+  constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<0>&)
   {
-    return begin()+i;
+    assert(begin <= end);
+    return Imp::DynamicIntegralRange<End>(begin, end);
   }
 
-private:
-  T begin_;
-  T end_;
-};
-
-template<class Begin, class End,
-  std::enable_if_t<IsIntegralConstant<Begin>::value and IsIntegralConstant<End>::value, int> = 0>
-constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<1>&)
-{
-  static_assert(Begin::value <= End::value, "You cannot create an integralRange where end<begin");
-  return Imp::StaticIntegralRange<Begin,End>();
-}
-
-template<class Begin, class End>
-constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<0>&)
-{
-  assert(begin <= end);
-  return Imp::DynamicIntegralRange<End>(begin, end);
-}
-
 } // namespace Imp
 
 
@@ -248,40 +246,59 @@ constexpr auto integralRange(const End& end)
 
 namespace Imp {
 
-template<class Range, class F, class Index, Index... i>
-constexpr void forEachIndex(Range&& range, F&& f, std::integer_sequence<Index, i...>)
-{
-  std::initializer_list<int>{(f(HybridAlgorithm::elementAt(range, std::integral_constant<Index,i>())), 0)...};
-}
+  template<class Range, class F, class Index, Index... i>
+  constexpr void forEachIndex(Range&& range, F&& f, std::integer_sequence<Index, i...>)
+  {
+    std::initializer_list<int>{(f(HybridAlgorithm::elementAt(range, std::integral_constant<Index,i>())), 0)...};
+  }
 
-template<class Range, class F,
-  std::enable_if_t<IsIntegralConstant<decltype(HybridAlgorithm::size(std::declval<Range>()))>::value, int> = 0>
-constexpr void forEach(Range&& range, F&& f, PriorityTag<1>)
-{
-  auto size = HybridAlgorithm::size(range);
-  auto indices = std::make_index_sequence<size>();
-  forEachIndex(std::forward<Range>(range), std::forward<F>(f), indices);
-}
+  template<class Range, class F,
+    std::enable_if_t<IsIntegralConstant<decltype(HybridAlgorithm::size(std::declval<Range>()))>::value, int> = 0>
+  constexpr void forEach(Range&& range, F&& f, PriorityTag<1>)
+  {
+    auto size = HybridAlgorithm::size(range);
+    auto indices = std::make_index_sequence<size>();
+    forEachIndex(std::forward<Range>(range), std::forward<F>(f), indices);
+  }
 
-template<class Range, class F>
-constexpr void forEach(Range&& range, F&& f, PriorityTag<0>)
-{
-  for(std::size_t i=0; i<range.size(); ++i)
-    f(range[i]);
-//  for(auto e : range)
-//    f(e);
-}
+  template<class Range, class F>
+  constexpr void forEach(Range&& range, F&& f, PriorityTag<0>)
+  {
+    for(std::size_t i=0; i<range.size(); ++i)
+      f(range[i]);
+  // \ToDo Why does the following not compile?
+  //  for(auto e : range)
+  //    f(e);
+  }
 
 } // namespace Imp
 
 
 
+/**
+ * \brief Range based for loop
+ *
+ * \tparam Range Type of given range
+ * \tparam F Type of given predicate
+ *
+ * \param range The range to loop over
+ * \param f A predicate that will be called with each entry of the range
+ *
+ * This supports looping over the following ranges
+ * * ranges obtained from integralRange()
+ * * all ranges that provide HybridAlgorithm::size() and HybridAlgorithm::elementAt()
+ *
+ * This especially included instances of std::integer_sequence,
+ * std::tuple, Dune::TupleVector, and Dune::MultiTypeBlockVector.
+ */
 template<class Range, class F>
 constexpr void forEach(Range&& range, F&& f)
 {
   HybridAlgorithm::Imp::forEach(std::forward<Range>(range), std::forward<F>(f), PriorityTag<42>());
 }
 
+
+
 } // namespace HybridAlgorithm