diff --git a/dune/solvers/common/algorithm.hh b/dune/solvers/common/algorithm.hh
index 9066ffc6042fec7f0571ffdc74f4b7f910ce110e..d40eab96cbb340fc7e432c7239a2518cc8db9a93 100644
--- a/dune/solvers/common/algorithm.hh
+++ b/dune/solvers/common/algorithm.hh
@@ -13,6 +13,281 @@
 
 namespace Dune {
 namespace Solvers {
+namespace HybridAlgorithm {
+
+namespace Imp {
+
+  // Try if tuple_size is implemented for class
+  template<class T, int i>
+  constexpr auto size(const Dune::FieldVector<T, i>*, const PriorityTag<5>&)
+    -> decltype(std::integral_constant<std::size_t,i>())
+  {
+    return {};
+  }
+
+  // Try if we have an instance of std::integer_sequence
+  template<class T, T... t, class Index>
+  constexpr auto size(std::integer_sequence<T, t...>, PriorityTag<4>)
+  {
+    using sizeAsType = std::tuple_size<decltype(std::make_tuple(t...))>;
+    return std::integral_constant<std::size_t, sizeAsType::value>();
+  }
+
+  // Try if tuple_size is implemented for class
+  template<class T>
+  constexpr auto size(const T*, const PriorityTag<3>&)
+    -> decltype(std::integral_constant<std::size_t,std::tuple_size<T>::value>())
+  {
+    return {};
+  }
+
+  // Try if there's a static constexpr size()
+  template<class T>
+  constexpr auto size(const T*, const PriorityTag<1>&)
+    -> decltype(std::integral_constant<std::size_t,T::size()>())
+  {
+    return {};
+  }
+
+  // As a last resort try if there's a static constexpr size()
+  template<class T>
+  constexpr auto size(const T* t, const PriorityTag<0>&)
+  {
+    return t->size();
+  }
+
+} // namespace Imp
+
+
+
+/**
+ * \brief Size query
+ *
+ * \tparam T Type of container whose size is queried
+ *
+ * \param t Container whose size is queried
+ *
+ * \return Size of t
+ *
+ * If the size of t is known at compile type the size is
+ * returned as std::integral_constant<std::size_t, size>.
+ * Otherwise the result of t.size() is returned.
+ *
+ * Supported types for deriving the size at compile time are:
+ * * instances of std::integer_sequence
+ * * all types std::tuple_size is implemented for
+ * * all typed that have a static method ::size()
+ * * instances of Dune::FieldVector
+ */
+template<class T>
+constexpr auto size(const T& t)
+{
+  return HybridAlgorithm::Imp::size(&t, PriorityTag<42>());
+}
+
+
+
+namespace Imp {
+
+  template<class Container, class Index,
+    std::enable_if_t<IsTuple<std::decay_t<Container>>::value, int> = 0>
+  constexpr auto elementAt(Container&& c, Index&&, PriorityTag<2>)
+  {
+    return std::get<Index::value>(c);
+  }
+
+  template<class T, T... t, class Index>
+  constexpr auto elementAt(std::integer_sequence<T, t...> c, Index&&, PriorityTag<1>)
+  {
+    return std::get<Index::value>(std::make_tuple(std::integral_constant<T, t>()...));
+  }
+
+  template<class Container, class Index>
+  constexpr auto elementAt(Container&& c, Index&& i, PriorityTag<0>)
+  {
+    return c[i];
+  }
+
+} // namespace Imp
+
+
+
+/**
+ * \brief Get element at given position from container
+ *
+ * \tparam Container Type of given container
+ * \tparam Index Type of index
+ *
+ * \param c Given container
+ * \param i Index of element to obtain
+ *
+ * \return The element at position i, i.e. c[i]
+ *
+ * If this returns the i-th entry of c. It supports the following
+ * containers
+ * * Containers providing dynamic access via operator[]
+ * * Heterogenous containers providing access via operator[](integral_constant<...>)
+ * * std::tuple<...>
+ * * std::integer_sequence
+ */
+template<class Container, class Index>
+constexpr auto elementAt(Container&& c, Index&& i)
+{
+  return HybridAlgorithm::Imp::elementAt(std::forward<Container>(c), std::forward<Index>(i), PriorityTag<42>());
+}
+
+
+
+namespace Imp {
+
+template<class Begin, class End>
+class StaticIntegralRange
+{
+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>();
+  }
+};
+
+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& end() const
+  { return end_; }
+
+  constexpr auto size() const
+  {
+    return end() - begin();
+  }
+
+  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>&)
+{
+  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
+
+
+
+/**
+ * \brief Create an integral range
+ *
+ * \tparam Begin Type of begin entry of the range
+ * \tparam End Type of end entry of the range
+ *
+ * \param begin First entry of the range
+ * \param end One past the last entry of the range
+ *
+ * \returns An object encoding the given range
+ *
+ * If Begin and End are both instances of type
+ * std::integral_constant, the returnes range
+ * encodes begin and end statically.
+ */
+template<class Begin, class End>
+constexpr auto integralRange(const Begin& begin, const End& end)
+{
+  return HybridAlgorithm::Imp::integralRange(begin, end, PriorityTag<42>());
+}
+
+/**
+ * \brief Create an integral range starting from 0
+ *
+ * \tparam End Type of end entry of the range
+ *
+ * \param end One past the last entry of the range
+ *
+ * \returns An object encoding the given range
+ *
+ * This is a short cut for integralRange(_0, end).
+ */
+template<class End>
+constexpr auto integralRange(const End& end)
+{
+  return HybridAlgorithm::Imp::integralRange(Dune::Indices::_0, end, PriorityTag<42>());
+}
+
+
+
+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,
+  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);
+}
+
+} // namespace Imp
+
+
+
+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
+
+
+
+
+
 
 
 // Implementation of integralRangeFor