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
104 if (type == 2 || type == 9 || type == 21 || type == 23 || type == 25) // tri
105 {
106 assert(cells_cols == -1 || cells_cols == 3);
107 cells_cols = 3;
108 num_els += e.num_elements_in_block;
109 }
110 else if (type == 3 || type == 10) // quad
111 {
112 assert(cells_cols == -1 || cells_cols == 4);
113 cells_cols = 4;
114 num_els += e.num_elements_in_block;
115 }
116 else if (type == 4 || type == 11 || type == 29 || type == 30 || type == 31) // tet
117 {
118 assert(cells_cols == -1 || cells_cols == 4);
119 cells_cols = 4;
120 num_els += e.num_elements_in_block;
121 }
122 else if (type == 5 || type == 12) // hex
123 {
124 assert(cells_cols == -1 || cells_cols == 8);
125 cells_cols = 8;
126 num_els += e.num_elements_in_block;
127 }
128 else if (type == 6) // prism
129 {
130 assert(cells_cols == -1 || cells_cols == 6);
131 cells_cols = 6;
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 body_ids.resize(num_els);
145 elements.resize(num_els);
146 weights.resize(num_els);
147 int cell_index = 0;
148 for (const auto &e : els.entity_blocks)
149 {
150 if (e.entity_dim != dim)
151 continue;
152 const int type = e.element_type;
153 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)
154 {
155 const size_t n_nodes = mshio::nodes_per_element(type);
156 for (int i = 0; i < e.data.size(); i += (n_nodes + 1))
157 {
158 int index = 0;
159 for (int j = i + 1; j <= i + cells_cols; ++j)
160 {
161 const int v_index = tag_to_index[e.data[j]];
162 assert(v_index < n_vertices);
163 cells(cell_index, index++) = v_index;
164 }
165
166 for (int j = i + 1; j < i + 1 + n_nodes; ++j)
167 {
168 const int v_index = tag_to_index[e.data[j]];
169 assert(v_index < n_vertices);
170 elements[cell_index].push_back(v_index);
171 }
172
173 const auto &it = entity_tag_to_physical_tag.find(e.entity_tag);
174 body_ids[cell_index] =
175 it != entity_tag_to_physical_tag.end() ? it->second : 0;
176
177 ++cell_index;
178 }
179 }
180 }
181
182 node_data.resize(spec.node_data.size());
183 int i = 0;
184 for (const auto &data : spec.node_data)
185 {
186 for (const auto &str : data.header.string_tags)
187 node_data_name.push_back(str);
188
189 for (const auto &entry : data.entries)
190 for (const auto &d : entry.data)
191 node_data[i].push_back(d);
192
193 i++;
194 }
195
196 // std::ifstream infile(path.c_str());
197
198 // std::string line;
199
200 // int phase = -1;
201 // int line_number = -1;
202 // bool size_read = false;
203
204 // int n_triangles = 0;
205 // int n_tets = 0;
206
207 // std::vector<std::vector<double>> all_elements;
208
209 // while (std::getline(infile, line))
210 // {
211 // line = StringUtils::trim(line);
212 // ++line_number;
213
214 // if (line.empty())
215 // continue;
216
217 // if (line[0] == '$')
218 // {
219 // if (line.substr(1, 3) == "End")
220 // phase = -1;
221 // else
222 // {
223 // const auto header = line.substr(1);
224
225 // if (header.find("MeshFormat") == 0)
226 // phase = 0;
227 // else if (header.find("Nodes") == 0)
228 // phase = 1;
229 // else if (header.find("Elements") == 0)
230 // phase = 2;
231 // else
232 // {
233 // logger().debug("{}: [Warning] ignoring {}", line_number, header);
234 // phase = -1;
235 // }
236 // }
237
238 // size_read = false;
239
240 // continue;
241 // }
242
243 // if (phase == -1)
244 // continue;
245
246 // std::istringstream iss(line);
247 // //header
248 // if (phase == 0)
249 // {
250 // double version_number;
251 // int file_type;
252 // int data_size;
253
254 // iss >> version_number >> file_type >> data_size;
255
256 // assert(version_number == 2.2);
257 // assert(file_type == 0);
258 // assert(data_size == 8);
259 // }
260 // //coordiantes
261 // else if (phase == 1)
262 // {
263 // if (!size_read)
264 // {
265 // int n_vertices;
266 // iss >> n_vertices;
267 // vertices.resize(n_vertices, 3);
268 // size_read = true;
269 // }
270 // else
271 // {
272 // int node_number;
273 // double x_coord, y_coord, z_coord;
274
275 // iss >> node_number >> x_coord >> y_coord >> z_coord;
276 // //node_numbers starts with 1
277 // vertices.row(node_number - 1) << x_coord, y_coord, z_coord;
278 // }
279 // }
280 // //elements
281 // else if (phase == 2)
282 // {
283 // if (!size_read)
284 // {
285 // int number_of_elements;
286 // iss >> number_of_elements;
287 // all_elements.resize(number_of_elements);
288 // size_read = true;
289 // }
290 // else
291 // {
292 // int elm_number, elm_type, number_of_tags;
293
294 // iss >> elm_number >> elm_type >> number_of_tags;
295
296 // //9-node third order incomplete triangle
297 // assert(elm_type != 20);
298
299 // //12-node fourth order incomplete triangle
300 // assert(elm_type != 22);
301
302 // //15-node fifth order incomplete triangle
303 // assert(elm_type != 24);
304
305 // //21-node fifth order complete triangle
306 // assert(elm_type != 25);
307
308 // //56-node fifth order tetrahedron
309 // assert(elm_type != 31);
310
311 // //60 is the new rational element
312 // if (elm_type == 2 || elm_type == 9 || elm_type == 21 || elm_type == 23 || elm_type == 60)
313 // ++n_triangles;
314 // else if (elm_type == 4 || elm_type == 11 || elm_type == 29 || elm_type == 30)
315 // ++n_tets;
316
317 // //skipping tags
318 // for (int i = 0; i < number_of_tags; ++i)
319 // {
320 // int tmp;
321 // iss >> tmp;
322 // }
323
324 // auto &node_list = all_elements[elm_number - 1];
325 // node_list.push_back(elm_type);
326
327 // while (iss.good())
328 // {
329 // double tmp;
330 // iss >> tmp;
331 // node_list.push_back(tmp);
332 // }
333 // }
334 // }
335 // else
336 // {
337 // assert(false);
338 // }
339 // }
340
341 // int index = 0;
342 // if (n_tets == 0)
343 // {
344 // elements.resize(n_triangles);
345 // weights.resize(n_triangles);
346 // cells.resize(n_triangles, 3);
347
348 // for (const auto &els : all_elements)
349 // {
350 // const int elm_type = els[0];
351 // if (elm_type != 2 && elm_type != 9 && elm_type != 21 && elm_type != 23 && elm_type != 60)
352 // continue;
353
354 // auto &el = elements[index];
355 // auto &wh = weights[index];
356 // for (size_t i = 1; i < (elm_type == 60 ? 7 : els.size()); ++i)
357 // el.push_back(int(els[i]) - 1);
358 // if (elm_type == 60)
359 // {
360 // for (size_t i = 7; i < els.size(); ++i)
361 // wh.push_back(els[i]);
362
363 // assert(wh.size() == el.size());
364 // assert(wh.size() == 6);
365 // }
366
367 // cells.row(index) << el[0], el[1], el[2];
368
369 // ++index;
370 // }
371 // }
372 // else
373 // {
374 // elements.resize(n_tets);
375 // weights.resize(n_tets);
376 // cells.resize(n_tets, 4);
377
378 // for (const auto &els : all_elements)
379 // {
380 // const int elm_type = els[0];
381 // if (elm_type != 4 && elm_type != 11 && elm_type != 29 && elm_type != 30)
382 // continue;
383
384 // auto &el = elements[index];
385 // auto &wh = weights[index];
386 // for (size_t i = 1; i < els.size(); ++i)
387 // el.push_back(els[i] - 1);
388
389 // cells.row(index) << el[0], el[1], el[2], el[3];
390 // ++index;
391 // }
392 // }
393
394 return true;
395 }
396} // 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