diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index b1ae2a4d0afbd06f2a31f5c75f44405df175a7d8..1788ed4d6bb8d72acaa4176b4a12633c31d24a52 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -17,10 +17,8 @@
 #include <rang.hpp>
 
 #include <charconv>
-#include <fstream>
 #include <iostream>
 #include <map>
-#include <regex>
 #include <set>
 #include <sstream>
 #include <unordered_map>
@@ -37,7 +35,7 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
 {
   ConnectivityDescriptor descriptor;
 
-  descriptor.setNodeNumberVector(gmsh_data.__verticesNumbers);
+  descriptor.setNodeNumberVector(gmsh_data.node_number);
   Array<CellType> cell_type_vector(nb_cells);
   Array<int> cell_number_vector(nb_cells);
 
@@ -48,7 +46,7 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
   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) {
     for (int i_node = 0; i_node < 2; ++i_node) {
-      cell_to_node_list[2 * cell_id + i_node] = gmsh_data.__edges[cell_id][i_node];
+      cell_to_node_list[2 * cell_id + i_node] = gmsh_data.edge_entity[cell_id][i_node];
     }
   }
   descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
@@ -57,14 +55,14 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
   descriptor.setCellTypeVector(cell_type_vector);
 
   for (size_t j = 0; j < nb_cells; ++j) {
-    cell_number_vector[j] = gmsh_data.__edges_number[j];
+    cell_number_vector[j] = gmsh_data.edge_entity_number[j];
   }
   descriptor.setCellNumberVector(cell_number_vector);
 
   std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
-    const unsigned int point_number = gmsh_data.__points[r];
-    ref_points_map[gmsh_data.__points_ref[r]].push_back(point_number);
+  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);
   }
 
   Array<size_t> node_nb_cell(descriptor.nodeNumberVector().size());
