1
2 // This file is part of libigl, a simple c++ geometry processing library.
3 //
4 // Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla Public License
7 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
8 // obtain one at http://mozilla.org/MPL/2.0/.
9 //
10 #include "readMSH.h"
11 #include <iostream>
12 #include <sstream>
13 #include <fstream>
14 #include <vector>
15 #include <map>
16
17 template <
18 typename DerivedV,
19 typename DerivedT>
readMSH(const std::string & filename,Eigen::PlainObjectBase<DerivedV> & V,Eigen::PlainObjectBase<DerivedT> & T)20 IGL_INLINE bool igl::readMSH(
21 const std::string & filename,
22 Eigen::PlainObjectBase<DerivedV> & V,
23 Eigen::PlainObjectBase<DerivedT> & T)
24 {
25 // https://github.com/Yixin-Hu/TetWild/blob/master/pymesh/MshSaver.cpp
26 // Original copyright: /* This file is part of PyMesh. Copyright (c) 2015 by Qingnan Zhou */
27 typedef typename DerivedV::Scalar Float;
28 typedef Eigen::Matrix<Float,Eigen::Dynamic,1> VectorF;
29 typedef Eigen::Matrix<int,Eigen::Dynamic,1> VectorI;
30 typedef std::map<std::string, VectorF> FieldMap;
31 typedef std::vector<std::string> FieldNames;
32 VectorF m_nodes;
33 VectorI m_elements;
34 FieldMap m_node_fields;
35 FieldMap m_element_fields;
36
37 bool m_binary;
38 size_t m_data_size;
39 size_t m_nodes_per_element;
40 size_t m_element_type;
41 std::ifstream fin(filename.c_str(), std::ios::in | std::ios::binary);
42 if (!fin.is_open())
43 {
44 std::stringstream err_msg;
45 err_msg << "failed to open file \"" << filename << "\"";
46 return false;
47 }
48 // Parse header
49 std::string buf;
50 double version;
51 int type;
52 fin >> buf;
53 const auto invalid_format = []()->bool
54 {
55 assert(false && "Invalid format");
56 return false;
57 };
58 const auto not_implemented = []()->bool
59 {
60 assert(false && "Not implemented");
61 return false;
62 };
63 if (buf != "$MeshFormat") { return invalid_format(); }
64
65 fin >> version >> type >> m_data_size;
66 m_binary = (type == 1);
67
68 // Some sanity check.
69 if (m_data_size != 8) {
70 std::cerr << "Error: data size must be 8 bytes." << std::endl;
71 return not_implemented();
72 }
73 if (sizeof(int) != 4) {
74 std::cerr << "Error: code must be compiled with int size 4 bytes." << std::endl;
75 return not_implemented();
76 }
77 const auto eat_white_space = [](std::ifstream& fin)
78 {
79 char next = fin.peek();
80 while (next == '\n' || next == ' ' || next == '\t' || next == '\r')
81 {
82 fin.get();
83 next = fin.peek();
84 }
85 };
86
87 // Read in extra info from binary header.
88 if (m_binary) {
89 int one;
90 eat_white_space(fin);
91 fin.read(reinterpret_cast<char*>(&one), sizeof(int));
92 if (one != 1) {
93 std::cerr << "Warning: binary msh file " << filename
94 << " is saved with different endianness than this machine."
95 << std::endl;
96 return not_implemented();
97 }
98 }
99
100 fin >> buf;
101 if (buf != "$EndMeshFormat") { return not_implemented(); }
102
103 const auto num_nodes_per_elem_type = [](int elem_type)->int
104 {
105 size_t nodes_per_element = 0;
106 switch (elem_type) {
107 case 2:
108 nodes_per_element = 3; // Triangle
109 break;
110 case 3:
111 nodes_per_element = 4; // Quad
112 break;
113 case 4:
114 nodes_per_element = 4; // Tet
115 break;
116 case 5:
117 nodes_per_element = 8; // hexahedron
118 break;
119 default:
120 assert(false && "not implemented");
121 nodes_per_element = -1;
122 break;
123 }
124 return nodes_per_element;
125 };
126
127 const auto parse_nodes = [&](std::ifstream& fin)
128 {
129 size_t num_nodes;
130 fin >> num_nodes;
131 m_nodes.resize(num_nodes*3);
132
133 if (m_binary) {
134 size_t num_bytes = (4+3*m_data_size) * num_nodes;
135 char* data = new char[num_bytes];
136 eat_white_space(fin);
137 fin.read(data, num_bytes);
138
139 for (size_t i=0; i<num_nodes; i++) {
140 int node_idx = *reinterpret_cast<int*> (&data[i*(4+3*m_data_size)]) - 1;
141 m_nodes[node_idx*3] = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4]);
142 m_nodes[node_idx*3+1] = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4 + m_data_size]);
143 m_nodes[node_idx*3+2] = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4 + 2*m_data_size]);
144 }
145
146 delete [] data;
147 } else {
148 int node_idx;
149 for (size_t i=0; i<num_nodes; i++) {
150 fin >> node_idx;
151 node_idx -= 1;
152 fin >> m_nodes[node_idx*3]
153 >> m_nodes[node_idx*3+1]
154 >> m_nodes[node_idx*3+2];
155 }
156 }
157 };
158
159 const auto parse_elements = [&](std::ifstream& fin)
160 {
161 size_t num_elements;
162 fin >> num_elements;
163
164 // Tmp storage of elements;
165 std::vector<int> triangle_element_idx;
166 std::vector<int> triangle_elements;
167 std::vector<int> quad_element_idx;
168 std::vector<int> quad_elements;
169 std::vector<int> tet_element_idx;
170 std::vector<int> tet_elements;
171 std::vector<int> hex_element_idx;
172 std::vector<int> hex_elements;
173
174 auto get_element_storage = [&](int elem_type) -> std::vector<int>* {
175 switch (elem_type) {
176 default:
177 assert(false && "Unsupported element type encountered");
178 case 2:
179 return &triangle_elements;
180 case 3:
181 return &quad_elements;
182 case 4:
183 return &tet_elements;
184 case 5:
185 return &hex_elements;
186 };
187 };
188
189 auto get_element_idx_storage = [&](int elem_type) -> std::vector<int>* {
190 switch (elem_type) {
191 default:
192 assert(false && "Unsupported element type encountered");
193 case 2:
194 return &triangle_element_idx;
195 case 3:
196 return &quad_element_idx;
197 case 4:
198 return &tet_element_idx;
199 case 5:
200 return &hex_element_idx;
201 };
202 };
203
204 size_t nodes_per_element;
205 int glob_elem_type = -1;
206
207
208 if (m_binary)
209 {
210 eat_white_space(fin);
211 int elem_read = 0;
212 while (elem_read < num_elements) {
213 // Parse element header.
214 int elem_type, num_elems, num_tags;
215 fin.read((char*)&elem_type, sizeof(int));
216 fin.read((char*)&num_elems, sizeof(int));
217 fin.read((char*)&num_tags, sizeof(int));
218 nodes_per_element = num_nodes_per_elem_type(elem_type);
219 std::vector<int>& elements = *get_element_storage(elem_type);
220 std::vector<int>& element_idx = *get_element_idx_storage(elem_type);
221
222 for (size_t i=0; i<num_elems; i++) {
223 int elem_idx;
224 fin.read((char*)&elem_idx, sizeof(int));
225 elem_idx -= 1;
226 element_idx.push_back(elem_idx);
227
228 // Eat up tags.
229 for (size_t j=0; j<num_tags; j++) {
230 int tag;
231 fin.read((char*)&tag, sizeof(int));
232 }
233
234 // Element values.
235 for (size_t j=0; j<nodes_per_element; j++) {
236 int idx;
237 fin.read((char*)&idx, sizeof(int));
238 elements.push_back(idx-1);
239 }
240 }
241
242 elem_read += num_elems;
243 }
244 } else
245 {
246 for (size_t i=0; i<num_elements; i++) {
247 // Parse per element header
248 int elem_num, elem_type, num_tags;
249 fin >> elem_num >> elem_type >> num_tags;
250 for (size_t j=0; j<num_tags; j++) {
251 int tag;
252 fin >> tag;
253 }
254 nodes_per_element = num_nodes_per_elem_type(elem_type);
255 std::vector<int>& elements = *get_element_storage(elem_type);
256 std::vector<int>& element_idx = *get_element_idx_storage(elem_type);
257
258 elem_num -= 1;
259 element_idx.push_back(elem_num);
260
261 // Parse node idx.
262 for (size_t j=0; j<nodes_per_element; j++) {
263 int idx;
264 fin >> idx;
265 elements.push_back(idx-1); // msh index starts from 1.
266 }
267 }
268 }
269
270 auto copy_to_array = [&](
271 const std::vector<int>& elements,
272 const int nodes_per_element) {
273 const size_t num_elements = elements.size() / nodes_per_element;
274 if (elements.size() % nodes_per_element != 0) {
275 assert(false && "parsing element failed");
276 return;
277 }
278 m_elements.resize(elements.size());
279 std::copy(elements.begin(), elements.end(), m_elements.data());
280 m_nodes_per_element = nodes_per_element;
281 };
282
283 if (!tet_elements.empty()) {
284 copy_to_array(tet_elements, 4);
285 m_element_type = 4;
286 } else if (!hex_elements.empty()) {
287 copy_to_array(hex_elements, 8);
288 m_element_type = 5;
289 } else if (!triangle_elements.empty()) {
290 copy_to_array(triangle_elements, 3);
291 m_element_type = 2;
292 } else if (!quad_elements.empty()) {
293 copy_to_array(quad_elements, 4);
294 m_element_type = 3;
295 } else {
296 // 0 elements, use triangle by default.
297 m_element_type = 2;
298 }
299 };
300 const auto parse_element_field = [&](std::ifstream& fin)
301 {
302 size_t num_string_tags;
303 size_t num_real_tags;
304 size_t num_int_tags;
305
306 fin >> num_string_tags;
307 std::string* str_tags = new std::string[num_string_tags];
308 for (size_t i=0; i<num_string_tags; i++) {
309 eat_white_space(fin);
310 if (fin.peek() == '\"') {
311 // Handle field name between quoates.
312 char buf[128];
313 fin.get(); // remove the quote at the beginning.
314 fin.getline(buf, 128, '\"');
315 str_tags[i] = std::string(buf);
316 } else {
317 fin >> str_tags[i];
318 }
319 }
320
321 fin >> num_real_tags;
322 Float* real_tags = new Float[num_real_tags];
323 for (size_t i=0; i<num_real_tags; i++)
324 fin >> real_tags[i];
325
326 fin >> num_int_tags;
327 int* int_tags = new int[num_int_tags];
328 for (size_t i=0; i<num_int_tags; i++)
329 fin >> int_tags[i];
330
331 if (num_string_tags <= 0 || num_int_tags <= 2) { assert(false && "invalid format"); return; }
332 std::string fieldname = str_tags[0];
333 int num_components = int_tags[1];
334 int num_entries = int_tags[2];
335 VectorF field(num_entries * num_components);
336
337 delete [] str_tags;
338 delete [] real_tags;
339 delete [] int_tags;
340
341 if (m_binary) {
342 size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
343 char* data = new char[num_bytes];
344 eat_white_space(fin);
345 fin.read(data, num_bytes);
346 for (size_t i=0; i<num_entries; i++) {
347 int elem_idx = *reinterpret_cast<int*>(&data[i*(4+num_components*m_data_size)]);
348 elem_idx -= 1;
349 size_t base_idx = i*(4+num_components*m_data_size) + 4;
350 for (size_t j=0; j<num_components; j++) {
351 field[elem_idx * num_components + j] = *reinterpret_cast<Float*>(&data[base_idx+j*m_data_size]);
352 }
353 }
354 delete [] data;
355 } else {
356 int elem_idx;
357 for (size_t i=0; i<num_entries; i++) {
358 fin >> elem_idx;
359 elem_idx -= 1;
360 for (size_t j=0; j<num_components; j++) {
361 fin >> field[elem_idx * num_components + j];
362 }
363 }
364 }
365
366 m_element_fields[fieldname] = field;
367 };
368
369 const auto parse_node_field = [&](std::ifstream& fin)
370 {
371 size_t num_string_tags;
372 size_t num_real_tags;
373 size_t num_int_tags;
374
375 fin >> num_string_tags;
376 std::string* str_tags = new std::string[num_string_tags];
377 for (size_t i=0; i<num_string_tags; i++) {
378 eat_white_space(fin);
379 if (fin.peek() == '\"') {
380 // Handle field name between quoates.
381 char buf[128];
382 fin.get(); // remove the quote at the beginning.
383 fin.getline(buf, 128, '\"');
384 str_tags[i] = std::string(buf);
385 } else {
386 fin >> str_tags[i];
387 }
388 }
389
390 fin >> num_real_tags;
391 Float* real_tags = new Float[num_real_tags];
392 for (size_t i=0; i<num_real_tags; i++)
393 fin >> real_tags[i];
394
395 fin >> num_int_tags;
396 int* int_tags = new int[num_int_tags];
397 for (size_t i=0; i<num_int_tags; i++)
398 fin >> int_tags[i];
399
400 if (num_string_tags <= 0 || num_int_tags <= 2) { assert(false && "invalid format"); return; }
401 std::string fieldname = str_tags[0];
402 int num_components = int_tags[1];
403 int num_entries = int_tags[2];
404 VectorF field(num_entries * num_components);
405
406 delete [] str_tags;
407 delete [] real_tags;
408 delete [] int_tags;
409
410 if (m_binary) {
411 size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
412 char* data = new char[num_bytes];
413 eat_white_space(fin);
414 fin.read(data, num_bytes);
415 for (size_t i=0; i<num_entries; i++) {
416 int node_idx = *reinterpret_cast<int*>(&data[i*(4+num_components*m_data_size)]);
417 node_idx -= 1;
418 size_t base_idx = i*(4+num_components*m_data_size) + 4;
419 for (size_t j=0; j<num_components; j++) {
420 field[node_idx * num_components + j] = *reinterpret_cast<Float*>(&data[base_idx+j*m_data_size]);
421 }
422 }
423 delete [] data;
424 } else {
425 int node_idx;
426 for (size_t i=0; i<num_entries; i++) {
427 fin >> node_idx;
428 node_idx -= 1;
429 for (size_t j=0; j<num_components; j++) {
430 fin >> field[node_idx * num_components + j];
431 }
432 }
433 }
434
435 m_node_fields[fieldname] = field;
436 };
437 const auto parse_unknown_field = [](std::ifstream& fin,
438 const std::string& fieldname)
439 {
440 std::cerr << "Warning: \"" << fieldname << "\" not supported yet. Ignored." << std::endl;
441 std::string endmark = fieldname.substr(0,1) + "End"
442 + fieldname.substr(1,fieldname.size()-1);
443
444 std::string buf("");
445 while (buf != endmark && !fin.eof()) {
446 fin >> buf;
447 }
448 };
449
450
451 while (!fin.eof()) {
452 buf.clear();
453 fin >> buf;
454 if (buf == "$Nodes") {
455 parse_nodes(fin);
456 fin >> buf;
457 if (buf != "$EndNodes") { return invalid_format(); }
458 } else if (buf == "$Elements") {
459 parse_elements(fin);
460 fin >> buf;
461 if (buf != "$EndElements") { return invalid_format(); }
462 } else if (buf == "$NodeData") {
463 parse_node_field(fin);
464 fin >> buf;
465 if (buf != "$EndNodeData") { return invalid_format(); }
466 } else if (buf == "$ElementData") {
467 parse_element_field(fin);
468 fin >> buf;
469 if (buf != "$EndElementData") { return invalid_format(); }
470 } else if (fin.eof()) {
471 break;
472 } else {
473 parse_unknown_field(fin, buf);
474 }
475 }
476 fin.close();
477 V.resize(m_nodes.rows()/3,3);
478 for (int i = 0; i < m_nodes.rows() / 3; i++)
479 {
480 for (int j = 0; j < 3; j++)
481 {
482 V(i,j) = m_nodes(i * 3 + j);
483 }
484 }
485 int ss = num_nodes_per_elem_type(m_element_type);
486 T.resize(m_elements.rows()/ss,ss);
487 for (int i = 0; i < m_elements.rows() / ss; i++)
488 {
489 for (int j = 0; j < ss; j++)
490 {
491 T(i, j) = m_elements(i * ss + j);
492 }
493 }
494 return true;
495 }
496
497
498 #ifdef IGL_STATIC_LIBRARY
499 // Explicit template instantiation
500 #endif
501