PolyFEM
Loading...
Searching...
No Matches
MshReader.cpp
Go to the documentation of this file.
1#include "MshReader.hpp"
2
3#include <polyfem/utils/Logger.hpp>
4#include <polyfem/utils/StringUtils.hpp>
5
6#include <mshio/mshio.h>
7
8#include <fstream>
9#include <string>
10#include <iostream>
11#include <vector>
12
13#include <filesystem> // filesystem
14
15namespace polyfem::io
16{
17 template <typename Entity>
18 void map_entity_tag_to_physical_tag(const std::vector<Entity> &entities, std::unordered_map<int, int> &entity_tag_to_physical_tag)
19 {
20 for (int i = 0; i < entities.size(); i++)
21 {
22 entity_tag_to_physical_tag[entities[i].tag] =
23 entities[i].physical_group_tags.size() > 0
24 ? entities[i].physical_group_tags.front()
25 : 0;
26 }
27 }
28
29 bool MshReader::load(const std::string &path, Eigen::MatrixXd &vertices, Eigen::MatrixXi &cells, std::vector<std::vector<int>> &elements, std::vector<std::vector<double>> &weights, std::vector<int> &body_ids)
30 {
31 std::vector<std::string> node_data_name;
32 std::vector<std::vector<double>> node_data;
33
34 return load(path, vertices, cells, elements, weights, body_ids, node_data_name, node_data);
35 }
36
37 bool MshReader::load(const std::string &path, Eigen::MatrixXd &vertices, Eigen::MatrixXi &cells, std::vector<std::vector<int>> &elements, std::vector<std::vector<double>> &weights, std::vector<int> &body_ids, std::vector<std::string> &node_data_name, std::vector<std::vector<double>> &node_data)
38 {
39 if (!std::filesystem::exists(path))
40 {
41 logger().error("Msh file does not exist: {}", path);
42 return false;
43 }
44
45 mshio::MshSpec spec;
46 try
47 {
48 spec = mshio::load_msh(path);
49 }
50 catch (const std::exception &err)
51 {
52 logger().error("{}", err.what());
53 return false;
54 }
55 catch (...)
56 {
57 logger().error("Unknown error while reading MSH file: {}", path);
58 return false;
59 }
60
61 const auto &nodes = spec.nodes;
62 const auto &els = spec.elements;
63 const int n_vertices = nodes.num_nodes;
64 const int max_tag = nodes.max_node_tag;
65 int dim = -1;
66
67 assert(els.entity_blocks.size() > 0);
68 for (const auto &e : els.entity_blocks)
69 {
70 dim = std::max(dim, e.entity_dim);
71 }
72 assert(dim == 2 || dim == 3);
73
74 vertices.resize(n_vertices, dim);
75 std::vector<int> tag_to_index = std::vector<int>(max_tag + 1, -1);
76 if (n_vertices != max_tag)
77 logger().warn("MSH file contains more node tags than nodes, condensing nodes which will break input node ordering.");
78
79 int index = 0;
80 for (const auto &n : nodes.entity_blocks)
81 {
82 for (int i = 0; i < n.num_nodes_in_block * 3; i += 3)
83 {
84 const int node_id = n_vertices != max_tag ? (index++) : (n.tags[i / 3] - 1);
85
86 if (dim == 2)
87 vertices.row(node_id) << n.data[i], n.data[i + 1];
88 if (dim == 3)
89 vertices.row(node_id) << n.data[i], n.data[i + 1], n.data[i + 2];
90
91 assert(n.tags[i / 3] < tag_to_index.size());
92 tag_to_index[n.tags[i / 3]] = node_id;
93 }
94 }
95
96 int cells_cols = -1;
97 int num_els = 0;
98 for (const auto &e : els.entity_blocks)
99 {
100 if (e.entity_dim != dim)
101 continue;
102 const int type = e.element_type;
103 // https://shipengcheng1230.github.io/GmshTools.jl/stable/element_types/
104 if (type == 2 || type == 9 || type == 21 || type == 23 || type == 25) // tri
105 {
106 cells_cols = std::max(cells_cols, 3);
107 num_els += e.num_elements_in_block;
108 }
109 else if (type == 3 || type == 10) // quad
110 {
111 cells_cols = std::max(cells_cols, 4);
112 num_els += e.num_elements_in_block;
113 }
114 else if (type == 4 || type == 11 || type == 29 || type == 30 || type == 31) // tet
115 {
116 cells_cols = std::max(cells_cols, 4);
117 num_els += e.num_elements_in_block;
118 }
119 else if (type == 5 || type == 12) // hex
120 {
121 cells_cols = std::max(cells_cols, 8);
122 num_els += e.num_elements_in_block;
123 }
124 else if (type == 6) // prism
125 {
126 cells_cols = std::max(cells_cols, 6);
127 num_els += e.num_elements_in_block;
128 }
129 else if (type == 7) // pyramid
130 {
131 cells_cols = std::max(cells_cols, 5);
132 num_els += e.num_elements_in_block;
133 }
134 }
135 assert(cells_cols > 0);
136
137 std::unordered_map<int, int> entity_tag_to_physical_tag;
138 if (dim == 2)
139 map_entity_tag_to_physical_tag(spec.entities.surfaces, entity_tag_to_physical_tag);
140 else
141 map_entity_tag_to_physical_tag(spec.entities.volumes, entity_tag_to_physical_tag);
142
143 cells.resize(num_els, cells_cols);
144 cells.setConstant(-1);
145 body_ids.resize(num_els);
146 elements.resize(num_els);
147 weights.resize(num_els);
148 int cell_index = 0;
149 for (const auto &e : els.entity_blocks)
150 {
151 if (e.entity_dim != dim)
152 continue;
153 const int type = e.element_type;
154 if (type == 2 || type == 9 || type == 21 || type == 23 || type == 25 || type == 3 || type == 10 || type == 4 || type == 11 || type == 29 || type == 30 || type == 31 || type == 5 || type == 12 || type == 6 || type == 7)
155 {
156 const size_t n_nodes = mshio::nodes_per_element(type);
157 int local_cells_cols = -1;
158
159 if (type == 2 || type == 9 || type == 21 || type == 23 || type == 25) // tri
160 local_cells_cols = 3;
161 else if (type == 3 || type == 10) // quad
162 local_cells_cols = 4;
163 else if (type == 4 || type == 11 || type == 29 || type == 30 || type == 31) // tet
164 local_cells_cols = 4;
165 else if (type == 5 || type == 12) // hex
166 local_cells_cols = 8;
167 else if (type == 6) // prism
168 local_cells_cols = 6;
169 else if (type == 7) // pyramid
170 local_cells_cols = 5;
171
172 for (int i = 0; i < e.data.size(); i += (n_nodes + 1))
173 {
174 int index = 0;
175 for (int j = i + 1; j <= i + local_cells_cols; ++j)
176 {
177 const int v_index = tag_to_index[e.data[j]];
178 assert(v_index < n_vertices);
179 cells(cell_index, index++) = v_index;
180 }
181
182 for (int j = i + 1; j < i + 1 + n_nodes; ++j)
183 {
184 const int v_index = tag_to_index[e.data[j]];
185 assert(v_index < n_vertices);
186 elements[cell_index].push_back(v_index);
187 }
188
189 const auto &it = entity_tag_to_physical_tag.find(e.entity_tag);
190 body_ids[cell_index] =
191 it != entity_tag_to_physical_tag.end() ? it->second : 0;
192
193 ++cell_index;
194 }
195 }
196 }
197
198 node_data.resize(spec.node_data.size());
199 int i = 0;
200 for (const auto &data : spec.node_data)
201 {
202 for (const auto &str : data.header.string_tags)
203 node_data_name.push_back(str);
204
205 for (const auto &entry : data.entries)
206 for (const auto &d : entry.data)
207 node_data[i].push_back(d);
208
209 i++;
210 }
211
212 // std::ifstream infile(path.c_str());
213
214 // std::string line;
215
216 // int phase = -1;
217 // int line_number = -1;
218 // bool size_read = false;
219
220 // int n_triangles = 0;
221 // int n_tets = 0;
222
223 // std::vector<std::vector<double>> all_elements;
224
225 // while (std::getline(infile, line))
226 // {
227 // line = StringUtils::trim(line);
228 // ++line_number;
229
230 // if (line.empty())
231 // continue;
232
233 // if (line[0] == '$')
234 // {
235 // if (line.substr(1, 3) == "End")
236 // phase = -1;
237 // else
238 // {
239 // const auto header = line.substr(1);
240
241 // if (header.find("MeshFormat") == 0)
242 // phase = 0;
243 // else if (header.find("Nodes") == 0)
244 // phase = 1;
245 // else if (header.find("Elements") == 0)
246 // phase = 2;
247 // else
248 // {
249 // logger().debug("{}: [Warning] ignoring {}", line_number, header);
250 // phase = -1;
251 // }
252 // }
253
254 // size_read = false;
255
256 // continue;
257 // }
258
259 // if (phase == -1)
260 // continue;
261
262 // std::istringstream iss(line);
263 // //header
264 // if (phase == 0)
265 // {
266 // double version_number;
267 // int file_type;
268 // int data_size;
269
270 // iss >> version_number >> file_type >> data_size;
271
272 // assert(version_number == 2.2);
273 // assert(file_type == 0);
274 // assert(data_size == 8);
275 // }
276 // //coordiantes
277 // else if (phase == 1)
278 // {
279 // if (!size_read)
280 // {
281 // int n_vertices;
282 // iss >> n_vertices;
283 // vertices.resize(n_vertices, 3);
284 // size_read = true;
285 // }
286 // else
287 // {
288 // int node_number;
289 // double x_coord, y_coord, z_coord;
290
291 // iss >> node_number >> x_coord >> y_coord >> z_coord;
292 // //node_numbers starts with 1
293 // vertices.row(node_number - 1) << x_coord, y_coord, z_coord;
294 // }
295 // }
296 // //elements
297 // else if (phase == 2)
298 // {
299 // if (!size_read)
300 // {
301 // int number_of_elements;
302 // iss >> number_of_elements;
303 // all_elements.resize(number_of_elements);
304 // size_read = true;
305 // }
306 // else
307 // {
308 // int elm_number, elm_type, number_of_tags;
309
310 // iss >> elm_number >> elm_type >> number_of_tags;
311
312 // //9-node third order incomplete triangle
313 // assert(elm_type != 20);
314
315 // //12-node fourth order incomplete triangle
316 // assert(elm_type != 22);
317
318 // //15-node fifth order incomplete triangle
319 // assert(elm_type != 24);
320
321 // //21-node fifth order complete triangle
322 // assert(elm_type != 25);
323
324 // //56-node fifth order tetrahedron
325 // assert(elm_type != 31);
326
327 // //60 is the new rational element
328 // if (elm_type == 2 || elm_type == 9 || elm_type == 21 || elm_type == 23 || elm_type == 60)
329 // ++n_triangles;
330 // else if (elm_type == 4 || elm_type == 11 || elm_type == 29 || elm_type == 30)
331 // ++n_tets;
332
333 // //skipping tags
334 // for (int i = 0; i < number_of_tags; ++i)
335 // {
336 // int tmp;
337 // iss >> tmp;
338 // }
339
340 // auto &node_list = all_elements[elm_number - 1];
341 // node_list.push_back(elm_type);
342
343 // while (iss.good())
344 // {
345 // double tmp;
346 // iss >> tmp;
347 // node_list.push_back(tmp);
348 // }
349 // }
350 // }
351 // else
352 // {
353 // assert(false);
354 // }
355 // }
356
357 // int index = 0;
358 // if (n_tets == 0)
359 // {
360 // elements.resize(n_triangles);
361 // weights.resize(n_triangles);
362 // cells.resize(n_triangles, 3);
363
364 // for (const auto &els : all_elements)
365 // {
366 // const int elm_type = els[0];
367 // if (elm_type != 2 && elm_type != 9 && elm_type != 21 && elm_type != 23 && elm_type != 60)
368 // continue;
369
370 // auto &el = elements[index];
371 // auto &wh = weights[index];
372 // for (size_t i = 1; i < (elm_type == 60 ? 7 : els.size()); ++i)
373 // el.push_back(int(els[i]) - 1);
374 // if (elm_type == 60)
375 // {
376 // for (size_t i = 7; i < els.size(); ++i)
377 // wh.push_back(els[i]);
378
379 // assert(wh.size() == el.size());
380 // assert(wh.size() == 6);
381 // }
382
383 // cells.row(index) << el[0], el[1], el[2];
384
385 // ++index;
386 // }
387 // }
388 // else
389 // {
390 // elements.resize(n_tets);
391 // weights.resize(n_tets);
392 // cells.resize(n_tets, 4);
393
394 // for (const auto &els : all_elements)
395 // {
396 // const int elm_type = els[0];
397 // if (elm_type != 4 && elm_type != 11 && elm_type != 29 && elm_type != 30)
398 // continue;
399
400 // auto &el = elements[index];
401 // auto &wh = weights[index];
402 // for (size_t i = 1; i < els.size(); ++i)
403 // el.push_back(els[i] - 1);
404
405 // cells.row(index) << el[0], el[1], el[2], el[3];
406 // ++index;
407 // }
408 // }
409
410 return true;
411 }
412} // namespace polyfem::io
polyfem::io::MshReader::load
static bool load(const std::string &path, Eigen::MatrixXd &vertices, Eigen::MatrixXi &cells, std::vector< std::vector< int > > &elements, std::vector< std::vector< double > > &weights, std::vector< int > &body_ids)
Definition MshReader.cpp:29
polyfem::io::map_entity_tag_to_physical_tag
void map_entity_tag_to_physical_tag(const std::vector< Entity > &entities, std::unordered_map< int, int > &entity_tag_to_physical_tag)
Definition MshReader.cpp:18
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44