diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 87b6d9706..2654a5863 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -621,8 +621,10 @@ if (Damaris_FOUND AND MPI_FOUND)
   list (APPEND PUBLIC_HEADER_FILES ebos/damaris_properties.hh)
   list (APPEND PUBLIC_HEADER_FILES ebos/damariswriter.hh)
   list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/DamarisVar.hpp)
-  list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/GridDataOutput.hpp)
-  list(APPEND MAIN_SOURCE_FILES opm/simulators/utils/DamarisVar.cpp)
+  list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/GridDataOutput.hpp
+                                   opm/simulators/utils/GridDataOutput_impl.hpp)
+  list(APPEND MAIN_SOURCE_FILES opm/simulators/utils/DamarisVar.cpp
+                                opm/simulators/utils/GridDataOutput.cpp)
 endif()
 
 if(HDF5_FOUND)
diff --git a/flow/flow_blackoil_alugrid.cpp b/flow/flow_blackoil_alugrid.cpp
index f7a46a075..4514670b8 100644
--- a/flow/flow_blackoil_alugrid.cpp
+++ b/flow/flow_blackoil_alugrid.cpp
@@ -35,6 +35,7 @@
 #include <ebos/eclgenericwriter_impl.hh>
 #include <ebos/ecltransmissibility_impl.hh>
 #include <ebos/equil/initstateequil_impl.hh>
+#include <opm/simulators/utils/GridDataOutput_impl.hpp>
 
 namespace Opm {
 namespace Properties {
diff --git a/flow/flow_blackoil_polyhedralgrid.cpp b/flow/flow_blackoil_polyhedralgrid.cpp
index ca1b5dd59..9eb77060f 100644
--- a/flow/flow_blackoil_polyhedralgrid.cpp
+++ b/flow/flow_blackoil_polyhedralgrid.cpp
@@ -33,6 +33,7 @@
 #include <ebos/ecltransmissibility_impl.hh>
 #include <ebos/eclgenericwriter_impl.hh>
 #include <ebos/equil/initstateequil_impl.hh>
+#include <opm/simulators/utils/GridDataOutput_impl.hpp>
 
 namespace Opm {
 namespace Properties {
diff --git a/opm/simulators/utils/GridDataOutput.cpp b/opm/simulators/utils/GridDataOutput.cpp
new file mode 100644
index 000000000..0b66606c3
--- /dev/null
+++ b/opm/simulators/utils/GridDataOutput.cpp
@@ -0,0 +1,54 @@
+/*
+  Copyright 2023 Inria, Bretagne–Atlantique Research Center
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include <config.h>
+#include <opm/simulators/utils/GridDataOutput_impl.hpp>
+
+#include <opm/grid/CpGrid.hpp>
+
+#include <opm/simulators/utils/DamarisVar.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/gridpart/common/gridpart2gridview.hh>
+#include <ebos/femcpgridcompat.hh>
+#endif // HAVE_DUNE_FEM
+
+namespace Opm::GridDataOutput {
+
+template<class T> using DV = DamarisOutput::DamarisVar<T>;
+
+#define INSTANCE(part, ...) \
+    template class SimMeshDataAccessor<__VA_ARGS__, part>; \
+    template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
+        writeGridPoints<DV<double>>(DV<double>&, DV<double>&, DV<double>&) const; \
+    template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
+        writeConnectivity<DV<int>>(DV<int>&, ConnectivityVertexOrder) const; \
+    template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
+        writeOffsetsCells<DV<int>>(DV<int>&) const; \
+    template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
+        writeCellTypes<DV<char>>(DV<char>&) const;
+
+INSTANCE(1, Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>)
+
+#if HAVE_DUNE_FEM
+INSTANCE(1, Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, (Dune::PartitionIteratorType)4, false> >)
+INSTANCE(1, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>)
+#endif
+
+} // namespace Opm::GridDataOutput
diff --git a/opm/simulators/utils/GridDataOutput.hpp b/opm/simulators/utils/GridDataOutput.hpp
index a5410e3a2..c24c7f758 100644
--- a/opm/simulators/utils/GridDataOutput.hpp
+++ b/opm/simulators/utils/GridDataOutput.hpp
@@ -22,11 +22,10 @@
 
 #include <opm/common/ErrorMacros.hpp>
 
-#include <dune/grid/common/rangegenerators.hh>
-#include <dune/grid/io/file/vtk/common.hh>
+#include <dune/grid/common/partitionset.hh>
 
-#include <cstddef>
-#include <sstream>
+#include <iosfwd>
+#include <string>
 
 /** @file
     @brief Allows model geometry data to be passed to external code - via a copy
@@ -129,14 +128,8 @@ public:
      * polyhedralCellPresent() to test if polyhedral cells are present and decide
      * what they want to do before copying data using the data accessor methods.
      */
-    explicit SimMeshDataAccessor(const GridView& gridView, Dune::PartitionSet<partitions> dunePartition)
-        : gridView_(gridView)
-        , dunePartition_(dunePartition)
-    {
-        dimw_ = GridView::dimension; // this is an enum
-        partition_value_ = dunePartition.value;
-        countEntities();
-    }
+    explicit SimMeshDataAccessor(const GridView& gridView,
+                                 Dune::PartitionSet<partitions> dunePartition);
 
     /**
       Checks for cells that have polyhedral type within the current partition of
@@ -146,16 +139,7 @@ public:
       partition is not going to be available for visualisation as this class does
       not yet handle polyhedral cells.
    */
-    bool polyhedralCellPresent()
-    {
-        for (const auto& cit : elements(gridView_, dunePartition_)) {
-            auto corner_geom = cit.geometry();
-            if (Dune::VTK::geometryType(corner_geom.type()) == Dune::VTK::polyhedron) {
-                return true;
-            }
-        }
-        return false;
-    }
+    bool polyhedralCellPresent() const;
 
     /**
         Count the vertices, cells and corners.
@@ -164,22 +148,7 @@ public:
         do not need to renumber the vertices as all the subsets use references to
         the full set.
     */
-    void countEntities()
-    {
-        // We include all the vertices for this ranks partition
-        const auto& vert_partition = vertices(gridView_, Dune::Partitions::all);
-        nvertices_ = std::distance(vert_partition.begin(), vert_partition.end());
-
-        const auto& cell_partition = elements(gridView_, dunePartition_);
-        ncells_ = 0;
-        ncorners_ = 0;
-        for (const auto& cit : cell_partition) {
-            auto corner_geom = cit.geometry();
-            ncorners_ += corner_geom.corners();
-            ++ncells_;
-        }
-    }
-
+    void countEntities();
 
     /**
       Write the positions of vertices - directly to the pointers given in
@@ -194,38 +163,7 @@ public:
       Returns the number of vertices written
    */
     template <typename T>
-    long writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size = 0)
-    {
-        if (max_size < nvertices_) {
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints( T&  x_inout,  T&  "
-                      "y_inout, T& z_inout )  "
-                          + " Input objects max_size (" + std::to_string(max_size)
-                          + ") is not sufficient to fit the nvertices_ values (" + std::to_string(nvertices_) + ")");
-        }
-
-        long i = 0;
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
-                x_inout[i] = static_cast<T>(xyz_local[0]);
-                y_inout[i] = static_cast<T>(xyz_local[1]);
-                z_inout[i] = static_cast<T>(xyz_local[2]);
-                i++;
-            }
-        } else if (dimw_ == 2) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
-                x_inout[i] = static_cast<T>(xyz_local[0]);
-                y_inout[i] = static_cast<T>(xyz_local[1]);
-                z_inout[i] = static_cast<T>(0.0);
-                i++;
-            }
-        }
-        assert(i == nvertices_); // As we are templated on the
-                                 // Dune::PartitionSet<partitions>, this cannot change
-        return i;
-    }
+    long writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size = 0) const;
 
     /**
       Write the positions of vertices - directly to the pointers given in
@@ -241,49 +179,7 @@ public:
       Returns the number of vertices written
    */
     template <typename VectType>
