Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "meshio.hpp"
5 :
6 : #include <unordered_map>
7 : #include <vector>
8 : #include <mfem.hpp>
9 :
10 : #define GMSH_BIN // Use binary Gmsh format
11 :
12 : namespace palace
13 : {
14 :
15 : namespace
16 : {
17 :
18 0 : inline int ElemTypeComsol(const std::string &type)
19 : {
20 0 : if (!type.compare("tri")) // 3-node triangle
21 : {
22 : return 2;
23 : }
24 0 : if (!type.compare("quad")) // 4-node quadrangle
25 : {
26 : return 3;
27 : }
28 0 : if (!type.compare("tet")) // 4-node tetrahedron
29 : {
30 : return 4;
31 : }
32 0 : if (!type.compare("hex")) // 8-node hexahedron
33 : {
34 : return 5;
35 : }
36 0 : if (!type.compare("prism")) // 6-node prism
37 : {
38 : return 6;
39 : }
40 0 : if (!type.compare("pyr")) // 5-node pyramid
41 : {
42 : return 7;
43 : }
44 0 : if (!type.compare("tri2")) // 6-node triangle
45 : {
46 : return 9;
47 : }
48 0 : if (!type.compare("quad2")) // 9-node quadrangle
49 : {
50 : return 10;
51 : }
52 0 : if (!type.compare("tet2")) // 10-node tetrahedron
53 : {
54 : return 11;
55 : }
56 0 : if (!type.compare("hex2")) // 27-node hexahedron
57 : {
58 : return 12;
59 : }
60 0 : if (!type.compare("prism2")) // 18-node prism
61 : {
62 : return 13;
63 : }
64 0 : if (!type.compare("pyr2")) // 14-node pyramid
65 : {
66 0 : return 14;
67 : }
68 : return 0; // Skip this element type
69 : }
70 :
71 0 : inline int ElemTypeNastran(const std::string &type)
72 : {
73 : // Returns only the low-order type for a given keyword.
74 0 : if (!type.compare(0, 5, "CTRIA")) // 3-node triangle
75 : {
76 : return 2;
77 : }
78 0 : if (!type.compare(0, 5, "CQUAD")) // 4-node quadrangle
79 : {
80 : return 3;
81 : }
82 0 : if (!type.compare(0, 6, "CTETRA")) // 4-node tetrahedron
83 : {
84 : return 4;
85 : }
86 0 : if (!type.compare(0, 5, "CHEXA")) // 8-node hexahedron
87 : {
88 : return 5;
89 : }
90 0 : if (!type.compare(0, 6, "CPENTA")) // 6-node prism
91 : {
92 : return 6;
93 : }
94 0 : if (!type.compare(0, 6, "CPYRAM")) // 5-node pyramid
95 : {
96 0 : return 7;
97 : }
98 : return 0; // Skip this element type
99 : }
100 :
101 0 : inline int HOElemTypeNastran(const int lo_type, const int num_nodes)
102 : {
103 : // Get high-order element type for corresponding low-order type.
104 0 : if (lo_type == 2 && num_nodes > 3) // 6-node triangle
105 : {
106 0 : MFEM_VERIFY(num_nodes == 6, "Invalid high-order Nastran element!");
107 : return 9;
108 : }
109 0 : if (lo_type == 3)
110 : {
111 0 : if (num_nodes == 9) // 9-node quadrangle
112 : {
113 : return 10;
114 : }
115 0 : if (num_nodes == 8) // 8-node quadrangle
116 : {
117 : return 16;
118 : }
119 0 : MFEM_VERIFY(num_nodes == 4, "Invalid high-order Nastran element!");
120 : return lo_type;
121 : }
122 0 : if (lo_type == 4 && num_nodes > 4) // 10-node tetrahedron
123 : {
124 0 : MFEM_VERIFY(num_nodes == 10, "Invalid high-order Nastran element!");
125 : return 11;
126 : }
127 0 : if (lo_type == 5 && num_nodes > 8) // 20-node hexahedron
128 : {
129 0 : MFEM_VERIFY(num_nodes == 20, "Invalid high-order Nastran element!");
130 : return 17;
131 : }
132 0 : if (lo_type == 6 && num_nodes > 6) // 15-node prism
133 : {
134 0 : MFEM_VERIFY(num_nodes == 15, "Invalid high-order Nastran element!");
135 : return 18;
136 : }
137 0 : if (lo_type == 7 && num_nodes > 5) // 13-node pyramid
138 : {
139 0 : MFEM_VERIFY(num_nodes == 13, "Invalid high-order Nastran element!");
140 : return 19;
141 : }
142 : return lo_type;
143 : }
144 :
145 0 : inline int LOElemTypeGmsh(int ho_type)
146 : {
147 : if (ho_type == 9) // 6-node triangle
148 : {
149 0 : return 2;
150 : }
151 : if (ho_type == 10 || ho_type == 16) // 9- or 8-node quadrangle
152 : {
153 0 : return 3;
154 : }
155 : if (ho_type == 11) // 10-node tetrahedron
156 : {
157 0 : return 4;
158 : }
159 : if (ho_type == 12 || ho_type == 17) // 27- or 20-node hexahedron
160 : {
161 0 : return 5;
162 : }
163 : if (ho_type == 13 || ho_type == 18) // 18- or 15-node prism
164 : {
165 0 : return 6;
166 : }
167 : if (ho_type == 14 || ho_type == 19) // 14- or 13-node pyramid
168 : {
169 0 : return 7;
170 : }
171 : return ho_type;
172 : }
173 :
174 : constexpr int ElemNumNodes[] = {-1, // 2-node edge
175 : 3, 4, 4, 8, 6, 5,
176 : -1, // 3-node edge
177 : 6, 9, 10, 27, 18, 14,
178 : -1, // 1-node node
179 : 8, 20, 15, 13};
180 :
181 : // From COMSOL or Nastran to Gmsh ordering. See:
182 : // - https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
183 : // - https://tinyurl.com/yezswzfv
184 : // - https://tinyurl.com/4d32zxtn
185 : constexpr int SkipElem[] = {-1};
186 : constexpr int Msh3[] = {0, 1, 2};
187 : constexpr int Msh4[] = {0, 1, 2, 3};
188 : constexpr int Msh5[] = {0, 1, 2, 3, 4};
189 : constexpr int Msh6[] = {0, 1, 2, 3, 4, 5};
190 : constexpr int Msh8[] = {0, 1, 2, 3, 4, 5, 6, 7};
191 : constexpr int Msh9[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
192 :
193 : constexpr int MphQuad4[] = {0, 1, 3, 2};
194 : constexpr int MphHex8[] = {0, 1, 3, 2, 4, 5, 7, 6};
195 : constexpr int MphPyr5[] = {0, 1, 3, 2, 4};
196 : constexpr int MphTri6[] = {0, 1, 2, 3, 5, 4};
197 : constexpr int MphQuad9[] = {0, 1, 3, 2, 4, 7, 8, 5, 6};
198 : constexpr int MphTet10[] = {0, 1, 2, 3, 4, 6, 5, 7, 9, 8};
199 : constexpr int MphHex27[] = {0, 1, 3, 2, 4, 5, 7, 6, 8, 9, 20, 11, 13, 10,
200 : 21, 12, 22, 26, 23, 15, 24, 14, 16, 17, 25, 18, 19};
201 : constexpr int MphWdg18[] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8, 15, 10, 16, 17, 11, 12, 13, 14};
202 : constexpr int MphPyr14[] = {0, 1, 3, 2, 4, 5, 6, 13, 8, 10, 7, 9, 12, 11};
203 :
204 : constexpr int NasTet10[] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8};
205 : constexpr int NasHex20[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11,
206 : 13, 9, 10, 12, 14, 15, 16, 18, 19, 17};
207 : constexpr int NasWdg15[] = {0, 1, 2, 3, 4, 5, 6, 9, 7, 8, 10, 11, 12, 14, 13};
208 : constexpr int NasPyr13[] = {0, 1, 2, 3, 4, 5, 8, 10, 6, 7, 9, 11, 12};
209 :
210 : constexpr const int *ElemNodesComsol[] = {SkipElem, Msh3, MphQuad4, Msh4, MphHex8,
211 : Msh6, MphPyr5, SkipElem, MphTri6, MphQuad9,
212 : MphTet10, MphHex27, MphWdg18, MphPyr14, SkipElem,
213 : SkipElem, SkipElem, SkipElem, SkipElem};
214 : constexpr const int *ElemNodesNastran[] = {SkipElem, Msh3, Msh4, Msh4, Msh8,
215 : Msh6, Msh5, SkipElem, Msh6, Msh9,
216 : NasTet10, SkipElem, SkipElem, SkipElem, SkipElem,
217 : Msh8, NasHex20, NasWdg15, NasPyr13};
218 :
219 : // Get line, strip comments, leading/trailing whitespace. Should not be called if end of
220 : // file is expected.
221 0 : inline std::string GetLineComsol(std::ifstream &input)
222 : {
223 : std::string str;
224 0 : std::getline(input, str);
225 0 : MFEM_VERIFY(input, "Unexpected read failure parsing mesh file!");
226 : const auto pos = str.find_first_of('#');
227 0 : if (pos != std::string::npos)
228 : {
229 0 : str.erase(pos);
230 : }
231 : const auto start = str.find_first_not_of(" \t");
232 0 : if (start == std::string::npos)
233 : {
234 0 : return "";
235 : }
236 : const auto stop = str.find_last_not_of(" \t");
237 0 : return str.substr(start, stop - start + 1);
238 : }
239 :
240 0 : inline std::string GetLineNastran(std::ifstream &input)
241 : {
242 : std::string str;
243 0 : std::getline(input, str);
244 0 : MFEM_VERIFY(input.good(), "Unexpected read failure parsing mesh file!");
245 0 : str.erase(std::remove(str.begin(), str.end(), '\r'), str.end());
246 0 : return str[0] == '$' ? "" : str;
247 : }
248 :
249 : // COMSOL strings are parsed as an integer length followed by array of integers for the
250 : // string characters.
251 0 : inline std::string ReadStringComsol(std::istream &input)
252 : {
253 : int n;
254 : std::string str;
255 0 : input >> n >> str;
256 0 : return str;
257 : }
258 :
259 0 : inline std::string ReadStringComsolBinary(std::istream &input)
260 : {
261 : int n;
262 0 : input.read(reinterpret_cast<char *>(&n), sizeof(int));
263 0 : std::vector<int> vstr(n);
264 0 : input.read(reinterpret_cast<char *>(vstr.data()), (std::streamsize)(n * sizeof(int)));
265 0 : return std::string(vstr.begin(), vstr.end());
266 : }
267 :
268 : // Nastran has a special floating point format: "-7.-1" instead of "-7.E-01" or "2.3+2"
269 : // instead of "2.3E+02".
270 0 : inline double ConvertDoubleNastran(const std::string &str)
271 : {
272 : double d;
273 : try
274 : {
275 : d = std::stod(str);
276 : }
277 0 : catch (const std::invalid_argument &ia)
278 : {
279 0 : const std::size_t start = str.find_first_not_of(' ');
280 0 : MFEM_VERIFY(start != std::string::npos,
281 : "Invalid number conversion parsing Nastran mesh!")
282 0 : std::string fstr = str.substr(start);
283 0 : std::size_t pos = fstr.find('+', 1); // Skip leading +/- sign
284 0 : if (pos != std::string::npos)
285 : {
286 0 : fstr.replace(pos, 1, "E+");
287 : }
288 0 : else if ((pos = fstr.find('-', 1)) != std::string::npos)
289 : {
290 0 : fstr.replace(pos, 1, "E-");
291 : }
292 : d = std::stod(fstr);
293 0 : }
294 0 : return d;
295 : }
296 :
297 : inline void WriteNode(std::ostream &buffer, const int tag, const double *coord)
298 : {
299 : #if defined(GMSH_BIN)
300 0 : buffer.write(reinterpret_cast<const char *>(&tag), sizeof(int));
301 0 : buffer.write(reinterpret_cast<const char *>(coord), 3 * sizeof(double));
302 : // No newline for binary data.
303 : #else
304 : // Always 3D coordinates (user sets floating point format/precision on buffer).
305 : buffer << tag << ' ' << coord[0] << ' ' << coord[1] << ' ' << coord[2] << '\n';
306 : #endif
307 : }
308 :
309 0 : inline void WriteElement(std::ostream &buffer, const int tag, const int type,
310 : const int geom, const int nodes[])
311 : {
312 : #if defined(GMSH_BIN)
313 0 : const int data[3] = {tag, geom, geom};
314 0 : buffer.write(reinterpret_cast<const char *>(data), 3 * sizeof(int));
315 0 : buffer.write(reinterpret_cast<const char *>(nodes),
316 0 : (std::streamsize)(ElemNumNodes[type - 1] * sizeof(int)));
317 : // No newline for binary data.
318 : #else
319 : buffer << tag << ' ' << type << " 2 " << geom << ' ' << geom;
320 : for (int i = 0; i < ElemNumNodes[type - 1]; i++)
321 : {
322 : buffer << ' ' << nodes[i];
323 : }
324 : buffer << '\n';
325 : #endif
326 0 : }
327 :
328 0 : void WriteGmsh(std::ostream &buffer, const std::vector<double> &node_coords,
329 : const std::vector<int> &node_tags,
330 : const std::unordered_map<int, std::vector<int>> &elem_nodes,
331 : const bool use_lo_type)
332 : {
333 : // Write the Gmsh file header (version 2.2).
334 : buffer << "$MeshFormat\n2.2 "
335 : <<
336 : #if defined(GMSH_BIN)
337 : "1 " <<
338 : #else
339 : "0 " <<
340 : #endif
341 0 : sizeof(double) << '\n';
342 : #if defined(GMSH_BIN)
343 0 : const int one = 1;
344 0 : buffer.write(reinterpret_cast<const char *>(&one), sizeof(int));
345 0 : buffer << '\n';
346 : #endif
347 0 : buffer << "$EndMeshFormat\n";
348 :
349 : // Write mesh nodes.
350 0 : const int num_nodes = (int)node_coords.size() / 3;
351 0 : MFEM_VERIFY(num_nodes > 0 && node_coords.size() % 3 == 0,
352 : "Gmsh nodes should always be in 3D space!");
353 0 : buffer << "$Nodes\n" << num_nodes << '\n';
354 : {
355 0 : if (!node_tags.empty())
356 : {
357 : // Use input node tags which should be positive but don't need to be contiguous.
358 0 : MFEM_VERIFY(node_tags.size() == (std::size_t)num_nodes,
359 : "Invalid size for node tags!");
360 0 : for (int i = 0; i < num_nodes; i++)
361 : {
362 0 : WriteNode(buffer, node_tags[i], &node_coords[3 * i]);
363 : }
364 : }
365 : else
366 : {
367 : // Label nodes as contiguous starting at 1.
368 0 : for (int i = 0; i < num_nodes; i++)
369 : {
370 0 : WriteNode(buffer, i + 1, &node_coords[3 * i]);
371 : }
372 : }
373 : }
374 : #if defined(GMSH_BIN)
375 0 : buffer << '\n';
376 : #endif
377 0 : buffer << "$EndNodes\n";
378 :
379 : // Write mesh elements.
380 : int tot_num_elem = 0;
381 0 : for (const auto &[elem_type, nodes] : elem_nodes)
382 : {
383 0 : MFEM_VERIFY(elem_type > 0, "Invalid element type writing Gmsh elements!");
384 0 : const int &num_elem_nodes = ElemNumNodes[elem_type - 1];
385 0 : tot_num_elem += ((int)nodes.size()) / (num_elem_nodes + 1);
386 0 : MFEM_VERIFY(nodes.size() % (num_elem_nodes + 1) == 0,
387 : "Unexpected data size when writing elements!");
388 : }
389 0 : MFEM_VERIFY(tot_num_elem > 0, "No mesh elements parsed from COMSOL mesh file!");
390 0 : buffer << "$Elements\n" << tot_num_elem << '\n';
391 : {
392 : int tag = 1; // Global element tag
393 0 : for (const auto &[elem_type, nodes] : elem_nodes)
394 : {
395 0 : const int elem_type_w = use_lo_type ? LOElemTypeGmsh(elem_type) : elem_type;
396 0 : const int &num_elem_nodes = ElemNumNodes[elem_type - 1];
397 0 : const int num_elem = (int)nodes.size() / (num_elem_nodes + 1);
398 : #if defined(GMSH_BIN)
399 : // For binary output, write the element header for each type. Always have 2 tags
400 : // (physical + geometry).
401 0 : const int header[3] = {elem_type_w, num_elem, 2};
402 0 : buffer.write(reinterpret_cast<const char *>(header), 3 * sizeof(int));
403 : #endif
404 0 : for (int i = 0; i < num_elem; i++)
405 : {
406 0 : WriteElement(buffer, tag++, elem_type_w,
407 0 : nodes[i * (num_elem_nodes + 1)], // Geometry tag
408 0 : &nodes[i * (num_elem_nodes + 1) + 1]); // Element nodes
409 : }
410 : }
411 : }
412 : #if defined(GMSH_BIN)
413 0 : buffer << '\n';
414 : #endif
415 0 : buffer << "$EndElements\n";
416 0 : }
417 :
418 : } // namespace
419 :
420 : namespace mesh
421 : {
422 :
423 0 : void ConvertMeshComsol(const std::string &filename, std::ostream &buffer,
424 : bool remove_curvature)
425 : {
426 : // Read a COMSOL format mesh.
427 0 : const int comsol_bin = !filename.compare(filename.length() - 7, 7, ".mphbin") ||
428 0 : !filename.compare(filename.length() - 7, 7, ".MPHBIN");
429 0 : MFEM_VERIFY(!filename.compare(filename.length() - 7, 7, ".mphtxt") ||
430 : !filename.compare(filename.length() - 7, 7, ".MPHTXT") || comsol_bin,
431 : "Invalid file extension for COMSOL mesh format conversion!");
432 0 : std::ifstream input(filename);
433 0 : if (!input.is_open())
434 : {
435 0 : MFEM_ABORT("Unable to open mesh file \"" << filename << "\"!");
436 : }
437 :
438 : // Parse COMSOL header. COMSOL encodes strings as integer-string pairs where the integer
439 : // is the string length. It also allows for blank lines and other whitespace wherever in
440 : // the file.
441 : {
442 0 : int version[2] = {-1, -1};
443 0 : int num_tags = -1;
444 0 : int num_types = -1;
445 0 : if (!comsol_bin)
446 : {
447 0 : while (num_types < 0)
448 : {
449 0 : auto line = GetLineComsol(input);
450 0 : if (!line.empty())
451 : {
452 0 : std::istringstream sline(line);
453 0 : if (version[0] < 0)
454 : {
455 0 : sline >> version[0] >> version[1];
456 : }
457 0 : else if (num_tags < 0)
458 : {
459 0 : sline >> num_tags;
460 : int i = 0;
461 0 : while (i < num_tags)
462 : {
463 0 : if (!GetLineComsol(input).empty())
464 : {
465 0 : i++;
466 : }
467 : }
468 : }
469 0 : else if (num_types < 0)
470 : {
471 0 : sline >> num_types;
472 : int i = 0;
473 0 : while (i < num_types)
474 : {
475 0 : if (!GetLineComsol(input).empty())
476 : {
477 0 : i++;
478 : }
479 : }
480 : }
481 0 : }
482 : }
483 : }
484 : else
485 : {
486 0 : input.read(reinterpret_cast<char *>(version), 2 * sizeof(int));
487 0 : input.read(reinterpret_cast<char *>(&num_tags), sizeof(int));
488 : {
489 : int i = 0;
490 0 : while (i < num_tags)
491 : {
492 0 : ReadStringComsolBinary(input);
493 0 : i++;
494 : }
495 : }
496 0 : input.read(reinterpret_cast<char *>(&num_types), sizeof(int));
497 : {
498 : int i = 0;
499 0 : while (i < num_types)
500 : {
501 0 : ReadStringComsolBinary(input);
502 0 : i++;
503 : }
504 : }
505 : }
506 0 : MFEM_VERIFY(version[0] == 0 && version[1] == 1, "Invalid COMSOL file version!");
507 : }
508 :
509 : // Parse mesh objects until we get to the mesh. Currently only supports a single mesh
510 : // object in the file, and selections are ignored.
511 : while (true)
512 : {
513 0 : int object[3] = {-1, -1, -1};
514 : std::string object_class;
515 0 : if (!comsol_bin)
516 : {
517 0 : while (object_class.empty())
518 : {
519 0 : auto line = GetLineComsol(input);
520 0 : if (!line.empty())
521 : {
522 0 : std::istringstream sline(line);
523 0 : if (object[0] < 0)
524 : {
525 0 : sline >> object[0] >> object[1] >> object[2];
526 : }
527 0 : else if (object_class.empty())
528 : {
529 0 : object_class = ReadStringComsol(sline);
530 : }
531 0 : }
532 : }
533 : }
534 : else
535 : {
536 0 : input.read(reinterpret_cast<char *>(object), 3 * sizeof(int));
537 0 : object_class = ReadStringComsolBinary(input);
538 : }
539 0 : MFEM_VERIFY(object[0] == 0 && object[1] == 0 && object[2] == 1,
540 : "Invalid COMSOL object version!");
541 :
542 : // If yes, then ready to parse the mesh.
543 0 : if (!object_class.compare(0, 4, "Mesh"))
544 : {
545 : break;
546 : }
547 :
548 : // Otherwise, parse over the selection to the next object.
549 0 : MFEM_VERIFY(!object_class.compare(0, 9, "Selection"),
550 : "COMSOL mesh file only supports Mesh and Selection objects!");
551 0 : int version = -1;
552 : std::string label_str;
553 : std::string tag_str;
554 0 : int sdim = -1;
555 0 : int num_ent = -1;
556 0 : if (!comsol_bin)
557 : {
558 0 : while (num_ent < 0)
559 : {
560 0 : auto line = GetLineComsol(input);
561 0 : if (!line.empty())
562 : {
563 0 : std::istringstream sline(line);
564 0 : if (version < 0)
565 : {
566 0 : sline >> version;
567 : }
568 0 : else if (label_str.empty())
569 : {
570 0 : label_str = ReadStringComsol(sline);
571 : }
572 0 : else if (tag_str.empty())
573 : {
574 0 : tag_str = ReadStringComsol(sline);
575 : }
576 0 : else if (sdim < 0)
577 : {
578 0 : sline >> sdim;
579 : }
580 0 : else if (num_ent < 0)
581 : {
582 0 : sline >> num_ent;
583 : }
584 0 : }
585 : }
586 : }
587 : else
588 : {
589 0 : input.read(reinterpret_cast<char *>(&version), sizeof(int));
590 0 : label_str = ReadStringComsolBinary(input);
591 0 : tag_str = ReadStringComsolBinary(input);
592 0 : input.read(reinterpret_cast<char *>(&sdim), sizeof(int));
593 0 : input.read(reinterpret_cast<char *>(&num_ent), sizeof(int));
594 : }
595 :
596 : // Parse over the entities in the selection.
597 : int i = 0;
598 0 : if (!comsol_bin)
599 : {
600 0 : while (i < num_ent)
601 : {
602 0 : if (!GetLineComsol(input).empty())
603 : {
604 0 : i++;
605 : }
606 : }
607 : }
608 : else
609 : {
610 0 : while (i < num_ent)
611 : {
612 : int dummy;
613 0 : input.read(reinterpret_cast<char *>(&dummy), sizeof(int));
614 0 : i++;
615 : }
616 : }
617 : } // Repeat until Mesh is found
618 :
619 : // Parse the mesh object header.
620 0 : int sdim = -1;
621 0 : int num_nodes = -1;
622 0 : int nodes_start = -1;
623 : {
624 0 : int version = -1;
625 0 : if (!comsol_bin)
626 : {
627 0 : while (nodes_start < 0)
628 : {
629 0 : auto line = GetLineComsol(input);
630 0 : if (!line.empty())
631 : {
632 0 : std::istringstream sline(line);
633 0 : if (version < 0)
634 : {
635 0 : sline >> version;
636 : }
637 0 : else if (sdim < 0)
638 : {
639 0 : sline >> sdim;
640 : }
641 0 : else if (num_nodes < 0)
642 : {
643 0 : sline >> num_nodes;
644 : }
645 0 : else if (nodes_start < 0)
646 : {
647 0 : sline >> nodes_start;
648 : }
649 0 : }
650 : }
651 : }
652 : else
653 : {
654 0 : input.read(reinterpret_cast<char *>(&version), sizeof(int));
655 0 : input.read(reinterpret_cast<char *>(&sdim), sizeof(int));
656 0 : input.read(reinterpret_cast<char *>(&num_nodes), sizeof(int));
657 0 : input.read(reinterpret_cast<char *>(&nodes_start), sizeof(int));
658 : }
659 0 : MFEM_VERIFY(version == 4, "Only COMSOL files with Mesh version 4 are supported!");
660 0 : MFEM_VERIFY(sdim == 2 || sdim == 3,
661 : "COMSOL mesh nodes are required to be in 2D or 3D space!");
662 0 : MFEM_VERIFY(num_nodes > 0, "COMSOL mesh file contains no nodes!");
663 0 : MFEM_VERIFY(nodes_start >= 0, "COMSOL mesh nodes have a negative starting tag!");
664 : }
665 :
666 : // Parse mesh nodes.
667 : std::vector<double> node_coords;
668 : {
669 : // Gmsh nodes are always 3D, so initialize to 0.0 in case z-coordinate isn't set.
670 0 : node_coords.resize(3 * num_nodes, 0.0);
671 : int i = 0;
672 0 : if (!comsol_bin)
673 : {
674 0 : while (i < num_nodes)
675 : {
676 0 : auto line = GetLineComsol(input);
677 0 : if (!line.empty())
678 : {
679 0 : std::istringstream sline(line);
680 0 : for (int j = 0; j < sdim; j++)
681 : {
682 0 : sline >> node_coords[3 * i + j];
683 : }
684 0 : i++;
685 0 : }
686 : }
687 : }
688 : else
689 : {
690 : // Don't read as a single block in case sdim < 3.
691 0 : while (i < num_nodes)
692 : {
693 0 : input.read(reinterpret_cast<char *>(node_coords.data() + 3 * i),
694 0 : (std::streamsize)(sdim * sizeof(double)));
695 0 : i++;
696 : }
697 : }
698 : }
699 :
700 : // Parse mesh elements. Store for each element of each type: [geometry tag, [node tags]].
701 : std::unordered_map<int, std::vector<int>> elem_nodes;
702 : {
703 0 : int num_elem_types = -1;
704 0 : if (!comsol_bin)
705 : {
706 0 : while (num_elem_types < 0)
707 : {
708 0 : auto line = GetLineComsol(input);
709 0 : if (!line.empty())
710 : {
711 0 : std::istringstream sline(line);
712 0 : if (num_elem_types < 0)
713 : {
714 0 : sline >> num_elem_types;
715 : }
716 0 : }
717 : }
718 : }
719 : else
720 : {
721 0 : input.read(reinterpret_cast<char *>(&num_elem_types), sizeof(int));
722 : }
723 0 : MFEM_VERIFY(num_elem_types > 0, "COMSOL mesh file contains no elements!");
724 :
725 : int parsed_types = 0; // COMSOL groups elements by type in file
726 0 : int elem_type = -1;
727 0 : int num_elem_nodes = -1;
728 0 : int num_elem = -1;
729 0 : int num_elem_geom = -1;
730 : bool skip_type = false;
731 0 : while (parsed_types < num_elem_types)
732 : {
733 0 : if (!comsol_bin)
734 : {
735 0 : auto line = GetLineComsol(input);
736 0 : if (!line.empty())
737 : {
738 0 : std::istringstream sline(line);
739 0 : if (elem_type < 0)
740 : {
741 0 : auto elem_str = ReadStringComsol(sline);
742 0 : MFEM_VERIFY(!elem_str.empty(),
743 : "Unexpected empty element type found in COMSOL mesh file!");
744 0 : elem_type = ElemTypeComsol(elem_str);
745 0 : skip_type = (elem_type == 0);
746 0 : MFEM_VERIFY(skip_type || elem_nodes.find(elem_type) == elem_nodes.end(),
747 : "Duplicate element types found in COMSOL mesh file!");
748 : }
749 0 : else if (num_elem_nodes < 0)
750 : {
751 0 : sline >> num_elem_nodes;
752 0 : MFEM_VERIFY(num_elem_nodes > 0,
753 : "COMSOL element type " << elem_type << " has no nodes!");
754 0 : MFEM_VERIFY(skip_type || num_elem_nodes == ElemNumNodes[elem_type - 1],
755 : "Mismatch between COMSOL and Gmsh element types!");
756 : }
757 0 : else if (num_elem < 0)
758 : {
759 0 : sline >> num_elem;
760 0 : MFEM_VERIFY(num_elem > 0,
761 : "COMSOL mesh file has no elements of type " << elem_type << "!");
762 : std::vector<int> *data = nullptr;
763 0 : if (!skip_type)
764 : {
765 : data = &elem_nodes[elem_type];
766 0 : data->resize(num_elem * (num_elem_nodes + 1)); // Node tags + geometry tag
767 : }
768 :
769 : // Parse all element nodes.
770 : int i = 0;
771 0 : while (i < num_elem)
772 : {
773 0 : line = GetLineComsol(input);
774 0 : if (!line.empty())
775 : {
776 0 : if (!skip_type)
777 : {
778 0 : std::istringstream isline(line);
779 0 : for (int j = 0; j < num_elem_nodes; j++)
780 : {
781 : // Permute and reset to 1-based node tags.
782 0 : const int &p = ElemNodesComsol[elem_type - 1][j];
783 0 : isline >> (*data)[i * (num_elem_nodes + 1) + 1 + p];
784 0 : (*data)[i * (num_elem_nodes + 1) + 1 + p] += (1 - nodes_start);
785 : }
786 0 : }
787 0 : i++;
788 : }
789 : }
790 : }
791 0 : else if (num_elem_geom < 0)
792 : {
793 0 : sline >> num_elem_geom;
794 0 : MFEM_VERIFY(num_elem_geom == num_elem,
795 : "COMSOL mesh file should have geometry tags for all elements!");
796 : std::vector<int> *data = nullptr;
797 0 : if (!skip_type)
798 : {
799 0 : MFEM_VERIFY(elem_nodes.find(elem_type) != elem_nodes.end(),
800 : "Can't find expected element type!");
801 : data = &elem_nodes[elem_type];
802 0 : MFEM_VERIFY(data->size() == (std::size_t)num_elem * (num_elem_nodes + 1),
803 : "Unexpected element data size!");
804 : }
805 :
806 : // Parse all element geometry tags (stored at beginning of element nodes). For
807 : // geometric entites in < 3D, the exported COMSOL tags are 0-based and need
808 : // correcting to 1-based for Gmsh.
809 : int i = 0;
810 : const int geom_start =
811 0 : (elem_type < 4 || (elem_type > 7 && elem_type < 11)) ? 1 : 0;
812 0 : while (i < num_elem)
813 : {
814 0 : line = GetLineComsol(input);
815 0 : if (!line.empty())
816 : {
817 0 : if (!skip_type)
818 : {
819 0 : std::istringstream ssline(line);
820 0 : ssline >> (*data)[i * (num_elem_nodes + 1)];
821 0 : (*data)[i * (num_elem_nodes + 1)] += geom_start;
822 0 : }
823 0 : i++;
824 : }
825 : }
826 :
827 : // Debug
828 : // std::cout << "Finished parsing " << num_elem
829 : // << " elements with type " << elem_type
830 : // << " (parsed types " << parsed_types + 1 << ")\n";
831 :
832 : // Finished with this element type, on to the next.
833 0 : parsed_types++;
834 0 : elem_type = num_elem_nodes = num_elem = num_elem_geom = -1;
835 : skip_type = false;
836 : }
837 0 : }
838 : }
839 : else
840 : {
841 0 : auto elem_str = ReadStringComsolBinary(input);
842 0 : MFEM_VERIFY(!elem_str.empty(),
843 : "Unexpected empty element type found in COMSOL mesh file!");
844 0 : elem_type = ElemTypeComsol(elem_str);
845 : skip_type = (elem_type == 0);
846 0 : MFEM_VERIFY(skip_type || elem_nodes.find(elem_type) == elem_nodes.end(),
847 : "Duplicate element types found in COMSOL mesh file!");
848 0 : input.read(reinterpret_cast<char *>(&num_elem_nodes), sizeof(int));
849 0 : MFEM_VERIFY(num_elem_nodes > 0,
850 : "COMSOL element type " << elem_type << " has no nodes!");
851 0 : MFEM_VERIFY(skip_type || num_elem_nodes == ElemNumNodes[elem_type - 1],
852 : "Mismatch between COMSOL and Gmsh element types!");
853 :
854 : // Parse all element nodes.
855 0 : input.read(reinterpret_cast<char *>(&num_elem), sizeof(int));
856 0 : MFEM_VERIFY(num_elem > 0,
857 : "COMSOL mesh file has no elements of type " << elem_type << "!");
858 : std::vector<int> *data = nullptr;
859 0 : if (!skip_type)
860 : {
861 : data = &elem_nodes[elem_type];
862 0 : data->resize(num_elem * (num_elem_nodes + 1)); // Node tags + geometry tag
863 : }
864 : int i = 0;
865 0 : std::vector<int> nodes(num_elem_nodes);
866 0 : while (i < num_elem)
867 : {
868 0 : input.read(reinterpret_cast<char *>(nodes.data()),
869 0 : (std::streamsize)(num_elem_nodes * sizeof(int)));
870 0 : if (!skip_type)
871 : {
872 0 : for (int j = 0; j < num_elem_nodes; j++)
873 : {
874 : // Permute and reset to 1-based node tags.
875 0 : const int &p = ElemNodesComsol[elem_type - 1][j];
876 0 : (*data)[i * (num_elem_nodes + 1) + 1 + p] = nodes[j] + (1 - nodes_start);
877 : }
878 : }
879 0 : i++;
880 : }
881 :
882 : // Parse element geometry tags.
883 0 : input.read(reinterpret_cast<char *>(&num_elem_geom), sizeof(int));
884 0 : MFEM_VERIFY(num_elem_geom == num_elem,
885 : "COMSOL mesh file should have geometry tags for all elements!");
886 :
887 : i = 0;
888 0 : const int geom_start = (elem_type < 4 || (elem_type > 7 && elem_type < 11)) ? 1 : 0;
889 : int geom_tag;
890 0 : while (i < num_elem)
891 : {
892 0 : input.read(reinterpret_cast<char *>(&geom_tag), sizeof(int));
893 0 : if (!skip_type)
894 : {
895 0 : (*data)[i * (num_elem_nodes + 1)] = geom_tag + geom_start;
896 : }
897 0 : i++;
898 : }
899 :
900 : // Debug
901 : // std::cout << "Finished parsing " << num_elem
902 : // << " elements with type " << elem_type
903 : // << " (parsed types " << parsed_types + 1 << ")\n";
904 :
905 : // Finished with this element type, on to the next.
906 0 : parsed_types++;
907 0 : elem_type = num_elem_nodes = num_elem = num_elem_geom = -1;
908 : skip_type = false;
909 : }
910 : }
911 : }
912 :
913 : // Finalize input, write the Gmsh mesh.
914 0 : input.close();
915 : std::vector<int> dummy;
916 0 : WriteGmsh(buffer, node_coords, dummy, elem_nodes, remove_curvature);
917 0 : }
918 :
919 0 : void ConvertMeshNastran(const std::string &filename, std::ostream &buffer,
920 : bool remove_curvature)
921 : {
922 : // Read a Nastran/BDF format mesh.
923 0 : MFEM_VERIFY(!filename.compare(filename.length() - 4, 4, ".nas") ||
924 : !filename.compare(filename.length() - 4, 4, ".NAS") ||
925 : !filename.compare(filename.length() - 4, 4, ".bdf") ||
926 : !filename.compare(filename.length() - 4, 4, ".BDF"),
927 : "Invalid file extension for Nastran mesh format conversion!");
928 0 : std::ifstream input(filename);
929 0 : if (!input.is_open())
930 : {
931 0 : MFEM_ABORT("Unable to open mesh file \"" << filename << "\"!");
932 : }
933 : const int NASTRAN_CHUNK = 8; // NASTRAN divides row into 10 columns of 8 spaces
934 : const int MAX_CHUNK = 9; // Never read the 10-th chunk
935 :
936 : // Parse until bulk data starts.
937 : while (true)
938 : {
939 0 : auto line = GetLineNastran(input);
940 0 : if (line.length() > 0)
941 : {
942 0 : if (!line.compare(0, 10, "BEGIN BULK"))
943 : {
944 : break;
945 : }
946 : }
947 : }
948 :
949 : // Parse mesh nodes and elements. It is expected that node tags start at 1 and are
950 : // contiguous. Store for each element of each type: [geometry tag, [node tags]].
951 : std::vector<double> node_coords;
952 : std::vector<int> node_tags;
953 : std::unordered_map<int, std::vector<int>> elem_nodes;
954 : int elem_type;
955 : while (true)
956 : {
957 0 : auto line = GetLineNastran(input);
958 0 : if (line.length() > 0 && !input.eof())
959 : {
960 0 : if (!line.compare(0, 7, "ENDDATA"))
961 : {
962 : break; // Done parsing file
963 : }
964 0 : else if (!line.compare(0, 5, "GRID*"))
965 : {
966 : // Coordinates in long field format (8 + 16 * 4 + 8).
967 0 : auto next = GetLineNastran(input);
968 0 : MFEM_VERIFY(!next.empty(), "Unexpected empty line parsing Nastran!");
969 :
970 0 : node_tags.push_back(std::stoi(line.substr(1 * NASTRAN_CHUNK, 2 * NASTRAN_CHUNK)));
971 0 : node_coords.insert(
972 : node_coords.end(),
973 0 : {ConvertDoubleNastran(line.substr(5 * NASTRAN_CHUNK, 2 * NASTRAN_CHUNK)),
974 0 : ConvertDoubleNastran(line.substr(7 * NASTRAN_CHUNK, 2 * NASTRAN_CHUNK)),
975 0 : ConvertDoubleNastran(next.substr(1 * NASTRAN_CHUNK, 2 * NASTRAN_CHUNK))});
976 : }
977 0 : else if (!line.compare(0, 4, "GRID"))
978 : {
979 0 : if (line.find_first_of(',') != std::string::npos)
980 : {
981 : // Free field format (comma separated).
982 0 : std::istringstream sline(line);
983 :
984 : std::string word;
985 0 : std::getline(sline, word, ','); // Discard "GRID"
986 :
987 0 : std::getline(sline, word, ',');
988 0 : node_tags.push_back(std::stoi(word));
989 :
990 0 : std::getline(sline, word, ','); // Discard coordinate system
991 :
992 0 : std::getline(sline, word, ',');
993 0 : double x = ConvertDoubleNastran(word);
994 0 : std::getline(sline, word, ',');
995 0 : double y = ConvertDoubleNastran(word);
996 0 : std::getline(sline, word, ',');
997 0 : double z = ConvertDoubleNastran(word);
998 0 : node_coords.insert(node_coords.end(), {x, y, z});
999 0 : }
1000 : else
1001 : {
1002 : // Short format (10 * 8).
1003 0 : node_tags.push_back(std::stoi(line.substr(1 * NASTRAN_CHUNK, NASTRAN_CHUNK)));
1004 0 : node_coords.insert(
1005 : node_coords.end(),
1006 0 : {ConvertDoubleNastran(line.substr(3 * NASTRAN_CHUNK, NASTRAN_CHUNK)),
1007 0 : ConvertDoubleNastran(line.substr(4 * NASTRAN_CHUNK, NASTRAN_CHUNK)),
1008 0 : ConvertDoubleNastran(line.substr(5 * NASTRAN_CHUNK, NASTRAN_CHUNK))});
1009 : }
1010 : }
1011 0 : else if ((elem_type = ElemTypeNastran(line)))
1012 : {
1013 : // Prepare to parse the element ID and nodes.
1014 : const bool free = (line.find_first_of(',') != std::string::npos);
1015 :
1016 : // Get the element type, tag, and geometry attribute. Then get the element nodes on
1017 : // this line.
1018 : std::string elem_str;
1019 : // int elem_tag;
1020 : int geom_tag;
1021 : std::vector<int> nodes;
1022 : std::string word;
1023 0 : if (!free)
1024 : {
1025 0 : elem_str = line.substr(0 * NASTRAN_CHUNK, NASTRAN_CHUNK);
1026 0 : const std::size_t stop = elem_str.find_last_not_of(' ');
1027 0 : MFEM_VERIFY(stop != std::string::npos, "Invalid element type string!");
1028 0 : elem_str.resize(stop + 1);
1029 : // elem_tag = std::stoi(line.substr(1*NASTRAN_CHUNK, NASTRAN_CHUNK));
1030 0 : geom_tag = std::stoi(line.substr(2 * NASTRAN_CHUNK, NASTRAN_CHUNK));
1031 :
1032 : int i = 3;
1033 0 : while (i < MAX_CHUNK)
1034 : {
1035 0 : word = line.substr((i++) * NASTRAN_CHUNK, NASTRAN_CHUNK);
1036 0 : if (word.find_first_not_of(' ') == std::string::npos)
1037 : {
1038 : break;
1039 : }
1040 0 : nodes.push_back(std::stoi(word));
1041 : }
1042 : }
1043 : else
1044 : {
1045 0 : std::istringstream sline(line);
1046 0 : std::getline(sline, elem_str, ',');
1047 0 : std::getline(sline, word, ',');
1048 : // elem_tag = std::stoi(word);
1049 0 : std::getline(sline, word, ',');
1050 : geom_tag = std::stoi(word);
1051 :
1052 : int i = 3;
1053 0 : while (i < MAX_CHUNK)
1054 : {
1055 0 : std::getline(sline, word, ',');
1056 0 : if (word.find_first_not_of(' ') == std::string::npos)
1057 : {
1058 : break;
1059 : }
1060 0 : nodes.push_back(std::stoi(word));
1061 0 : i++;
1062 : }
1063 0 : }
1064 :
1065 : // Handle line continuation.
1066 0 : while (input.peek() == '+')
1067 : {
1068 0 : auto next = GetLineNastran(input);
1069 0 : MFEM_VERIFY(!next.empty(), "Unexpected empty line parsing Nastran!");
1070 :
1071 0 : if (!free)
1072 : {
1073 : int i = 1;
1074 0 : while (i < MAX_CHUNK)
1075 : {
1076 0 : word = next.substr((i++) * NASTRAN_CHUNK, NASTRAN_CHUNK);
1077 0 : if (word.find_first_not_of(' ') == std::string::npos)
1078 : {
1079 : break;
1080 : }
1081 0 : nodes.push_back(std::stoi(word));
1082 : }
1083 : }
1084 : else
1085 : {
1086 0 : std::istringstream snext(next);
1087 : int i = 1;
1088 0 : while (i < MAX_CHUNK)
1089 : {
1090 0 : std::getline(snext, word, ',');
1091 0 : if (word.find_first_not_of(' ') == std::string::npos)
1092 : {
1093 : break;
1094 : }
1095 0 : nodes.push_back(std::stoi(word));
1096 0 : i++;
1097 : }
1098 0 : }
1099 : }
1100 :
1101 : // Save the element and its geometry tag.
1102 0 : elem_type = HOElemTypeNastran(elem_type, (int)nodes.size());
1103 0 : const int &num_elem_nodes = ElemNumNodes[elem_type - 1];
1104 0 : MFEM_VERIFY((std::size_t)num_elem_nodes == nodes.size(),
1105 : "Mismatch between Nastran and Gmsh element types!");
1106 : std::vector<int> &data = elem_nodes[elem_type];
1107 0 : const int i = (int)data.size();
1108 0 : data.resize(i + 1 + num_elem_nodes);
1109 0 : data[i] = geom_tag;
1110 0 : for (int j = 0; j < num_elem_nodes; j++)
1111 : {
1112 : // Permute back to Gmsh ordering.
1113 0 : const int &p = ElemNodesNastran[elem_type - 1][j];
1114 0 : data[i + 1 + p] = nodes[j];
1115 : }
1116 : }
1117 : }
1118 : }
1119 :
1120 : // Finalize input, write the Gmsh mesh.
1121 0 : input.close();
1122 0 : WriteGmsh(buffer, node_coords, node_tags, elem_nodes, remove_curvature);
1123 0 : }
1124 :
1125 : } // namespace mesh
1126 :
1127 : } // namespace palace
|
