diff --git a/src/mesh/CellType.hpp b/src/mesh/CellType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1f444323b5f5073b9aa09d020cc1fc962b1182a5
--- /dev/null
+++ b/src/mesh/CellType.hpp
@@ -0,0 +1,17 @@
+#ifndef CELL_TYPE_HPP
+#define  CELL_TYPE_HPP
+
+enum class CellType : unsigned short
+{
+  Line,
+
+  Triangle,
+  Quadrangle,
+
+  Tetrahedron,
+  Pyramid,
+  Prism,
+  Hexahedron
+};
+
+#endif // CELL_TYPE_HPP
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 93354d6de0cdb2714fa156e7f809aa6b750050d9..5ea41281219fd8a4e8643c1c3974a8e7ddadf9ae 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -215,14 +215,24 @@ void Connectivity<2>::_computeFaceCellConnectivities()
 
 template<size_t Dimension>
 Connectivity<Dimension>::
-Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+             const std::vector<CellType>& cell_type_vector)
 {
+  Assert(cell_by_node_vector.size() == cell_type_vector.size());
+
   auto& cell_to_node_matrix
       = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
   cell_to_node_matrix = cell_by_node_vector;
 
   Assert(this->numberOfCells()>0);
 
+  {
+    Kokkos::View<CellType*> cell_type("cell_type", this->numberOfCells());
+    Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+        cell_type[j] = cell_type_vector[j];
+      });
+    m_cell_type = cell_type;
+  }
   {
     Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
     Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -239,10 +249,13 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
 
 
 template Connectivity1D::
-Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+             const std::vector<CellType>& cell_type_vector);
 
 template Connectivity2D::
-Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+             const std::vector<CellType>& cell_type_vector);
 
 template Connectivity3D::
-Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+             const std::vector<CellType>& cell_type_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 96433d7d796f1ff66827435dade86b202be724a7..65ac93594c3c5566025c79b568e6c64b6f6783b0 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -16,6 +16,8 @@
 #include <unordered_map>
 #include <algorithm>
 
+#include <CellType.hpp>
+
 #include <RefId.hpp>
 #include <ItemType.hpp>
 #include <RefNodeList.hpp>
@@ -232,6 +234,7 @@ class Connectivity final
 
  private:
   ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
+  Kokkos::View<const CellType*> m_cell_type;
 
   FaceValuePerCell<const bool> m_cell_face_is_reversed;
 
@@ -479,7 +482,8 @@ class Connectivity final
 
   Connectivity(const Connectivity&) = delete;
 
-  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector);
+  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+               const std::vector<CellType>& cell_type_vector);
 
   ~Connectivity()
   {
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 43074d6f9c37936187aa4c5ee880fb21fde53f71..f16ab2fa02e06aad3715eeb7c84028476c7565ae 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -5,6 +5,7 @@
 #include <set>
 #include <rang.hpp>
 
+#include <CellType.hpp>
 #include <Connectivity.hpp>
 
 #include <Mesh.hpp>
@@ -794,6 +795,8 @@ GmshReader::__proceedData()
   if ((dimension3_mask, elementNumber)>0) {
     const size_t nb_cells = (dimension3_mask, elementNumber);
 
+    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);
     for (size_t j=0; j<nb_tetrahedra; ++j) {
@@ -801,6 +804,7 @@ GmshReader::__proceedData()
       for (int r=0; r<4; ++r) {
         cell_by_node_vector[j][r] = __tetrahedra[j][r];
       }
+      cell_type_vector[j] = CellType::Tetrahedron;
     }
     const size_t nb_hexahedra = __hexahedra.extent(0);
     for (size_t j=0; j<nb_hexahedra; ++j) {
@@ -809,9 +813,11 @@ GmshReader::__proceedData()
       for (int r=0; r<8; ++r) {
         cell_by_node_vector[jh][r] = __hexahedra[j][r];
       }
+      cell_type_vector[jh] = CellType::Hexahedron;
     }
 
-    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector));
+    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector,
+                                                                      cell_type_vector));
     Connectivity3D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -851,6 +857,8 @@ GmshReader::__proceedData()
   } else if ((dimension2_mask, elementNumber)>0) {
     const size_t nb_cells = (dimension2_mask, elementNumber);
 
+    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);
     for (size_t j=0; j<nb_triangles; ++j) {
@@ -858,6 +866,7 @@ GmshReader::__proceedData()
       for (int r=0; r<3; ++r) {
         cell_by_node_vector[j][r] = __triangles[j][r];
       }
+      cell_type_vector[j] = CellType::Triangle;
     }
 
     const size_t nb_quadrangles = __quadrangles.extent(0);
@@ -867,9 +876,11 @@ GmshReader::__proceedData()
       for (int r=0; r<4; ++r) {
         cell_by_node_vector[jq][r] = __quadrangles[j][r];
       }
+      cell_type_vector[jq] = CellType::Quadrangle;
     }
 
-    std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector));
+    std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector,
+                                                                      cell_type_vector));
     Connectivity2D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -915,18 +926,21 @@ GmshReader::__proceedData()
     m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
 
   } else if ((dimension1_mask, elementNumber)>0) {
-
     const size_t nb_cells = (dimension1_mask, elementNumber);
 
+    std::vector<CellType> cell_type_vector(nb_cells);
+
     std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
     for (size_t j=0; j<nb_cells; ++j) {
       cell_by_node_vector[j].resize(2);
       for (int r=0; r<2; ++r) {
         cell_by_node_vector[j][r] = __edges[j][r];
       }
+      cell_type_vector[j] = CellType::Line;
     }
 
-    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector));
+    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector,
+                                                                      cell_type_vector));
     Connectivity1D& connectivity = *p_connectivity;
 
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;