-    long writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout)
-    {
-        const std::size_t check_size_x = x_inout.size();
-        const std::size_t check_size_y = y_inout.size();
-        const std::size_t check_size_z = z_inout.size();
-
-        using VT = decltype(x_inout.data()[0]);
-
-        if ((check_size_x < static_cast<std::size_t>(nvertices_)) ||
-            (check_size_y < static_cast<std::size_t>(nvertices_)) ||
-            (check_size_z < static_cast<std::size_t>(nvertices_))) {
-            // assert(check_size >= nvertices_);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints( VectType&  x_inout,  VectType&  "
-                      "y_inout, VectType& z_inout )  At least one of the inputs"
-                          + "  object x size " + std::to_string(check_size_x) + "  object y size "
-                          + std::to_string(check_size_y) + "  object z size " + std::to_string(check_size_z)
-                          + " is not sufficient to fit the nvertices_ values( " + std::to_string(nvertices_) + " )");
-        }
-
-        long i = 0;
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
-                x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
-                y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
-                z_inout.data()[i] = static_cast<VT>(xyz_local[2]);
-                i++;
-            }
-        } else if (dimw_ == 2) {
-            double td = 0.0;
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
-                x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
-                y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
-                z_inout.data()[i] = static_cast<VT>(td);
-                i++;
-            }
-        }
-        assert(i == nvertices_); // As we are templated on the
-                                 // Dune::PartitionSet<partitions>, this cannot change
-        return i;
-    }
+    long writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout) const;
 
     /**
       Write the positions of vertices - directly to the pointers given in
@@ -296,35 +192,7 @@ public:
       Returns the number of vertices written
    */
     template <typename T>
