diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index c774c6a82f039cdf89cc11cc03362bdcee4cd1c6..d777c66750bad343a970dc9e8148437b4bd1491a 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -39,7 +39,8 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
   Array<int> cell_number_vector(nb_cells);
 
   Array<unsigned int> cell_to_node_row(nb_cells + 1);
-  parallel_for(cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
+  parallel_for(
+    cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
 
   Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
   for (CellId cell_id = 0; cell_id < nb_cells; ++cell_id) {
@@ -62,7 +63,10 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (unsigned int r = 0; r < gmsh_data.point_entity.size(); ++r) {
     const unsigned int point_number = gmsh_data.point_entity[r];
-    ref_points_map[gmsh_data.point_entity_ref[r]].push_back(point_number);
+    auto ref_range                  = point_to_physical_tag.equal_range(gmsh_data.point_entity_ref[r]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_points_map[i_ref->second].push_back(point_number);
+    }
   }
 
   Array<size_t> node_nb_cell(descriptor.nodeNumberVector().size());
@@ -125,7 +129,10 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (unsigned int j = 0; j < gmsh_data.edge_entity_ref.size(); ++j) {
     const unsigned int elem_number = j;
-    ref_cells_map[gmsh_data.edge_entity_ref[j]].push_back(elem_number);
+    auto ref_range                 = cell_to_physical_tag.equal_range(gmsh_data.edge_entity_ref[j]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_cells_map[i_ref->second].push_back(elem_number);
+    }
   }
 
   for (const auto& ref_cell_list : ref_cells_map) {
@@ -295,7 +302,10 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   for (unsigned int e = 0; e < gmsh_data.edge_entity.size(); ++e) {
     const unsigned int edge_id = find_face(gmsh_data.edge_entity[e][0], gmsh_data.edge_entity[e][1]);
 
-    ref_faces_map[gmsh_data.edge_entity_ref[e]].push_back(edge_id);
+    auto ref_range = face_to_physical_tag.equal_range(gmsh_data.edge_entity_ref[e]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_faces_map[i_ref->second].push_back(edge_id);
+    }
 
     if (face_number_vector[edge_id] != gmsh_data.edge_entity_number[e]) {
       if (auto i_face = face_number_id_map.find(gmsh_data.edge_entity_number[e]); i_face != face_number_id_map.end()) {
@@ -388,7 +398,10 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (unsigned int r = 0; r < gmsh_data.point_entity.size(); ++r) {
     const unsigned int point_number = gmsh_data.point_entity[r];
-    ref_points_map[gmsh_data.point_entity_ref[r]].push_back(point_number);
+    auto ref_range                  = point_to_physical_tag.equal_range(gmsh_data.point_entity_ref[r]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_points_map[i_ref->second].push_back(point_number);
+    }
   }
 
   for (const auto& ref_point_list : ref_points_map) {
@@ -548,22 +561,34 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (unsigned int j = 0; j < gmsh_data.tetrahedron_entity_ref.size(); ++j) {
     const unsigned int elem_number = j;
-    ref_cells_map[gmsh_data.tetrahedron_entity_ref[j]].push_back(elem_number);
+    auto ref_range                 = cell_to_physical_tag.equal_range(gmsh_data.tetrahedron_entity_ref[j]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_cells_map[i_ref->second].push_back(elem_number);
+    }
   }
 
   for (unsigned int j = 0; j < gmsh_data.hexahedron_entity_ref.size(); ++j) {
     const size_t elem_number = nb_tetrahedra + j;
-    ref_cells_map[gmsh_data.hexahedron_entity_ref[j]].push_back(elem_number);
+    auto ref_range           = cell_to_physical_tag.equal_range(gmsh_data.hexahedron_entity_ref[j]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_cells_map[i_ref->second].push_back(elem_number);
+    }
   }
 
   for (unsigned int j = 0; j < gmsh_data.prism_entity_ref.size(); ++j) {
     const size_t elem_number = nb_tetrahedra + nb_hexahedra + j;
-    ref_cells_map[gmsh_data.prism_entity_ref[j]].push_back(elem_number);
+    auto ref_range           = cell_to_physical_tag.equal_range(gmsh_data.prism_entity_ref[j]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_cells_map[i_ref->second].push_back(elem_number);
+    }
   }
 
   for (unsigned int j = 0; j < gmsh_data.pyramid_entity_ref.size(); ++j) {
     const size_t elem_number = nb_tetrahedra + nb_hexahedra + nb_prisms + j;
-    ref_cells_map[gmsh_data.pyramid_entity_ref[j]].push_back(elem_number);
+    auto ref_range           = cell_to_physical_tag.equal_range(gmsh_data.pyramid_entity_ref[j]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_cells_map[i_ref->second].push_back(elem_number);
+    }
   }
 
   for (const auto& ref_cell_list : ref_cells_map) {
@@ -676,7 +701,10 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       const unsigned int face_id =
         find_face({gmsh_data.triangle_entity[f][0], gmsh_data.triangle_entity[f][1], gmsh_data.triangle_entity[f][2]});
 
-      ref_faces_map[gmsh_data.triangle_entity_ref[f]].push_back(face_id);
+      auto ref_range = face_to_physical_tag.equal_range(gmsh_data.triangle_entity_ref[f]);
+      for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+        ref_faces_map[i_ref->second].push_back(face_id);
+      }
 
       if (face_number_vector[face_id] != gmsh_data.triangle_entity_number[f]) {
         if (auto i_face = face_number_id_map.find(gmsh_data.triangle_entity_number[f]);
@@ -701,7 +729,10 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       const unsigned int face_id = find_face({gmsh_data.quadrangle_entity[f][0], gmsh_data.quadrangle_entity[f][1],
                                               gmsh_data.quadrangle_entity[f][2], gmsh_data.quadrangle_entity[f][3]});
 
-      ref_faces_map[gmsh_data.quadrangle_entity_ref[f]].push_back(face_id);
+      auto ref_range = face_to_physical_tag.equal_range(gmsh_data.quadrangle_entity_ref[f]);
+      for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+        ref_faces_map[i_ref->second].push_back(face_id);
+      }
 
       if (face_number_vector[face_id] != gmsh_data.quadrangle_entity_number[f]) {
         if (auto i_face = face_number_id_map.find(gmsh_data.quadrangle_entity_number[f]);
@@ -824,7 +855,10 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (unsigned int e = 0; e < gmsh_data.edge_entity.size(); ++e) {
       const unsigned int edge_id = find_edge(gmsh_data.edge_entity[e][0], gmsh_data.edge_entity[e][1]);
 
-      ref_edges_map[gmsh_data.edge_entity_ref[e]].push_back(edge_id);
+      auto ref_range = edge_to_physical_tag.equal_range(gmsh_data.edge_entity_ref[e]);
+      for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+        ref_edges_map[i_ref->second].push_back(edge_id);
+      }
 
       if (edge_number_vector[edge_id] != gmsh_data.edge_entity_number[e]) {
         if (auto i_edge = edge_number_id_map.find(gmsh_data.edge_entity_number[e]);
@@ -899,7 +933,11 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   for (unsigned int r = 0; r < gmsh_data.point_entity.size(); ++r) {
     const unsigned int point_number = gmsh_data.point_entity[r];
-    ref_points_map[gmsh_data.point_entity_ref[r]].push_back(point_number);
+
+    auto ref_range = point_to_physical_tag.equal_range(gmsh_data.point_entity_ref[r]);
+    for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+      ref_points_map[i_ref->second].push_back(point_number);
+    }
   }
 
   for (const auto& ref_point_list : ref_points_map) {
@@ -1080,7 +1118,8 @@ inline auto remove_embedding_spaces = [](std::string_view line) -> std::string_v
 
 inline constexpr int number_of_gmsh_entities = 93;
 
-inline constexpr auto primitive_number_of_nodes = []() constexpr {
+inline constexpr auto primitive_number_of_nodes = []() constexpr
+{
   std::array<short, number_of_gmsh_entities> number_of_nodes;
   number_of_nodes.fill(-1);
 
@@ -1119,9 +1158,11 @@ inline constexpr auto primitive_number_of_nodes = []() constexpr {
   number_of_nodes[92] = 125;   // 125-node fourth order hexahedron
 
   return number_of_nodes;
-}();
+}
+();
 
-inline constexpr auto entity_name = []() constexpr {
+inline constexpr auto entity_name = []() constexpr
+{
   std::array<std::string_view, number_of_gmsh_entities> name;
   name.fill("unknown");
   name[0]  = "edges";
@@ -1159,9 +1200,11 @@ inline constexpr auto entity_name = []() constexpr {
   name[92] = "125-node fourth order hexahedron";
 
   return name;
-}();
+}
+();
 
-inline constexpr auto supported_entity = []() constexpr {
+inline constexpr auto supported_entity = []() constexpr
+{
   std::array<bool, number_of_gmsh_entities> supported;
   supported.fill(false);
   supported[0]  = true;   // edge
@@ -1174,7 +1217,8 @@ inline constexpr auto supported_entity = []() constexpr {
   supported[14] = true;   // point
 
   return supported;
-}();
+}
+();
 
 inline auto tokenize = [](std::string_view& line, auto& tokens) {
   size_t token_idx = 0;