diff --git a/dumux/io/vtkfunction.hh b/dumux/io/vtkfunction.hh
index 6da7340f7d0cf273c6e75e3ed4be27f1179f8961..20ca6bd5f38e7c7e12703bcfa396778e24ae6151 100644
--- a/dumux/io/vtkfunction.hh
+++ b/dumux/io/vtkfunction.hh
@@ -36,6 +36,24 @@
 
 namespace Dumux::Vtk {
 
+namespace Impl {
+
+template<class Field>
+double accessEntry(const Field& f, int mycomp, int i)
+{
+    if constexpr (IsIndexable<std::decay_t<decltype(std::declval<Field>()[0])>>{})
+    {
+        if constexpr (IsIndexable<std::decay_t<decltype(std::declval<Field>()[0][0])>>{})
+            DUNE_THROW(Dune::InvalidStateException, "Invalid field type");
+        else
+            return f[i][mycomp];
+    }
+    else
+        return f[i];
+}
+
+} // end namespace Impl
+
 /*!
  * \ingroup InputOutput
  * \brief a VTK function that supports both scalar and vector values for each element
@@ -54,18 +72,18 @@ struct VectorP0VTKFunction : Dune::VTKFunction<GridView>
 public:
 
     //! return number of components
-    virtual int ncomps() const { return nComps_; }
+    int ncomps() const final { return nComps_; }
 
     //! get name
-    virtual std::string name() const { return name_; }
+    std::string name() const final { return name_; }
 
     //! evaluate
-    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>&) const
-    { return accessChooser_(mycomp, mapper_.index(e), IsIndexable<decltype(field_[0])>()); }
+    double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>&) const final
+    { return Impl::accessEntry(field_, mycomp, mapper_.index(e)); }
 
 #if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
     //! get output precision for the field
-    Dumux::Vtk::Precision precision() const override
+    Dumux::Vtk::Precision precision() const final
     { return precision_; }
 #endif
 
@@ -84,20 +102,6 @@ public:
     }
 
 private:
-
-    //! access for vectorial fields
-    double accessChooser_(int mycomp, int i, std::true_type) const
-    { return vectorFieldAccess_(mycomp, i, IsIndexable<decltype(field_[0][0])>()); }
-
-    //! access for scalar fields
-    double accessChooser_(int mycomp, int i, std::false_type) const { return field_[i]; }
-
-    //! access to permissive vectorial fields
-    double vectorFieldAccess_(int mycomp, int i, std::false_type) const { return field_[i][mycomp]; }
-
-    //! if the field is indexable more than two times, throw error
-    double vectorFieldAccess_(int mycomp, int i, std::true_type) const { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
-
     const F& field_;
     const std::string name_;
     int nComps_;
@@ -123,20 +127,20 @@ struct VectorP1VTKFunction : Dune::VTKFunction<GridView>
 public:
 
     //! return number of components
-    virtual int ncomps() const { return nComps_; }
+    int ncomps() const final { return nComps_; }
 
     //! get name
-    virtual std::string name() const { return name_; }
+    std::string name() const final { return name_; }
 
     //! evaluate
-    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const
+    double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const final
     {
         const unsigned int dim = Element::mydimension;
         const unsigned int nVertices = e.subEntities(dim);
 
         std::vector<Dune::FieldVector<ctype, 1>> cornerValues(nVertices);
         for (unsigned i = 0; i < nVertices; ++i)
-            cornerValues[i] = accessChooser_(mycomp, mapper_.subIndex(e, i, dim), IsIndexable<decltype(field_[0])>());
+            cornerValues[i] = Impl::accessEntry(field_, mycomp, mapper_.subIndex(e, i, dim));
 
         // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
         const Dune::MultiLinearGeometry<ctype, dim, 1> interpolation(e.type(), std::move(cornerValues));
@@ -145,7 +149,7 @@ public:
 
 #if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
     //! get output precision for the field
-    Dumux::Vtk::Precision precision() const override
+    Dumux::Vtk::Precision precision() const final
     { return precision_; }
 #endif
 
@@ -164,20 +168,6 @@ public:
                                       << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")");
     }
 private:
-
-    //! access for vectorial fields
-    double accessChooser_(int mycomp, int i, std::true_type) const
-    { return vectorFieldAccess_(mycomp, i, IsIndexable<decltype(field_[0][0])>()); }
-
-    //! access for scalar fields
-    double accessChooser_(int mycomp, int i, std::false_type) const { return field_[i]; }
-
-    //! access to permissive vectorial fields
-    double vectorFieldAccess_(int mycomp, int i, std::false_type) const { return field_[i][mycomp]; }
-
-    //! if the field is indexable more than two times, throw error
-    double vectorFieldAccess_(int mycomp, int i, std::true_type) const { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
-
     const F& field_;
     const std::string name_;
     int nComps_;
@@ -207,20 +197,20 @@ struct VectorP1NonConformingVTKFunction : Dune::VTKFunction<GridView>
 public:
 
     //! return number of components
-    virtual int ncomps() const { return nComps_; }
+    int ncomps() const final { return nComps_; }
 
     //! get name
-    virtual std::string name() const { return name_; }
+    std::string name() const final { return name_; }
 
     //! evaluate
-    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const
+    double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const final
     {
         const unsigned int dim = Element::mydimension;
         const unsigned int nVertices = e.subEntities(dim);
 
         std::vector<Dune::FieldVector<ctype, 1>> cornerValues(nVertices);
         for (unsigned i = 0; i < nVertices; ++i)
-            cornerValues[i] = accessChooser_(mycomp, mapper_.index(e), i, IsIndexable<decltype(field_[0])>());
+            cornerValues[i] = accessEntry_(mycomp, mapper_.index(e), i);
 
         // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
         const Dune::MultiLinearGeometry<ctype, dim, 1> interpolation(e.type(), std::move(cornerValues));
@@ -229,7 +219,7 @@ public:
 
 #if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
     //! get output precision for the field
-    Dumux::Vtk::Precision precision() const override
+    Dumux::Vtk::Precision precision() const final
     { return precision_; }
 #endif
 
@@ -247,22 +237,19 @@ public:
                                       << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")");
     }
 private:
-
     //! access to the field
-    double accessChooser_(int mycomp, int eIdx, int cornerIdx, std::true_type) const
-    { return fieldAccess_(mycomp, eIdx, cornerIdx, IsIndexable<decltype(field_[0][0])>()); }
-
-    //! fields have to be indexable at least twice
-    double accessChooser_(int mycomp, int eIdx, int cornerIdx, std::false_type) const
-    { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
-
-    //! scalar field access
-    double fieldAccess_(int mycomp, int eIdx, int cornerIdx, std::false_type) const
-    { return field_[eIdx][cornerIdx]; }
-
-    //! vectorial field access
-    double fieldAccess_(int mycomp, int eIdx, int cornerIdx, std::true_type) const
-    { return field_[eIdx][cornerIdx][mycomp]; }
+    double accessEntry_(int mycomp, int eIdx, int cornerIdx) const
+    {
+        if constexpr (IsIndexable<decltype(std::declval<F>()[0])>{})
+        {
+            if constexpr (IsIndexable<decltype(std::declval<F>()[0][0])>{})
+                return field_[eIdx][cornerIdx][mycomp];
+            else
+                return field_[eIdx][cornerIdx];
+        }
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Invalid field type");
+    }
 
     const F& field_;
     const std::string name_;
@@ -307,16 +294,14 @@ public:
             DUNE_THROW(Dune::NotImplemented, "Only element or vertex quantities allowed.");
     }
 
-    virtual ~Field() {}
-
     //! return the name of this field
-    virtual std::string name () const { return field_->name(); }
+    std::string name () const { return field_->name(); }
 
     //! return the number of components of this field
-    virtual int ncomps() const { return field_->ncomps(); }
+    int ncomps() const { return field_->ncomps(); }
 
     //! return the precision of this field
-    virtual Dumux::Vtk::Precision precision() const
+    Dumux::Vtk::Precision precision() const
     {
 #if DUNE_VERSION_LT(DUNE_GRID, 2, 7)
         return Dumux::Vtk::Precision::float32;
@@ -329,9 +314,9 @@ public:
     int codim() const { return codim_; }
 
     //! element-local evaluation of the field
-    virtual double evaluate(int mycomp,
-                            const Element &element,
-                            const Dune::FieldVector< ctype, dim > &xi) const
+    double evaluate(int mycomp,
+                    const Element &element,
+                    const Dune::FieldVector< ctype, dim > &xi) const
     { return field_->evaluate(mycomp, element, xi); }
 
     //! returns the underlying vtk function
diff --git a/dumux/io/vtkmultiwriter.hh b/dumux/io/vtkmultiwriter.hh
index fe87d1bfd22a8af28b7e38b3e6544eea36eebca5..573326cdb602b7a897aae6dc7fb9f7dc8f64f0bb 100644
--- a/dumux/io/vtkmultiwriter.hh
+++ b/dumux/io/vtkmultiwriter.hh
@@ -32,12 +32,11 @@
 #include <memory>
 #include <string>
 
-#include "vtknestedfunction.hh"
-
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
 
 #include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/io/file/vtk/function.hh>
 
 #include <dumux/common/valgrind.hh>
 
@@ -47,6 +46,98 @@
 #endif
 
 namespace Dumux {
+
+/*!
+ * \ingroup InputOutput
+ * \brief Provides a vector-valued function using Dune::FieldVectors
+ *        as elements.
+ * DEPRECATED will be removed once this header is removed
+ */
+template <class GridView, class Mapper, class Buffer>
+class VtkNestedFunction : public Dune::VTKFunction<GridView>
+{
+    enum { dim = GridView::dimension };
+    using ctype = typename GridView::ctype;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using ReferenceElements = Dune::ReferenceElements<ctype, dim>;
+
+public:
+    VtkNestedFunction(std::string name,
+                      const GridView &gridView,
+                      const Mapper &mapper,
+                      const Buffer &buf,
+                      int codim,
+                      int numComp)
+        : name_(name)
+        , gridView_(gridView)
+        , mapper_(mapper)
+        , buf_(buf)
+        , codim_(codim)
+        , numComp_(numComp)
+    {
+        assert(std::size_t(buf_.size()) == std::size_t(mapper_.size()));
+    }
+
+    virtual std::string name () const
+    { return name_; }
+
+    virtual int ncomps() const
+    { return numComp_; }
+
+    virtual double evaluate(int mycomp,
+                            const Element &element,
+                            const Dune::FieldVector< ctype, dim > &xi) const
+    {
+        int idx;
+        if (codim_ == 0) {
+            // cells. map element to the index
+            idx = mapper_.index(element);
+        }
+        else if (codim_ == dim) {
+            // find vertex which is closest to xi in local
+            // coordinates. This code is based on Dune::P1VTKFunction
+            double min=1e100;
+            int imin=-1;
+            Dune::GeometryType geomType = element.type();
+            int n = element.subEntities(dim);
+
+            for (int i=0; i < n; ++i)
+            {
+                Dune::FieldVector<ctype,dim> local =
+                    ReferenceElements::general(geomType).position(i,dim);
+                local -= xi;
+                if (local.infinity_norm()<min)
+                {
+                    min = local.infinity_norm();
+                    imin = i;
+                }
+            }
+
+            // map vertex to an index
+            idx = mapper_.subIndex(element, imin, codim_);
+        }
+        else
+            DUNE_THROW(Dune::InvalidStateException,
+                       "Only element and vertex based vector "
+                       " fields are supported so far.");
+
+        double val = buf_[idx][mycomp];
+        using std::abs;
+        if (abs(val) < std::numeric_limits<float>::min())
+            val = 0;
+
+        return val;
+    }
+
+private:
+    const std::string name_;
+    const GridView gridView_;
+    const Mapper &mapper_;
+    const Buffer &buf_;
+    int codim_;
+    int numComp_;
+};
+
 /*!
  * \ingroup InputOutput
  * \brief Simplifies writing multi-file VTK datasets.
diff --git a/dumux/io/vtknestedfunction.hh b/dumux/io/vtknestedfunction.hh
index 1fa34c2f8226e7a1bc4c0aab84b741ad44c4365e..cb35bec5780c1c2247712631080286846a2bbb3b 100644
--- a/dumux/io/vtknestedfunction.hh
+++ b/dumux/io/vtknestedfunction.hh
@@ -25,104 +25,9 @@
 #ifndef VTK_NESTED_FUNCTION_HH
 #define VTK_NESTED_FUNCTION_HH
 
-#include <string>
-
-#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
-#include <dune/grid/io/file/vtk/function.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup InputOutput
- * \brief Provides a vector-valued function using Dune::FieldVectors
- *        as elements.
- */
-template <class GridView, class Mapper, class Buffer>
-class VtkNestedFunction : public Dune::VTKFunction<GridView>
-{
-    enum { dim = GridView::dimension };
-    using ctype = typename GridView::ctype;
-    using Element = typename GridView::template Codim<0>::Entity;
-    using ReferenceElements = Dune::ReferenceElements<ctype, dim>;
-
-public:
-    VtkNestedFunction(std::string name,
-                      const GridView &gridView,
-                      const Mapper &mapper,
-                      const Buffer &buf,
-                      int codim,
-                      int numComp)
-        : name_(name)
-        , gridView_(gridView)
-        , mapper_(mapper)
-        , buf_(buf)
-        , codim_(codim)
-        , numComp_(numComp)
-    {
-        assert(std::size_t(buf_.size()) == std::size_t(mapper_.size()));
-    }
-
-    virtual std::string name () const
-    { return name_; }
-
-    virtual int ncomps() const
-    { return numComp_; }
-
-    virtual double evaluate(int mycomp,
-                            const Element &element,
-                            const Dune::FieldVector< ctype, dim > &xi) const
-    {
-        int idx;
-        if (codim_ == 0) {
-            // cells. map element to the index
-            idx = mapper_.index(element);
-        }
-        else if (codim_ == dim) {
-            // find vertex which is closest to xi in local
-            // coordinates. This code is based on Dune::P1VTKFunction
-            double min=1e100;
-            int imin=-1;
-            Dune::GeometryType geomType = element.type();
-            int n = element.subEntities(dim);
-
-            for (int i=0; i < n; ++i)
-            {
-                Dune::FieldVector<ctype,dim> local =
-                    ReferenceElements::general(geomType).position(i,dim);
-                local -= xi;
-                if (local.infinity_norm()<min)
-                {
-                    min = local.infinity_norm();
-                    imin = i;
-                }
-            }
-
-            // map vertex to an index
-            idx = mapper_.subIndex(element, imin, codim_);
-        }
-        else
-            DUNE_THROW(Dune::InvalidStateException,
-                       "Only element and vertex based vector "
-                       " fields are supported so far.");
-
-        double val = buf_[idx][mycomp];
-        using std::abs;
-        if (abs(val) < std::numeric_limits<float>::min())
-            val = 0;
-
-        return val;
-    }
-
-private:
-    const std::string name_;
-    const GridView gridView_;
-    const Mapper &mapper_;
-    const Buffer &buf_;
-    int codim_;
-    int numComp_;
-};
-
-} // end namespace Dumux
+// This header is deprecated. The class implementation has been moved
+// to the deprecated header dumux/io/vtkmultiwriter.hh to avoid double deprecation warnings.
+// Remove this header once dumux/io/vtkmultiwriter.hh is removed.
+#include <dumux/io/vtkmultiwriter.hh>
 
 #endif