diff --git a/src/main.cpp b/src/main.cpp
index 12fe62ea61829affcfea2b2a1de7aad57fb02041..6f34e713ecba2de04440bb6146db0a3dc7aef6c4 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -223,10 +223,11 @@ int main(int argc, char *argv[])
           MeshDataType mesh_data(mesh);
 
           std::vector<BoundaryConditionHandler> bc_list;
-          { // quite dirty!
-            for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
-              ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
-              const RefId& ref = nodes_boundary.first;
+          {
+            for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfBoundaries(); ++i_boundary) {
+              const Boundary& boundary = mesh.connectivity().boundary(i_boundary);
+              const RefId& ref = boundary.refId();
+
               TinyVector<2> normal(0,0);
               if ((ref.tagName()== std::string("XMIN")) or (ref.tagName()=="XMAX")) {
                 normal = TinyVector<2>(1,0);
@@ -237,7 +238,7 @@ int main(int argc, char *argv[])
                 break;
               }
 
-              const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
+              const Kokkos::View<const unsigned int*> nodes_ids = boundary.nodeList();
               std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
               for (size_t r=0; r<nodes_ids.extent(0); ++r) {
                 node_boundary_vector[r] = nodes_ids[r];
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index 8b53668de644cefdd9f4419a14a48842a5f2682a..0aa15f6515ab375493d0158abb8c0d77d768813f 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -11,6 +11,7 @@
 #include <algorithm>
 
 #include <RefId.hpp>
+#include <Boundary.hpp>
 
 class EntityListManager
 {
@@ -90,6 +91,8 @@ public:
   typedef std::vector<NodesBoundary> NodesBoundaryList;
 
   EntityListManager m_entity_list_manager;
+  std::vector<Boundary> m_boundary_list;
+
 private:
   const size_t  m_number_of_cells;
   size_t m_number_of_faces;
@@ -252,6 +255,21 @@ private:
   }
 
  public:
+  void addBoundary(const Boundary& boundary)
+  {
+    m_boundary_list.push_back(boundary);
+  }
+
+  size_t numberOfBoundaries() const
+  {
+    return m_boundary_list.size();
+  }
+
+  const Boundary& boundary(const size_t& i) const
+  {
+    return m_boundary_list[i];
+  }
+
   const size_t& numberOfNodes() const
   {
     return m_number_of_nodes;
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index c695118da58212adea5c7cf301e72d975c6ccd7c..816ab2cefd4c0e3456cb9f291fe8ba4761a9b39b 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -8,6 +8,9 @@
 #include <Connectivity1D.hpp>
 #include <Connectivity2D.hpp>
 #include <Mesh.hpp>
+
+#include <Boundary.hpp>
+
 #include <map>
 #include <regex>
 
@@ -836,6 +839,8 @@ GmshReader::__proceedData()
       }
       const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
       connectivity.entityListManager().addFaceIdArray(physical_ref_id.refId(), face_list);
+      Boundary boundary(connectivity, physical_ref_id.refId(), face_list);
+      connectivity.addBoundary(boundary);
     }
 
     const Kokkos::View<const unsigned int**> face_nodes = connectivity.faceNodes();