@@ -123,9 +121,9 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
-  for (unsigned int j = 0; j < gmsh_data.__edges_ref.size(); ++j) {
+  for (unsigned int j = 0; j < gmsh_data.edge_entity_ref.size(); ++j) {
     const unsigned int elem_number = j;
-    ref_cells_map[gmsh_data.__edges_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.edge_entity_ref[j]].push_back(elem_number);
   }
 
   for (const auto& ref_cell_list : ref_cells_map) {
@@ -166,18 +164,18 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
 {
   ConnectivityDescriptor descriptor;
 
-  descriptor.setNodeNumberVector(gmsh_data.__verticesNumbers);
+  descriptor.setNodeNumberVector(gmsh_data.node_number);
   Array<CellType> cell_type_vector(nb_cells);
   Array<int> cell_number_vector(nb_cells);
 
-  Array<unsigned int> cell_to_node_row(gmsh_data.__triangles.size() + gmsh_data.__quadrangles.size() + 1);
+  Array<unsigned int> cell_to_node_row(gmsh_data.triangle_entity.size() + gmsh_data.quadrangle_entity.size() + 1);
   {
     cell_to_node_row[0] = 0;
     size_t i_cell       = 0;
-    for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle, ++i_cell) {
+    for (size_t i_triangle = 0; i_triangle < gmsh_data.triangle_entity.size(); ++i_triangle, ++i_cell) {
       cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 3;
     }
-    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle, ++i_cell) {
+    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.quadrangle_entity.size(); ++i_quadrangle, ++i_cell) {
       cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 4;
     }
   }
@@ -185,46 +183,46 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
   {
     size_t i_cell_node = 0;
-    for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][2];
+    for (size_t i_triangle = 0; i_triangle < gmsh_data.triangle_entity.size(); ++i_triangle) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.triangle_entity[i_triangle][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.triangle_entity[i_triangle][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.triangle_entity[i_triangle][2];
     }
-    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][2];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][3];
+    for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.quadrangle_entity.size(); ++i_quadrangle) {
+      cell_to_node_list[i_cell_node++] = gmsh_data.quadrangle_entity[i_quadrangle][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.quadrangle_entity[i_quadrangle][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.quadrangle_entity[i_quadrangle][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.quadrangle_entity[i_quadrangle][3];
     }
   }
   descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
 
-  const size_t nb_triangles = gmsh_data.__triangles.size();
+  const size_t nb_triangles = gmsh_data.triangle_entity.size();
   for (size_t i_triangle = 0; i_triangle < nb_triangles; ++i_triangle) {
     cell_type_vector[i_triangle]   = CellType::Triangle;
-    cell_number_vector[i_triangle] = gmsh_data.__triangles_number[i_triangle];
+    cell_number_vector[i_triangle] = gmsh_data.triangle_entity_number[i_triangle];
   }
 
-  const size_t nb_quadrangles = gmsh_data.__quadrangles.size();
+  const size_t nb_quadrangles = gmsh_data.quadrangle_entity.size();
   for (size_t i_quadrangle = 0; i_quadrangle < nb_quadrangles; ++i_quadrangle) {
     const size_t i_cell = i_quadrangle + nb_triangles;
 
     cell_type_vector[i_cell]   = CellType::Quadrangle;
-    cell_number_vector[i_cell] = gmsh_data.__quadrangles_number[i_quadrangle];
+    cell_number_vector[i_cell] = gmsh_data.quadrangle_entity_number[i_quadrangle];
   }
 
   descriptor.setCellNumberVector(cell_number_vector);
   descriptor.setCellTypeVector(cell_type_vector);
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
-  for (unsigned int j = 0; j < gmsh_data.__triangles_ref.size(); ++j) {
+  for (unsigned int j = 0; j < gmsh_data.triangle_entity_ref.size(); ++j) {
     const unsigned int elem_number = j;
-    ref_cells_map[gmsh_data.__triangles_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.triangle_entity_ref[j]].push_back(elem_number);
   }
 
-  for (unsigned int j = 0; j < gmsh_data.__quadrangles_ref.size(); ++j) {
+  for (unsigned int j = 0; j < gmsh_data.quadrangle_entity_ref.size(); ++j) {
     const size_t elem_number = nb_triangles + j;
-    ref_cells_map[gmsh_data.__quadrangles_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.quadrangle_entity_ref[j]].push_back(elem_number);
   }
 
   for (const auto& ref_cell_list : ref_cells_map) {
@@ -257,12 +255,12 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     return mutable_face_number_id_map;
   }();
 
-  const auto& node_number_vector  = descriptor.nodeNumberVector();
+  const auto& node_number         = descriptor.nodeNumberVector();
   const auto& node_to_face_matrix = descriptor.nodeToFaceMatrix();
   const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
 
   const auto find_face = [&](uint32_t node0, uint32_t node1) {
-    if (node_number_vector[node0] > node_number_vector[node1]) {
+    if (node_number[node0] > node_number[node1]) {
       std::swap(node0, node1);
     }
     const auto& face_id_vector = node_to_face_matrix[node0];
@@ -281,13 +279,13 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   Array<int> face_number_vector = copy(descriptor.faceNumberVector());
   std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-  for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
-    const unsigned int edge_id = find_face(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+  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.__edges_ref[e]].push_back(edge_id);
+    ref_faces_map[gmsh_data.edge_entity_ref[e]].push_back(edge_id);
 
-    if (face_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
-      if (auto i_face = face_number_id_map.find(gmsh_data.__edges_number[e]); i_face != face_number_id_map.end()) {
+    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()) {
         const int other_edge_id = i_face->second;
         std::swap(face_number_vector[edge_id], face_number_vector[other_edge_id]);
 
@@ -298,7 +296,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
         face_number_id_map[face_number_vector[other_edge_id]] = other_edge_id;
       } else {
         face_number_id_map.erase(face_number_vector[edge_id]);
-        face_number_vector[edge_id]                     = gmsh_data.__edges_number[e];
+        face_number_vector[edge_id]                     = gmsh_data.edge_entity_number[e];
         face_number_id_map[face_number_vector[edge_id]] = edge_id;
       }
     }
@@ -360,7 +358,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     descriptor.addRefItemList(RefFaceList{ref_id, face_list, list_type});
   }
 
-  Array<bool> is_boundary_node(node_number_vector.size());
+  Array<bool> is_boundary_node(node_number.size());
   is_boundary_node.fill(false);
 
   for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
@@ -373,9 +371,9 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
-    const unsigned int point_number = gmsh_data.__points[r];
-    ref_points_map[gmsh_data.__points_ref[r]].push_back(point_number);
+  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);
   }
 
   for (const auto& ref_point_list : ref_points_map) {
@@ -424,7 +422,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
   }());
 
   descriptor.setNodeOwnerVector([&] {
-    Array<int> node_owner_vector(node_number_vector.size());
+    Array<int> node_owner_vector(node_number.size());
     node_owner_vector.fill(parallel::rank());
     return node_owner_vector;
   }());
@@ -437,14 +435,14 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 {
   ConnectivityDescriptor descriptor;
 
-  descriptor.setNodeNumberVector(gmsh_data.__verticesNumbers);
+  descriptor.setNodeNumberVector(gmsh_data.node_number);
   Array<CellType> cell_type_vector(nb_cells);
   Array<int> cell_number_vector(nb_cells);
 
-  const size_t nb_tetrahedra = gmsh_data.__tetrahedra.size();
-  const size_t nb_hexahedra  = gmsh_data.__hexahedra.size();
-  const size_t nb_prisms     = gmsh_data.__prisms.size();
-  const size_t nb_pyramids   = gmsh_data.__pyramids.size();
+  const size_t nb_tetrahedra = gmsh_data.tetrahedron_entity.size();
+  const size_t nb_hexahedra  = gmsh_data.hexahedron_entity.size();
+  const size_t nb_prisms     = gmsh_data.prism_entity.size();
+  const size_t nb_pyramids   = gmsh_data.pyramid_entity.size();
 
   Array<unsigned int> cell_to_node_row(nb_cells + 1);
   {
@@ -468,62 +466,62 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   {
     size_t i_cell_node = 0;
     for (size_t i_tetrahedron = 0; i_tetrahedron < nb_tetrahedra; ++i_tetrahedron) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][2];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.tetrahedron_entity[i_tetrahedron][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.tetrahedron_entity[i_tetrahedron][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.tetrahedron_entity[i_tetrahedron][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.tetrahedron_entity[i_tetrahedron][3];
     }
     for (size_t i_hexahedron = 0; i_hexahedron < nb_hexahedra; ++i_hexahedron) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][2];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][3];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][4];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][5];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][6];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][7];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][4];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][5];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][6];
+      cell_to_node_list[i_cell_node++] = gmsh_data.hexahedron_entity[i_hexahedron][7];
     }
     for (size_t i_prism = 0; i_prism < nb_prisms; ++i_prism) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][2];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][3];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][4];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][5];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][4];
