Merge branch 'issue/allow_containers_in_discretefunction' into 'master'

Allow std containers as coefficients in DiscreteFunction and allow DefaultGlobalBasis

See merge request !248

amdis/ProblemInstat.hpp

@@ -62,12 +62,16 @@ namespace AMDiS
     }
 
     /// Return a const view to a oldSolution component
-    template <class... Indices>
+    /**
+     * \tparam Range  The range type return by evaluating the view in coordinates. If not specified,
+     *                it is automatically selected using \ref RangeType_t template.
+     **/
+    template <class Range = void, class... Indices>
     auto oldSolution(Indices... ii) const
     {
       test_exit_dbg(bool(oldSolution_),
         "OldSolution need to be created. Call initialize with INIT_UH_OLD.");
-      return valueOf(*oldSolution_, ii...);
+      return valueOf<Range>(*oldSolution_, ii...);
     }
 
     /// Implementation of \ref ProblemTimeInterface::transferInitialSolution().

amdis/ProblemStat.hpp

@@ -393,19 +393,27 @@ namespace AMDiS
 
 
     /// Return a mutable view to a solution component
-    template <class... Indices>
+    /**
+     * \tparam Range  The range type return by evaluating the view in coordinates. If not specified,
+     *                it is automatically selected using \ref RangeType_t template.
+     **/
+    template <class Range = void, class... Indices>
     auto solution(Indices... ii)
     {
       assert(bool(solution_) && "You have to call initialize() before.");
-      return valueOf(*solution_, ii...);
+      return valueOf<Range>(*solution_, ii...);
     }
 
     /// Return a const view to a solution component
-    template <class... Indices>
+    /**
+     * \tparam Range  The range type return by evaluating the view in coordinates. If not specified,
+     *                it is automatically selected using \ref RangeType_t template.
+     **/
+    template <class Range = void, class... Indices>
     auto solution(Indices... ii) const
     {
       assert(bool(solution_) && "You have to call initialize() before.");
-      return valueOf(*solution_, ii...);
+      return valueOf<Range>(*solution_, ii...);
     }
 
 

amdis/gridfunctions/DiscreteFunction.hpp

@@ -3,14 +3,16 @@
 #include <optional>
 #include <vector>
 
+#include <dune/common/ftraits.hh>
+#include <dune/common/typeutilities.hh>
 #include <dune/functions/common/defaultderivativetraits.hh>
-#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
 #include <dune/functions/functionspacebases/flatvectorview.hh>
 #include <dune/functions/gridfunctions/gridviewentityset.hh>
 #include <dune/typetree/childextraction.hh>
 
 #include <amdis/common/Tags.hpp>
 #include <amdis/functions/NodeCache.hpp>
+#include <amdis/linearalgebra/Access.hpp>
 #include <amdis/typetree/FiniteElementType.hpp>
 #include <amdis/typetree/RangeType.hpp>
 #include <amdis/typetree/TreePath.hpp>
@@ -22,28 +24,28 @@ namespace AMDiS
   /**
   * \ingroup GridFunctions
   *
-  * \tparam Coeff     Const or mutable coefficient type of the DOFVector
-  * \tparam GB        Thy type of the global basis
+  * \tparam Coeff     Const or mutable coefficient vector
+  * \tparam GB        The type of the global basis
   * \tparam TreePath  A realization of \ref Dune::TypeTree::HybridTreePath
+  * \tparam Range     The range type for th evaluation of the discrete function
   *
   * **Requirements:**
   * - GB models \ref Concepts::GlobalBasis
   **/
-  template <class Coeff, class GB, class TreePath>
+  template <class Coeff, class GB, class TreePath, class Range = void>
   class DiscreteFunction;
 
 
   /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
-  template <class Coeff, class GB, class TreePath>
+  template <class Coeff, class GB, class TreePath, class R>
   class DiscreteFunction
