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