+      cell_to_node_list[i_cell_node++] = gmsh_data.prism_entity[i_prism][5];
     }
     for (size_t i_pyramid = 0; i_pyramid < nb_pyramids; ++i_pyramid) {
-      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][0];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][1];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][2];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][3];
-      cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][4];
+      cell_to_node_list[i_cell_node++] = gmsh_data.pyramid_entity[i_pyramid][0];
+      cell_to_node_list[i_cell_node++] = gmsh_data.pyramid_entity[i_pyramid][1];
+      cell_to_node_list[i_cell_node++] = gmsh_data.pyramid_entity[i_pyramid][2];
+      cell_to_node_list[i_cell_node++] = gmsh_data.pyramid_entity[i_pyramid][3];
+      cell_to_node_list[i_cell_node++] = gmsh_data.pyramid_entity[i_pyramid][4];
     }
   }
 
   for (size_t j = 0; j < nb_tetrahedra; ++j) {
     cell_type_vector[j]   = CellType::Tetrahedron;
-    cell_number_vector[j] = gmsh_data.__tetrahedra_number[j];
+    cell_number_vector[j] = gmsh_data.tetrahedron_entity_number[j];
   }
 
   for (size_t j = 0; j < nb_hexahedra; ++j) {
     const size_t jh = nb_tetrahedra + j;
 
     cell_type_vector[jh]   = CellType::Hexahedron;
-    cell_number_vector[jh] = gmsh_data.__hexahedra_number[j];
+    cell_number_vector[jh] = gmsh_data.hexahedron_entity_number[j];
   }
 
   for (size_t j = 0; j < nb_prisms; ++j) {
     const size_t jp = nb_tetrahedra + nb_hexahedra + j;
 
     cell_type_vector[jp]   = CellType::Prism;
-    cell_number_vector[jp] = gmsh_data.__prisms_number[j];
+    cell_number_vector[jp] = gmsh_data.prism_entity_number[j];
   }
 
   for (size_t j = 0; j < nb_pyramids; ++j) {
     const size_t jh = nb_tetrahedra + nb_hexahedra + nb_prisms + j;
 
     cell_type_vector[jh]   = CellType::Pyramid;
-    cell_number_vector[jh] = gmsh_data.__pyramids_number[j];
+    cell_number_vector[jh] = gmsh_data.pyramid_entity_number[j];
   }
 
   descriptor.setCellNumberVector(cell_number_vector);
@@ -531,24 +529,24 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
-  for (unsigned int j = 0; j < gmsh_data.__tetrahedra_ref.size(); ++j) {
+  for (unsigned int j = 0; j < gmsh_data.tetrahedron_entity_ref.size(); ++j) {
     const unsigned int elem_number = j;
-    ref_cells_map[gmsh_data.__tetrahedra_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.tetrahedron_entity_ref[j]].push_back(elem_number);
   }
 
-  for (unsigned int j = 0; j < gmsh_data.__hexahedra_ref.size(); ++j) {
+  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.__hexahedra_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.hexahedron_entity_ref[j]].push_back(elem_number);
   }
 
-  for (unsigned int j = 0; j < gmsh_data.__prisms_ref.size(); ++j) {
+  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.__prisms_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.prism_entity_ref[j]].push_back(elem_number);
   }
 
-  for (unsigned int j = 0; j < gmsh_data.__pyramids_ref.size(); ++j) {
+  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.__pyramids_ref[j]].push_back(elem_number);
+    ref_cells_map[gmsh_data.pyramid_entity_ref[j]].push_back(elem_number);
   }
 
   for (const auto& ref_cell_list : ref_cells_map) {
@@ -571,7 +569,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
 
-  const auto& node_number_vector = descriptor.nodeNumberVector();
+  const auto& node_number = descriptor.nodeNumberVector();
 
   Array<size_t> face_nb_cell(descriptor.faceNumberVector().size());
   face_nb_cell.fill(0);
@@ -599,10 +597,10 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     const auto& node_to_face_matrix = descriptor.nodeToFaceMatrix();
     const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
 
-    const auto find_face = [&](std::vector<uint32_t> node_list) {
+    const auto find_face = [&](std::vector<unsigned int> node_list) {
       size_t i_node_smallest_number = 0;
       for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
-        if (node_number_vector[node_list[i_node]] < node_number_vector[node_list[i_node_smallest_number]]) {
+        if (node_number[node_list[i_node]] < node_number[node_list[i_node_smallest_number]]) {
           i_node_smallest_number = i_node;
         }
       }
@@ -621,7 +619,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
         }
       }
 
-      if (node_number_vector[node_list[1]] > node_number_vector[node_list[node_list.size() - 1]]) {
+      if (node_number[node_list[1]] > node_number[node_list[node_list.size() - 1]]) {
         for (size_t i_node = 1; i_node <= (node_list.size() + 1) / 2 - 1; ++i_node) {
           std::swap(node_list[i_node], node_list[node_list.size() - i_node]);
         }
@@ -654,14 +652,14 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
     Array<int> face_number_vector = copy(descriptor.faceNumberVector());
     std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
-    for (unsigned int f = 0; f < gmsh_data.__triangles.size(); ++f) {
+    for (unsigned int f = 0; f < gmsh_data.triangle_entity.size(); ++f) {
       const unsigned int face_id =
-        find_face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]});
+        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.__triangles_ref[f]].push_back(face_id);
+      ref_faces_map[gmsh_data.triangle_entity_ref[f]].push_back(face_id);
 
-      if (face_number_vector[face_id] != gmsh_data.__triangles_number[f]) {
-        if (auto i_face = face_number_id_map.find(gmsh_data.__triangles_number[f]);
+      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]);
             i_face != face_number_id_map.end()) {
           const int other_face_id = i_face->second;
           std::swap(face_number_vector[face_id], face_number_vector[other_face_id]);
@@ -673,20 +671,20 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
           face_number_id_map[face_number_vector[other_face_id]] = other_face_id;
         } else {
           face_number_id_map.erase(face_number_vector[face_id]);
-          face_number_vector[face_id]                     = gmsh_data.__triangles_number[f];
+          face_number_vector[face_id]                     = gmsh_data.triangle_entity_number[f];
           face_number_id_map[face_number_vector[face_id]] = face_id;
         }
       }
     }
 
