libMesh
stl_io.C
Go to the documentation of this file.
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 void test_and_add(std::unique_ptr<libMesh::Elem> e,
49 {
50  const libMesh::Real volume = e->volume();
52  libmesh_warning
53  ("Warning: STL file contained sliver element with volume " <<
54  volume << "!\n");
55  mesh.add_elem(std::move(e));
56 }
57 
58 }
59 
60 
61 namespace libMesh
62 {
63 
66  _subdivide_second_order (true)
67 {
68 }
69 
70 
71 
75  _subdivide_second_order (true)
76 {
77 }
78 
79 
80 
81 void STLIO::write (const std::string & fname)
82 {
83  LOG_SCOPE("write()", "STLIO");
84 
85  // Get a reference to the mesh
86  const MeshBase * mesh_ptr = &MeshOutput<MeshBase>::mesh();
87 
88  libmesh_parallel_only(mesh_ptr->comm());
89 
90  // If necessary, serialize to proc 0, which will do all the writing
91  std::unique_ptr<DistributedMesh> mesh_copy;
92  std::unique_ptr<MeshSerializer> serializer;
93  if (!mesh_ptr->is_serial())
94  {
95  mesh_copy = std::make_unique<DistributedMesh>(*mesh_ptr);
96  mesh_ptr = mesh_copy.get();
97  serializer = std::make_unique<MeshSerializer>(*mesh_copy, true, true);
98  }
99 
100  const MeshBase & mesh = *mesh_ptr;
101 
102  if (mesh.processor_id() != 0)
103  return;
104 
105  // Open the output file stream
106  std::ofstream out_stream (fname.c_str());
107 
108  out_stream << std::setprecision(this->ascii_precision());
109 
110  // Make sure it opened correctly
111  if (!out_stream.good())
112  libmesh_file_error(fname.c_str());
113 
114  // Write the header
115  out_stream << "solid " << this->_name << '\n';
116 
117  // Write all our triangles
118  // scream if we see a non-triangle
119  for (auto elem : mesh.element_ptr_range())
120  {
121  // Get to this later
122  if (elem->type() == TRI6 || elem->type() == TRI7)
123  {
124  if (this->_subdivide_second_order)
125  libmesh_not_implemented();
126 
127  libmesh_error_msg("Tried to write a non-linear triangle to an STL file");
128  }
129 
130  libmesh_error_msg_if(elem->type() != TRI3,
131  "Tried to write a non-triangle to an STL file");
132 
133  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  if (auto length = n.norm())
137  n /= length;
138 
139  out_stream << "facet normal " <<
140  n(0) << ' ' <<
141  n(1) << ' ' <<
142  n(2) << "\n"
143  " outer loop\n"
144  " vertex " <<
145  elem->point(0)(0) << ' ' <<
146  elem->point(0)(1) << ' ' <<
147  elem->point(0)(2) << "\n"
148  " vertex " <<
149  elem->point(1)(0) << ' ' <<
150  elem->point(1)(1) << ' ' <<
151  elem->point(1)(2) << "\n"
152  " vertex " <<
153  elem->point(2)(0) << ' ' <<
154  elem->point(2)(1) << ' ' <<
155  elem->point(2)(2) << "\n"
156  " endloop\n"
157  "endfacet\n";
158  }
159 
160  out_stream << "endsolid" << std::endl;
161 }
162 
163 
164 
165 void STLIO::read (const std::string & filename)
166 {
167  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  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.
176  mesh.clear();
178 
179  std::unique_ptr<std::istream> fstream =
180  this->open_file(filename);
181 
182  char c;
183  const char * expected_header = "solid!";
184  bool is_ascii_stl = false,
185  is_binary_stl = false;
186  while (fstream->get(c))
187  {
188  if (c == ' ')
189  continue;
190  if (c == *expected_header)
191  {
192  ++expected_header;
193  if (*expected_header == '!')
194  {
195  is_ascii_stl = true;
196  break;
197  }
198  }
199  else
200  {
201  is_binary_stl = true; // probably
202  break;
203  }
204  }
205 
206  if (is_ascii_stl)
207  {
208  fstream->seekg(0);
209  this->read_ascii(*fstream);
210  }
211  else if (is_binary_stl)
212  {
213  fstream->seekg(0, std::ios_base::end);
214  std::size_t length = fstream->tellg();
215  fstream->seekg(0);
216  this->read_binary(*fstream, length);
217  }
218  else
219  libmesh_error_msg("Failed to read an STL header in " << filename);
220 }
221 
222 
223 std::unique_ptr<std::istream> STLIO::open_file (const std::string & filename)
224 {
225  const bool gzipped_file = Utility::ends_with(filename, ".gz");
226 
227  std::unique_ptr<std::istream> file;
228 
229  if (gzipped_file)
230  {
231 #ifdef LIBMESH_HAVE_GZSTREAM
232  auto inf = std::make_unique<igzstream>();
233  libmesh_assert(inf);
234  inf->open(filename.c_str(), std::ios::in);
235  file = std::move(inf);
236 #else
237  libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
238 #endif
239  }
240  else
241  {
242  auto inf = std::make_unique<std::ifstream>();
243  libmesh_assert(inf);
244 
245  std::string new_name = Utility::unzip_file(filename);
246 
247  inf->open(new_name.c_str(), std::ios::in);
248  file = std::move(inf);
249  }
250 
251  return file;
252 }
253 
254 
255 
256 void STLIO::read_ascii (std::istream & file)
257 {
258  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 
266 
267  const std::regex all_expected_chars
268  ("^[\\w\\d\\-\\.\\^\\$]*$");
269 
270  const std::regex start_regex
271  ("^\\s*solid");
272  const std::regex start_with_name_regex
273  ("^\\s*solid\\s+(\\w+)");
274 
275  const std::regex start_facet_regex
276  ("^\\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  ("^\\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  "\\s+[+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?))");
295 
296  const std::regex end_facet_regex
297  ("^\\s*end\\s*facet");
298 
299  const std::regex end_loop_regex
300  ("^\\s*end\\s*loop");
301 
302  bool have_started = false;
303  bool in_facet = false;
304  bool in_vertex_loop = false;
305  std::unique_ptr<Tri3> triangle;
306  int next_vertex = 0;
307  int line_num = 0;
308 
309  std::unordered_map<Point, Node *> mesh_points;
310 
311  for (std::string line; std::getline(file, line);)
312  {
313  ++line_num;
314  std::smatch sm;
315 
316  for (char & c : line)
317  {
318  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  if (std::regex_search(line, sm, start_regex))
326  {
327  if (std::regex_search(line, sm, start_with_name_regex))
328  this->_name = sm[1];
329  libmesh_error_msg_if
330  (have_started,
331  "Found two 'solid' lines starting the same STL file?");
332  have_started = true;
333  }
334 
335  else if (std::regex_search(line, sm, start_facet_regex))
336  {
337  libmesh_error_msg_if
338  (in_facet,
339  "Found a repeated 'facet' line with no 'endfacet' between them.");
340  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  else if (std::regex_search(line, end_facet_regex))
358  {
359  libmesh_error_msg_if
360  (!in_facet,
361  "Found an 'endfacet' line with no matching 'facet'.");
362  in_facet = false;
363  }
364 
365  else if (std::regex_search(line, start_loop_regex))
366  {
367  libmesh_error_msg_if
368  (!in_facet,
369  "Found an 'outer loop' line with no matching 'facet'.");
370  libmesh_error_msg_if
371  (in_vertex_loop,
372  "Found a repeated 'outer loop' line with no 'endloop' between them.");
373  in_vertex_loop = true;
374  triangle = std::make_unique<Tri3>();
375  next_vertex = 0;
376  }
377 
378  else if (std::regex_search(line, sm, vertex_regex))
379  {
380  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  std::stringstream ss(normalvec);
385  Point p;
386  ss >> p(0);
387  ss >> p(1);
388  ss >> p(2);
389 
390  Node * node;
391  if (auto it = mesh_points.find(p); it != mesh_points.end())
392  {
393  node = it->second;
394  }
395  else
396  {
397  node = mesh.add_point(p);
398  mesh_points[p] = node;
399  }
400 
401  libmesh_error_msg_if
402  (next_vertex > 2,
403  "Found more than 3 vertices in a loop; STLIO only supports Tri3.");
404  triangle->set_node(next_vertex++, node);
405  }
406 
407  else if (std::regex_search(line, end_loop_regex))
408  {
409  libmesh_error_msg_if
410  (next_vertex != 3,
411  "Found an 'endloop' line after only seeing " << next_vertex << " vertices in 'loop'.");
412  libmesh_error_msg_if
413  (!in_vertex_loop,
414  "Found an 'endloop' line with no matching 'loop'.");
415  in_vertex_loop = false;
416 
417  test_and_add(std::move(triangle), mesh);
418  }
419  }
420 
421  libmesh_error_msg_if
422  (in_facet,
423  "File ended without ending a facet first.");
424 
425  libmesh_error_msg_if
426  (in_vertex_loop,
427  "File ended without ending an outer loop first.");
428 #endif // LIBMESH_HAVE_CXX11_REGEX
429 }
430 
431 
432 
433 void STLIO::read_binary (std::istream & file,
434  std::size_t input_size)
435 {
436  LOG_SCOPE("read_binary()", "STLIO");
437 
439 
440  char header_buffer[80];
441 
442  // 80-character header which is generally ignored - Wikipedia
443  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  uint32_t test_int = 0x87654321;
449  const bool big_endian = ((*reinterpret_cast<char *>(&test_int)) != 0x21);
450  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  uint32_t n_elem = 0;
459 
460  file.read(reinterpret_cast<char*>(&n_elem), 4);
461  endian_fix(n_elem);
462 
463  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  std::unique_ptr<Tri3> triangle;
470  std::unordered_map<Point, Node *> mesh_points;
471 
472  for (unsigned int e : make_range(n_elem))
473  {
474  libmesh_ignore(e);
475 
476  triangle = std::make_unique<Tri3>();
477 
478  // We'll ignore facets' normals for now
479  char ignored_buffer[12];
480  file.read(ignored_buffer, 12);
481 
482  // Read vertex locations
483  for (int i : make_range(3))
484  {
485  float point_buffer[3];
486  file.read(reinterpret_cast<char*>(point_buffer), 12);
487  endian_fix(point_buffer[0]);
488  endian_fix(point_buffer[1]);
489  endian_fix(point_buffer[2]);
490  const Point p {point_buffer[0],
491  point_buffer[1],
492  point_buffer[2]};
493 
494  Node * node;
495  if (auto it = mesh_points.find(p); it != mesh_points.end())
496  {
497  node = it->second;
498  }
499  else
500  {
501  node = mesh.add_point(p);
502  mesh_points[p] = node;
503  }
504 
505  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  file.read(ignored_buffer, 2);
511 
512  test_and_add(std::move(triangle), mesh);
513  }
514 }
515 
516 
517 } // namespace libMesh
const MT & mesh() const
Definition: mesh_output.h:259
This Functor simply takes an object and reverses its byte representation.
Definition: utility.h:449
std::unique_ptr< std::istream > open_file(const std::string &filename)
Helper to open possibly-zipped files.
Definition: stl_io.C:223
A Node is like a Point, but with more information.
Definition: node.h:52
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
bool ends_with(std::string_view superstring, std::string_view suffix)
Look for a substring at the very end of a string.
Definition: utility.C:213
static constexpr Real TOLERANCE
virtual void read(const std::string &mesh_file) override
This method implements reading a mesh from a specified file.
Definition: stl_io.C:165
MeshBase & mesh
unsigned int & ascii_precision()
Return/set the precision to use when writing ASCII files.
Definition: mesh_output.h:269
const Parallel::Communicator & comm() const
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
The libMesh namespace provides an interface to certain functionality in the library.
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
This is the MeshBase class.
Definition: mesh_base.h:75
virtual void read_ascii(std::istream &input)
This method implements reading a mesh from a specified ASCII input stream.
Definition: stl_io.C:256
virtual bool is_serial() const
Definition: mesh_base.h:211
void libmesh_ignore(const Args &...)
virtual void write(const std::string &) override
This method implements writing a mesh to a specified file.
Definition: stl_io.C:81
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:57
bool _subdivide_second_order
Flag to subdivide second order elements.
Definition: stl_io.h:120
std::string unzip_file(std::string_view name)
Create an unzipped copy of a bz2 or xz file, returning the name of the now-unzipped file that can be ...
Definition: utility.C:164
libmesh_assert(ctx)
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
Find the total volume of a mesh (interpreting that as area for dim = 2, or total arc length for dim =...
Definition: mesh_tools.C:985
virtual void read_binary(std::istream &input, std::size_t input_size=0)
This method implements reading a mesh from a specified binary input stream.
Definition: stl_io.C:433
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::string _name
Definition: stl_io.h:122
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
STLIO(const MeshBase &)
Constructor.
Definition: stl_io.C:64
processor_id_type processor_id() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39