From cd009e481ea6f421d1e0f558239717c156a1b413 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 10 Jul 2025 19:25:34 +0200
Subject: [PATCH] Begin entity tag to physical tags handling in mesh
construction
---
src/mesh/GmshReader.cpp | 223 ++++++++++++++++++++++++++++++++++------
src/mesh/GmshReader.hpp | 4 +-
2 files changed, 192 insertions(+), 35 deletions(-)
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f2b57eb72..c89836dad 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) {
@@ -107,8 +108,8 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
}();
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_point_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_point_list.first};
@@ -131,8 +132,8 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_cell_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_cell_list.first};
@@ -212,10 +213,14 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
descriptor.setCellNumberVector(cell_number_vector);
descriptor.setCellTypeVector(cell_type_vector);
+ auto& cell_to_physical_tag = gmsh_data.entity_to_physical_tag_map[2];
std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
for (unsigned int j = 0; j < gmsh_data.triangle_entity_ref.size(); ++j) {
const unsigned int elem_number = j;
- ref_cells_map[gmsh_data.triangle_entity_ref[j]].push_back(elem_number);
+ auto ref_range = cell_to_physical_tag.equal_range(gmsh_data.triangle_entity_ref[j]);
+ for (auto i_ref = ref_range.first; i_ref != ref_range.second; ++i_ref) {
+ ref_cells_map[i_ref->second].push_back(elem_number);
+ }
}
for (unsigned int j = 0; j < gmsh_data.quadrangle_entity_ref.size(); ++j) {
@@ -230,8 +235,8 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 2)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_cell_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 2)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_cell_list.first};
@@ -319,8 +324,8 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_face_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_face_list.first};
@@ -381,8 +386,8 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_point_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_point_list.first};
@@ -554,8 +559,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 3)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_cell_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 3)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_cell_list.first};
@@ -708,8 +713,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 2)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_face_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 2)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_face_list.first};
@@ -828,8 +833,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_edge_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_edge_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 1)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_edge_list.first};
@@ -882,8 +887,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
auto ref_id = [&] {
- auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
- if ((i_physical_ref != gmsh_data.m_physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
+ auto i_physical_ref = gmsh_data.physical_ref_map.find(ref_point_list.first);
+ if ((i_physical_ref != gmsh_data.physical_ref_map.end()) and (i_physical_ref->second.dimension() == 0)) {
return i_physical_ref->second.refId();
} else {
return RefId{ref_point_list.first};
@@ -1053,7 +1058,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);
@@ -1092,9 +1098,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";
@@ -1132,9 +1140,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
@@ -1147,7 +1157,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;
@@ -1340,14 +1351,14 @@ struct actions<physical_section>
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()) {
+ if (auto i_searched_physical_ref_id = gmsh_data.physical_ref_map.find(physical_number);
+ i_searched_physical_ref_id != gmsh_data.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;
+ gmsh_data.physical_ref_map[physical_number] = physical_ref_id;
}
}
catch (std::runtime_error& e) {
@@ -1356,11 +1367,6 @@ struct actions<physical_section>
<< e.what();
throw NormalError(error_msg.str());
}
-
- std::cout << "===== PhysicalNames =====\n";
- for (auto i : gmsh_data.m_physical_ref_map) {
- std::cout << i.first << " -> " << i.second << '\n';
- }
} else {
throw UnexpectedError("invalid file format version: " + std::to_string(gmsh_data.format.version));
}
@@ -1379,6 +1385,145 @@ struct actions<entities_section>
} else {
state.has_entities_section = true;
}
+
+ if (gmsh_data.format.version == 2.2) {
+ throw NormalError("unexpected $Entities section in 2.2 file format");
+ }
+
+ 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, $Entities section is empty");
+ }
+
+ body = body.substr(pos + 1);
+
+ std::vector<std::string_view> tokens;
+ tokens.reserve(20);
+
+ try {
+ ++line_number;
+ auto header_line = remove_embedding_spaces(get_next_line(body));
+
+ if (header_line.empty()) {
+ throw std::runtime_error("empty entry");
+ }
+
+ std::array<std::string_view, 4> header_tokens;
+ const size_t number_of_header_tokens = tokenize(header_line, header_tokens);
+
+ if (number_of_header_tokens != 4) {
+ throw std::runtime_error{"expecting 4 entries"};
+ }
+
+ const size_t number_of_points = parse_size_t(header_tokens[0]);
+ const size_t number_of_curves = parse_size_t(header_tokens[1]);
+ const size_t number_of_surfaces = parse_size_t(header_tokens[2]);
+ const size_t number_of_volumes = parse_size_t(header_tokens[3]);
+
+ auto& entity_to_physical_tag_map = gmsh_data.entity_to_physical_tag_map;
+
+ for (size_t i_point = 0; i_point < number_of_points; ++i_point) {
+ ++line_number;
+ auto point_line = remove_embedding_spaces(get_next_line(body));
+ tokenize_variable(point_line, tokens);
+ if (tokens.size() < 5) {
+ throw std::runtime_error("expecting at least 5 entries");
+ }
+ const int point_tag = parse_int(tokens[0]);
+ parse_int(tokens[1]); // x
+ parse_int(tokens[2]); // y
+ parse_int(tokens[3]); // z
+ const size_t nb_physical_tags = parse_size_t(tokens[4]);
+ if (tokens.size() < 5 + nb_physical_tags) {
+ throw std::runtime_error("missing physical tags");
+ }
+ for (size_t i_tag = 0; i_tag < nb_physical_tags; ++i_tag) {
+ const int physical_tag = parse_int(tokens[5 + i_tag]);
+ entity_to_physical_tag_map[0].insert({point_tag, physical_tag});
+ }
+ }
+
+ for (size_t i_curve = 0; i_curve < number_of_curves; ++i_curve) {
+ ++line_number;
+ auto curve_line = remove_embedding_spaces(get_next_line(body));
+ tokenize_variable(curve_line, tokens);
+ if (tokens.size() < 9) {
+ throw std::runtime_error("expecting at least 9 entries");
+ }
+ const int curve_tag = parse_int(tokens[0]);
+ parse_int(tokens[1]); // x0
+ parse_int(tokens[2]); // y0
+ parse_int(tokens[3]); // z0
+ parse_int(tokens[4]); // x1
+ parse_int(tokens[5]); // y1
+ parse_int(tokens[6]); // z1
+ const size_t nb_physical_tags = parse_size_t(tokens[7]);
+ if (tokens.size() < 9 + nb_physical_tags) {
+ throw std::runtime_error("missing physical tags");
+ }
+ for (size_t i_tag = 0; i_tag < nb_physical_tags; ++i_tag) {
+ const int physical_tag = parse_int(tokens[8 + i_tag]);
+ entity_to_physical_tag_map[1].insert({curve_tag, physical_tag});
+ }
+ }
+
+ for (size_t i_surface = 0; i_surface < number_of_surfaces; ++i_surface) {
+ ++line_number;
+ auto surface_line = remove_embedding_spaces(get_next_line(body));
+ tokenize_variable(surface_line, tokens);
+ if (tokens.size() < 9) {
+ throw std::runtime_error("expecting at least 9 entries");
+ }
+ const int surface_tag = parse_int(tokens[0]);
+ parse_int(tokens[1]); // x0
+ parse_int(tokens[2]); // y0
+ parse_int(tokens[3]); // z0
+ parse_int(tokens[4]); // x1
+ parse_int(tokens[5]); // y1
+ parse_int(tokens[6]); // z1
+ const size_t nb_physical_tags = parse_size_t(tokens[7]);
+ if (tokens.size() < 9 + nb_physical_tags) {
+ throw std::runtime_error("missing physical tags");
+ }
+ for (size_t i_tag = 0; i_tag < nb_physical_tags; ++i_tag) {
+ const int physical_tag = parse_int(tokens[8 + i_tag]);
+ entity_to_physical_tag_map[2].insert({surface_tag, physical_tag});
+ }
+ }
+
+ for (size_t i_volume = 0; i_volume < number_of_volumes; ++i_volume) {
+ ++line_number;
+ auto volume_line = remove_embedding_spaces(get_next_line(body));
+ tokenize_variable(volume_line, tokens);
+ if (tokens.size() < 9) {
+ throw std::runtime_error("expecting at least 9 entries");
+ }
+ const int volume_tag = parse_int(tokens[0]);
+ parse_int(tokens[1]); // x0
+ parse_int(tokens[2]); // y0
+ parse_int(tokens[3]); // z0
+ parse_int(tokens[4]); // x1
+ parse_int(tokens[5]); // y1
+ parse_int(tokens[6]); // z1
+ const size_t nb_physical_tags = parse_size_t(tokens[7]);
+ if (tokens.size() < 9 + nb_physical_tags) {
+ throw std::runtime_error("missing physical tags");
+ }
+ for (size_t i_tag = 0; i_tag < nb_physical_tags; ++i_tag) {
+ const int physical_tag = parse_int(tokens[8 + i_tag]);
+ entity_to_physical_tag_map[3].insert({volume_tag, physical_tag});
+ }
+ }
+ }
+ catch (std::runtime_error& e) {
+ std::stringstream error_msg;
+ error_msg << "reading " << in.position().source << ':' << line_number << ": " << e.what();
+ throw NormalError(error_msg.str());
+ }
}
};
@@ -1886,6 +2031,16 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
if (not state.has_elements_section) {
throw NormalError("malformed Gmsh file, $Elements section missing");
}
+
+ if (m_mesh_data.format.version == 2.2) {
+ // populate fake entity_to_physical_tag_map for simpler ref
+ // handling: entity tag == physical tag
+ for (auto i_ref : m_mesh_data.physical_ref_map) {
+ const size_t key = i_ref.second.refId().tagNumber();
+ const int value = key;
+ m_mesh_data.entity_to_physical_tag_map[i_ref.second.dimension()].insert({key, value});
+ }
+ }
}
catch (std::filesystem::filesystem_error& e) {
throw NormalError("cannot open file '" + m_filename + "'");
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 4a5479aac..e8a546081 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -148,7 +148,9 @@ class GmshReader : public MeshBuilderBase
*/
std::vector<std::vector<int>> element_vertices;
- std::map<unsigned int, PhysicalRefId> m_physical_ref_map;
+ std::map<unsigned int, PhysicalRefId> physical_ref_map;
+
+ std::array<std::unordered_multimap<size_t, int>, 4> entity_to_physical_tag_map;
};
private:
--
GitLab