diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index cda9cb9b5f60c9f1ca74d2346120e4198e958d52..0ff14c610c5bae3e9fc71507c54026ddafb42bc7 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -39,8 +39,7 @@ 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) {
@@ -1019,6 +1018,15 @@ inline auto parse_int = [](std::string_view sv) {
   return integer;
 };
 
+inline auto parse_size_t = [](std::string_view sv) {
+  size_t size;
+  auto [p1, ec1] = std::from_chars(sv.data(), sv.data() + sv.size(), size);
+  if ((ec1 != std::errc()) or (p1 != sv.data() + sv.size())) {
+    throw std::runtime_error("invalid unsigned integer: '" + std::string{sv} + "'");
+  }
+  return size;
+};
+
 inline auto get_next_line = [](std::string_view& body) -> std::string_view {
   std::string_view line;
   auto nl = body.find('\n');
@@ -1045,8 +1053,7 @@ 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);
 
@@ -1085,11 +1092,9 @@ 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";
@@ -1127,11 +1132,9 @@ 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
@@ -1144,8 +1147,7 @@ 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;
@@ -1392,19 +1394,7 @@ struct actions<nodes_section>
 
     body = body.substr(pos + 1);
 
-    int nb_nodes;
-    try {
-      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;
+    const size_t number_of_nodes = parse_size_t(remove_embedding_spaces(get_next_line(body)));
     std::cout << "- Number Of Vertices: " << number_of_nodes << '\n';
 
     Array<int> vertex_number_list{number_of_nodes};
@@ -1453,6 +1443,7 @@ struct actions<nodes_section>
   static void
   apply_4_1(const Input& in, GmshReader::GmshData& gmsh_data)
   {
+    size_t line_number    = in.position().line;
     std::string_view body = in.string_view();
 
     auto pos = body.find('\n');
@@ -1463,8 +1454,10 @@ struct actions<nodes_section>
 
     body = body.substr(pos + 1);
 
-    int nb_nodes;
+    size_t nb_nodes;
+    size_t nb_entity_blocks;
     try {
+      ++line_number;
       auto line = remove_embedding_spaces(get_next_line(body));
 
       std::array<std::string_view, 4> tokens{};
@@ -1474,60 +1467,74 @@ struct actions<nodes_section>
         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]);
+      nb_entity_blocks = parse_size_t(tokens[0]);
+      nb_nodes         = parse_size_t(tokens[1]);
+      parse_size_t(tokens[2]);   // min node tag
+      parse_size_t(tokens[3]);   // max node tag
     }
     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;
+    const size_t number_of_entity_blocks = nb_entity_blocks;
+    const size_t number_of_nodes         = nb_nodes;
     std::cout << "- Number Of Vertices: " << number_of_nodes << '\n';
 
     Array<int> vertex_number_list{number_of_nodes};
     Array<TinyVector<3>> vertex_list{number_of_nodes};
 
-    size_t i_node = 0;
+    size_t i_block = 0;
+    size_t i_node  = 0;
     try {
-      for (; i_node < number_of_nodes; ++i_node) {
+      for (; i_block < number_of_entity_blocks; ++i_block) {
         if (body.empty()) {
           std::stringstream error_msg;
-          error_msg << "$Nodes section contains " << i_node << " entries, expecting " << number_of_nodes;
+          error_msg << "$Nodes section contains " << i_block << " block entries, expecting " << number_of_entity_blocks;
           throw std::runtime_error(error_msg.str());
         }
+        ++line_number;
         auto line = remove_embedding_spaces(get_next_line(body));
 
         if (line.empty()) {
           throw std::runtime_error("empty entry");
         }
 
-        std::array<std::string_view, 4> tokens{};
-        size_t token_number = tokenize(line, tokens);
-
-        if (token_number != 4) {
+        std::array<std::string_view, 4> block_tokens{};
+        if (tokenize(line, block_tokens) != 4) {
           std::stringstream error_msg;
-          error_msg << "malformed Gmsh file, $Nodes section contains " << i_node << " entries, expecting "
-                    << number_of_nodes;
+          error_msg << "malformed Gmsh file, $Nodes section contains " << i_block << " entries, expecting "
+                    << number_of_entity_blocks;
           throw std::runtime_error(error_msg.str());
         }
 
-        vertex_number_list[i_node] = parse_int(tokens[0]);
+        parse_int(block_tokens[0]);   // entity dim
+        parse_int(block_tokens[1]);   // entity tag
+        parse_int(block_tokens[2]);   // parametric
+        const size_t number_of_block_nodes = parse_size_t(block_tokens[3]);
 
-        vertex_list[i_node] = TinyVector<3>{parse_double(tokens[1]), parse_double(tokens[2]), parse_double(tokens[3])};
+        for (size_t i_block_node = 0; i_block_node < number_of_block_nodes; ++i_block_node) {
+          ++line_number;
+          vertex_number_list[i_node + i_block_node] = parse_size_t(remove_embedding_spaces(get_next_line(body)));
+        }
+        for (size_t i_block_node = 0; i_block_node < number_of_block_nodes; ++i_block_node) {
+          ++line_number;
+          auto coord_line = remove_embedding_spaces(get_next_line(body));
+          std::array<std::string_view, 3> coord_tokens{};
+          if (tokenize(coord_line, coord_tokens) != 3) {
+            std::cout << "coord_line = " << coord_line << '\n';
+            throw std::runtime_error("expecting 3 coordinates");
+          }
+
+          vertex_list[i_node + i_block_node] =
+            TinyVector<3>{parse_double(coord_tokens[0]), parse_double(coord_tokens[1]), parse_double(coord_tokens[2])};
+        }
+        i_node += number_of_block_nodes;
       }
     }
     catch (std::runtime_error& e) {
       std::stringstream error_msg;
-      error_msg << "reading " << in.position().source << ':' << in.position().line + i_node + 2 << ": " << e.what();
+      error_msg << "reading " << in.position().source << ':' << line_number << ": " << e.what();
       throw NormalError(error_msg.str());
     }
 
@@ -1560,14 +1567,8 @@ struct actions<elements_section>
 {
   template <typename Input>
   static void
-  apply(const Input& in, State& state, GmshReader::GmshData& gmsh_data)
+  apply_2_2(const Input& in, GmshReader::GmshData& gmsh_data)
   {
-    if (state.has_elements_section) {
-      throw NormalError("malformed Gmsh file, $Elements section should appear only once");
-    } else {
-      state.has_elements_section = true;
-    }
-
     std::string_view body = in.string_view();
 
     auto pos = body.find('\n');
@@ -1578,20 +1579,7 @@ struct actions<elements_section>
 
     body = body.substr(pos + 1);
 
-    int nb_elements;
-    try {
-      nb_elements = parse_int(remove_embedding_spaces(get_next_line(body)));
-      if (nb_elements < 0) {
-        throw std::runtime_error("number of elements 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_elements = nb_elements;
+    const size_t number_of_elements = parse_size_t(remove_embedding_spaces(get_next_line(body)));
     std::cout << "- Number Of Elements: " << number_of_elements << '\n' << std::flush;
 
     gmsh_data.entity_number.resize(number_of_elements);
@@ -1632,12 +1620,7 @@ struct actions<elements_section>
         }
         gmsh_data.entity_type[i_element] = entity_type;
 
-        int nb_tags = parse_int(tokens[2]);
-        if (nb_tags < 0) {
-          throw std::runtime_error(std::string{entity_name[entity_type]} + " number of tags is negative");
-        }
-
-        const size_t number_of_tags           = nb_tags;
+        const size_t number_of_tags           = parse_size_t(tokens[2]);
         gmsh_data.entity_reference[i_element] = parse_int(tokens[3]);
         // In Gmsh file format 2.2 the first tag IS the physical tag
         // other are spurious construction element tags
@@ -1657,6 +1640,133 @@ struct actions<elements_section>
       throw NormalError(error_msg.str());
     }
   }
+
+  template <typename Input>
+  static void
+  apply_4_1(const Input& in, GmshReader::GmshData& gmsh_data)
+  {
+    size_t line_number = in.position().line;
+
+    std::string_view body = in.string_view();
+
+    auto pos = body.find('\n');
+
+    if (pos == std::string_view::npos) {
+      throw NormalError("malformed Gmsh file, $Elements section is empty");
+    }
+
+    body = body.substr(pos + 1);
+
+    size_t nb_elements;
+    size_t nb_entity_blocks;
+    try {
+      ++line_number;
+      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");
+      }
+
+      nb_entity_blocks = parse_size_t(tokens[0]);
+      nb_elements      = parse_size_t(tokens[1]);
+      parse_size_t(tokens[2]);   // min element tag
+      parse_size_t(tokens[3]);   // max element tag
+    }
+    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_entity_blocks = nb_entity_blocks;
+    const size_t number_of_elements      = nb_elements;
+
+    std::cout << "- Number Of Elements: " << number_of_elements << '\n' << std::flush;
+
+    gmsh_data.entity_number.resize(number_of_elements);
+    gmsh_data.entity_type.resize(number_of_elements);
+    gmsh_data.__elementVertices.resize(number_of_elements);
+
+    std::vector<std::string_view> element_tokens;
+    element_tokens.reserve(20);
+
+    try {
+      size_t i_element = 0;
+      for (size_t i_block_entity = 0; i_block_entity < number_of_entity_blocks; ++i_block_entity) {
+        if (body.empty()) {
+          std::stringstream error_msg;
+          error_msg << "$Elements section contains " << i_block_entity << " entries, expecting " << number_of_elements;
+          throw std::runtime_error(error_msg.str());
+        }
+
+        ++line_number;
+        auto line = remove_embedding_spaces(get_next_line(body));
+
+        if (line.empty()) {
+          throw std::runtime_error("empty entry");
+        }
+
+        std::array<std::string_view, 4> block_tokens{};
+        if (tokenize(line, block_tokens) != 4) {
+          std::stringstream error_msg;
+          error_msg << "malformed Gmsh file, $Elements section contains " << i_block_entity << " entries, expecting "
+                    << number_of_entity_blocks;
+          throw std::runtime_error(error_msg.str());
+        }
+
+        parse_int(block_tokens[0]);   // entity dim
+        parse_int(block_tokens[1]);   // entity tag
+        const size_t element_type             = parse_int(block_tokens[2]);
+        const size_t number_of_block_elements = parse_size_t(block_tokens[3]);
+
+        for (size_t i_block_element = 0; i_block_element < number_of_block_elements; ++i_block_element) {
+          gmsh_data.entity_type[i_element + i_block_element] = element_type;
+        }
+
+        const size_t number_of_nodes = primitive_number_of_nodes[element_type];
+        for (size_t i_block_element = 0; i_block_element < number_of_block_elements; ++i_block_element) {
+          gmsh_data.__elementVertices[i_element + i_block_element].resize(number_of_nodes);
+
+          ++line_number;
+          auto element_line = remove_embedding_spaces(get_next_line(body));
+          tokenize_variable(element_line, element_tokens);
+
+          gmsh_data.entity_number[i_element + i_block_element] = parse_int(element_tokens[0]);
+          for (size_t i_node = 0; i_node < number_of_nodes; ++i_node) {
+            gmsh_data.__elementVertices[i_element + i_block_element][i_node] = parse_int(element_tokens[1 + i_node]);
+          }
+        }
+
+        i_element += number_of_block_elements;
+      }
+    }
+    catch (std::runtime_error& e) {
+      std::stringstream error_msg;
+      error_msg << "reading " << in.position().source << ':' << line_number << ": " << e.what();
+      throw NormalError(error_msg.str());
+    }
+  }
+
+  template <typename Input>
+  static void
+  apply(const Input& in, State& state, GmshReader::GmshData& gmsh_data)
+  {
+    if (state.has_elements_section) {
+      throw NormalError("malformed Gmsh file, $Elements section should appear only once");
+    } else {
+      state.has_elements_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 <>
@@ -1911,7 +2021,6 @@ GmshReader::__proceedData()
     throw NotImplementedError("reading file '" + m_filename + "': negative vertices number");
   }
 
-#warning consider the case of very sparse node lists (max-min>3*nb nodes). Uses unordered_map?
   // A value of -1 means that the vertex is unknown
   m_mesh_data.gmsh_node_to_node_id_correspondance.resize(max_node_number + 1, -1);