-    long writeGridPoints_AOS(T* xyz_inout, long max_size = 0)
-    {
-        if (max_size < nvertices_ * 3) {
-            assert(max_size >= nvertices_ * 3);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints_AOS( T*  xyz_inout )  " + " Input objects max_size ("
-                          + std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
-                          + std::to_string(nvertices_ * 3) + ")");
-        }
-
-        long i = 0;
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout[i++] = static_cast<T>(xyz_local[0]);
-                xyz_inout[i++] = static_cast<T>(xyz_local[1]);
-                xyz_inout[i++] = static_cast<T>(xyz_local[2]);
-            }
-        } else if (dimw_ == 2) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout[i++] = static_cast<T>(xyz_local[0]);
-                xyz_inout[i++] = static_cast<T>(xyz_local[1]);
-                xyz_inout[i++] = static_cast<T>(0.0);
-            }
-        }
-        return ((i) / 3);
-    }
-
+    long writeGridPoints_AOS(T* xyz_inout, long max_size = 0) const;
 
     /**
       Write the positions of vertices - directly to the pointers given in
@@ -336,40 +204,7 @@ public:
       Returns the number of vertices written
    */
     template <typename VectType>
-    long writeGridPoints_AOS(VectType& xyz_inout)
-    {
-        const std::size_t check_size = xyz_inout.size();
-
-        using VT = decltype(xyz_inout.data()[0]);
-
-        if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
-            assert(check_size >= nvertices_ * 3);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints_AOS( VectType&  xyz_inout )  "
-                          + " Input objects check_size (" + std::to_string(check_size)
-                          + ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
-                          + ")");
-        }
-
-        long i = 0;
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout.data()[i++] = static_cast<VT>(xyz_local[0]);
-                xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
-                xyz_inout[i++] = static_cast<VT>(xyz_local[2]);
-            }
-        } else if (dimw_ == 2) {
-            double td = 0.0;
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout[i++] = static_cast<VT>(xyz_local[0]);
-                xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
-                xyz_inout[i++] = static_cast<VT>(td);
-            }
-        }
-        return ((i) / 3);
-    }
+    long writeGridPoints_AOS(VectType& xyz_inout) const;
 
     /**
         Write the positions of vertices - directly to the pointers given in
@@ -382,40 +217,7 @@ public:
         Returns the number of vertices written
      */
     template <typename T>
-    long writeGridPoints_SOA(T* xyz_inout, long max_size = 0)
-    {
-        if (max_size < nvertices_ * 3) {
-            // assert(max_size >= nvertices_ * 3);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints_SOA( T&  xyz_inout )  " + " Input objects max_size ("
-                          + std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
-                          + std::to_string(nvertices_ * 3) + ")");
-        }
-
-        long i = 0;
-        // Get offsets into structure
-        T* xyz_inout_y = xyz_inout + nvertices_;
-        T* xyz_inout_z = xyz_inout + (2 * nvertices_);
-
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout[i] = static_cast<T>(xyz_local[0]);
-                xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
-                xyz_inout_z[i] = static_cast<T>(xyz_local[2]);
-                i++;
-            }
-        } else if (dimw_ == 2) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout[i] = static_cast<T>(xyz_local[0]);
-                xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
-                xyz_inout_z[i] = static_cast<T>(0.0);
-                i++;
-            }
-        }
-        return (i);
-    }
+    long writeGridPoints_SOA(T* xyz_inout, long max_size = 0) const;
 
     /**
       Write the positions of vertices - directly to the pointers given in
@@ -427,46 +229,7 @@ public:
       Returns the number of vertices written
    */
     template <typename VectType>
-    long writeGridPoints_SOA(VectType& xyz_inout)
-    {
-        const std::size_t check_size = xyz_inout.size();
-
-        if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
-            // assert(check_size >= nvertices_ * 3);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeGridPoints_SOA( VectType&  xyz_inout )  "
-                          + " Input objects check_size (" + std::to_string(check_size)
-                          + ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
-                          + ")");
-        }
-
-        using VT = decltype(xyz_inout.data()[0]);
-
-        long i = 0;
-        // Get offsets into structure
-        VT* xyz_inout_y = xyz_inout.data() + nvertices_;
-        VT* xyz_inout_z = xyz_inout.data() + (2 * nvertices_);
-
-        if (dimw_ == 3) {
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
-                xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
-                xyz_inout_z[i] = static_cast<VT>(xyz_local[2]);
-                i++;
-            }
-        } else if (dimw_ == 2) {
-            double td = 0.0;
-            for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
-                auto xyz_local = vit.geometry().corner(0);
-                xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
-                xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
-                xyz_inout_z[i] = static_cast<VT>(td);
-                i++;
-            }
-        }
-        return (i);
-    }
+    long writeGridPoints_SOA(VectType& xyz_inout) const;
 
     /**
     * Write the connectivity array - directly to the pointer given in parameter 1
@@ -480,42 +243,8 @@ public:
       Returns the number of corner indices written.
     */
     template <typename Integer>
