diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 53feb3375a9de19370d03a002cf869fe4f039c0f..18b39a67a8839e1a893bb7bf2cec50aa60c407ee 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -2,6 +2,7 @@
 #define CONNECTIVITY_MATRIX_HPP
 
 #include <PastisUtils.hpp>
+#include <Array.hpp>
 #include <Kokkos_StaticCrsGraph.hpp>
 
 class ConnectivityMatrix
@@ -20,9 +21,11 @@ class ConnectivityMatrix
     return m_is_built;
   }
 
-  const auto& entries() const
+  auto entries() const
   {
-    return m_host_matrix.entries;
+    return encapsulate(m_host_matrix.entries);
+    // using DataType = typename decltype(m_host_matrix.entries)::value_type;
+    // return Array<DataType>(m_host_matrix.entries);
   }
 
   const auto& rowsMap() const
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index bf7e86ebca7d279a03586a6ad183b24887457b23..d10aa5ac11e33c95aa54012828944c45a34024c1 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -67,7 +67,7 @@ class ItemToItemMatrix
     return m_connectivity_matrix.numEntries();
   }
 
-  const auto& entries() const
+  auto entries() const
   {
     return m_connectivity_matrix.entries();
   }
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index bd1a00f7fe6095e31b4b40f91dec974a9d15f60f..dc3d1bb12018184c6b75a3997a1f889d1f3cc3d3 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -279,14 +279,14 @@ class VTKWriter
         using Rd = TinyVector<MeshType::dimension>;
         const NodeValue<const Rd>& xr = mesh.xr();
         Array<TinyVector<3>> positions(mesh.numberOfNodes());
-        for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
-          for (unsigned short i=0; i<MeshType::dimension; ++i) {
-            positions[r][i] = xr[r][i];
-          }
-          for (unsigned short i=MeshType::dimension; i<3; ++i) {
-            positions[r][i] = 0;
-          }
-        }
+        parallel_for(mesh.numberOfNodes(), PASTIS_LAMBDA(NodeId r) {
+            for (unsigned short i=0; i<MeshType::dimension; ++i) {
+              positions[r][i] = xr[r][i];
+            }
+            for (unsigned short i=MeshType::dimension; i<3; ++i) {
+              positions[r][i] = 0;
+            }
+          });
         _write_array(fout, "Positions", positions);
       }
       fout << "</Points>\n";
@@ -295,12 +295,8 @@ class VTKWriter
       {
         const auto& cell_to_node_matrix
             = mesh.connectivity().cellToNodeMatrix();
-        Array<int> connectivity(cell_to_node_matrix.numberOfEntries());
 
-        for (size_t i=0; i<cell_to_node_matrix.numberOfEntries(); ++i) {
-          connectivity[i] = cell_to_node_matrix.entries()[i];
-        }
-        _write_array(fout, "connectivity", connectivity);
+        _write_array(fout, "connectivity", cell_to_node_matrix.entries());
       }
 
       {
@@ -317,38 +313,45 @@ class VTKWriter
       }
 
       {
-        const auto& cell_to_node_matrix
-            = mesh.connectivity().cellToNodeMatrix();
         Array<int8_t> types(mesh.numberOfCells());
-        for (CellId j=0; j<mesh.numberOfCells(); ++j) {
-          const auto& cell_nodes = cell_to_node_matrix[j];
-          switch (cell_nodes.size()) {
-            case 2: {
-              types[j]=3;
-              break;
-            }
-            case 3: {
-              types[j]=5;
-              break;
-            }
-            case 4: {
-              if (mesh.meshDimension() == 3) {
-                types[j]=10;
-              } else {
+        const auto& cell_type
+            = mesh.connectivity().cellType();
+        parallel_for(mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+            switch (cell_type[j]) {
+              case CellType::Line: {
+                types[j]=3;
+                break;
+              }
+              case CellType::Triangle: {
+                types[j]=5;
+                break;
+              }
+              case CellType::Quadrangle: {
                 types[j]=9;
+                break;
+              }
+              case CellType::Tetrahedron: {
+                types[j]=10;
+                break;
+              }
+              case CellType::Pyramid: {
+                types[j]=14;
+                break;
+              }
+              case CellType::Prism: {
+                types[j]=13;
+                break;
+              }
+              case CellType::Hexahedron: {
+                types[j]=12;
+                break;
+              }
+              default: {
+                std::cerr << __FILE__ << ':' << __LINE__ << ": unknown cell type\n";
+                std::exit(1);
               }
-              break;
-            }
-            case 8: {
-              types[j]=12;
-              break;
-            }
-            default: {
-              types[j]=7;
-              break;
             }
-          }
-        }
+          });
         _write_array(fout, "types", types);
       }
 
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 489ae32254ef0641506181d633de1d481eb1e6b6..a563c6462ecb71bbdc3b3aecf8f9a4c605b56a71 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -38,6 +38,10 @@ class Array
     return image;
   }
 
+  template <typename DataType2, typename ... RT>
+  friend PASTIS_INLINE
+  Array<DataType2> encapsulate(const Kokkos::View<DataType2*, RT...>& values);
+
   PASTIS_INLINE
   DataType& operator[](const index_type& i) const noexcept(NO_ASSERT)
   {
@@ -106,6 +110,15 @@ class Array
   ~Array() = default;
 };
 
+template <typename DataType, typename ... RT>
+PASTIS_INLINE
+Array<DataType> encapsulate(const Kokkos::View<DataType*, RT...>& values)
+{
+  Array<DataType> array;
+  array.m_values = values;
+  return array;
+}
+
 // map, multimap, unordered_map and stack cannot be copied this way
 template <typename Container>
 PASTIS_INLINE