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 :
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 2 : 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 1 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
42 :
43 2 : std::string name_node, name_ele, dummy;
44 :
45 : // tetgen only works in 3D
46 2 : 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 2 : if (Utility::contains(name, ".node"))
56 : {
57 0 : name_node = name;
58 0 : dummy = name;
59 0 : std::size_t position = dummy.rfind(".node");
60 0 : name_ele = dummy.replace(position, 5, ".ele");
61 : }
62 2 : else if (Utility::contains(name, ".ele"))
63 : {
64 1 : name_ele = name;
65 1 : dummy = name;
66 1 : std::size_t position = dummy.rfind(".ele");
67 1 : name_node = dummy.replace(position, 4, ".node");
68 : }
69 : else
70 0 : libmesh_error_msg("ERROR: Unrecognized file name: " << name);
71 :
72 :
73 :
74 : // Set the streams from which to read in
75 4 : std::ifstream node_stream (name_node.c_str());
76 4 : std::ifstream ele_stream (name_ele.c_str());
77 :
78 2 : 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 1 : libMesh::out<< "TetGenIO found the tetgen files to read " <<std::endl;
85 :
86 : // Skip the comment lines at the beginning
87 2 : this->skip_comment_lines (node_stream, '#');
88 2 : this->skip_comment_lines (ele_stream, '#');
89 :
90 : // Read the nodes and elements from the streams
91 2 : this->read_nodes_and_elem (node_stream, ele_stream);
92 1 : libMesh::out<< "TetGenIO read in nodes and elements " <<std::endl;
93 2 : }
94 :
95 :
96 :
97 2 : void TetGenIO::read_nodes_and_elem (std::istream & node_stream,
98 : std::istream & ele_stream)
99 : {
100 2 : _num_nodes = 0;
101 2 : _num_elements = 0;
102 :
103 : // Read all the datasets.
104 2 : this->node_in (node_stream);
105 2 : this->element_in (ele_stream);
106 :
107 : // some more clean-up
108 1 : _assign_nodes.clear();
109 2 : }
110 :
111 :
112 :
113 : //----------------------------------------------------------------------
114 : // Function to read in the node table.
115 2 : void TetGenIO::node_in (std::istream & node_stream)
116 : {
117 : // Check input buffer
118 1 : libmesh_assert (node_stream.good());
119 :
120 : // Get a reference to the mesh
121 2 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
122 :
123 2 : unsigned int dimension=0, nAttributes=0, BoundaryMarkers=0;
124 :
125 2 : node_stream >> _num_nodes // Read the number of nodes from the stream
126 1 : >> dimension // Read the dimension from the stream
127 1 : >> nAttributes // Read the number of attributes from stream
128 1 : >> 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 2 : unsigned int node_lab=0;
132 : Real dummy;
133 :
134 : // If present, make room for node attributes to be stored.
135 2 : this->node_attributes.resize(nAttributes);
136 2 : for (unsigned i=0; i<nAttributes; ++i)
137 0 : this->node_attributes[i].resize(_num_nodes);
138 :
139 :
140 22 : for (unsigned int i=0; i<_num_nodes; i++)
141 : {
142 : // Check input buffer
143 10 : libmesh_assert (node_stream.good());
144 :
145 : std::array<Real, 3> xyz;
146 :
147 10 : node_stream >> node_lab // node number
148 10 : >> xyz[0] // x-coordinate value
149 10 : >> xyz[1] // y-coordinate value
150 10 : >> xyz[2]; // z-coordinate value
151 :
152 : // Read and store attributes from the stream.
153 20 : for (unsigned int j=0; j<nAttributes; j++)
154 0 : node_stream >> node_attributes[j][i];
155 :
156 : // Read (and discard) boundary marker if BoundaryMarker=1.
157 : // TODO: should we store this somehow?
158 20 : if (BoundaryMarkers == 1)
159 0 : node_stream >> dummy;
160 :
161 : // Store the new position of the node under its label.
162 : //_assign_nodes.emplace(node_lab,i);
163 20 : _assign_nodes[node_lab] = i;
164 :
165 20 : Point p(xyz[0]);
166 : #if LIBMESH_DIM > 1
167 20 : p(1) = xyz[1];
168 : #else
169 : libmesh_assert_equal_to(xyz[1], 0);
170 : #endif
171 : #if LIBMESH_DIM > 2
172 20 : 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 20 : mesh.add_point(p, i);
179 : }
180 2 : }
181 :
182 :
183 :
184 : //----------------------------------------------------------------------
185 : // Function to read in the element table.
186 2 : void TetGenIO::element_in (std::istream & ele_stream)
187 : {
188 : // Check input buffer
189 1 : libmesh_assert (ele_stream.good());
190 :
191 : // Get a reference to the mesh
192 2 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
193 :
194 : // Read the elements from the ele_stream (*.ele file).
195 2 : unsigned int element_lab=0, n_nodes=0, region_attribute=0;
196 :
197 2 : ele_stream >> _num_elements // Read the number of tetrahedrons from the stream.
198 1 : >> n_nodes // Read the number of nodes per tetrahedron from the stream (defaults to 4).
199 1 : >> 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 2 : 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 4 : for (dof_id_type i=0; i<_num_elements; i++)
216 : {
217 1 : libmesh_assert (ele_stream.good());
218 :
219 : // TetGen only supports Tet4 and Tet10 elements.
220 1 : Elem * elem = nullptr;
221 :
222 2 : if (n_nodes==4)
223 0 : elem = mesh.add_elem(Elem::build_with_id(TET4, i));
224 :
225 2 : else if (n_nodes==10)
226 3 : elem = mesh.add_elem(Elem::build_with_id(TET10, i));
227 :
228 : else
229 0 : libmesh_error_msg("Elements with " << n_nodes <<
230 : " nodes are not supported in the LibMesh tetgen module.");
231 :
232 1 : libmesh_assert(elem);
233 1 : 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 1 : ele_stream >> element_lab;
240 :
241 : // Read node labels
242 22 : for (dof_id_type j=0; j<n_nodes; j++)
243 : {
244 : dof_id_type node_label;
245 10 : ele_stream >> node_label;
246 :
247 : // Assign node to element
248 20 : elem->set_node(assign_elm_nodes[j],
249 20 : 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 2 : if (region_attribute)
254 : {
255 : unsigned int region;
256 0 : 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 0 : elem->subdomain_id() = cast_int<subdomain_id_type>(region);
261 : }
262 : }
263 2 : }
264 :
265 :
266 : /**
267 : * This method implements writing a mesh to a specified ".poly" file.
268 : * ".poly" files defines so called Piecewise Linear Complex (PLC).
269 : */
270 0 : void TetGenIO::write (const std::string & fname)
271 : {
272 : // libmesh_assert three dimensions (should be extended later)
273 0 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().mesh_dimension(), 3);
274 :
275 0 : libmesh_error_msg_if(!Utility::contains(fname, ".poly"),
276 : "ERROR: Unrecognized file name: " << fname);
277 :
278 : // Open the output file stream
279 0 : std::ofstream out_stream (fname.c_str());
280 :
281 : // Make sure it opened correctly
282 0 : if (!out_stream.good())
283 0 : libmesh_file_error(fname.c_str());
284 :
285 : // Get a reference to the mesh
286 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
287 :
288 : // Begin interfacing with the .poly file
289 : {
290 : // header:
291 0 : out_stream << "# poly file output generated by libmesh\n"
292 0 : << mesh.n_nodes() << " 3 0 0\n";
293 :
294 : // write the nodes:
295 0 : for (auto v : make_range(mesh.n_nodes()))
296 0 : out_stream << v << " "
297 0 : << mesh.point(v)(0) << " "
298 0 : << mesh.point(v)(1) << " "
299 0 : << mesh.point(v)(2) << "\n";
300 : }
301 :
302 : {
303 : // write the connectivity:
304 0 : out_stream << "# Facets:\n"
305 0 : << mesh.n_elem() << " 0\n";
306 :
307 0 : for (const auto & elem : mesh.active_element_ptr_range())
308 0 : out_stream << "1\n3 " // no. of facet polygons
309 : // << elem->n_nodes() << " "
310 0 : << elem->node_id(0) << " "
311 0 : << elem->node_id(1) << " "
312 0 : << elem->node_id(2) << "\n";
313 : }
314 :
315 : // end of the file
316 0 : out_stream << "0\n"; // no holes output!
317 0 : out_stream << "\n\n# end of file\n";
318 0 : }
319 :
320 : } // namespace libMesh
|