-    long writeConnectivity(Integer* connectivity_inout, ConnectivityVertexOrder whichOrder, long max_size = 0)
-    {
-        if (max_size < ncorners_) {
-
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeConnectivity( T*  connectivity_inout,... )  "
-                          + " Input max_size value (" + std::to_string(max_size)
-                          + ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
-        }
-
-        long i = 0;
-        if (whichOrder == DUNE) {
-            // DUNE order
-            for (const auto& cit : elements(gridView_, dunePartition_)) {
-                auto cell_corners = cit.geometry().corners();
-                for (auto vx = 0; vx < cell_corners; ++vx) {
-                    const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
-                    connectivity_inout[i + vx] = vxIdx;
-                }
-                i += cell_corners;
-            }
-        } else {
-            // VTK order
-            for (const auto& cit : elements(gridView_, dunePartition_)) {
-                auto cell_corners = cit.geometry().corners();
-                for (auto vx = 0; vx < cell_corners; ++vx) {
-                    const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
-                    int vtkOrder;
-                    vtkOrder = Dune::VTK::renumber(cit.type(), vx);
-                    connectivity_inout[i + vtkOrder] = vxIdx;
-                }
-                i += cell_corners;
-            }
-        }
-        return (i);
-    }
+    long writeConnectivity(Integer* connectivity_inout,
+                           ConnectivityVertexOrder whichOrder, long max_size = 0) const;
 
     /**
     * Write the connectivity array - directly to a VectType object given in parameter 1
@@ -530,46 +259,8 @@ public:
       Returns the number of corner indices written.
     */
     template <typename VectType>
-    long writeConnectivity(VectType& connectivity_inout, ConnectivityVertexOrder whichOrder)
-    {
-        const std::size_t check_size = connectivity_inout.size();
-
-        if (check_size < static_cast<std::size_t>(ncorners_)) {
-            // assert(check_size >= ncorners_);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeConnectivity( VectType&  "
-                      "connectivity_inout )  "
-                          + " Input objects size (" + std::to_string(check_size)
-                          + ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
-        }
-
-        long i = 0;
-        if (whichOrder == DUNE) {
-            // DUNE order
-            for (const auto& cit : elements(gridView_, dunePartition_)) {
-                auto cell_corners = cit.geometry().corners();
-                for (auto vx = 0; vx < cell_corners; ++vx) {
-                    const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
-                    connectivity_inout.data()[i + vx] = vxIdx;
-                }
-                i += cell_corners;
-            }
-        } else {
-            // VTK order
-            for (const auto& cit : elements(gridView_, dunePartition_)) {
-                auto cell_corners = cit.geometry().corners();
-                for (auto vx = 0; vx < cell_corners; ++vx) {
-                    const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
-                    int vtkOrder;
-                    vtkOrder = Dune::VTK::renumber(cit.type(), vx);
-                    connectivity_inout.data()[i + vtkOrder] = vxIdx;
-                }
-                i += cell_corners;
-            }
-        }
-        return (i);
-    }
-
+    long writeConnectivity(VectType& connectivity_inout,
+                           ConnectivityVertexOrder whichOrder) const;
 
     /**
      * Write the offsets values  - directly to the pointer given in parameter 1
@@ -583,24 +274,7 @@ public:
        Returns number of offset values written + 1
    */
     template <typename Integer>
-    long writeOffsetsCells(Integer* offsets_inout, long max_size = 0)
-    {
-        if (max_size < ncells_) {
-            // assert(max_size >= ncells_);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeOffsetsCells( T*  offsets_inout )  " + " Input objects max_size ("
-                          + std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
-                          + std::to_string(ncells_) + ")");
-        }
-        long i = 1;
-        offsets_inout[0] = 0;
-        for (const auto& cit : elements(gridView_, dunePartition_)) {
-            auto cell_corners = cit.geometry().corners();
-            offsets_inout[i] = offsets_inout[i - 1] + cell_corners;
-            i++;
-        }
-        return (i); // This should be 1 greater than ncells_
-    }
+    long writeOffsetsCells(Integer* offsets_inout, long max_size = 0) const;
 
     /**
       * Write the offsets values  -  directly to a VectType object given in parameter 1
@@ -613,29 +287,7 @@ public:
         Returns number of offset values written + 1
     */
     template <typename VectType>
-    long writeOffsetsCells(VectType& offsets_inout)
-    {
-        const std::size_t check_size = offsets_inout.size();
-        if (check_size < static_cast<std::size_t>(ncells_)) {
-            // assert(check_size >= ncells_);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeOffsetsCells( VectType& "
-                      "offsets_inout )  "
-                          + " Input objects size (" + std::to_string(offsets_inout.size())
-                          + ") is not sufficient to fit the ncells_ values (" + std::to_string(ncells_) + ")");
-        }
-
-        // using VT = decltype(offsets_inout.data()[0]);
-
-        long i = 1;
-        offsets_inout.data()[0] = 0;
-        for (const auto& cit : elements(gridView_, dunePartition_)) {
-            auto cell_corners = cit.geometry().corners();
-            offsets_inout.data()[i] = offsets_inout.data()[i - 1] + cell_corners;
-            i++;
-        }
-        return (i); // This should be 1 greater than ncells_
-    }
+    long writeOffsetsCells(VectType& offsets_inout) const;
 
     /**
      * Write the cell types values  - directly to the pointer given in parameter 1
@@ -648,22 +300,7 @@ public:
         Returns number of cells type values written
     */
     template <typename Integer>
-    long writeCellTypes(Integer* types_inout, long max_size = 0)
-    {
-        if (max_size < ncells_) {
-            // assert(max_size >= ncells_);
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeCellTypes( T*  types_inout )  " + " Input objects max_size ("
-                          + std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
-                          + std::to_string(ncells_) + ")");
-        }
-        int i = 0;
-        for (const auto& cit : elements(gridView_, dunePartition_)) {
-            Integer vtktype = static_cast<Integer>(Dune::VTK::geometryType(cit.type()));
-            types_inout[i++] = vtktype;
-        }
-        return (i);
-    }
+    long writeCellTypes(Integer* types_inout, long max_size = 0) const;
 
     /**
      * Write the cell types values  - directly to the VectType object given in parameter 1
@@ -674,91 +311,45 @@ public:
         Returns number of cells type values written
     */
     template <typename VectType>
-    long writeCellTypes(VectType& types_inout)
-    {
-        const std::size_t check_size = types_inout.size();
+    long writeCellTypes(VectType& types_inout) const;
 
-        if (check_size < static_cast<std::size_t>(ncells_)) {
-            OPM_THROW(std::runtime_error,
-                      "Opm::GridDataOutput::writeCellTypes( VectType&  types_inout )  " + " Input objects check_size ("
-                          + std::to_string(check_size) + ") is not sufficient to fit the ncells_ values ("
-                          + std::to_string(ncells_) + ")");
-        }
-
-        int i = 0;
-        for (const auto& cit : elements(gridView_, dunePartition_)) {
-            int vtktype = static_cast<int>(Dune::VTK::geometryType(cit.type()));
-            types_inout.data()[i++] = vtktype;
-        }
-        return (i);
-    }
-
-    std::string getPartitionTypeString()
-    {
-        if (this->dunePartition_ == Dune::Partitions::all)
-            return (std::string("Dune::Partitions::all"));
-        if (this->dunePartition_ == Dune::Partitions::interior)
-            return (std::string("Dune::Partitions::interior"));
-        if (this->dunePartition_ == Dune::Partitions::interiorBorder)
-            return (std::string("Dune::Partitions::interiorBorder"));
-        if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlap)
-            return (std::string("Dune::Partitions::interiorBorderOverlap"));
-        if (this->dunePartition_ == Dune::Partitions::front)
-            return (std::string("Dune::Partitions::front"));
-        if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlapFront)
-            return (std::string("Dune::Partitions::InteriorBorderOverlapFront"));
-        if (this->dunePartition_ == Dune::Partitions::border)
-            return (std::string("Dune::Partitions::border"));
-        if (this->dunePartition_ == Dune::Partitions::ghost)
-            return (std::string("Dune::Partitions::ghost"));
-
-        return (std::string("Unknown Dune::PartitionSet<>"));
-    }
+    std::string getPartitionTypeString() const;
 
     Dune::PartitionSet<partitions> getPartition(void)
     {
-        return (this->dunePartition_);
+        return this->dunePartition_;
     }
 
-    void printGridDetails(std::ostream& outstr)
+    void printGridDetails(std::ostream& outstr) const;
+
+    int getNCells() const
     {
-        outstr << "Dune Partition = " << partition_value_ << ", " << getPartitionTypeString() << std::endl;
-        outstr << "ncells_: " << getNCells() << std::endl;
-        outstr << "nvertices_: " << getNVertices() << std::endl;
-        outstr << "ncorners_: " << getNCorners() << std::endl;
+        return ncells_;
     }
 
-    int getNCells()
+    int getNVertices() const
     {
-        return (ncells_);
+        return nvertices_;
     }
 
-    int getNVertices()
+    int getNCorners() const
     {
-        return (nvertices_);
+        return ncorners_;
     }
 
-    int getNCorners()
+    std::string getError() const
     {
-        return (ncorners_);
-    }
-
-    std::string getError()
-    {
-        return error_strm_.str();
+        return error_;
     }
 
     void clearError()
     {
-        error_strm_.str("");
+        error_.clear();
     }
 
-    bool hasError()
+    bool hasError() const
     {
-        if (error_strm_.str().length() > 0)
-            return true;
-        else
-            return false;
+        return !error_.empty();
     }
 
 protected:
@@ -783,7 +374,7 @@ protected:
     int dimw_; // dimensions of the input grid
 
 private:
-    std::stringstream error_strm_;
+    std::string error_;
 };
 
 } // namespace Opm::GridDataOutput
