Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // Local includes
19 : #include "libmesh/stl_io.h"
20 :
21 : #include "libmesh/distributed_mesh.h"
22 : #include "libmesh/enum_elem_type.h"
23 : #include "libmesh/face_tri3.h"
24 : #include "libmesh/libmesh_config.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/mesh_base.h"
27 : #include "libmesh/point.h"
28 : #include "libmesh/utility.h"
29 :
30 : // gzstream for reading compressed files as a stream
31 : #ifdef LIBMESH_HAVE_GZSTREAM
32 : # include "gzstream.h"
33 : #endif
34 :
35 : // C++ includes
36 : #include <iomanip>
37 : #include <fstream>
38 : #include <vector>
39 : #include <cctype> // isspace
40 :
41 : #ifdef LIBMESH_HAVE_CXX11_REGEX
42 : #include <regex>
43 : #endif
44 :
45 : namespace {
46 :
47 5592 : void test_and_add(std::unique_ptr<libMesh::Elem> e,
48 : libMesh::MeshBase & mesh)
49 : {
50 5592 : const libMesh::Real volume = e->volume();
51 : if (volume < libMesh::TOLERANCE * libMesh::TOLERANCE)
52 : libmesh_warning
53 : ("Warning: STL file contained sliver element with volume " <<
54 : volume << "!\n");
55 6058 : mesh.add_elem(std::move(e));
56 5592 : }
57 :
58 : }
59 :
60 :
61 : namespace libMesh
62 : {
63 :
64 0 : STLIO::STLIO (const MeshBase & mesh) :
65 : MeshOutput<MeshBase> (mesh),
66 0 : _subdivide_second_order (true)
67 : {
68 0 : }
69 :
70 :
71 :
72 142 : STLIO::STLIO (MeshBase & mesh) :
73 : MeshInput<MeshBase> (mesh),
74 : MeshOutput<MeshBase>(mesh),
75 142 : _subdivide_second_order (true)
76 : {
77 142 : }
78 :
79 :
80 :
81 0 : void STLIO::write (const std::string & fname)
82 : {
83 0 : LOG_SCOPE("write()", "STLIO");
84 :
85 : // Get a reference to the mesh
86 0 : const MeshBase * mesh_ptr = &MeshOutput<MeshBase>::mesh();
87 :
88 0 : libmesh_parallel_only(mesh_ptr->comm());
89 :
90 : // If necessary, serialize to proc 0, which will do all the writing
91 0 : std::unique_ptr<DistributedMesh> mesh_copy;
92 0 : std::unique_ptr<MeshSerializer> serializer;
93 0 : if (!mesh_ptr->is_serial())
94 : {
95 0 : mesh_copy = std::make_unique<DistributedMesh>(*mesh_ptr);
96 0 : mesh_ptr = mesh_copy.get();
97 0 : serializer = std::make_unique<MeshSerializer>(*mesh_copy, true, true);
98 : }
99 :
100 0 : const MeshBase & mesh = *mesh_ptr;
101 :
102 0 : if (mesh.processor_id() != 0)
103 0 : return;
104 :
105 : // Open the output file stream
106 0 : std::ofstream out_stream (fname.c_str());
107 :
108 0 : out_stream << std::setprecision(this->ascii_precision());
109 :
110 : // Make sure it opened correctly
111 0 : if (!out_stream.good())
112 0 : libmesh_file_error(fname.c_str());
113 :
114 : // Write the header
115 0 : out_stream << "solid " << this->_name << '\n';
116 :
117 : // Write all our triangles
118 : // scream if we see a non-triangle
119 0 : for (auto elem : mesh.element_ptr_range())
120 : {
121 : // Get to this later
122 0 : if (elem->type() == TRI6 || elem->type() == TRI7)
123 : {
124 0 : if (this->_subdivide_second_order)
125 0 : libmesh_not_implemented();
126 :
127 0 : libmesh_error_msg("Tried to write a non-linear triangle to an STL file");
128 : }
129 :
130 0 : libmesh_error_msg_if(elem->type() != TRI3,
131 : "Tried to write a non-triangle to an STL file");
132 :
133 0 : auto n = (elem->point(1)-elem->point(0)).cross(elem->point(2)-elem->point(0));
134 :
135 : // Other STL files have slivers, I guess ours can too
136 0 : if (auto length = n.norm())
137 0 : n /= length;
138 :
139 0 : out_stream << "facet normal " <<
140 0 : n(0) << ' ' <<
141 0 : n(1) << ' ' <<
142 0 : n(2) << "\n"
143 : " outer loop\n"
144 0 : " vertex " <<
145 0 : elem->point(0)(0) << ' ' <<
146 0 : elem->point(0)(1) << ' ' <<
147 0 : elem->point(0)(2) << "\n"
148 0 : " vertex " <<
149 0 : elem->point(1)(0) << ' ' <<
150 0 : elem->point(1)(1) << ' ' <<
151 0 : elem->point(1)(2) << "\n"
152 0 : " vertex " <<
153 0 : elem->point(2)(0) << ' ' <<
154 0 : elem->point(2)(1) << ' ' <<
155 0 : elem->point(2)(2) << "\n"
156 : " endloop\n"
157 0 : "endfacet\n";
158 0 : }
159 :
160 0 : out_stream << "endsolid" << std::endl;
161 0 : }
162 :
163 :
164 :
165 24 : void STLIO::read (const std::string & filename)
166 : {
167 4 : LOG_SCOPE("read()", "STLIO");
168 :
169 : // This is a serial-only process for now;
170 : // the Mesh should be read on processor 0 and
171 : // broadcast later
172 2 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
173 :
174 : // Clear the mesh so we are sure to start from a pristine state.
175 24 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
176 24 : mesh.clear();
177 24 : mesh.set_mesh_dimension(2);
178 :
179 : std::unique_ptr<std::istream> fstream =
180 26 : this->open_file(filename);
181 :
182 : char c;
183 2 : const char * expected_header = "solid!";
184 2 : bool is_ascii_stl = false,
185 2 : is_binary_stl = false;
186 72 : while (fstream->get(c))
187 : {
188 72 : if (c == ' ')
189 0 : continue;
190 72 : if (c == *expected_header)
191 : {
192 60 : ++expected_header;
193 60 : if (*expected_header == '!')
194 : {
195 1 : is_ascii_stl = true;
196 1 : break;
197 : }
198 : }
199 : else
200 : {
201 1 : is_binary_stl = true; // probably
202 1 : break;
203 : }
204 : }
205 :
206 24 : if (is_ascii_stl)
207 : {
208 12 : fstream->seekg(0);
209 13 : this->read_ascii(*fstream);
210 : }
211 12 : else if (is_binary_stl)
212 : {
213 12 : fstream->seekg(0, std::ios_base::end);
214 12 : std::size_t length = fstream->tellg();
215 12 : fstream->seekg(0);
216 13 : this->read_binary(*fstream, length);
217 : }
218 : else
219 0 : libmesh_error_msg("Failed to read an STL header in " << filename);
220 24 : }
221 :
222 :
223 24 : std::unique_ptr<std::istream> STLIO::open_file (const std::string & filename)
224 : {
225 24 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
226 :
227 24 : std::unique_ptr<std::istream> file;
228 :
229 24 : if (gzipped_file)
230 : {
231 : #ifdef LIBMESH_HAVE_GZSTREAM
232 0 : auto inf = std::make_unique<igzstream>();
233 0 : libmesh_assert(inf);
234 0 : inf->open(filename.c_str(), std::ios::in);
235 0 : file = std::move(inf);
236 : #else
237 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
238 : #endif
239 0 : }
240 : else
241 : {
242 26 : auto inf = std::make_unique<std::ifstream>();
243 2 : libmesh_assert(inf);
244 :
245 26 : std::string new_name = Utility::unzip_file(filename);
246 :
247 24 : inf->open(new_name.c_str(), std::ios::in);
248 22 : file = std::move(inf);
249 20 : }
250 :
251 24 : return file;
252 0 : }
253 :
254 :
255 :
256 12 : void STLIO::read_ascii (std::istream & file)
257 : {
258 2 : LOG_SCOPE("read_ascii()", "STLIO");
259 :
260 : #ifndef LIBMESH_HAVE_CXX11_REGEX
261 : libmesh_not_implemented(); // What is your compiler?!? Email us!
262 : libmesh_ignore(file);
263 : #else
264 :
265 12 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
266 :
267 : const std::regex all_expected_chars
268 14 : ("^[\\w\\d\\-\\.\\^\\$]*$");
269 :
270 : const std::regex start_regex
271 14 : ("^\\s*solid");
272 : const std::regex start_with_name_regex
273 14 : ("^\\s*solid\\s+(\\w+)");
274 :
275 : const std::regex start_facet_regex
276 14 : ("^\\s*facet");
277 :
278 : // We'll ignore facets' normals for now
279 : /*
280 : const std::regex facet_with_normal_regex
281 : ("^\\s*facet\\s+normal"
282 : "\\s+([+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?)"
283 : "\\s+[+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?)"
284 : "\\s+[+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?))");
285 : */
286 :
287 : const std::regex start_loop_regex
288 14 : ("^\\s*outer\\s+loop");
289 :
290 : const std::regex vertex_regex
291 : ("^\\s*vertex"
292 : "\\s+([+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?)"
293 : "\\s+[+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?)"
294 14 : "\\s+[+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?))");
295 :
296 : const std::regex end_facet_regex
297 14 : ("^\\s*end\\s*facet");
298 :
299 : const std::regex end_loop_regex
300 14 : ("^\\s*end\\s*loop");
301 :
302 1 : bool have_started = false;
303 1 : bool in_facet = false;
304 1 : bool in_vertex_loop = false;
305 11 : std::unique_ptr<Tri3> triangle;
306 1 : int next_vertex = 0;
307 1 : int line_num = 0;
308 :
309 2 : std::unordered_map<Point, Node *> mesh_points;
310 :
311 3396 : for (std::string line; std::getline(file, line);)
312 : {
313 3384 : ++line_num;
314 564 : std::smatch sm;
315 :
316 64596 : for (char & c : line)
317 : {
318 61212 : libmesh_error_msg_if
319 : (c != '\n' && c != '\r' &&
320 : (c < ' ' || c > '~'),
321 : "Found non-ASCII character " << int(c) << " on line " <<
322 : line_num << "\nSTLIO does not yet support binary files.");
323 : }
324 :
325 3384 : if (std::regex_search(line, sm, start_regex))
326 : {
327 12 : if (std::regex_search(line, sm, start_with_name_regex))
328 23 : this->_name = sm[1];
329 12 : libmesh_error_msg_if
330 : (have_started,
331 : "Found two 'solid' lines starting the same STL file?");
332 1 : have_started = true;
333 : }
334 :
335 3372 : else if (std::regex_search(line, sm, start_facet_regex))
336 : {
337 480 : libmesh_error_msg_if
338 : (in_facet,
339 : "Found a repeated 'facet' line with no 'endfacet' between them.");
340 40 : in_facet = true;
341 :
342 : /*
343 : if (std::regex_search(line, sm, facet_with_normal_regex))
344 : {
345 : const std::string normalvec = sm[1];
346 :
347 : // Try to be compatible with higher-than-double T; maybe
348 : // *someone* out there is doing STL in 128 bits? ...
349 : std::stringstream ss(normalvec);
350 : ss >> normal(0);
351 : ss >> normal(1);
352 : ss >> normal(2);
353 : }
354 : */
355 : }
356 :
357 2892 : else if (std::regex_search(line, end_facet_regex))
358 : {
359 480 : libmesh_error_msg_if
360 : (!in_facet,
361 : "Found an 'endfacet' line with no matching 'facet'.");
362 40 : in_facet = false;
363 : }
364 :
365 2412 : else if (std::regex_search(line, start_loop_regex))
366 : {
367 480 : libmesh_error_msg_if
368 : (!in_facet,
369 : "Found an 'outer loop' line with no matching 'facet'.");
370 480 : libmesh_error_msg_if
371 : (in_vertex_loop,
372 : "Found a repeated 'outer loop' line with no 'endloop' between them.");
373 40 : in_vertex_loop = true;
374 920 : triangle = std::make_unique<Tri3>();
375 40 : next_vertex = 0;
376 : }
377 :
378 1932 : else if (std::regex_search(line, sm, vertex_regex))
379 : {
380 1560 : const std::string normalvec = sm[1];
381 :
382 : // Try to be compatible with higher-than-double-precision T;
383 : // maybe *someone* out there is doing STL in 128 bits? ...
384 1680 : std::stringstream ss(normalvec);
385 120 : Point p;
386 120 : ss >> p(0);
387 120 : ss >> p(1);
388 120 : ss >> p(2);
389 :
390 : Node * node;
391 1440 : if (auto it = mesh_points.find(p); it != mesh_points.end())
392 : {
393 720 : node = it->second;
394 : }
395 : else
396 : {
397 720 : node = mesh.add_point(p);
398 720 : mesh_points[p] = node;
399 : }
400 :
401 1440 : libmesh_error_msg_if
402 : (next_vertex > 2,
403 : "Found more than 3 vertices in a loop; STLIO only supports Tri3.");
404 1440 : triangle->set_node(next_vertex++, node);
405 1200 : }
406 :
407 492 : else if (std::regex_search(line, end_loop_regex))
408 : {
409 480 : libmesh_error_msg_if
410 : (next_vertex != 3,
411 : "Found an 'endloop' line after only seeing " << next_vertex << " vertices in 'loop'.");
412 480 : libmesh_error_msg_if
413 : (!in_vertex_loop,
414 : "Found an 'endloop' line with no matching 'loop'.");
415 40 : in_vertex_loop = false;
416 :
417 520 : test_and_add(std::move(triangle), mesh);
418 : }
419 : }
420 :
421 12 : libmesh_error_msg_if
422 : (in_facet,
423 : "File ended without ending a facet first.");
424 :
425 12 : libmesh_error_msg_if
426 : (in_vertex_loop,
427 : "File ended without ending an outer loop first.");
428 : #endif // LIBMESH_HAVE_CXX11_REGEX
429 12 : }
430 :
431 :
432 :
433 12 : void STLIO::read_binary (std::istream & file,
434 : std::size_t input_size)
435 : {
436 2 : LOG_SCOPE("read_binary()", "STLIO");
437 :
438 12 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
439 :
440 : char header_buffer[80];
441 :
442 : // 80-character header which is generally ignored - Wikipedia
443 12 : file.read(header_buffer, 80);
444 :
445 : // What's our endianness here? Input binary files are specified to
446 : // be little endian, and we might need to do conversions.
447 :
448 1 : uint32_t test_int = 0x87654321;
449 1 : const bool big_endian = ((*reinterpret_cast<char *>(&test_int)) != 0x21);
450 1 : const Utility::ReverseBytes endian_fix{big_endian};
451 :
452 : // C++ doesn't specify a size for float, but we really need 4-byte
453 : // floats here. Fortunately basically every implementation ever
454 : // uses exactly 4 bytes for float.
455 : if constexpr (sizeof(float) != 4)
456 : libmesh_error_msg("Trying to read 4 byte floats without 4-byte float?");
457 :
458 12 : uint32_t n_elem = 0;
459 :
460 12 : file.read(reinterpret_cast<char*>(&n_elem), 4);
461 11 : endian_fix(n_elem);
462 :
463 12 : libmesh_error_msg_if
464 : (input_size &&
465 : (input_size < 84+n_elem*50),
466 : "Not enough data for " << n_elem << " STL triangles in " <<
467 : input_size << " uncompressed bytes.");
468 :
469 11 : std::unique_ptr<Tri3> triangle;
470 2 : std::unordered_map<Point, Node *> mesh_points;
471 :
472 5124 : for (unsigned int e : make_range(n_elem))
473 : {
474 426 : libmesh_ignore(e);
475 :
476 9372 : triangle = std::make_unique<Tri3>();
477 :
478 : // We'll ignore facets' normals for now
479 : char ignored_buffer[12];
480 5112 : file.read(ignored_buffer, 12);
481 :
482 : // Read vertex locations
483 20448 : for (int i : make_range(3))
484 : {
485 : float point_buffer[3];
486 15336 : file.read(reinterpret_cast<char*>(point_buffer), 12);
487 14058 : endian_fix(point_buffer[0]);
488 14058 : endian_fix(point_buffer[1]);
489 14058 : endian_fix(point_buffer[2]);
490 15336 : const Point p {point_buffer[0],
491 15336 : point_buffer[1],
492 15336 : point_buffer[2]};
493 :
494 : Node * node;
495 15336 : if (auto it = mesh_points.find(p); it != mesh_points.end())
496 : {
497 12756 : node = it->second;
498 : }
499 : else
500 : {
501 2580 : node = mesh.add_point(p);
502 2580 : mesh_points[p] = node;
503 : }
504 :
505 15336 : triangle->set_node(i, node);
506 : }
507 :
508 : // The 2-byte "attribute byte count" is unstandardized; typically
509 : // 0, or sometimes triangle color. Ignore it.
510 5112 : file.read(ignored_buffer, 2);
511 :
512 5538 : test_and_add(std::move(triangle), mesh);
513 : }
514 12 : }
515 :
516 :
517 : } // namespace libMesh
|