-    for (unsigned int f = 0; f < gmsh_data.__quadrangles.size(); ++f) {
-      const unsigned int face_id = find_face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
-                                              gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]});
+    for (unsigned int f = 0; f < gmsh_data.quadrangle_entity.size(); ++f) {
+      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.__quadrangles_ref[f]].push_back(face_id);
+      ref_faces_map[gmsh_data.quadrangle_entity_ref[f]].push_back(face_id);
 
-      if (face_number_vector[face_id] != gmsh_data.__quadrangles_number[f]) {
-        if (auto i_face = face_number_id_map.find(gmsh_data.__quadrangles_number[f]);
+      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]);
             i_face != face_number_id_map.end()) {
           const int other_face_id = i_face->second;
           std::swap(face_number_vector[face_id], face_number_vector[other_face_id]);
@@ -698,7 +696,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
           face_number_id_map[face_number_vector[other_face_id]] = other_face_id;
         } else {
           face_number_id_map.erase(face_number_vector[face_id]);
-          face_number_vector[face_id]                     = gmsh_data.__quadrangles_number[f];
+          face_number_vector[face_id]                     = gmsh_data.quadrangle_entity_number[f];
           face_number_id_map[face_number_vector[face_id]] = face_id;
         }
       }
@@ -771,7 +769,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     const auto& edge_to_node_matrix = descriptor.edgeToNodeMatrix();
 
     const auto find_edge = [&](uint32_t node0, uint32_t node1) {
-      if (node_number_vector[node0] > node_number_vector[node1]) {
+      if (node_number[node0] > node_number[node1]) {
         std::swap(node0, node1);
       }
       const auto& edge_id_vector = node_to_edge_matrix[node0];
@@ -800,13 +798,14 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
 
     Array<int> edge_number_vector = copy(descriptor.edgeNumberVector());
     std::map<unsigned int, std::vector<unsigned int>> ref_edges_map;
-    for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
-      const unsigned int edge_id = find_edge(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+    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.__edges_ref[e]].push_back(edge_id);
+      ref_edges_map[gmsh_data.edge_entity_ref[e]].push_back(edge_id);
 
-      if (edge_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
-        if (auto i_edge = edge_number_id_map.find(gmsh_data.__edges_number[e]); i_edge != edge_number_id_map.end()) {
+      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]);
+            i_edge != edge_number_id_map.end()) {
           const int other_edge_id = i_edge->second;
           std::swap(edge_number_vector[edge_id], edge_number_vector[other_edge_id]);
 
@@ -817,7 +816,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
           edge_number_id_map[edge_number_vector[other_edge_id]] = other_edge_id;
         } else {
           edge_number_id_map.erase(edge_number_vector[edge_id]);
-          edge_number_vector[edge_id]                     = gmsh_data.__edges_number[e];
+          edge_number_vector[edge_id]                     = gmsh_data.edge_entity_number[e];
           edge_number_id_map[edge_number_vector[edge_id]] = edge_id;
         }
       }
@@ -857,7 +856,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     }
   }
 
-  Array<bool> is_boundary_node(node_number_vector.size());
+  Array<bool> is_boundary_node(node_number.size());
   is_boundary_node.fill(false);
 
   const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
@@ -873,9 +872,9 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-  for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
-    const unsigned int point_number = gmsh_data.__points[r];
-    ref_points_map[gmsh_data.__points_ref[r]].push_back(point_number);
+  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);
   }
 
   for (const auto& ref_point_list : ref_points_map) {
@@ -929,7 +928,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
   }());
 
   descriptor.setNodeOwnerVector([&] {
-    Array<int> node_owner_vector(node_number_vector.size());
+    Array<int> node_owner_vector(node_number.size());
     node_owner_vector.fill(parallel::rank());
     return node_owner_vector;
   }());
@@ -1316,8 +1315,8 @@ struct actions<nodes_section>
       throw NormalError(error_msg.str());
     }
 
-    gmsh_data.__verticesNumbers = vertex_number_list;
-    gmsh_data.__vertices        = vertex_list;
+    gmsh_data.node_number     = vertex_number_list;
+    gmsh_data.node_coordinate = vertex_list;
   }
 };
 
@@ -1360,9 +1359,9 @@ struct actions<elements_section>
     const size_t number_of_elements = nb_elements;
     std::cout << "- Number Of Elements: " << number_of_elements << '\n' << std::flush;
 
-    gmsh_data.__elementNumber.resize(number_of_elements);
-    gmsh_data.__elementType.resize(number_of_elements);
-    gmsh_data.__references.resize(number_of_elements);
+    gmsh_data.entity_number.resize(number_of_elements);
+    gmsh_data.entity_type.resize(number_of_elements);
+    gmsh_data.entity_reference.resize(number_of_elements);
     gmsh_data.__elementVertices.resize(number_of_elements);
 
     std::vector<std::string_view> tokens;
@@ -1387,29 +1386,29 @@ struct actions<elements_section>
           throw std::runtime_error("empty element line");
         }
 
-        gmsh_data.__elementNumber[i_element] = parse_int(tokens[0]);
+        gmsh_data.entity_number[i_element] = parse_int(tokens[0]);
 
-        const int element_type = parse_int(tokens[1]) - 1;
-        if ((element_type < 0) or (element_type > 15)) {
+        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 (not supported_element[element_type]) {
-          throw std::runtime_error(std::string{element_name[element_type]} + " is not supported");
+        if (not supported_element[entity_type]) {
+          throw std::runtime_error(std::string{element_name[entity_type]} + " is not supported");
         }
-        gmsh_data.__elementType[i_element] = element_type;
+        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[element_type]} + " number of tags is negative");
+          throw std::runtime_error(std::string{element_name[entity_type]} + " number of tags is negative");
         }
 
-        const size_t number_of_tags       = nb_tags;
-        gmsh_data.__references[i_element] = parse_int(tokens[3]);
+        const size_t number_of_tags           = nb_tags;
+        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
 
         const size_t first_node_index = 3 + number_of_tags;
-        const size_t number_of_nodes  = primitive_number_of_nodes[element_type];
+        const size_t number_of_nodes  = primitive_number_of_nodes[entity_type];
 
         gmsh_data.__elementVertices[i_element].resize(number_of_nodes);
         for (size_t i_node = 0; i_node < number_of_nodes; ++i_node) {
@@ -1519,25 +1518,25 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 void
 GmshReader::__proceedData()
 {
-  TinyVector<15, int> elementNumber = zero;
-  std::vector<std::set<size_t>> elementReferences(15);
+  TinyVector<15, int> element_number = zero;
+  std::vector<std::set<size_t>> element_reference(15);
 
-  for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
-    const int elementType = m_mesh_data.__elementType[i];
-    elementNumber[elementType]++;
-    elementReferences[elementType].insert(m_mesh_data.__references[i]);
+  for (size_t i = 0; i < m_mesh_data.entity_type.size(); ++i) {
+    const int elementType = m_mesh_data.entity_type[i];
+    element_number[elementType]++;
+    element_reference[elementType].insert(m_mesh_data.entity_reference[i]);
   }
 
-  for (size_t i = 0; i < elementNumber.dimension(); ++i) {
-    if (elementNumber[i] > 0) {
-      std::cout << "  - Number of " << gmsh_parser::element_name[i] << ": " << elementNumber[i];
+  for (size_t i = 0; i < element_number.dimension(); ++i) {
+    if (element_number[i] > 0) {
+      std::cout << "  - Number of " << gmsh_parser::element_name[i] << ": " << element_number[i];
       if (not(gmsh_parser::supported_element[i])) {
         std::cout << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
       } else {
         std::cout << " references: ";
-        for (std::set<size_t>::const_iterator iref = elementReferences[i].begin(); iref != elementReferences[i].end();
+        for (std::set<size_t>::const_iterator iref = element_reference[i].begin(); iref != element_reference[i].end();
              ++iref) {
-          if (iref != elementReferences[i].begin())
+          if (iref != element_reference[i].begin())
             std::cout << ',';
           std::cout << *iref;
         }
@@ -1545,51 +1544,51 @@ GmshReader::__proceedData()
         switch (i) {
         // Supported entities
         case 0: {   // edges
-          m_mesh_data.__edges = Array<GmshData::Edge>(elementNumber[i]);
-          m_mesh_data.__edges_ref.resize(elementNumber[i]);
-          m_mesh_data.__edges_number.resize(elementNumber[i]);
+          m_mesh_data.edge_entity        = Array<GmshData::Edge>(element_number[i]);
+          m_mesh_data.edge_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.edge_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 1: {   // triangles
-          m_mesh_data.__triangles = Array<GmshData::Triangle>(elementNumber[i]);
-          m_mesh_data.__triangles_ref.resize(elementNumber[i]);
-          m_mesh_data.__triangles_number.resize(elementNumber[i]);
+          m_mesh_data.triangle_entity        = Array<GmshData::Triangle>(element_number[i]);
+          m_mesh_data.triangle_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.triangle_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 2: {   // quadrangles
-          m_mesh_data.__quadrangles = Array<GmshData::Quadrangle>(elementNumber[i]);
-          m_mesh_data.__quadrangles_ref.resize(elementNumber[i]);
-          m_mesh_data.__quadrangles_number.resize(elementNumber[i]);
+          m_mesh_data.quadrangle_entity        = Array<GmshData::Quadrangle>(element_number[i]);
+          m_mesh_data.quadrangle_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.quadrangle_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 3: {   // tetrahedra
-          m_mesh_data.__tetrahedra = Array<GmshData::Tetrahedron>(elementNumber[i]);
-          m_mesh_data.__tetrahedra_ref.resize(elementNumber[i]);
-          m_mesh_data.__tetrahedra_number.resize(elementNumber[i]);
+          m_mesh_data.tetrahedron_entity        = Array<GmshData::Tetrahedron>(element_number[i]);
+          m_mesh_data.tetrahedron_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.tetrahedron_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 4: {   // hexahedra
-          m_mesh_data.__hexahedra = Array<GmshData::Hexahedron>(elementNumber[i]);
-          m_mesh_data.__hexahedra_ref.resize(elementNumber[i]);
-          m_mesh_data.__hexahedra_number.resize(elementNumber[i]);
+          m_mesh_data.hexahedron_entity        = Array<GmshData::Hexahedron>(element_number[i]);
+          m_mesh_data.hexahedron_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.hexahedron_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 5: {   // prism
-          m_mesh_data.__prisms = Array<GmshData::Prism>(elementNumber[i]);
-          m_mesh_data.__prisms_ref.resize(elementNumber[i]);
-          m_mesh_data.__prisms_number.resize(elementNumber[i]);
+          m_mesh_data.prism_entity        = Array<GmshData::Prism>(element_number[i]);
+          m_mesh_data.prism_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.prism_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 6: {   // pyramid
-          m_mesh_data.__pyramids = Array<GmshData::Pyramid>(elementNumber[i]);
-          m_mesh_data.__pyramids_ref.resize(elementNumber[i]);
-          m_mesh_data.__pyramids_number.resize(elementNumber[i]);
+          m_mesh_data.pyramid_entity        = Array<GmshData::Pyramid>(element_number[i]);
+          m_mesh_data.pyramid_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.pyramid_entity_number = Array<int>(element_number[i]);
           break;
         }
         case 14: {   // point
-          m_mesh_data.__points = Array<GmshData::Point>(elementNumber[i]);
-          m_mesh_data.__points_ref.resize(elementNumber[i]);
-          m_mesh_data.__points_number.resize(elementNumber[i]);
+          m_mesh_data.point_entity        = Array<GmshData::Point>(element_number[i]);
+          m_mesh_data.point_entity_ref    = Array<int>(element_number[i]);
+          m_mesh_data.point_entity_number = Array<int>(element_number[i]);
           break;
         }
           // Unsupported entities
@@ -1609,156 +1608,135 @@ GmshReader::__proceedData()
     }
   }
 
-  const int minNumber = min(m_mesh_data.__verticesNumbers);
-  const int maxNumber = max(m_mesh_data.__verticesNumbers);
+  const int min_node_number = min(m_mesh_data.node_number);
+  const int max_node_number = max(m_mesh_data.node_number);
 
-  if (minNumber < 0) {
+  if (min_node_number < 0) {
     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.__verticesCorrepondance.resize(maxNumber + 1, -1);
+  m_mesh_data.gmsh_node_to_node_id_correspondance.resize(max_node_number + 1, -1);
 
-  for (size_t i = 0; i < m_mesh_data.__verticesNumbers.size(); ++i) {
-    m_mesh_data.__verticesCorrepondance[m_mesh_data.__verticesNumbers[i]] = i;
+  for (size_t i = 0; i < m_mesh_data.node_number.size(); ++i) {
+    m_mesh_data.gmsh_node_to_node_id_correspondance[m_mesh_data.node_number[i]] = i;
   }
 
   // reset element number to count them while filling structures
-  elementNumber = zero;
+  element_number = zero;
 
   // Element structures filling
-  for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
+  for (size_t i = 0; i < m_mesh_data.entity_type.size(); ++i) {
     std::vector<int>& elementVertices = m_mesh_data.__elementVertices[i];
-    switch (m_mesh_data.__elementType[i]) {
+    switch (m_mesh_data.entity_type[i]) {
     // Supported entities
     case 0: {   // edge
-      int& edgeNumber = elementNumber[0];
-      const int a     = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b     = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      if ((a < 0) or (b < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__edges[edgeNumber]        = GmshData::Edge(a, b);
-      m_mesh_data.__edges_ref[edgeNumber]    = m_mesh_data.__references[i];
-      m_mesh_data.__edges_number[edgeNumber] = m_mesh_data.__elementNumber[i];
-      edgeNumber++;   // one more edge
+      int& edge_number = element_number[0];
+      const auto a     = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b     = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+
+      m_mesh_data.edge_entity[edge_number]        = GmshData::Edge{a, b};
+      m_mesh_data.edge_entity_ref[edge_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.edge_entity_number[edge_number] = m_mesh_data.entity_number[i];
+      edge_number++;   // one more edge
       break;
     }
     case 1: {   // triangles
-      int& triangleNumber = elementNumber[1];
-
-      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      if ((a < 0) or (b < 0) or (c < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__triangles[triangleNumber]        = GmshData::Triangle(a, b, c);
-      m_mesh_data.__triangles_ref[triangleNumber]    = m_mesh_data.__references[i];
-      m_mesh_data.__triangles_number[triangleNumber] = m_mesh_data.__elementNumber[i];
-      triangleNumber++;   // one more triangle
+      int& triangle_number = element_number[1];
+
+      const auto a = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+
+      m_mesh_data.triangle_entity[triangle_number]        = GmshData::Triangle{a, b, c};
+      m_mesh_data.triangle_entity_ref[triangle_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.triangle_entity_number[triangle_number] = m_mesh_data.entity_number[i];
+      triangle_number++;   // one more triangle
       break;
     }
     case 2: {   // quadrangle
-      int& quadrilateralNumber = elementNumber[2];
-
-      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
-      if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__quadrangles[quadrilateralNumber]        = GmshData::Quadrangle(a, b, c, d);
-      m_mesh_data.__quadrangles_ref[quadrilateralNumber]    = m_mesh_data.__references[i];
-      m_mesh_data.__quadrangles_number[quadrilateralNumber] = m_mesh_data.__elementNumber[i];
-      quadrilateralNumber++;   // one more quadrangle
+      int& quadrilateral_number = element_number[2];
+
+      const auto a = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+      const auto d = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[3]];
+
+      m_mesh_data.quadrangle_entity[quadrilateral_number]        = GmshData::Quadrangle{a, b, c, d};
+      m_mesh_data.quadrangle_entity_ref[quadrilateral_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.quadrangle_entity_number[quadrilateral_number] = m_mesh_data.entity_number[i];
+      quadrilateral_number++;   // one more quadrangle
       break;
     }
     case 3: {   // tetrahedra
-      int& tetrahedronNumber = elementNumber[3];
-
-      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
-      if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__tetrahedra[tetrahedronNumber]        = GmshData::Tetrahedron(a, b, c, d);
-      m_mesh_data.__tetrahedra_ref[tetrahedronNumber]    = m_mesh_data.__references[i];
-      m_mesh_data.__tetrahedra_number[tetrahedronNumber] = m_mesh_data.__elementNumber[i];
-      tetrahedronNumber++;   // one more tetrahedron
+      int& tetrahedron_number = element_number[3];
+
+      const auto a = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+      const auto d = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[3]];
+
+      m_mesh_data.tetrahedron_entity[tetrahedron_number]        = GmshData::Tetrahedron{a, b, c, d};
+      m_mesh_data.tetrahedron_entity_ref[tetrahedron_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.tetrahedron_entity_number[tetrahedron_number] = m_mesh_data.entity_number[i];
+      tetrahedron_number++;   // one more tetrahedron
       break;
     }
     case 4: {   // hexaredron
-      int& hexahedronNumber = elementNumber[4];
-
-      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
-      const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
-      const int f = m_mesh_data.__verticesCorrepondance[elementVertices[5]];
-      const int g = m_mesh_data.__verticesCorrepondance[elementVertices[6]];
-      const int h = m_mesh_data.__verticesCorrepondance[elementVertices[7]];
-      if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0) or (g < 0) or (h < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__hexahedra[hexahedronNumber]        = GmshData::Hexahedron(a, b, c, d, e, f, g, h);
-      m_mesh_data.__hexahedra_ref[hexahedronNumber]    = m_mesh_data.__references[i];
-      m_mesh_data.__hexahedra_number[hexahedronNumber] = m_mesh_data.__elementNumber[i];
-      hexahedronNumber++;   // one more hexahedron
+      int& hexahedron_number = element_number[4];
+
+      const auto a = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+      const auto d = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[3]];
+      const auto e = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[4]];
+      const auto f = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[5]];
+      const auto g = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[6]];
+      const auto h = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[7]];
+
+      m_mesh_data.hexahedron_entity[hexahedron_number]        = GmshData::Hexahedron{a, b, c, d, e, f, g, h};
+      m_mesh_data.hexahedron_entity_ref[hexahedron_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.hexahedron_entity_number[hexahedron_number] = m_mesh_data.entity_number[i];
+      hexahedron_number++;   // one more hexahedron
       break;
     }
     case 5: {   // prism
-      int& prism_number = elementNumber[5];
-      const int a       = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b       = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c       = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      const int d       = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
-      const int e       = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
-      const int f       = m_mesh_data.__verticesCorrepondance[elementVertices[5]];
-      if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__prisms[prism_number]        = GmshData::Prism(a, b, c, d, e, f);
-      m_mesh_data.__prisms_ref[prism_number]    = m_mesh_data.__references[i];
-      m_mesh_data.__prisms_number[prism_number] = m_mesh_data.__elementNumber[i];
+      int& prism_number = element_number[5];
+      const auto a      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+      const auto d      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[3]];
+      const auto e      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[4]];
+      const auto f      = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[5]];
+
+      m_mesh_data.prism_entity[prism_number]        = GmshData::Prism{a, b, c, d, e, f};
+      m_mesh_data.prism_entity_ref[prism_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.prism_entity_number[prism_number] = m_mesh_data.entity_number[i];
       prism_number++;   // one more prism
       break;
     }
     case 6: {   // pyramid
-      int& pyramid_number = elementNumber[6];
-
-      const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
-      const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
-      const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
-      const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
-      if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
-                          " [bad vertices definition]");
-      }
-      m_mesh_data.__pyramids[pyramid_number]        = GmshData::Pyramid(a, b, c, d, e);
-      m_mesh_data.__pyramids_ref[pyramid_number]    = m_mesh_data.__references[i];
-      m_mesh_data.__pyramids_number[pyramid_number] = m_mesh_data.__elementNumber[i];
+      int& pyramid_number = element_number[6];
+
+      const auto a = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      const auto b = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[1]];
+      const auto c = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[2]];
+      const auto d = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[3]];
+      const auto e = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[4]];
+
+      m_mesh_data.pyramid_entity[pyramid_number]        = GmshData::Pyramid{a, b, c, d, e};
+      m_mesh_data.pyramid_entity_ref[pyramid_number]    = m_mesh_data.entity_reference[i];
+      m_mesh_data.pyramid_entity_number[pyramid_number] = m_mesh_data.entity_number[i];
       pyramid_number++;   // one more hexahedron
       break;
     }
       // Unsupported entities
     case 14: {   // point
-      int& point_number                         = elementNumber[14];
-      const int a                               = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
-      m_mesh_data.__points[point_number]        = a;
-      m_mesh_data.__points_ref[point_number]    = m_mesh_data.__references[i];
-      m_mesh_data.__points_number[point_number] = m_mesh_data.__elementNumber[i];
+      int& point_number                          = element_number[14];
+      m_mesh_data.point_entity[point_number]     = m_mesh_data.gmsh_node_to_node_id_correspondance[elementVertices[0]];
+      m_mesh_data.point_entity_ref[point_number] = m_mesh_data.entity_reference[i];
+      m_mesh_data.point_entity_number[point_number] = m_mesh_data.entity_number[i];
       point_number++;
       break;
     }
@@ -1771,7 +1749,7 @@ GmshReader::__proceedData()
     case 13: {   // second order pyramid
     }
     default: {
-      throw NormalError("reading file '" + m_filename + "': ff3d cannot read this kind of element");
+      throw NormalError("reading file '" + m_filename + "': cannot read this kind of element");
     }
     }
   }
@@ -1796,8 +1774,8 @@ GmshReader::__proceedData()
   dimension3_mask[12]                 = 1;
   dimension3_mask[13]                 = 1;
 
-  if (dot(dimension3_mask, elementNumber) > 0) {
-    const size_t nb_cells = dot(dimension3_mask, elementNumber);
+  if (dot(dimension3_mask, element_number) > 0) {
+    const size_t nb_cells = dot(dimension3_mask, element_number);
 
     GmshConnectivityBuilder<3> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1808,17 +1786,12 @@ GmshReader::__proceedData()
     using MeshType = Mesh<3>;
     using Rd       = TinyVector<3, double>;
 
-    NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
-      xr[i][0] = m_mesh_data.__vertices[i][0];
-      xr[i][1] = m_mesh_data.__vertices[i][1];
-      xr[i][2] = m_mesh_data.__vertices[i][2];
-    }
+    NodeValue<Rd> xr(connectivity, m_mesh_data.node_coordinate);
     m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(p_connectivity, xr));
     this->_checkMesh<3>();
 
-  } else if (dot(dimension2_mask, elementNumber) > 0) {
-    const size_t nb_cells = dot(dimension2_mask, elementNumber);
+  } else if (dot(dimension2_mask, element_number) > 0) {
+    const size_t nb_cells = dot(dimension2_mask, element_number);
 
     GmshConnectivityBuilder<2> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1830,15 +1803,15 @@ GmshReader::__proceedData()
     using Rd       = TinyVector<2, double>;
 
     NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
-      xr[i][0] = m_mesh_data.__vertices[i][0];
-      xr[i][1] = m_mesh_data.__vertices[i][1];
+    for (NodeId i = 0; i < m_mesh_data.node_coordinate.size(); ++i) {
+      xr[i][0] = m_mesh_data.node_coordinate[i][0];
+      xr[i][1] = m_mesh_data.node_coordinate[i][1];
     }
     m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(p_connectivity, xr));
     this->_checkMesh<2>();
 
-  } else if (dot(dimension1_mask, elementNumber) > 0) {
-    const size_t nb_cells = dot(dimension1_mask, elementNumber);
+  } else if (dot(dimension1_mask, element_number) > 0) {
+    const size_t nb_cells = dot(dimension1_mask, element_number);
 
     GmshConnectivityBuilder<1> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1850,8 +1823,8 @@ GmshReader::__proceedData()
     using Rd       = TinyVector<1, double>;
 
     NodeValue<Rd> xr(connectivity);
-    for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
-      xr[i][0] = m_mesh_data.__vertices[i][0];
+    for (NodeId i = 0; i < m_mesh_data.node_coordinate.size(); ++i) {
+      xr[i][0] = m_mesh_data.node_coordinate[i][0];
     }
 
     m_mesh = std::make_shared<MeshVariant>(std::make_shared<const MeshType>(p_connectivity, xr));
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 325622bdbe1c245db96ac47f07164729c7fd6ccb..24682243937853556e89764e2ca906d4505407af 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -78,70 +78,70 @@ class GmshReader : public MeshBuilderBase
      * Gmsh format provides a numbered, none ordrered and none dense
      * vertices list, this stores the number of read vertices.
      */
-    Array<int> __verticesNumbers;
+    Array<int> node_number;
 
-    Array<R3> __vertices;
+    Array<R3> node_coordinate;
 
     using Point = unsigned int;
-    Array<Point> __points;
-    std::vector<int> __points_ref;
-    std::vector<int> __points_number;
-
-    using Edge = TinyVector<2, unsigned int>;
-    Array<Edge> __edges;
-    std::vector<int> __edges_ref;
-    std::vector<int> __edges_number;
-
-    using Triangle = TinyVector<3, unsigned int>;
-    Array<Triangle> __triangles;
-    std::vector<int> __triangles_ref;
-    std::vector<int> __triangles_number;
-
-    using Quadrangle = TinyVector<4, unsigned int>;
-    Array<Quadrangle> __quadrangles;
-    std::vector<int> __quadrangles_ref;
-    std::vector<int> __quadrangles_number;
-
-    using Tetrahedron = TinyVector<4, unsigned int>;
-    Array<Tetrahedron> __tetrahedra;
-    std::vector<int> __tetrahedra_ref;
-    std::vector<int> __tetrahedra_number;
-
-    using Pyramid = TinyVector<5, unsigned int>;
-    Array<Pyramid> __pyramids;
-    std::vector<int> __pyramids_ref;
-    std::vector<int> __pyramids_number;
-
-    using Prism = TinyVector<6, unsigned int>;
-    Array<Prism> __prisms;
-    std::vector<int> __prisms_ref;
-    std::vector<int> __prisms_number;
-
-    using Hexahedron = TinyVector<8, unsigned int>;
-    Array<Hexahedron> __hexahedra;
-    std::vector<int> __hexahedra_ref;
-    std::vector<int> __hexahedra_number;
+    Array<Point> point_entity;
+    Array<int> point_entity_ref;
+    Array<int> point_entity_number;
+
+    using Edge = std::array<unsigned int, 2>;
+    Array<Edge> edge_entity;
+    Array<int> edge_entity_ref;
+    Array<int> edge_entity_number;
+
+    using Triangle = std::array<unsigned int, 3>;
+    Array<Triangle> triangle_entity;
+    Array<int> triangle_entity_ref;
+    Array<int> triangle_entity_number;
+
+    using Quadrangle = std::array<unsigned int, 4>;
+    Array<Quadrangle> quadrangle_entity;
+    Array<int> quadrangle_entity_ref;
+    Array<int> quadrangle_entity_number;
+
+    using Tetrahedron = std::array<unsigned int, 4>;
+    Array<Tetrahedron> tetrahedron_entity;
+    Array<int> tetrahedron_entity_ref;
+    Array<int> tetrahedron_entity_number;
+
+    using Pyramid = std::array<unsigned int, 5>;
+    Array<Pyramid> pyramid_entity;
+    Array<int> pyramid_entity_ref;
+    Array<int> pyramid_entity_number;
+
+    using Prism = std::array<unsigned int, 6>;
+    Array<Prism> prism_entity;
+    Array<int> prism_entity_ref;
+    Array<int> prism_entity_number;
+
+    using Hexahedron = std::array<unsigned int, 8>;
+    Array<Hexahedron> hexahedron_entity;
+    Array<int> hexahedron_entity_ref;
+    Array<int> hexahedron_entity_number;
 
     /**
      * Gmsh format provides a numbered, none ordrered and none dense
      * vertices list, this provides vertices renumbering correspondence
      */
-    std::vector<int> __verticesCorrepondance;
+    std::vector<unsigned int> gmsh_node_to_node_id_correspondance;
 
     /**
      * elements types
      */
-    std::vector<int> __elementNumber;
+    std::vector<int> entity_number;
 
     /**
      * elements types
      */
-    std::vector<short> __elementType;
+    std::vector<short> entity_type;
 
     /**
      * References
      */
-    std::vector<int> __references;
+    std::vector<int> entity_reference;
 
     /**
      * References