libMesh
tetgen_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 
19 // Local includes
20 #include "libmesh/tetgen_io.h"
21 
22 #include "libmesh/cell_tet4.h"
23 #include "libmesh/cell_tet10.h"
24 #include "libmesh/mesh_base.h"
25 #include "libmesh/utility.h"
26 
27 // C++ includes
28 #include <array>
29 #include <fstream>
30 
31 namespace libMesh
32 {
33 
34 // ------------------------------------------------------------
35 // TetgenIO class members
36 void TetGenIO::read (const std::string & name)
37 {
38  // This is a serial-only process for now;
39  // the Mesh should be read on processor 0 and
40  // broadcast later
41  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
42 
43  std::string name_node, name_ele, dummy;
44 
45  // tetgen only works in 3D
46  MeshInput<MeshBase>::mesh().set_mesh_dimension(3);
47 
48 #if LIBMESH_DIM < 3
49  libmesh_error_msg("Cannot open dimension 3 mesh file when configured without 3D support.");
50 #endif
51 
52  // Check name for *.node or *.ele extension.
53  // Set std::istream for node_stream and ele_stream.
54  //
55  if (Utility::contains(name, ".node"))
56  {
57  name_node = name;
58  dummy = name;
59  std::size_t position = dummy.rfind(".node");
60  name_ele = dummy.replace(position, 5, ".ele");
61  }
62  else if (Utility::contains(name, ".ele"))
63  {
64  name_ele = name;
65  dummy = name;
66  std::size_t position = dummy.rfind(".ele");
67  name_node = dummy.replace(position, 4, ".node");
68  }
69  else
70  libmesh_error_msg("ERROR: Unrecognized file name: " << name);
71 
72 
73 
74  // Set the streams from which to read in
75  std::ifstream node_stream (name_node.c_str());
76  std::ifstream ele_stream (name_ele.c_str());
77 
78  libmesh_error_msg_if(!node_stream.good() || !ele_stream.good(),
79  "Error while opening either "
80  << name_node
81  << " or "
82  << name_ele);
83 
84  libMesh::out<< "TetGenIO found the tetgen files to read " <<std::endl;
85 
86  // Skip the comment lines at the beginning
87  this->skip_comment_lines (node_stream, '#');
88  this->skip_comment_lines (ele_stream, '#');
89 
90  // Read the nodes and elements from the streams
91  this->read_nodes_and_elem (node_stream, ele_stream);
92  libMesh::out<< "TetGenIO read in nodes and elements " <<std::endl;
93 }
94 
95 
96 
97 void TetGenIO::read_nodes_and_elem (std::istream & node_stream,
98  std::istream & ele_stream)
99 {
100  _num_nodes = 0;
101  _num_elements = 0;
102 
103  // Read all the datasets.
104  this->node_in (node_stream);
105  this->element_in (ele_stream);
106 
107  // some more clean-up
108  _assign_nodes.clear();
109 }
110 
111 
112 
113 //----------------------------------------------------------------------
114 // Function to read in the node table.
115 void TetGenIO::node_in (std::istream & node_stream)
116 {
117  // Check input buffer
118  libmesh_assert (node_stream.good());
119 
120  // Get a reference to the mesh
122 
123  unsigned int dimension=0, nAttributes=0, BoundaryMarkers=0;
124 
125  node_stream >> _num_nodes // Read the number of nodes from the stream
126  >> dimension // Read the dimension from the stream
127  >> nAttributes // Read the number of attributes from stream
128  >> BoundaryMarkers; // Read if or not boundary markers are included in *.node (0 or 1)
129 
130  // Read the nodal coordinates from the node_stream (*.node file).
131  unsigned int node_lab=0;
132  Real dummy;
133 
134  // If present, make room for node attributes to be stored.
135  this->node_attributes.resize(nAttributes);
136  for (unsigned i=0; i<nAttributes; ++i)
137  this->node_attributes[i].resize(_num_nodes);
138 
139 
140  for (unsigned int i=0; i<_num_nodes; i++)
141  {
142  // Check input buffer
143  libmesh_assert (node_stream.good());
144 
145  std::array<Real, 3> xyz;
146 
147  node_stream >> node_lab // node number
148  >> xyz[0] // x-coordinate value
149  >> xyz[1] // y-coordinate value
150  >> xyz[2]; // z-coordinate value
151 
152  // Read and store attributes from the stream.
153  for (unsigned int j=0; j<nAttributes; j++)
154  node_stream >> node_attributes[j][i];
155 
156  // Read (and discard) boundary marker if BoundaryMarker=1.
157  // TODO: should we store this somehow?
158  if (BoundaryMarkers == 1)
159  node_stream >> dummy;
160 
161  // Store the new position of the node under its label.
162  //_assign_nodes.emplace(node_lab,i);
163  _assign_nodes[node_lab] = i;
164 
165  Point p(xyz[0]);
166 #if LIBMESH_DIM > 1
167  p(1) = xyz[1];
168 #else
169  libmesh_assert_equal_to(xyz[1], 0);
170 #endif
171 #if LIBMESH_DIM > 2
172  p(2) = xyz[2];
173 #else
174  libmesh_assert_equal_to(xyz[2], 0);
175 #endif
176 
177  // Add this point to the Mesh.
178  mesh.add_point(p, i);
179  }
180 }
181 
182 
183 
184 //----------------------------------------------------------------------
185 // Function to read in the element table.
186 void TetGenIO::element_in (std::istream & ele_stream)
187 {
188  // Check input buffer
189  libmesh_assert (ele_stream.good());
190 
191  // Get a reference to the mesh
193 
194  // Read the elements from the ele_stream (*.ele file).
195  unsigned int element_lab=0, n_nodes=0, region_attribute=0;
196 
197  ele_stream >> _num_elements // Read the number of tetrahedrons from the stream.
198  >> n_nodes // Read the number of nodes per tetrahedron from the stream (defaults to 4).
199  >> region_attribute; // Read the number of attributes from stream.
200 
201  // According to the Tetgen docs for .ele files:
202  // http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual006.html#ff_ele
203  // region_attribute can either 0 or 1, and specifies whether, for
204  // each tetrahedron, there is an extra integer specifying which
205  // region it belongs to. Normally, this id matches a value in a
206  // corresponding .poly or .smesh file, but here we simply use it to
207  // set the subdomain_id of the element in question.
208  libmesh_error_msg_if(region_attribute > 1,
209  "Invalid region_attribute " << region_attribute << " specified in .ele file.");
210 
211  // Vector that maps Tetgen node numbering to libMesh node numbering. Tet4s are
212  // numbered identically, but Tet10s are not.
213  static const unsigned int assign_elm_nodes[] = {0, 1, 2, 3, 9, 7, 4, 5, 8, 6};
214 
215  for (dof_id_type i=0; i<_num_elements; i++)
216  {
217  libmesh_assert (ele_stream.good());
218 
219  // TetGen only supports Tet4 and Tet10 elements.
220  Elem * elem = nullptr;
221 
222  if (n_nodes==4)
224 
225  else if (n_nodes==10)
227 
228  else
229  libmesh_error_msg("Elements with " << n_nodes <<
230  " nodes are not supported in the LibMesh tetgen module.");
231 
232  libmesh_assert(elem);
233  libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
234 
235  // The first number on the line is the tetrahedron number. We
236  // have previously ignored this, preferring to set our own ids,
237  // but this could be changed to respect the Tetgen numbering if
238  // desired.
239  ele_stream >> element_lab;
240 
241  // Read node labels
242  for (dof_id_type j=0; j<n_nodes; j++)
243  {
244  dof_id_type node_label;
245  ele_stream >> node_label;
246 
247  // Assign node to element
248  elem->set_node(assign_elm_nodes[j],
249  mesh.node_ptr(_assign_nodes[node_label]));
250  }
251 
252  // Read the region attribute (if present) and use it to set the subdomain id.
253  if (region_attribute)
254  {
255  unsigned int region;
256  ele_stream >> region;
257 
258  // Make sure that the id we read can be successfully cast to
259  // an integral value of type subdomain_id_type.
260  elem->subdomain_id() = cast_int<subdomain_id_type>(region);
261  }
262  }
263 }
264 
265 
270 void TetGenIO::write (const std::string & fname)
271 {
272  // libmesh_assert three dimensions (should be extended later)
273  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().mesh_dimension(), 3);
274 
275  libmesh_error_msg_if(!Utility::contains(fname, ".poly"),
276  "ERROR: Unrecognized file name: " << fname);
277 
278  // Open the output file stream
279  std::ofstream out_stream (fname.c_str());
280 
281  // Make sure it opened correctly
282  if (!out_stream.good())
283  libmesh_file_error(fname.c_str());
284 
285  // Get a reference to the mesh
287 
288  // Begin interfacing with the .poly file
289  {
290  // header:
291  out_stream << "# poly file output generated by libmesh\n"
292  << mesh.n_nodes() << " 3 0 0\n";
293 
294  // write the nodes:
295  for (auto v : make_range(mesh.n_nodes()))
296  out_stream << v << " "
297  << mesh.point(v)(0) << " "
298  << mesh.point(v)(1) << " "
299  << mesh.point(v)(2) << "\n";
300  }
301 
302  {
303  // write the connectivity:
304  out_stream << "# Facets:\n"
305  << mesh.n_elem() << " 0\n";
306 
307  for (const auto & elem : mesh.active_element_ptr_range())
308  out_stream << "1\n3 " // no. of facet polygons
309  // << elem->n_nodes() << " "
310  << elem->node_id(0) << " "
311  << elem->node_id(1) << " "
312  << elem->node_id(2) << "\n";
313  }
314 
315  // end of the file
316  out_stream << "0\n"; // no holes output!
317  out_stream << "\n\n# end of file\n";
318 }
319 
320 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
dof_id_type _num_elements
total number of elements.
Definition: tetgen_io.h:130
const MT & mesh() const
Definition: mesh_output.h:259
std::map< dof_id_type, dof_id_type > _assign_nodes
stores new positions of nodes.
Definition: tetgen_io.h:120
virtual void read(const std::string &) override
This method implements reading a mesh from a specified file in TetGen format.
Definition: tetgen_io.C:36
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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.
void skip_comment_lines(std::istream &in, const char comment_start)
Reads input from in, skipping all the lines that start with the character comment_start.
Definition: mesh_input.h:187
This is the MeshBase class.
Definition: mesh_base.h:75
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libmesh_assert(ctx)
bool contains(std::string_view superstring, std::string_view substring)
Look for a substring within a string.
Definition: utility.C:205
dof_id_type _num_nodes
total number of nodes.
Definition: tetgen_io.h:125
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem&#39;s ...
Definition: elem.C:558
void read_nodes_and_elem(std::istream &node_stream, std::istream &ele_stream)
Reads a mesh (nodes & elements) from the file provided through node_stream and ele_stream.
Definition: tetgen_io.C:97
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
virtual void write(const std::string &) override
This method implements writing a mesh to a specified ".poly" file.
Definition: tetgen_io.C:270
virtual const Point & point(const dof_id_type i) const =0
void node_in(std::istream &node_stream)
Method reads nodes from node_stream and stores them in vector<Node *> nodes in the order they come in...
Definition: tetgen_io.C:115
void element_in(std::istream &ele_stream)
Method reads elements and stores them in vector<Elem *> elements in the same order as they come in...
Definition: tetgen_io.C:186
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
std::vector< std::vector< Real > > node_attributes
Data structure to hold node attributes read in from file.
Definition: tetgen_io.h:82
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual dof_id_type n_nodes() const =0
uint8_t dof_id_type
Definition: id_types.h:67