-      : public DiscreteFunction<Coeff const,GB,TreePath>
+      : public DiscreteFunction<Coeff const,GB,TreePath,R>
   {
-    using Self = DiscreteFunction<Coeff,GB,TreePath>;
-    using Super = DiscreteFunction<Coeff const,GB,TreePath>;
+    using Self = DiscreteFunction<Coeff,GB,TreePath,R>;
+    using Super = DiscreteFunction<Coeff const,GB,TreePath,R>;
 
     using Coefficients = Coeff;
     using GlobalBasis = GB;
-    using ValueType = typename Coeff::value_type;
 
   public:
     /// Constructor. Stores a pointer to the mutable `dofvector`.
@@ -54,16 +56,17 @@ namespace AMDiS
     {}
 
     template <class... Path>
-    DiscreteFunction(Coefficients& dofVector, std::shared_ptr<GlobalBasis const> const& basis, Path... path)
-      : DiscreteFunction(dofVector, *basis, path...)
+    DiscreteFunction(Coefficients& coefficients, std::shared_ptr<GlobalBasis const> const& basis, Path... path)
+      : DiscreteFunction(coefficients, *basis, path...)
     {}
 
     /// Construct a DiscreteFunction directly from a DOFVector
     template <class DV, class... Path,
+      Dune::disableCopyMove<DiscreteFunction,DV> = 0,
       class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
-      class GB_    = TYPEOF(*std::declval<DV>().basis())>
-    DiscreteFunction(DV&& dofVector, Path... path)
-      : DiscreteFunction(dofVector.coefficients(), *dofVector.basis(), path...)
+      class GB_    = TYPEOF(std::declval<DV>().basis())>
+    explicit DiscreteFunction(DV&& dofVector, Path... path)
+      : DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
     {}
 
   public:
@@ -125,11 +128,11 @@ namespace AMDiS
     /// Return the const DOFVector
     using Super::coefficients;
 
-    template <class... Indices>
+    template <class Range = void, class... Indices>
     auto child(Indices... ii)
     {
       auto tp = cat(this->treePath_, makeTreePath(ii...));
-      return DiscreteFunction<Coeff, GB, TYPEOF(makeTreePath(tp))>{*mutableCoeff_, this->basis(), tp};
+      return DiscreteFunction<Coeff, GB, TYPEOF(makeTreePath(tp)), Range>{*mutableCoeff_, this->basis(), tp};
     }
 
     using Super::child;
@@ -139,22 +142,20 @@ namespace AMDiS
   };
 
   /// A Const DiscreteFunction
-  template <class Coeff, class GB, class TreePath>
-  class DiscreteFunction<Coeff const,GB,TreePath>
+  template <class Coeff, class GB, class TreePath, class R>
+  class DiscreteFunction<Coeff const,GB,TreePath,R>
   {
   private:
     using Coefficients = std::remove_const_t<Coeff>;
     using GlobalBasis = GB;
     using LocalView = typename GlobalBasis::LocalView;
     using GridView = typename GlobalBasis::GridView;
-    using ValueType = typename Coefficients::value_type;
+
+    using Coefficient = TYPEOF(std::declval<Traits::Access>()(std::declval<Coeff const>(), std::declval<typename GB::MultiIndex>()));
 
     using Tree = typename GlobalBasis::LocalView::Tree;
     using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
 
-    using TreeCache = NodeCache_t<Tree>;
-    using SubTreeCache = typename Dune::TypeTree::ChildForTreePath<TreeCache, TreePath>;
-
   public:
     /// Set of entities the DiscreteFunction is defined on
     using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
@@ -162,14 +163,14 @@ namespace AMDiS
     /// Global coordinates of the EntitySet
     using Domain = typename EntitySet::GlobalCoordinate;
 
-    /// Range type of this DiscreteFunction
-    using Range = RangeType_t<SubTree,ValueType>;
+    /// Range type of this DiscreteFunction. If R=void deduce the Range automatically, using RangeType_t
+    using Range = std::conditional_t<std::is_same_v<R,void>, RangeType_t<SubTree,Coefficient>, R>;
 
     /// \brief This GridFunction has no derivative function, it can be created
     /// by \ref DiscreteGridFunction.
     enum { hasDerivative = false };
 
-  public:
+  private:
     /// A LocalFunction representing the localfunction and derivative of the DOFVector on a bound element
     template <class Type>
     class LocalFunction;
@@ -191,9 +192,10 @@ namespace AMDiS
 
     /// Construct a DiscreteFunction directly from a DOFVector
     template <class DV, class... Path,
+      Dune::disableCopyMove<DiscreteFunction,DV> = 0,
       class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
       class GB_    = TYPEOF(std::declval<DV>().basis())>
-    DiscreteFunction(DV const& dofVector, Path... path)
+    explicit DiscreteFunction(DV const& dofVector, Path... path)
       : DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
     {}
 
@@ -230,11 +232,11 @@ namespace AMDiS
       return *coefficients_;
     }
 
-    template <class... Indices>
+    template <class Range = void, class... Indices>
     auto child(Indices... ii) const
     {
       auto tp = cat(this->treePath_, makeTreePath(ii...));
-      return DiscreteFunction<Coeff const, GB, TYPEOF(makeTreePath(tp))>{*coefficients_, *basis_, tp};
+      return DiscreteFunction<Coeff const, GB, TYPEOF(makeTreePath(tp)), Range>{*coefficients_, *basis_, tp};
     }
 
   protected:
@@ -247,45 +249,47 @@ namespace AMDiS
 
   // deduction guides
 
-  template <class Coeff, class Basis, class... Path,
-    class TP = TYPEOF(makeTreePath(std::declval<Path>()...)),
+  template <class C, class Basis, class... Indices,
     class GB = Underlying_t<Basis>,
+    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...)),
     REQUIRES(Concepts::GlobalBasis<GB>)>
-  DiscreteFunction(Coeff&, Basis const&, Path...)
-    -> DiscreteFunction<Coeff,GB,TP>;
+  DiscreteFunction(C&, Basis const&, Indices...)
+    -> DiscreteFunction<C,GB,TP>;
 
-  template <class DV, class... Path,
-    class Coeff = decltype(std::declval<DV>().coefficients()),
-    class GB    = decltype(std::declval<DV>().basis()),
-    class TP    = TYPEOF(makeTreePath(std::declval<Path>()...))>
-  DiscreteFunction(DV&, Path...)
-    -> DiscreteFunction<std::remove_reference_t<Coeff>,Underlying_t<GB>,TP>;
+  template <class DV, class... Indices,
+    class C  = decltype(std::declval<DV>().coefficients()),
+    class GB = decltype(std::declval<DV>().basis()),
+    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
+  DiscreteFunction(DV&, Indices...)
+    -> DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP>;
 
 
   // grid functions representing the DOFVector
   // -----------------------------------------
 
   /// A Generator for the childs of a mutable \ref DiscreteFunction
-  template <class Coeff, class GB, class... Path, class... Indices>
-  auto valueOf(DiscreteFunction<Coeff,GB,Path...>& df, Indices... ii)
+  template <class Range = void, class C, class GB, class TP, class R, class... Indices>
+  auto valueOf(DiscreteFunction<C,GB,TP,R>& df, Indices... ii)
   {
-    return df.child(ii...);
+    return df.template child<Range>(ii...);
   }
 
   /// A Generator for the childs of a const \ref DiscreteFunction
-  template <class Coeff, class GB, class... Path, class... Indices>
-  auto valueOf(DiscreteFunction<Coeff,GB,Path...> const& df, Indices... ii)
+  template <class Range = void, class C, class GB, class TP, class R, class... Indices>
+  auto valueOf(DiscreteFunction<C,GB,TP,R> const& df, Indices... ii)
   {
-    return df.child(ii...);
+    return df.template child<Range>(ii...);
   }
 
   /// A Generator to transform a DOFVector into a \ref DiscreteFunction
-  template <class DV, class... Indices,
-    class = decltype(std::declval<DV>().coefficients()),
-    class = decltype(std::declval<DV>().basis())>
+  template <class Range = void, class DV, class... Indices,
+    class C  = decltype(std::declval<DV>().coefficients()),
+    class GB = decltype(std::declval<DV>().basis()),
+    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
   auto valueOf(DV& dofVec, Indices... ii)
   {
-    return DiscreteFunction{dofVec.coefficients(), dofVec.basis(), makeTreePath(ii...)};
+    using DF = DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP,Range>;
+    return DF{dofVec.coefficients(), dofVec.basis(), makeTreePath(ii...)};
   }
 
 } // end namespace AMDiS

amdis/gridfunctions/DiscreteFunction.inc.hpp

@@ -11,8 +11,8 @@
 namespace AMDiS {
 
 // Evaluate DiscreteFunction in global coordinates
-template <class Coeff, class GB, class TP>
-typename DiscreteFunction<Coeff const,GB,TP>::Range DiscreteFunction<Coeff const,GB,TP>::
+template <class C, class GB, class TP, class R>
+typename DiscreteFunction<C const,GB,TP,R>::Range DiscreteFunction<C const,GB,TP,R>::
   operator()(Domain const& x) const
 {
   using Grid = typename GlobalBasis::GridView::Grid;
@@ -30,9 +30,9 @@ typename DiscreteFunction<Coeff const,GB,TP>::Range DiscreteFunction<Coeff const
 
 
 // Interpolation of GridFunction to DOFVector
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Expr, class Tag>
-void DiscreteFunction<Coeff,GB,TP>::
+void DiscreteFunction<C,GB,TP,R>::
   interpolate_noalias(Expr&& expr, Tag strategy)
 {
   auto const& basis = this->basis();
@@ -41,7 +41,7 @@ void DiscreteFunction<Coeff,GB,TP>::
   auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
 
   if (std::is_same_v<Tag, tag::average>) {
-    VectorType_t<short,Coeff> counter(basis);
+    VectorType_t<short,Coefficients> counter(basis);
     AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
 
     coefficients().forEach([&counter](auto dof, auto& coeff)
@@ -55,13 +55,13 @@ void DiscreteFunction<Coeff,GB,TP>::
 
 
 // Interpolation of GridFunction to DOFVector
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Expr, class Tag>
-void DiscreteFunction<Coeff,GB,TP>::
+void DiscreteFunction<C,GB,TP,R>::
   interpolate(Expr&& expr, Tag strategy)
 {
   // create temporary copy of data
-  Coeff tmp(coefficients());
+  Coefficients tmp(coefficients());
 
   Self tmpView{tmp, this->basis(), this->treePath()};
   tmpView.interpolate_noalias(FWD(expr), strategy);

amdis/gridfunctions/DiscreteLocalFunction.inc.hpp

@@ -4,6 +4,8 @@
 #include <amdis/common/DerivativeTraits.hpp>
 #include <amdis/common/FieldMatVec.hpp>
 #include <amdis/common/Math.hpp>
+#include <amdis/common/RecursiveForEach.hpp>
+#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
 #include <amdis/functions/NodeIndices.hpp>
 #include <amdis/functions/Order.hpp>
 #include <amdis/typetree/FiniteElementType.hpp>
@@ -14,7 +16,7 @@
 namespace AMDiS {
 namespace Impl {
 
-// specialization of Coeff has gather method
+// specialization for containers with gather method
 template <class Coeff, class LocalView, class LocalCoeff,
   class = decltype(std::declval<Coeff>().gather(std::declval<LocalView>(), std::declval<LocalCoeff&>()))>
 void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<2>)
@@ -22,8 +24,9 @@ void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoe
   coeff.gather(localView, localCoeff);
 }
 
-// fallback implementation
-template <class Coeff, class LocalView, class LocalCoeff>
+// implementation for containers with multi-index access
+template <class Coeff, class LocalView, class LocalCoeff,
+  class = decltype(std::declval<Coeff>()[std::declval<typename LocalView::MultiIndex>()])>
 void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<1>)
 {
   localCoeff.resize(localView.size());
@@ -32,12 +35,24 @@ void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoe
     *it++ = coeff[idx];
 }
 
+// fallback implementation
+template <class Coeff, class LocalView, class LocalCoeff,
+  class = decltype(std::declval<Coeff>()[std::integral_constant<std::size_t,0>{}])>
+void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<0>)
+{
+  auto backend = Dune::Functions::istlVectorBackend(coeff);
+  localCoeff.resize(localView.size());
+  auto it = localCoeff.begin();
+  for (auto const& idx : nodeIndices(localView))
+    *it++ = backend[idx];
+}
+
 } // end namespace Impl
 
 
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Type>
-class DiscreteFunction<Coeff const,GB,TP>::LocalFunction
+class DiscreteFunction<C const,GB,TP,R>::LocalFunction
 {
   using DomainType = typename DiscreteFunction::Domain;
   using RangeType = typename DiscreteFunction::Range;
@@ -55,17 +70,16 @@ private:
   using LocalView = typename GlobalBasis::LocalView;
   using Element = typename EntitySet::Element;
   using Geometry = typename Element::Geometry;
-  using NodeToRangeMap = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
+  using NodeToRangeMap = HierarchicNodeToRangeMap;
 
 public:
   /// Constructor. Stores a copy of the DiscreteFunction.
   template <class LV>
-  LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, Coeff const& coefficients, Type type)
+  LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, C const& coefficients, Type type)
     : localView_(std::move(localView))
     , treePath_(treePath)
     , coefficients_(coefficients)
     , type_(type)
-    , nodeToRangeMap_(subTree())
   {}
 
   /// \brief Bind the LocalView to the element
@@ -103,7 +117,7 @@ public:
         return evaluateDivergence(local);
     }
 
-    return Range(0);
+    return Range{};
   }
 
   /// \brief Create a LocalFunction representing the derivative.
@@ -119,7 +133,7 @@ public:
     -> decltype(AMDiS::order(std::declval<SubTree>()))
   {
     assert(bound_);
-    return AMDiS::order(subTreeCache());
+    return AMDiS::order(subTree());
   }
 
   /// \brief Return the bound element
@@ -143,42 +157,62 @@ private:
   template <class LocalCoordinate>
   auto evaluatePartialDerivative (const LocalCoordinate& local) const;
 
+  // get the geometry of the bound element, that is constructed in the bind() method
   Geometry const& geometry() const
   {
     assert(bound_);
     return *geometry_;
   }
 
+  // return the sub-tree of the basis-tree associated to the tree-path
   SubTree const& subTree() const
   {
     return Dune::TypeTree::child(localView_->tree(), treePath_);
   }
 
-  SubTreeCache const& subTreeCache() const
+  // return the sub-tree of the basis-tree-cache associated to the tree-path
+  // NOTE: this is only return in case the local-view does provide a treeCache() function
+  template <class LV,
+    class = decltype(std::declval<LV>().treeCache())>
+  auto const& subTreeCache(Dune::PriorityTag<1>) const
   {
     return Dune::TypeTree::child(localView_->treeCache(), treePath_);
   }
 
+  // construct a new tree-cache that stores a reference to the sub-tree
+  template <class>
+  auto subTreeCache(Dune::PriorityTag<0>) const
+  {
+    return makeNodeCache(subTree());
+  }
+
+  decltype(auto) subTreeCache() const
+  {
+    return subTreeCache<LocalView>(Dune::PriorityTag<2>());
+  }
+
 private:
   std::shared_ptr<LocalView> localView_;
   TP treePath_;
-  Coeff const& coefficients_;
+  Coefficients const& coefficients_;
   Type type_;
   NodeToRangeMap nodeToRangeMap_;
 
   std::optional<Geometry> geometry_;
-  std::vector<ValueType> localCoefficients_;
+  std::vector<Coefficient> localCoefficients_;
   bool bound_ = false;
 };
 
 
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Type>
     template <class LocalCoordinate>
-auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
+auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
 ::evaluateFunction(LocalCoordinate const& local) const
 {
-  Range y(0);
+  Range y;
+  Recursive::forEach(y, [](auto& v) { v = 0; });
+
   Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
   {
     auto const& shapeFunctionValues = node.localBasisValuesAt(local);
@@ -209,13 +243,15 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
 }
 
 
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Type>
     template <class LocalCoordinate>
-auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
+auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
 ::evaluateJacobian(LocalCoordinate const& local) const
 {
-  Range dy(0);
+  Range dy;
+  Recursive::forEach(dy, [](auto& v) { v = 0; });
+
   Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
   {
     LocalToGlobalBasisAdapter localBasis(node, this->geometry());
@@ -246,15 +282,16 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
 }
 
 
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Type>
     template <class LocalCoordinate>
-auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
+auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
 ::evaluateDivergence(LocalCoordinate const& local) const
 {
   static_assert(static_size_v<Range> == 1, "Divergence implemented for scalar coefficients only.");
 
-  Range dy(0);
+  Range dy;
+  Recursive::forEach(dy, [](auto& v) { v = 0; });
 
   auto&& node = this->subTreeCache();
 
@@ -275,15 +312,17 @@ auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
 }
 
 
-template <class Coeff, class GB, class TP>
+template <class C, class GB, class TP, class R>
   template <class Type>
     template <class LocalCoordinate>
-auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
+auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
 ::evaluatePartialDerivative(LocalCoordinate const& local) const
 {
   std::size_t comp = this->type_.comp;
 
-  Range dy(0);
+  Range dy;
+  Recursive::forEach(dy, [](auto& v) { v = 0; });
+
   Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
   {
     LocalToGlobalBasisAdapter localBasis(node, this->geometry());

amdis/linearalgebra/Access.hpp

+#pragma once
+
+#include <dune/common/typeutilities.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+
+namespace AMDiS {
+namespace Traits {
+
+/// Unified accessor class
+struct Access
+{
+private:
+  // if C implements the vector facade interface
+  template <class C, class MI>
+  auto operator()(C const& c, MI const& mi, Dune::PriorityTag<2>) const
+    -> decltype(c.at(mi))
+  {
+    return c.at(mi);
+  }
+
+  // if C supports multi-index access
+  template <class C, class MI>
+  auto operator()(C const& c, MI const& mi, Dune::PriorityTag<1>) const
+    -> decltype(c[mi])
+  {
+    return c[mi];
+  }
+
+  // fallback implementation for simple vectors
+  template <class C, class MI>
+  auto operator()(C const& c, MI const& mi, Dune::PriorityTag<0>) const
+    -> decltype(Dune::Functions::istlVectorBackend(c)[mi])
+  {
+    return Dune::Functions::istlVectorBackend(c)[mi];
+  }
+
+public:
+  template <class C, class MI>
+  auto operator()(C const& c, MI const& mi) const