From c6891a465204bb20e5234a8a890cfbb9f9339b0e Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 8 Jul 2025 08:38:21 +0200
Subject: [PATCH] Prepare support for 4.1 file format

---
 src/mesh/GmshReader.cpp | 223 +++++++++++++++++++++++++++-------------
 1 file changed, 150 insertions(+), 73 deletions(-)

diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 904804c3a..01ad632cf 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -995,7 +995,7 @@ struct grammar : must<ignored, meshformat, star<data_section>>
 {
 };
 
-auto parse_double = [](std::string_view sv) {
+inline auto parse_double = [](std::string_view sv) {
   double value;
   auto [p1, ec1] = std::from_chars(sv.data(), sv.data() + sv.size(), value);
   if ((ec1 != std::errc()) or (p1 != sv.data() + sv.size())) {
@@ -1004,7 +1004,7 @@ auto parse_double = [](std::string_view sv) {
   return value;
 };
 
-auto parse_int = [](std::string_view sv) {
+inline auto parse_int = [](std::string_view sv) {
   int integer;
   auto [p1, ec1] = std::from_chars(sv.data(), sv.data() + sv.size(), integer);
   if ((ec1 != std::errc()) or (p1 != sv.data() + sv.size())) {
@@ -1013,7 +1013,7 @@ auto parse_int = [](std::string_view sv) {
   return integer;
 };
 
-auto get_next_line = [](std::string_view& body) -> std::string_view {
+inline auto get_next_line = [](std::string_view& body) -> std::string_view {
   std::string_view line;
   auto nl = body.find('\n');
   if (nl != std::string_view::npos) {
@@ -1026,7 +1026,7 @@ auto get_next_line = [](std::string_view& body) -> std::string_view {
   return line;
 };
 
-auto remove_embedding_spaces = [](std::string_view line) -> std::string_view {
+inline auto remove_embedding_spaces = [](std::string_view line) -> std::string_view {
   auto first = line.find_first_not_of(" \t\r");
   if (first == std::string_view::npos) {
     line = {};
@@ -1037,46 +1037,105 @@ auto remove_embedding_spaces = [](std::string_view line) -> std::string_view {
   return line;
 };
 
-constexpr std::array<size_t, 16> primitive_number_of_nodes = {
-  2,    // edge
-  3,    // triangle
-  4,    // quadrangle
-  4,    // tetrahedron
-  8,    // hexaredron
-  6,    // prism
-  5,    // pyramid
-  3,    // second order edge
-  6,    // second order triangle
-  9,    // second order quadrangle
-  10,   // second order tetrahedron
-  27,   // second order hexahedron
-  18,   // second order prism
-  14,   // second order pyramid
-  1     // point
-};
-
-constexpr std::array<std::string_view, 16> element_name = {
-  "edges",
-  "triangles",
-  "quadrangles",
-  "tetrahedra",
-  "hexahedra",
-  "prisms",
-  "pyramids",
-  "second order edges",
-  "second order triangles",
-  "second order quadrangles",
-  "second order tetrahedra",
-  "second order hexahedra",
-  "second order prisms",
-  "second order pyramids",
-  "points"   //
-};
-
-constexpr std::array<bool, 16> supported_element   //
-  = {true, true, true, true, true, true, true, false, false, false, false, false, false, false, true};
-
-auto tokenize = [](std::string_view& line, auto& tokens) {
+inline constexpr int number_of_gmsh_entities = 93;
+
+inline constexpr auto primitive_number_of_nodes = []() constexpr {
+  std::array<short, number_of_gmsh_entities> number_of_nodes;
+  number_of_nodes.fill(-1);
+
+  number_of_nodes[0]  = 2;     // edge
+  number_of_nodes[1]  = 3;     // triangle
+  number_of_nodes[2]  = 4;     // quadrangle
+  number_of_nodes[3]  = 4;     // tetrahedron
+  number_of_nodes[4]  = 8;     // hexaredron
+  number_of_nodes[5]  = 6;     // prism
+  number_of_nodes[6]  = 5;     // pyramid
+  number_of_nodes[7]  = 3;     // second order edge
+  number_of_nodes[8]  = 6;     // second order triangle
+  number_of_nodes[9]  = 9;     // second order quadrangle
+  number_of_nodes[10] = 10;    // second order tetrahedron
+  number_of_nodes[11] = 27;    // second order hexahedron
+  number_of_nodes[12] = 18;    // second order prism
+  number_of_nodes[13] = 14;    // second order pyramid
+  number_of_nodes[14] = 1;     // point
+  number_of_nodes[15] = 8;     // 8-node second order quadrangle
+  number_of_nodes[16] = 20;    // 20-node second order hexahedron
+  number_of_nodes[17] = 15;    // 15-node second order prism
+  number_of_nodes[18] = 13;    // 13-node second order pyramid
+  number_of_nodes[19] = 9;     // 9-node third order incomplete triangle
+  number_of_nodes[20] = 10;    // 10-node third order triangle
+  number_of_nodes[21] = 12;    // 12-node fourth order incomplete triangle
+  number_of_nodes[22] = 15;    // 15-node fourth order triangle
+  number_of_nodes[33] = 15;    // 15-node fifth order incomplete triangle
+  number_of_nodes[24] = 21;    // 21-node fifth order complete triangle
+  number_of_nodes[25] = 4;     // 4-node third order edge
+  number_of_nodes[26] = 5;     // 5-node fourth order edge
+  number_of_nodes[27] = 6;     // 6-node fifth order edge
+  number_of_nodes[28] = 20;    // 20-node third order tetrahedron
+  number_of_nodes[29] = 35;    // 35-node fourth order tetrahedron
+  number_of_nodes[30] = 56;    // 56-node fifth order tetrahedron
+  number_of_nodes[91] = 64;    // 64-node third order hexahedron
+  number_of_nodes[92] = 125;   // 125-node fourth order hexahedron
+
+  return number_of_nodes;
+}();
+
+inline constexpr auto entity_name = []() constexpr {
+  std::array<std::string_view, number_of_gmsh_entities> name;
+  name.fill("unknown");
+  name[0]  = "edges";
+  name[1]  = "triangles";
+  name[2]  = "quadrangles";
+  name[3]  = "tetrahedra";
+  name[4]  = "hexahedra";
+  name[5]  = "prisms";
+  name[6]  = "pyramids";
+  name[7]  = "second order edges";
+  name[8]  = "second order triangles";
+  name[9]  = "second order quadrangles";
+  name[10] = "second order tetrahedra";
+  name[11] = "second order hexahedra";
+  name[12] = "second order prisms";
+  name[13] = "second order pyramids";
+  name[14] = "points";
+  name[15] = "8-node second order quadrangle";
+  name[16] = "20-node second order hexahedron";
+  name[17] = "15-node second order prism";
+  name[18] = "13-node second order pyramid";
+  name[19] = "9-node third order incomplete triangle";
+  name[20] = "10-node third order triangle";
+  name[21] = "12-node fourth order incomplete triangle";
+  name[22] = "15-node fourth order triangle";
+  name[33] = "15-node fifth order incomplete triangle";
+  name[24] = "21-node fifth order complete triangle";
+  name[25] = "4-node third order edge";
+  name[26] = "5-node fourth order edge";
+  name[27] = "6-node fifth order edge";
+  name[28] = "20-node third order tetrahedron";
+  name[29] = "35-node fourth order tetrahedron";
+  name[30] = "56-node fifth order tetrahedron";
+  name[91] = "64-node third order hexahedron";
+  name[92] = "125-node fourth order hexahedron";
+
+  return name;
+}();
+
+inline constexpr auto supported_entity = []() constexpr {
+  std::array<bool, number_of_gmsh_entities> supported;
+  supported.fill(false);
+  supported[0]  = true;   // edge
+  supported[1]  = true;   // triangle
+  supported[2]  = true;   // quadrangle
+  supported[3]  = true;   // tetrahedron
+  supported[4]  = true;   // hexaredron
+  supported[5]  = true;   // prism
+  supported[6]  = true;   // pyramid
+  supported[14] = true;   // point
+
+  return supported;
+}();
+
+inline auto tokenize = [](std::string_view& line, auto& tokens) {
   size_t token_idx = 0;
   while (not line.empty() && token_idx < tokens.size()) {
     auto space = line.find_first_of(" \t\r");
@@ -1098,7 +1157,7 @@ auto tokenize = [](std::string_view& line, auto& tokens) {
   return token_idx;
 };
 
-auto tokenize_variable = [](std::string_view& line, std::vector<std::string_view>& tokens) {
+inline auto tokenize_variable = [](std::string_view& line, std::vector<std::string_view>& tokens) {
   tokens.clear();
   while (not line.empty()) {
     auto space = line.find_first_of(" \t\r");
@@ -1172,6 +1231,23 @@ struct actions<meshformat>
       }
       std::cout << "- Reading Gmsh format " << gmsh_data.format.version << '\n';
 
+      const std::set<double> supported_versions = {2.2};
+
+      if (not supported_versions.contains(gmsh_data.format.version)) {
+        std::stringstream error_msg;
+        error_msg << "Gmsh file format " << gmsh_data.format.version << " is not supported. Supported versions: ";
+        bool first = true;
+        for (auto version : supported_versions) {
+          if (not first) {
+            error_msg << ", ";
+          } else {
+            first = false;
+          }
+          error_msg << version;
+        }
+        throw NormalError(error_msg.str());
+      }
+
       if (gmsh_data.format.binary) {
         throw NotImplementedError("binary format is not supported yet");
       }
@@ -1417,17 +1493,17 @@ struct actions<elements_section>
         gmsh_data.entity_number[i_element] = parse_int(tokens[0]);
 
         const int entity_type = parse_int(tokens[1]) - 1;
-        if ((entity_type < 0) or (entity_type > 15)) {
-          throw std::runtime_error("unknown element type");
+        if ((entity_type < 0) or (entity_type >= number_of_gmsh_entities)) {
+          throw std::runtime_error("unknown entity type");
         }
-        if (not supported_element[entity_type]) {
-          throw std::runtime_error(std::string{element_name[entity_type]} + " is not supported");
+        if (not supported_entity[entity_type]) {
+          throw std::runtime_error(std::string{entity_name[entity_type]} + " is not supported");
         }
         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{element_name[entity_type]} + " number of tags is negative");
+          throw std::runtime_error(std::string{entity_name[entity_type]} + " number of tags is negative");
         }
 
         const size_t number_of_tags           = nb_tags;
@@ -1539,7 +1615,6 @@ struct actions<interpolation_scheme_section>
     std::cout << "  [" << rang::fg::yellow << "ignoring" << rang::fg::reset << " $InterpolationScheme section]\n";
   }
 };
-
 }   // namespace gmsh_parser
 
 GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
@@ -1617,8 +1692,10 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 void
 GmshReader::__proceedData()
 {
-  TinyVector<15, int> element_number = zero;
-  std::vector<std::set<size_t>> element_reference(15);
+  constexpr int number_of_supported_entities = gmsh_parser::number_of_gmsh_entities;
+
+  TinyVector<number_of_supported_entities, int> element_number = zero;
+  std::vector<std::set<size_t>> element_reference(number_of_supported_entities);
 
   for (size_t i = 0; i < m_mesh_data.entity_type.size(); ++i) {
     const int elementType = m_mesh_data.entity_type[i];
@@ -1628,7 +1705,7 @@ GmshReader::__proceedData()
 
   for (size_t i = 0; i < element_number.dimension(); ++i) {
     if (element_number[i] > 0) {
-      if (not(gmsh_parser::supported_element[i])) {
+      if (not(gmsh_parser::supported_entity[i])) {
       } else {
         switch (i) {
         // Supported entities
@@ -1842,25 +1919,25 @@ GmshReader::__proceedData()
     }
   }
 
-  TinyVector<15, int> dimension0_mask = zero;
-  dimension0_mask[14]                 = 1;
-  TinyVector<15, int> dimension1_mask = zero;
-  dimension1_mask[0]                  = 1;
-  dimension1_mask[7]                  = 1;
-  TinyVector<15, int> dimension2_mask = zero;
-  dimension2_mask[1]                  = 1;
-  dimension2_mask[2]                  = 1;
-  dimension2_mask[8]                  = 1;
-  dimension2_mask[9]                  = 1;
-  TinyVector<15, int> dimension3_mask = zero;
-  dimension3_mask[3]                  = 1;
-  dimension3_mask[4]                  = 1;
-  dimension3_mask[5]                  = 1;
-  dimension3_mask[6]                  = 1;
-  dimension3_mask[10]                 = 1;
-  dimension3_mask[11]                 = 1;
-  dimension3_mask[12]                 = 1;
-  dimension3_mask[13]                 = 1;
+  TinyVector<number_of_supported_entities, int> dimension0_mask = zero;
+  dimension0_mask[14]                                           = 1;
+  TinyVector<number_of_supported_entities, int> dimension1_mask = zero;
+  dimension1_mask[0]                                            = 1;
+  dimension1_mask[7]                                            = 1;
+  TinyVector<number_of_supported_entities, int> dimension2_mask = zero;
+  dimension2_mask[1]                                            = 1;
+  dimension2_mask[2]                                            = 1;
+  dimension2_mask[8]                                            = 1;
+  dimension2_mask[9]                                            = 1;
+  TinyVector<number_of_supported_entities, int> dimension3_mask = zero;
+  dimension3_mask[3]                                            = 1;
+  dimension3_mask[4]                                            = 1;
+  dimension3_mask[5]                                            = 1;
+  dimension3_mask[6]                                            = 1;
+  dimension3_mask[10]                                           = 1;
+  dimension3_mask[11]                                           = 1;
+  dimension3_mask[12]                                           = 1;
+  dimension3_mask[13]                                           = 1;
 
   if (dot(dimension3_mask, element_number) > 0) {
     const size_t nb_cells = dot(dimension3_mask, element_number);
-- 
GitLab