diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 01ad632cf2fbfd013d0da8bc7e39b7b63f7c4129..cda9cb9b5f60c9f1ca74d2346120e4198e958d52 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) {
@@ -949,6 +950,10 @@ struct physical_section : seq<TAO_PEGTL_STRING("$PhysicalNames"), until<TAO_PEGT
 {
 };
 
+struct entities_section : seq<TAO_PEGTL_STRING("$Entities"), until<TAO_PEGTL_STRING("$EndEntities")>, ignored>
+{
+};
+
 struct nodes_section : seq<TAO_PEGTL_STRING("$Nodes"), until<TAO_PEGTL_STRING("$EndNodes")>, ignored>
 {
 };
@@ -981,6 +986,7 @@ struct interpolation_scheme_section
 
 struct data_section
   : sor<physical_section,
+        entities_section,
         nodes_section,
         elements_section,
         periodic_section,
@@ -1039,7 +1045,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);
 
@@ -1078,9 +1085,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";
@@ -1118,9 +1127,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
@@ -1133,7 +1144,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;
@@ -1181,6 +1193,7 @@ inline auto tokenize_variable = [](std::string_view& line, std::vector<std::stri
 struct State
 {
   bool has_physiscal_section            = false;
+  bool has_entities_section             = false;
   bool has_nodes_section                = false;
   bool has_elements_section             = false;
   bool has_periodic_section             = false;
@@ -1231,7 +1244,7 @@ struct actions<meshformat>
       }
       std::cout << "- Reading Gmsh format " << gmsh_data.format.version << '\n';
 
-      const std::set<double> supported_versions = {2.2};
+      const std::set<double> supported_versions = {2.2, 4.1};
 
       if (not supported_versions.contains(gmsh_data.format.version)) {
         std::stringstream error_msg;
@@ -1270,32 +1283,139 @@ struct actions<physical_section>
       state.has_physiscal_section = true;
     }
 
+    if ((gmsh_data.format.version == 2.2) or (gmsh_data.format.version == 4.1)) {
+      std::string_view body = in.string_view();
+
+      auto pos = body.find('\n');
+
+      if (pos == std::string_view::npos) {
+        throw NormalError("malformed Gmsh file, $PhysicalNames section is empty");
+      }
+      body = body.substr(pos + 1);
+
+      size_t number_of_physical_names;
+      try {
+        number_of_physical_names = parse_int(remove_embedding_spaces(get_next_line(body)));
+      }
+      catch (std::runtime_error& e) {
+        std::stringstream error_msg;
+        error_msg << "reading " << in.position().source << ':' << in.position().line + 1 << ": " << e.what();
+        throw NormalError(error_msg.str());
+      }
+
+      size_t i_physical_name = 0;
+      try {
+        for (; i_physical_name < number_of_physical_names; ++i_physical_name) {
+          if (body.empty()) {
+            std::stringstream error_msg;
+            error_msg << "malformed Gmsh file, $PhysicalNames section contains " << i_physical_name
+                      << " entries, expecting " << number_of_physical_names;
+            throw std::runtime_error(error_msg.str());
+          }
+          auto line = remove_embedding_spaces(get_next_line(body));
+
+          if (line.empty()) {
+            throw std::runtime_error("empty entry");
+          }
+
+          std::array<std::string_view, 3> tokens{};
+          size_t number_of_tokens = tokenize(line, tokens);
+
+          if (number_of_tokens != 3) {
+            throw std::runtime_error("expecting 3 entries");
+          }
+
+          if ((tokens[2][0] != '"') or (tokens[2][tokens[2].size() - 1] != '"')) {
+            throw std::runtime_error("physical name must be enclosed by double quotes");
+          }
+          if (tokens[2].size() == 2) {
+            throw std::runtime_error("physical name cannot be an empty string");
+          }
+
+          const size_t physical_dimension = parse_int(tokens[0]);
+          const size_t physical_number    = parse_int(tokens[1]);
+          const std::string physical_name{tokens[2].substr(1, tokens[2].size() - 2)};
+
+          GmshReader::PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
+
+          if (auto i_searched_physical_ref_id = gmsh_data.m_physical_ref_map.find(physical_number);
+              i_searched_physical_ref_id != gmsh_data.m_physical_ref_map.end()) {
+            std::stringstream os;
+            os << "Physical reference id '" << physical_ref_id << "' already defined as '"
+               << i_searched_physical_ref_id->second << "'!";
+            throw std::runtime_error(os.str());
+          }
+          gmsh_data.m_physical_ref_map[physical_number] = physical_ref_id;
+        }
+      }
+      catch (std::runtime_error& e) {
+        std::stringstream error_msg;
+        error_msg << "reading " << in.position().source << ':' << in.position().line + i_physical_name + 2 << ": "
+                  << e.what();
+        throw NormalError(error_msg.str());
+      }
+    } else {
+      throw UnexpectedError("invalid file format version: " + std::to_string(gmsh_data.format.version));
+    }
+  }
+};
+
+template <>
+struct actions<entities_section>
+{
+  template <typename Input>
+  static void
+  apply([[maybe_unused]] const Input& in, State& state, [[maybe_unused]] GmshReader::GmshData& gmsh_data)
+  {
+    if (state.has_entities_section) {
+      throw NormalError("malformed Gmsh file, $Entities section should appear only once");
+    } else {
+      state.has_entities_section = true;
+    }
+  }
+};
+
+template <>
+struct actions<nodes_section>
+{
+  template <typename Input>
+  static void
+  apply_2_2(const Input& in, GmshReader::GmshData& gmsh_data)
+  {
     std::string_view body = in.string_view();
 
     auto pos = body.find('\n');
 
     if (pos == std::string_view::npos) {
-      throw NormalError("malformed Gmsh file, $PhysicalNames section is empty");
+      throw NormalError("malformed Gmsh file, $Nodes section is empty");
     }
+
     body = body.substr(pos + 1);
 
-    size_t number_of_physical_names;
+    int nb_nodes;
     try {
-      number_of_physical_names = parse_int(remove_embedding_spaces(get_next_line(body)));
+      nb_nodes = parse_int(remove_embedding_spaces(get_next_line(body)));
+      if (nb_nodes < 0) {
+        throw std::runtime_error("number of vertices is negative");
+      }
     }
     catch (std::runtime_error& e) {
       std::stringstream error_msg;
       error_msg << "reading " << in.position().source << ':' << in.position().line + 1 << ": " << e.what();
       throw NormalError(error_msg.str());
     }
+    const size_t number_of_nodes = nb_nodes;
+    std::cout << "- Number Of Vertices: " << number_of_nodes << '\n';
 
-    size_t i_physical_name = 0;
+    Array<int> vertex_number_list{number_of_nodes};
+    Array<TinyVector<3>> vertex_list{number_of_nodes};
+
+    size_t i_node = 0;
     try {
-      for (; i_physical_name < number_of_physical_names; ++i_physical_name) {
+      for (; i_node < number_of_nodes; ++i_node) {
         if (body.empty()) {
           std::stringstream error_msg;
-          error_msg << "malformed Gmsh file, $PhysicalNames section contains " << i_physical_name
-                    << " entries, expecting " << number_of_physical_names;
+          error_msg << "$Nodes section contains " << i_node << " entries, expecting " << number_of_nodes;
           throw std::runtime_error(error_msg.str());
         }
         auto line = remove_embedding_spaces(get_next_line(body));
@@ -1304,58 +1424,35 @@ struct actions<physical_section>
           throw std::runtime_error("empty entry");
         }
 
-        std::array<std::string_view, 3> tokens{};
-        size_t number_of_tokens = tokenize(line, tokens);
-
-        if (number_of_tokens != 3) {
-          throw std::runtime_error("expecting 3 entries");
-        }
+        std::array<std::string_view, 4> tokens{};
+        size_t token_number = tokenize(line, tokens);
 
-        if ((tokens[2][0] != '"') or (tokens[2][tokens[2].size() - 1] != '"')) {
-          throw std::runtime_error("physical name must be enclosed by double quotes");
-        }
-        if (tokens[2].size() == 2) {
-          throw std::runtime_error("physical name cannot be an empty string");
+        if (token_number != 4) {
+          std::stringstream error_msg;
+          error_msg << "malformed Gmsh file, $Nodes section contains " << i_node << " entries, expecting "
+                    << number_of_nodes;
+          throw std::runtime_error(error_msg.str());
         }
 
-        const size_t physical_dimension = parse_int(tokens[0]);
-        const size_t physical_number    = parse_int(tokens[1]);
-        const std::string physical_name{tokens[2].substr(1, tokens[2].size() - 2)};
-
-        GmshReader::PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
+        vertex_number_list[i_node] = parse_int(tokens[0]);
 
-        if (auto i_searched_physical_ref_id = gmsh_data.m_physical_ref_map.find(physical_number);
-            i_searched_physical_ref_id != gmsh_data.m_physical_ref_map.end()) {
-          std::stringstream os;
-          os << "Physical reference id '" << physical_ref_id << "' already defined as '"
-             << i_searched_physical_ref_id->second << "'!";
-          throw std::runtime_error(os.str());
-        }
-        gmsh_data.m_physical_ref_map[physical_number] = physical_ref_id;
+        vertex_list[i_node] = TinyVector<3>{parse_double(tokens[1]), parse_double(tokens[2]), parse_double(tokens[3])};
       }
     }
     catch (std::runtime_error& e) {
       std::stringstream error_msg;
-      error_msg << "reading " << in.position().source << ':' << in.position().line + i_physical_name + 2 << ": "
-                << e.what();
+      error_msg << "reading " << in.position().source << ':' << in.position().line + i_node + 2 << ": " << e.what();
       throw NormalError(error_msg.str());
     }
+
+    gmsh_data.node_number     = vertex_number_list;
+    gmsh_data.node_coordinate = vertex_list;
   }
-};
 
-template <>
-struct actions<nodes_section>
-{
   template <typename Input>
   static void
-  apply(const Input& in, State& state, GmshReader::GmshData& gmsh_data)
+  apply_4_1(const Input& in, GmshReader::GmshData& gmsh_data)
   {
-    if (state.has_nodes_section) {
-      throw NormalError("malformed Gmsh file, $Nodes section should appear only once");
-    } else {
-      state.has_nodes_section = true;
-    }
-
     std::string_view body = in.string_view();
 
     auto pos = body.find('\n');
@@ -1368,10 +1465,25 @@ struct actions<nodes_section>
 
     int nb_nodes;
     try {
-      nb_nodes = parse_int(remove_embedding_spaces(get_next_line(body)));
+      auto line = remove_embedding_spaces(get_next_line(body));
+
+      std::array<std::string_view, 4> tokens{};
+      size_t token_number = tokenize(line, tokens);
+
+      if (token_number != 4) {
+        throw std::runtime_error("malformed $Nodes descriptor expecting 4 values");
+      }
+
+      // size_t nb_entity_blocks =
+      parse_int(tokens[0]);
+      nb_nodes = parse_int(tokens[1]);
       if (nb_nodes < 0) {
         throw std::runtime_error("number of vertices is negative");
       }
+      // size_t min_node_tag =
+      parse_int(tokens[2]);
+      // size_t max_node_tag =
+      parse_int(tokens[3]);
     }
     catch (std::runtime_error& e) {
       std::stringstream error_msg;
@@ -1422,6 +1534,25 @@ struct actions<nodes_section>
     gmsh_data.node_number     = vertex_number_list;
     gmsh_data.node_coordinate = vertex_list;
   }
+
+  template <typename Input>
+  static void
+  apply(const Input& in, State& state, GmshReader::GmshData& gmsh_data)
+  {
+    if (state.has_nodes_section) {
+      throw NormalError("malformed Gmsh file, $Nodes section should appear only once");
+    } else {
+      state.has_nodes_section = true;
+    }
+
+    if (gmsh_data.format.version == 2.2) {
+      apply_2_2(in, gmsh_data);
+    } else if (gmsh_data.format.version == 4.1) {
+      apply_4_1(in, gmsh_data);
+    } else {
+      throw UnexpectedError("invalid file format " + std::to_string(gmsh_data.format.version));
+    }
+  }
 };
 
 template <>