tetgen_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 #include "libmesh/mesh_data.h"
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // TetgenIO class members
34 void TetGenIO::read (const std::string& name)
35 {
36  // This is a serial-only process for now;
37  // the Mesh should be read on processor 0 and
38  // broadcast later
39  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
40 
41  std::string name_node, name_ele, dummy;
42 
43  // tetgen only works in 3D
44  MeshInput<MeshBase>::mesh().set_mesh_dimension(3);
45 
46 #if LIBMESH_DIM < 3
47  libMesh::err << "Cannot open dimension 3 mesh file when configured without 3D support." <<
48  std::endl;
49  libmesh_error();
50 #endif
51 
52  // Check name for *.node or *.ele extension.
53  // Set std::istream for node_stream and ele_stream.
54  //
55  if (name.rfind(".node") < name.size())
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 (name.rfind(".ele") < name.size())
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  {
71  libMesh::err << "ERROR: Unrecognized file name: "
72  << name << std::endl;
73  libmesh_error();
74  }
75 
76 
77 
78  // Set the streams from which to read in
79  std::ifstream node_stream (name_node.c_str());
80  std::ifstream ele_stream (name_ele.c_str());
81 
82  if ( !node_stream.good() || !ele_stream.good() )
83  {
84  libMesh::err << "ERROR: One or both Input file(s) not good." << std::endl
85  << "Error checking files "
86  << name_node << " and "
87  << name_ele << std::endl;
88  libmesh_error();
89  }
90  libMesh::out<< "TetGenIO found the tetgen files to read " <<std::endl;
91 
92  // Skip the comment lines at the beginning
93  this->skip_comment_lines (node_stream, '#');
94  this->skip_comment_lines (ele_stream, '#');
95 
96  // Read the nodes and elements from the streams
97  this->read_nodes_and_elem (node_stream, ele_stream);
98  libMesh::out<< "TetGenIO read in nodes and elements " <<std::endl;
99 }
100 
101 
102 
103 void TetGenIO::read_nodes_and_elem (std::istream& node_stream,
104  std::istream& ele_stream)
105 {
106  _num_nodes = 0;
107  _num_elements = 0;
108 
109  // Read all the datasets.
110  this->node_in (node_stream);
111  this->element_in (ele_stream);
112 
113  // Tell the MeshData object that we are finished
114  // reading data.
115  if (this->_mesh_data != NULL)
117 
118  // some more clean-up
119  _assign_nodes.clear();
120 }
121 
122 
123 
124 //----------------------------------------------------------------------
125 // Function to read in the node table.
126 void TetGenIO::node_in (std::istream& node_stream)
127 {
128  // Check input buffer
129  libmesh_assert (node_stream.good());
130 
131  // Get a reference to the mesh
133 
134  unsigned int dimension=0, nAttributes=0, BoundaryMarkers=0;
135 
136  node_stream >> _num_nodes // Read the number of nodes from the stream
137  >> dimension // Read the dimension from the stream
138  >> nAttributes // Read the number of attributes from stream
139  >> BoundaryMarkers; // Read if or not boundary markers are included in *.node (0 or 1)
140 
141  // Read the nodal coordinates from the node_stream (*.node file).
142  unsigned int node_lab=0;
143  Point xyz;
144  Real dummy;
145 
146  // If present, make room for node attributes to be stored.
147  this->node_attributes.resize(nAttributes);
148  for (unsigned i=0; i<nAttributes; ++i)
150 
151 
152  for (unsigned int i=0; i<_num_nodes; i++)
153  {
154  // Check input buffer
155  libmesh_assert (node_stream.good());
156 
157  node_stream >> node_lab // node number
158  >> xyz(0) // x-coordinate value
159  >> xyz(1) // y-coordinate value
160  >> xyz(2); // z-coordinate value
161 
162  // Read and store attributes from the stream.
163  for (unsigned int j=0; j<nAttributes; j++)
164  node_stream >> node_attributes[j][i];
165 
166  // Read (and discard) boundary marker if BoundaryMarker=1.
167  // TODO: should we store this somehow?
168  if (BoundaryMarkers == 1)
169  node_stream >> dummy;
170 
171  // Store the new position of the node under its label.
172  //_assign_nodes.insert (std::make_pair(node_lab,i));
173  _assign_nodes[node_lab] = i;
174 
175  // do this irrespective whether MeshData exists
176  Node* newnode = mesh.add_point(xyz, i);
177 
178  // Add node to the nodes vector &
179  // tell the MeshData object the foreign node id.
180  if (this->_mesh_data != NULL)
181  this->_mesh_data->add_foreign_node_id (newnode, node_lab);
182  }
183 }
184 
185 
186 
187 //----------------------------------------------------------------------
188 // Function to read in the element table.
189 void TetGenIO::element_in (std::istream& ele_stream)
190 {
191  // Check input buffer
192  libmesh_assert (ele_stream.good());
193 
194  // Get a reference to the mesh
196 
197  // Read the elements from the ele_stream (*.ele file).
198  unsigned int element_lab=0, n_nodes=0, nAttri=0;
199 
200  ele_stream >> _num_elements // Read the number of tetrahedrons from the stream.
201  >> n_nodes // Read the number of nodes per tetrahedron from the stream (defaults to 4).
202  >> nAttri; // Read the number of attributes from stream.
203 
204  // Vector that assigns element nodes to their correct position.
205  // TetGen is normaly 0-based
206  // (right now this is strictly not necessary since it is the identity map,
207  // but in the future TetGen could change their numbering scheme.)
208  static const unsigned int assign_elm_nodes[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
209 
210  // If present, make room for element attributes to be stored.
211  this->element_attributes.resize(nAttri);
212  for (unsigned i=0; i<nAttri; ++i)
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;
221 
222  if (n_nodes==4)
223  elem = new Tet4;
224 
225  else if (n_nodes==10)
226  elem = new Tet10;
227 
228  else
229  {
230  libMesh::err << "Elements with " << n_nodes
231  << " nodes are not supported in the LibMesh tetgen module\n";
232  libmesh_error();
233  }
234  elem->set_id(i);
235 
236  mesh.add_elem (elem);
237 
238  libmesh_assert(elem);
239  libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
240 
241  // Read the element label
242  ele_stream >> element_lab;
243 
244  // Add the element to the mesh &
245  // tell the MeshData object the foreign element id
246  if (this->_mesh_data != NULL)
247  this->_mesh_data->add_foreign_elem_id (elem, element_lab);
248 
249  // Read node labels
250  for (dof_id_type j=0; j<n_nodes; j++)
251  {
252  dof_id_type node_label;
253  ele_stream >> node_label;
254 
255  // Assign node to element
256  elem->set_node(assign_elm_nodes[j]) =
257  mesh.node_ptr(_assign_nodes[node_label]);
258  }
259 
260  // Read and store attributes from the stream.
261  for (unsigned int j=0; j<nAttri; j++)
262  ele_stream >> this->element_attributes[j][i];
263  }
264 }
265 
266 
271 void TetGenIO::write (const std::string& fname)
272 {
273  // libmesh_assert three dimensions (should be extended later)
274  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().mesh_dimension(), 3);
275 
276  if (!(fname.rfind(".poly") < fname.size()))
277  {
278  libMesh::err << "ERROR: Unrecognized file name: "
279  << fname << std::endl;
280  libmesh_error();
281  }
282 
283  // Open the output file stream
284  std::ofstream out_stream (fname.c_str());
285 
286  // Make sure it opened correctly
287  if (!out_stream.good())
288  libmesh_file_error(fname.c_str());
289 
290  // Get a reference to the mesh
292 
293  // Begin interfacing with the .poly file
294  {
295  // header:
296  out_stream << "# poly file output generated by libmesh\n"
297  << mesh.n_nodes() << " 3 0 0\n";
298 
299  // write the nodes:
300  for (dof_id_type v=0; v<mesh.n_nodes(); v++)
301  out_stream << v << " "
302  << mesh.point(v)(0) << " "
303  << mesh.point(v)(1) << " "
304  << mesh.point(v)(2) << "\n";
305  }
306 
307  {
308  // write the connectivity:
309  out_stream << "# Facets:\n"
310  << mesh.n_elem() << " 0\n";
311 
312 // const_active_elem_iterator it (mesh.elements_begin());
313 // const const_active_elem_iterator end(mesh.elements_end());
314 
317 
318  for ( ; it != end; ++it)
319  out_stream << "1\n3 " // no. of facet polygons
320  // << (*it)->n_nodes() << " "
321  << (*it)->node(0) << " "
322  << (*it)->node(1) << " "
323  << (*it)->node(2) << "\n";
324  }
325 
326  // end of the file
327  out_stream << "0\n"; // no holes output!
328  out_stream << "\n\n# end of file\n";
329 }
330 
331 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo