diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index d483ee92421713aa1d9a38b7e99118c7235a6b1a..736bc94709c6e2d1ade3273a6f1b52541a41448e 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -310,7 +310,7 @@ void GmshReader::__readVertices()
   }
 
   __verticesNumbers.resize(numberOfVerices);
-  __vertices = Kokkos::View<R3*>("vertices",numberOfVerices);
+  __vertices = Array<R3>(numberOfVerices);
 
   if (not __binary) {
     for (int i=0; i<numberOfVerices; ++i) {
@@ -552,38 +552,32 @@ GmshReader::__proceedData()
         switch (i) {
           // Supported entities
           case 0: { // edges
-            __edges = Kokkos::View<Edge*>("edges",
-                                          elementNumber[i]);
+            __edges = Array<Edge>(elementNumber[i]);
             __edges_ref.resize(elementNumber[i]);
             break;
           }
           case 1: { // triangles
-            __triangles = Kokkos::View<Triangle*>("triangles",
-                                                  elementNumber[i]);
+            __triangles = Array<Triangle>(elementNumber[i]);
             __triangles_ref.resize(elementNumber[i]);
             break;
           }
           case  2: { // quadrangles
-            __quadrangles = Kokkos::View<Quadrangle*>("quadrangles",
-                                                      elementNumber[i]);
+            __quadrangles = Array<Quadrangle>(elementNumber[i]);
             __quadrangles_ref.resize(elementNumber[i]);
             break;
           }
           case 3: { // tetrahedra
-            __tetrahedra = Kokkos::View<Tetrahedron*>("tetrahedra",
-                                                      elementNumber[i]);
+            __tetrahedra = Array<Tetrahedron>(elementNumber[i]);
             __tetrahedra_ref.resize(elementNumber[i]);
             break;
           }
           case  4: {// hexahedra
-            __hexahedra = Kokkos::View<Hexahedron*>("hexahedra",
-                                                    elementNumber[i]);
+            __hexahedra = Array<Hexahedron>(elementNumber[i]);
             __hexahedra_ref.resize(elementNumber[i]);
           }
             // Ignored entities
           case 14: {// point
-            __points = Kokkos::View<Point*>("points",
-                                            elementNumber[i]);
+            __points = Array<Point>(elementNumber[i]);
             __points_ref.resize(elementNumber[i]);
             break;
           }
@@ -798,7 +792,7 @@ GmshReader::__proceedData()
     std::vector<CellType> cell_type_vector(nb_cells);
 
     std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
-    const size_t nb_tetrahedra = __tetrahedra.extent(0);
+    const size_t nb_tetrahedra = __tetrahedra.size();
     for (size_t j=0; j<nb_tetrahedra; ++j) {
       cell_by_node_vector[j].resize(4);
       for (int r=0; r<4; ++r) {
@@ -806,7 +800,7 @@ GmshReader::__proceedData()
       }
       cell_type_vector[j] = CellType::Tetrahedron;
     }
-    const size_t nb_hexahedra = __hexahedra.extent(0);
+    const size_t nb_hexahedra = __hexahedra.size();
     for (size_t j=0; j<nb_hexahedra; ++j) {
       const size_t jh = nb_tetrahedra+j;
       cell_by_node_vector[jh].resize(8);
@@ -821,13 +815,13 @@ GmshReader::__proceedData()
     Connectivity3D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-    for (unsigned int f=0; f<__triangles.extent(0); ++f) {
+    for (unsigned int f=0; f<__triangles.size(); ++f) {
       const unsigned int face_number
           = connectivity.getFaceNumber({__triangles[f][0], __triangles[f][1], __triangles[f][2]});
       const unsigned int& ref = __triangles_ref[f];
       ref_faces_map[ref].push_back(face_number);
     }
-    for (unsigned int f=0; f<__quadrangles.extent(0); ++f) {
+    for (unsigned int f=0; f<__quadrangles.size(); ++f) {
       const unsigned int face_number
           = connectivity.getFaceNumber({__quadrangles[f][0], __quadrangles[f][1],
                                         __quadrangles[f][2], __quadrangles[f][3]});
@@ -835,7 +829,7 @@ GmshReader::__proceedData()
       ref_faces_map[ref].push_back(face_number);
     }
     for (const auto& ref_face_list : ref_faces_map) {
-      Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
+      Array<unsigned int> face_list(ref_face_list.second.size());
       for (size_t j=0; j<ref_face_list.second.size(); ++j) {
         face_list[j]=ref_face_list.second[j];
       }
@@ -847,7 +841,7 @@ GmshReader::__proceedData()
     typedef TinyVector<3, double> Rd;
 
     NodeValue<Rd> xr(connectivity);
-    for (size_t i=0; i<__vertices.extent(0); ++i) {
+    for (size_t i=0; i<__vertices.size(); ++i) {
       xr[i][0] = __vertices[i][0];
       xr[i][1] = __vertices[i][1];
       xr[i][2] = __vertices[i][2];
@@ -860,7 +854,7 @@ GmshReader::__proceedData()
     std::vector<CellType> cell_type_vector(nb_cells);
 
     std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
-    const size_t nb_triangles = __triangles.extent(0);
+    const size_t nb_triangles = __triangles.size();
     for (size_t j=0; j<nb_triangles; ++j) {
       cell_by_node_vector[j].resize(3);
       for (int r=0; r<3; ++r) {
@@ -869,7 +863,7 @@ GmshReader::__proceedData()
       cell_type_vector[j] = CellType::Triangle;
     }
 
-    const size_t nb_quadrangles = __quadrangles.extent(0);
+    const size_t nb_quadrangles = __quadrangles.size();
     for (size_t j=0; j<nb_quadrangles; ++j) {
       const size_t jq = j+nb_triangles;
       cell_by_node_vector[jq].resize(4);
@@ -884,14 +878,14 @@ GmshReader::__proceedData()
     Connectivity2D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-    for (unsigned int e=0; e<__edges.extent(0); ++e) {
+    for (unsigned int e=0; e<__edges.size(); ++e) {
       const unsigned int edge_number = connectivity.getFaceNumber({__edges[e][0], __edges[e][1]});
       const unsigned int& ref = __edges_ref[e];
       ref_faces_map[ref].push_back(edge_number);
     }
 
     for (const auto& ref_face_list : ref_faces_map) {
-      Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
+      Array<unsigned int> face_list(ref_face_list.second.size());
       for (size_t j=0; j<ref_face_list.second.size(); ++j) {
         face_list[j]=ref_face_list.second[j];
       }
@@ -900,14 +894,14 @@ GmshReader::__proceedData()
     }
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-    for (unsigned int r=0; r<__points.extent(0); ++r) {
+    for (unsigned int r=0; r<__points.size(); ++r) {
       const unsigned int point_number = __points[r];
       const unsigned int& ref = __points_ref[r];
       ref_points_map[ref].push_back(point_number);
     }
 
     for (const auto& ref_point_list : ref_points_map) {
-      Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
+      Array<unsigned int> point_list(ref_point_list.second.size());
       for (size_t j=0; j<ref_point_list.second.size(); ++j) {
         point_list[j]=ref_point_list.second[j];
       }
@@ -919,7 +913,7 @@ GmshReader::__proceedData()
     typedef TinyVector<2, double> Rd;
 
     NodeValue<Rd> xr(connectivity);
-    for (size_t i=0; i<__vertices.extent(0); ++i) {
+    for (size_t i=0; i<__vertices.size(); ++i) {
       xr[i][0] = __vertices[i][0];
       xr[i][1] = __vertices[i][1];
     }
@@ -944,14 +938,14 @@ GmshReader::__proceedData()
     Connectivity1D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-    for (unsigned int r=0; r<__points.extent(0); ++r) {
+    for (unsigned int r=0; r<__points.size(); ++r) {
       const unsigned int point_number = __points[r];
       const unsigned int& ref = __points_ref[r];
       ref_points_map[ref].push_back(point_number);
     }
 
     for (const auto& ref_point_list : ref_points_map) {
-      Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
+      Array<unsigned int> point_list(ref_point_list.second.size());
       for (size_t j=0; j<ref_point_list.second.size(); ++j) {
         point_list[j]=ref_point_list.second[j];
       }
@@ -963,7 +957,7 @@ GmshReader::__proceedData()
     typedef TinyVector<1, double> Rd;
 
     NodeValue<Rd> xr(connectivity);
-    for (size_t i=0; i<__vertices.extent(0); ++i) {
+    for (size_t i=0; i<__vertices.size(); ++i) {
       xr[i][0] = __vertices[i][0];
     }
 
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 8825a67df3d6cf1593ff54bb6af8b7a96b079b5b..5a43e3f6116021a8a27ea36b3b36db43cf44babd 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -1,7 +1,7 @@
 #ifndef GMSH_READER_HPP
 #define GMSH_READER_HPP
 
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
 #include <TinyVector.hpp>
 
 #include <Mesh.hpp>
@@ -78,25 +78,25 @@ private:
    */
   std::vector<int> __verticesNumbers;
 
-  Kokkos::View<R3*> __vertices;
+  Array<R3> __vertices;
 
   typedef unsigned int Point;
-  Kokkos::View<Point*> __points;
+  Array<Point> __points;
   std::vector<int> __points_ref;
   typedef TinyVector<2,unsigned int> Edge;
-  Kokkos::View<Edge*> __edges;
+  Array<Edge> __edges;
   std::vector<int> __edges_ref;
   typedef TinyVector<3,unsigned int> Triangle;
-  Kokkos::View<Triangle*> __triangles;
+  Array<Triangle> __triangles;
   std::vector<int> __triangles_ref;
   typedef TinyVector<4,unsigned int> Quadrangle;
-  Kokkos::View<Quadrangle*> __quadrangles;
+  Array<Quadrangle> __quadrangles;
   std::vector<int> __quadrangles_ref;
   typedef TinyVector<4,unsigned int> Tetrahedron;
-  Kokkos::View<Tetrahedron*> __tetrahedra;
+  Array<Tetrahedron> __tetrahedra;
   std::vector<int> __tetrahedra_ref;
   typedef TinyVector<8,unsigned int> Hexahedron;
-  Kokkos::View<Hexahedron*> __hexahedra;
+  Array<Hexahedron> __hexahedra;
   std::vector<int> __hexahedra_ref;
 
   /**
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index 13ef536a962f80abf9aa812fe40ed155a73c479a..ef4b5310768f0016b7b3e8f118e9c5d29bc85be4 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -97,7 +97,8 @@ class ItemValue
       : m_is_built{true},
         m_values(connectivity.numberOf<item_type>())
   {
-    ;
+    static_assert(not std::is_const<DataType>(),
+                  "Cannot allocate ItemValue of const data: only view is supported"); ;
   }
 
   ~ItemValue() = default;
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 91a0943dea549e8cf101d6045659201236ecd3b2..8dff6d526d5342d5f3862c2fa84dbd3845c2b1af 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -1,7 +1,7 @@
 #ifndef MESH_NODE_BOUNDARY_HPP
 #define MESH_NODE_BOUNDARY_HPP
 
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
 #include <ItemValue.hpp>
 
 #include <Kokkos_Vector.hpp>
@@ -18,12 +18,12 @@ template <size_t dimension>
 class MeshNodeBoundary
 {
  protected:
-  Kokkos::View<const unsigned int*> m_node_list;
+  Array<const unsigned int> m_node_list;
  public:
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
   MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
-  const Kokkos::View<const unsigned int*>& nodeList() const
+  const Array<const unsigned int>& nodeList() const
   {
     return m_node_list;
   }
@@ -36,8 +36,8 @@ class MeshNodeBoundary
     const auto& face_to_cell_matrix
         = mesh.connectivity().getMatrix(ItemType::face, ItemType::cell);
 
-    const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
-    Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
+    const Array<const unsigned int>& face_list = ref_face_list.faceList();
+    Kokkos::parallel_for(face_list.size(), KOKKOS_LAMBDA(const int& l){
         const auto& face_cells = face_to_cell_matrix.rowConst(face_list[l]);
         if (face_cells.length>1) {
           std::cerr << "internal faces cannot be used to define mesh boundaries\n";
@@ -47,12 +47,12 @@ class MeshNodeBoundary
 
     Kokkos::vector<unsigned int> node_ids;
     // not enough but should reduce significantly the number of resizing
-    node_ids.reserve(dimension*face_list.extent(0));
+    node_ids.reserve(dimension*face_list.size());
     const auto& face_to_node_matrix
         = mesh.connectivity().getMatrix(ItemType::face,
                                         ItemType::node);
 
-    for (size_t l=0; l<face_list.extent(0); ++l) {
+    for (size_t l=0; l<face_list.size(); ++l) {
       const size_t face_number = face_list[l];
       const auto& face_nodes = face_to_node_matrix.rowConst(face_number);
 
@@ -64,7 +64,7 @@ class MeshNodeBoundary
     auto last = std::unique(node_ids.begin(), node_ids.end());
     node_ids.resize(std::distance(node_ids.begin(), last));
 
-    Kokkos::View<unsigned int*> node_list("node_list", node_ids.size());
+    Array<unsigned int> node_list(node_ids.size());
     Kokkos::parallel_for(node_ids.size(), KOKKOS_LAMBDA(const int& r){
         node_list[r] = node_ids[r];
       });
@@ -153,7 +153,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
 
   const NodeValue<const R2>& xr = mesh.xr();
 
-  Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) {
+  Kokkos::parallel_for(m_node_list.size(), KOKKOS_LAMBDA(const size_t& r) {
       const R2& x = xr[m_node_list[r]];
       if ((x-origin,normal)>1E-13*length) {
         std::cerr << "this FlatBoundary is not flat!\n";
@@ -171,7 +171,7 @@ _getNormal(const MeshType& mesh)
   static_assert(MeshType::dimension == 1);
   typedef TinyVector<1,double> R;
 
-  if (m_node_list.extent(0) != 1) {
+  if (m_node_list.size() != 1) {
     std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
     std::exit(1);
   }
@@ -196,7 +196,7 @@ _getNormal(const MeshType& mesh)
   R2 xmax(-std::numeric_limits<double>::max(),
           -std::numeric_limits<double>::max());
 
-  for (size_t r=0; r<m_node_list.extent(0); ++r) {
+  for (size_t r=0; r<m_node_list.size(); ++r) {
     const R2& x = xr[m_node_list[r]];
     if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
       xmin = x;
@@ -243,7 +243,7 @@ _getNormal(const MeshType& mesh)
 
   const NodeValue<const R3>& xr = mesh.xr();
 
-  for (size_t r=0; r<m_node_list.extent(0); ++r) {
+  for (size_t r=0; r<m_node_list.size(); ++r) {
     const R3& x = xr[m_node_list[r]];
     if (x[0]<xmin[0]) {
       xmin = x;
diff --git a/src/mesh/RefFaceList.hpp b/src/mesh/RefFaceList.hpp
index e891ea2581bc1dbb19a0e82f4178a2a172112574..da08dd63e8a64e303a601cdec79831877f6e5c5b 100644
--- a/src/mesh/RefFaceList.hpp
+++ b/src/mesh/RefFaceList.hpp
@@ -1,14 +1,14 @@
 #ifndef REF_FACE_LIST_HPP
 #define REF_FACE_LIST_HPP
 
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
 #include <RefId.hpp>
 
 class RefFaceList
 {
  private:
   const RefId m_ref_id;
-  const Kokkos::View<const unsigned int*> m_face_id_list;
+  const Array<const unsigned int> m_face_id_list;
 
  public:
   const RefId& refId() const
@@ -16,13 +16,13 @@ class RefFaceList
     return m_ref_id;
   }
 
-  const Kokkos::View<const unsigned int*>& faceList() const
+  const Array<const unsigned int>& faceList() const
   {
     return m_face_id_list;
   }
 
   RefFaceList(const RefId& ref_id,
-              const Kokkos::View<const unsigned int*>& face_id_list)
+              const Array<const unsigned int>& face_id_list)
       : m_ref_id(ref_id),
         m_face_id_list(face_id_list)
   {
diff --git a/src/mesh/RefNodeList.hpp b/src/mesh/RefNodeList.hpp
index 7bdc36fac84b1dd9cb713be5a5a64e9d992232c0..da65635511cc328c278bd03578cdff9ef05e31bc 100644
--- a/src/mesh/RefNodeList.hpp
+++ b/src/mesh/RefNodeList.hpp
@@ -1,14 +1,14 @@
 #ifndef REF_NODE_LIST_HPP
 #define REF_NODE_LIST_HPP
 
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
 #include <RefId.hpp>
 
 class RefNodeList
 {
  private:
   RefId m_ref_id;
-  Kokkos::View<const unsigned int*> m_node_id_list;
+  Array<const unsigned int> m_node_id_list;
 
  public:
   const RefId& refId() const
@@ -16,13 +16,13 @@ class RefNodeList
     return m_ref_id;
   }
 
-  const Kokkos::View<const unsigned int*>& nodeList() const
+  const Array<const unsigned int>& nodeList() const
   {
     return m_node_id_list;
   }
 
   RefNodeList(const RefId& ref_id,
-              const Kokkos::View<const unsigned int*>& node_id_list)
+              const Array<const unsigned int>& node_id_list)
       : m_ref_id(ref_id),
         m_node_id_list(node_id_list)
   {
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index be60e2634e7ae4ff380e36b52d15fc9e77492923..a5b96b56565b77a00a5874334c0fae9aa2545e70 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -2,6 +2,8 @@
 #define SUBITEM_VALUE_PER_ITEM_HPP
 
 #include <Kokkos_StaticCrsGraph.hpp>
+
+#include <Array.hpp>
 #include <ItemType.hpp>
 #include <PastisAssert.hpp>
 
@@ -29,7 +31,7 @@ class SubItemValuePerItem<DataType,
   bool m_is_built{false};
 
   ConnectivityMatrix::HostRowType m_host_row_map;
-  Kokkos::View<DataType*> m_values;
+  Array<DataType> m_values;
 
   // Allows const version to access our data
   friend SubItemValuePerItem<std::add_const_t<DataType>,
@@ -69,14 +71,14 @@ class SubItemValuePerItem<DataType,
     SubView(SubView&&) = default;
 
     KOKKOS_INLINE_FUNCTION
-    SubView(const Kokkos::View<DataType*>& values,
+    SubView(const Array<DataType>& values,
             const size_t& begin,
             const size_t& end)
         : m_sub_values(&(values[begin])),
           m_size(end-begin)
     {
       Assert(begin<=end);
-      Assert(end<=values.extent(0));
+      Assert(end<=values.size());
     }
   };
 
@@ -99,7 +101,7 @@ class SubItemValuePerItem<DataType,
   size_t numberOfValues() const
   {
     Assert(m_is_built);
-    return m_values.extent(0);
+    return m_values.size();
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
@@ -157,7 +159,7 @@ class SubItemValuePerItem<DataType,
     // ensures that const is not lost through copy
     static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
                    or not std::is_const<DataType2>()),
-                  "Cannot assign const  SubItemValuePerItem to a non const SubItemValuePerItem");
+                  "Cannot assign SubItemValuePerItem of const data to SubItemValuePerItem of non-const data");
     m_host_row_map = sub_item_value_per_item.m_host_row_map;
     m_values = sub_item_value_per_item.m_values;
 
@@ -178,11 +180,14 @@ class SubItemValuePerItem<DataType,
   SubItemValuePerItem(const IConnectivity& connectivity)
       : m_is_built{true}
   {
+    static_assert(not std::is_const<DataType>(),
+                  "Cannot allocate SubItemValuePerItem of const data: only view is supported"); ;
+
     ConnectivityMatrix connectivity_matrix
         = connectivity.getMatrix(item_type, sub_item_type);
 
     m_host_row_map = connectivity_matrix.rowsMap();
-    m_values = Kokkos::View<std::remove_const_t<DataType>*>("values", connectivity_matrix.numEntries());
+    m_values = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries());
   }
 
   ~SubItemValuePerItem() = default;
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
index c59eae3d347d6354a6772c3c2001bfbf412595d5..b1fe0d34e98b252296b4b505d82d651712a965c8 100644
--- a/src/scheme/BoundaryCondition.hpp
+++ b/src/scheme/BoundaryCondition.hpp
@@ -4,7 +4,7 @@
 #include <vector>
 #include <memory>
 
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
 
 #include <RefNodeList.hpp>
 #include <MeshNodeBoundary.hpp>
@@ -44,7 +44,7 @@ private:
   const double m_value;
 
   const size_t m_number_of_faces;
-  Kokkos::View<unsigned int*> m_face_list;
+  Array<const unsigned int> m_face_list;
 
 public:
   double value() const
@@ -57,7 +57,7 @@ public:
     return m_number_of_faces;
   }
 
-  const Kokkos::View<const unsigned int*> faceList() const
+  const Array<const unsigned int>& faceList() const
   {
     return m_face_list;
   }
@@ -66,12 +66,13 @@ public:
                             const std::vector<unsigned int>& faces)
     : BoundaryCondition(BoundaryCondition::pressure),
       m_value(value),
-      m_number_of_faces(faces.size()),
-      m_face_list("faces_list", m_number_of_faces)
+      m_number_of_faces(faces.size())
   {
+    Array<unsigned int> face_list(faces.size());
     Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const int& f){
-        m_face_list[f]=faces[f];
+        face_list[f]=faces[f];
       });
+    m_face_list = face_list;
   }
 
   ~PressureBoundaryCondition() = default;
@@ -95,10 +96,10 @@ public:
 
   size_t numberOfNodes() const
   {
-    return m_mesh_flat_node_boundary.nodeList().extent(0);
+    return m_mesh_flat_node_boundary.nodeList().size();
   }
 
-  const Kokkos::View<const unsigned int*> nodeList() const
+  const Array<const unsigned int>& nodeList() const
   {
     return m_mesh_flat_node_boundary.nodeList();
   }
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 8cb91fcd5b70aa8384c40a75104a183fc2b154fd..b702f61e738d888f51cef5b8d86b417bc330561f 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -48,7 +48,7 @@ class Array
       : m_values("anonymous", size)
   {
     static_assert(not std::is_const<DataType>(),
-                  "Cannot build Array of const data of given size");
+                  "Cannot allocate Array of const data: only view is supported");
   }
 
   template <typename DataType2>
@@ -75,6 +75,13 @@ class Array
   KOKKOS_INLINE_FUNCTION
   Array() = default;
 
+  template <typename DataType2>
+  KOKKOS_INLINE_FUNCTION
+  Array(const Array<DataType2>& array)
+  {
+    this->operator=(array);
+  }
+
   KOKKOS_INLINE_FUNCTION
   Array(const Array&) = default;