diff --git a/opm/simulators/utils/GridDataOutput_impl.hpp b/opm/simulators/utils/GridDataOutput_impl.hpp
new file mode 100644
index 000000000..3b4b62912
--- /dev/null
+++ b/opm/simulators/utils/GridDataOutput_impl.hpp
@@ -0,0 +1,526 @@
+/*
+  Copyright 2023 Inria, Bretagne–Atlantique Research Center
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include <opm/simulators/utils/GridDataOutput.hpp>
+
+#include <dune/grid/common/rangegenerators.hh>
+#include <dune/grid/io/file/vtk/common.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <cassert>
+#include <cstddef>
+#include <iterator>
+#include <ostream>
+
+namespace Opm::GridDataOutput {
+
+template <class GridView, unsigned int partitions>
+SimMeshDataAccessor<GridView,partitions>::
+ SimMeshDataAccessor(const GridView& gridView, Dune::PartitionSet<partitions> dunePartition)
+    : gridView_(gridView)
+    , dunePartition_(dunePartition)
+{
+    dimw_ = GridView::dimension; // this is an enum
+    partition_value_ = dunePartition.value;
+    countEntities();
+}
+
+template <class GridView, unsigned int partitions>
+bool SimMeshDataAccessor<GridView,partitions>::polyhedralCellPresent() const
+{
+    for (const auto& cit : elements(gridView_, dunePartition_)) {
+        auto corner_geom = cit.geometry();
+        if (Dune::VTK::geometryType(corner_geom.type()) == Dune::VTK::polyhedron) {
+            return true;
+        }
+    }
+    return false;
+}
+
+template <class GridView, unsigned int partitions>
+void SimMeshDataAccessor<GridView,partitions>::countEntities()
+{
+    // We include all the vertices for this ranks partition
+    const auto& vert_partition = vertices(gridView_, Dune::Partitions::all);
+    nvertices_ = std::distance(vert_partition.begin(), vert_partition.end());
+
+    const auto& cell_partition = elements(gridView_, dunePartition_);
+    ncells_ = 0;
+    ncorners_ = 0;
+    for (const auto& cit : cell_partition) {
+        auto corner_geom = cit.geometry();
+        ncorners_ += corner_geom.corners();
+        ++ncells_;
+    }
+}
+
+template <class GridView, unsigned int partitions>
+template <typename T>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size) const
+{
+    if (max_size < nvertices_) {
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints( T&  x_inout,  T&  "
+                  "y_inout, T& z_inout )  "
+                      + " Input objects max_size (" + std::to_string(max_size)
+                      + ") is not sufficient to fit the nvertices_ values (" + std::to_string(nvertices_) + ")");
+    }
+
+    long i = 0;
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
+            x_inout[i] = static_cast<T>(xyz_local[0]);
+            y_inout[i] = static_cast<T>(xyz_local[1]);
+            z_inout[i] = static_cast<T>(xyz_local[2]);
+            i++;
+        }
+    } else if (dimw_ == 2) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
+            x_inout[i] = static_cast<T>(xyz_local[0]);
+            y_inout[i] = static_cast<T>(xyz_local[1]);
+            z_inout[i] = static_cast<T>(0.0);
+            i++;
+        }
+    }
+    assert(i == nvertices_); // As we are templated on the
+                             // Dune::PartitionSet<partitions>, this cannot change
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout) const
+{
+    const std::size_t check_size_x = x_inout.size();
+    const std::size_t check_size_y = y_inout.size();
+    const std::size_t check_size_z = z_inout.size();
+
+    using VT = decltype(x_inout.data()[0]);
+
+    if ((check_size_x < static_cast<std::size_t>(nvertices_)) ||
+        (check_size_y < static_cast<std::size_t>(nvertices_)) ||
+        (check_size_z < static_cast<std::size_t>(nvertices_))) {
+        // assert(check_size >= nvertices_);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints( VectType&  x_inout,  VectType&  "
+                  "y_inout, VectType& z_inout )  At least one of the inputs"
+                      + "  object x size " + std::to_string(check_size_x) + "  object y size "
+                      + std::to_string(check_size_y) + "  object z size " + std::to_string(check_size_z)
+                      + " is not sufficient to fit the nvertices_ values( " + std::to_string(nvertices_) + " )");
+    }
+
+    long i = 0;
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
+            x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
+            y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
+            z_inout.data()[i] = static_cast<VT>(xyz_local[2]);
+            i++;
+        }
+    } else if (dimw_ == 2) {
+        double td = 0.0;
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
+            x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
+            y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
+            z_inout.data()[i] = static_cast<VT>(td);
+            i++;
+        }
+    }
+    assert(i == nvertices_); // As we are templated on the
+                             // Dune::PartitionSet<partitions>, this cannot change
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename T>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints_AOS(T* xyz_inout, long max_size) const
+{
+    if (max_size < nvertices_ * 3) {
+        assert(max_size >= nvertices_ * 3);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints_AOS( T*  xyz_inout )  " + " Input objects max_size ("
+                      + std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+                      + std::to_string(nvertices_ * 3) + ")");
+    }
+
+    long i = 0;
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout[i++] = static_cast<T>(xyz_local[0]);
+            xyz_inout[i++] = static_cast<T>(xyz_local[1]);
+            xyz_inout[i++] = static_cast<T>(xyz_local[2]);
+        }
+    } else if (dimw_ == 2) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout[i++] = static_cast<T>(xyz_local[0]);
+            xyz_inout[i++] = static_cast<T>(xyz_local[1]);
+            xyz_inout[i++] = static_cast<T>(0.0);
+        }
+    }
+    return i / 3;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints_AOS(VectType& xyz_inout) const
+{
+    const std::size_t check_size = xyz_inout.size();
+
+    using VT = decltype(xyz_inout.data()[0]);
+
+    if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
+        assert(check_size >= nvertices_ * 3);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints_AOS( VectType&  xyz_inout )  "
+                      + " Input objects check_size (" + std::to_string(check_size)
+                      + ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+                      + ")");
+    }
+
+    long i = 0;
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout.data()[i++] = static_cast<VT>(xyz_local[0]);
+            xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
+            xyz_inout[i++] = static_cast<VT>(xyz_local[2]);
+        }
+    } else if (dimw_ == 2) {
+        double td = 0.0;
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout[i++] = static_cast<VT>(xyz_local[0]);
+            xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
+            xyz_inout[i++] = static_cast<VT>(td);
+        }
+    }
+    return i / 3;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename T>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints_SOA(T* xyz_inout, long max_size) const
+{
+    if (max_size < nvertices_ * 3) {
+        // assert(max_size >= nvertices_ * 3);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints_SOA( T&  xyz_inout )  " + " Input objects max_size ("
+                      + std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+                      + std::to_string(nvertices_ * 3) + ")");
+    }
+
+    long i = 0;
+    // Get offsets into structure
+    T* xyz_inout_y = xyz_inout + nvertices_;
+    T* xyz_inout_z = xyz_inout + (2 * nvertices_);
+
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout[i] = static_cast<T>(xyz_local[0]);
+            xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
+            xyz_inout_z[i] = static_cast<T>(xyz_local[2]);
+            i++;
+        }
+    } else if (dimw_ == 2) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout[i] = static_cast<T>(xyz_local[0]);
+            xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
+            xyz_inout_z[i] = static_cast<T>(0.0);
+            i++;
+        }
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeGridPoints_SOA(VectType& xyz_inout) const
+{
+    const std::size_t check_size = xyz_inout.size();
+
+    if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
+        // assert(check_size >= nvertices_ * 3);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeGridPoints_SOA( VectType&  xyz_inout )  "
+                      + " Input objects check_size (" + std::to_string(check_size)
+                      + ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+                      + ")");
+    }
+
+    using VT = decltype(xyz_inout.data()[0]);
+
+    long i = 0;
+    // Get offsets into structure
+    VT* xyz_inout_y = xyz_inout.data() + nvertices_;
+    VT* xyz_inout_z = xyz_inout.data() + (2 * nvertices_);
+
+    if (dimw_ == 3) {
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
+            xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
+            xyz_inout_z[i] = static_cast<VT>(xyz_local[2]);
+            i++;
+        }
+    } else if (dimw_ == 2) {
+        double td = 0.0;
+        for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
+            auto xyz_local = vit.geometry().corner(0);
+            xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
+            xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
+            xyz_inout_z[i] = static_cast<VT>(td);
+            i++;
+        }
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename Integer>
+long SimMeshDataAccessor<GridView,partitions>::
+writeConnectivity(Integer* connectivity_inout,
+                  ConnectivityVertexOrder whichOrder, long max_size) const
+{
+    if (max_size < ncorners_) {
+
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeConnectivity( T*  connectivity_inout,... )  "
+                      + " Input max_size value (" + std::to_string(max_size)
+                      + ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
+    }
+
+    long i = 0;
+    if (whichOrder == DUNE) {
+        // DUNE order
+        for (const auto& cit : elements(gridView_, dunePartition_)) {
+            auto cell_corners = cit.geometry().corners();
+            for (auto vx = 0; vx < cell_corners; ++vx) {
+                const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
+                connectivity_inout[i + vx] = vxIdx;
+            }
+            i += cell_corners;
+        }
+    } else {
+        // VTK order
+        for (const auto& cit : elements(gridView_, dunePartition_)) {
+            auto cell_corners = cit.geometry().corners();
+            for (auto vx = 0; vx < cell_corners; ++vx) {
+                const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
+                int vtkOrder;
+                vtkOrder = Dune::VTK::renumber(cit.type(), vx);
+                connectivity_inout[i + vtkOrder] = vxIdx;
+            }
+            i += cell_corners;
+        }
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeConnectivity(VectType& connectivity_inout,
+                  ConnectivityVertexOrder whichOrder) const
+{
+    const std::size_t check_size = connectivity_inout.size();
+
+    if (check_size < static_cast<std::size_t>(ncorners_)) {
+        // assert(check_size >= ncorners_);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeConnectivity( VectType&  "
+                  "connectivity_inout )  "
+                      + " Input objects size (" + std::to_string(check_size)
+                      + ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
+    }
+
+    long i = 0;
+    if (whichOrder == DUNE) {
+        // DUNE order
+        for (const auto& cit : elements(gridView_, dunePartition_)) {
+            auto cell_corners = cit.geometry().corners();
+            for (auto vx = 0; vx < cell_corners; ++vx) {
+                const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
+                connectivity_inout.data()[i + vx] = vxIdx;
+            }
+            i += cell_corners;
+        }
+    } else {
+        // VTK order
+        for (const auto& cit : elements(gridView_, dunePartition_)) {
+            auto cell_corners = cit.geometry().corners();
+            for (auto vx = 0; vx < cell_corners; ++vx) {
+                const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
+                int vtkOrder;
+                vtkOrder = Dune::VTK::renumber(cit.type(), vx);
+                connectivity_inout.data()[i + vtkOrder] = vxIdx;
+            }
+            i += cell_corners;
+        }
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename Integer>
+long SimMeshDataAccessor<GridView,partitions>::
+writeOffsetsCells(Integer* offsets_inout, long max_size) const
+{
+    if (max_size < ncells_) {
+        // assert(max_size >= ncells_);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeOffsetsCells( T*  offsets_inout )  " + " Input objects max_size ("
+                      + std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+                      + std::to_string(ncells_) + ")");
+    }
+    long i = 1;
+    offsets_inout[0] = 0;
+    for (const auto& cit : elements(gridView_, dunePartition_)) {
+        auto cell_corners = cit.geometry().corners();
+        offsets_inout[i] = offsets_inout[i - 1] + cell_corners;
+        i++;
+    }
+    return i; // This should be 1 greater than ncells_
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeOffsetsCells(VectType& offsets_inout) const
+{
+    const std::size_t check_size = offsets_inout.size();
+    if (check_size < static_cast<std::size_t>(ncells_)) {
+        // assert(check_size >= ncells_);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeOffsetsCells( VectType& "
+                  "offsets_inout )  "
+                      + " Input objects size (" + std::to_string(offsets_inout.size())
+                      + ") is not sufficient to fit the ncells_ values (" + std::to_string(ncells_) + ")");
+    }
+
+    // using VT = decltype(offsets_inout.data()[0]);
+
+    long i = 1;
+    offsets_inout.data()[0] = 0;
+    for (const auto& cit : elements(gridView_, dunePartition_)) {
+        auto cell_corners = cit.geometry().corners();
+        offsets_inout.data()[i] = offsets_inout.data()[i - 1] + cell_corners;
+        i++;
+    }
+    return i; // This should be 1 greater than ncells_
+}
+
+template <class GridView, unsigned int partitions>
+template <typename Integer>
+long SimMeshDataAccessor<GridView,partitions>::
+writeCellTypes(Integer* types_inout, long max_size) const
+{
+    if (max_size < ncells_) {
+        // assert(max_size >= ncells_);
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeCellTypes( T*  types_inout )  " + " Input objects max_size ("
+                      + std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+                      + std::to_string(ncells_) + ")");
+    }
+    int i = 0;
+    for (const auto& cit : elements(gridView_, dunePartition_)) {
+        Integer vtktype = static_cast<Integer>(Dune::VTK::geometryType(cit.type()));
+        types_inout[i++] = vtktype;
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+template <typename VectType>
+long SimMeshDataAccessor<GridView,partitions>::
+writeCellTypes(VectType& types_inout) const
+{
+    const std::size_t check_size = types_inout.size();
+
+    if (check_size < static_cast<std::size_t>(ncells_)) {
+        OPM_THROW(std::runtime_error,
+                  "Opm::GridDataOutput::writeCellTypes( VectType&  types_inout )  " + " Input objects check_size ("
+                      + std::to_string(check_size) + ") is not sufficient to fit the ncells_ values ("
+                      + std::to_string(ncells_) + ")");
+    }
+
+    int i = 0;
+    for (const auto& cit : elements(gridView_, dunePartition_)) {
+        int vtktype = static_cast<int>(Dune::VTK::geometryType(cit.type()));
+        types_inout.data()[i++] = vtktype;
+    }
+    return i;
+}
+
+template <class GridView, unsigned int partitions>
+std::string SimMeshDataAccessor<GridView,partitions>::
+getPartitionTypeString() const
+{
+    if (this->dunePartition_ == Dune::Partitions::all) {
+        return "Dune::Partitions::all";
+    }
+    if (this->dunePartition_ == Dune::Partitions::interior) {
+        return "Dune::Partitions::interior";
+    }
+    if (this->dunePartition_ == Dune::Partitions::interiorBorder) {
+        return "Dune::Partitions::interiorBorder";
+    }
+    if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlap) {
+        return "Dune::Partitions::interiorBorderOverlap";
+    }
+    if (this->dunePartition_ == Dune::Partitions::front) {
+        return "Dune::Partitions::front";
+    }
+    if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlapFront) {
+        return "Dune::Partitions::InteriorBorderOverlapFront";
+    }
+    if (this->dunePartition_ == Dune::Partitions::border) {
+        return "Dune::Partitions::border";
+    }
+    if (this->dunePartition_ == Dune::Partitions::ghost) {
+        return "Dune::Partitions::ghost";
+    }
+
+    return "Unknown Dune::PartitionSet<>";
+}
+
+template <class GridView, unsigned int partitions>
+void SimMeshDataAccessor<GridView,partitions>::
+printGridDetails(std::ostream& outstr) const
+{
+    outstr << "Dune Partition = " << partition_value_ << ", " << getPartitionTypeString() << std::endl;
+    outstr << "ncells_: " << getNCells() << std::endl;
+    outstr << "nvertices_: " << getNVertices() << std::endl;
+    outstr << "ncorners_: " << getNCorners() << std::endl;
+}
+
+} // namespace Opm::GridDataOutput