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