libMesh
ucd_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/libmesh_config.h"
21 #include "libmesh/ucd_io.h"
22 #include "libmesh/mesh_base.h"
23 #include "libmesh/face_quad4.h"
24 #include "libmesh/face_tri3.h"
25 #include "libmesh/cell_tet4.h"
26 #include "libmesh/cell_hex8.h"
27 #include "libmesh/cell_prism6.h"
28 #include "libmesh/enum_io_package.h"
29 #include "libmesh/enum_elem_type.h"
30 #include "libmesh/int_range.h"
31 #include "libmesh/utility.h"
32 
33 #ifdef LIBMESH_HAVE_GZSTREAM
34 # include "gzstream.h" // For reading/writing compressed streams
35 #endif
36 
37 // C++ includes
38 #include <array>
39 #include <fstream>
40 
41 
42 namespace libMesh
43 {
44 
45 // Initialize the static data members by calling the static build functions.
46 std::map<ElemType, std::string> UCDIO::_writing_element_map = UCDIO::build_writing_element_map();
47 std::map<std::string, ElemType> UCDIO::_reading_element_map = UCDIO::build_reading_element_map();
48 
49 
50 
51 // Static function used to build the _writing_element_map.
52 std::map<ElemType, std::string> UCDIO::build_writing_element_map()
53 {
54  std::map<ElemType, std::string> ret;
55  ret[EDGE2] = "edge";
56  ret[TRI3] = "tri";
57  ret[QUAD4] = "quad";
58  ret[TET4] = "tet";
59  ret[HEX8] = "hex";
60  ret[PRISM6] = "prism";
61  ret[PYRAMID5] = "pyramid";
62  return ret;
63 }
64 
65 
66 
67 // Static function used to build the _reading_element_map.
68 std::map<std::string, ElemType> UCDIO::build_reading_element_map()
69 {
70  std::map<std::string, ElemType> ret;
71  ret["edge"] = EDGE2;
72  ret["tri"] = TRI3;
73  ret["quad"] = QUAD4;
74  ret["tet"] = TET4;
75  ret["hex"] = HEX8;
76  ret["prism"] = PRISM6;
77  ret["pyramid"] = PYRAMID5;
78  return ret;
79 }
80 
81 
82 void UCDIO::read (const std::string & file_name)
83 {
84  if (file_name.rfind(".gz") < file_name.size())
85  {
86 #ifdef LIBMESH_HAVE_GZSTREAM
87  igzstream in_stream (file_name.c_str());
88  this->read_implementation (in_stream);
89 #else
90  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
91 #endif
92  }
93 
94  else
95  {
96  std::ifstream in_stream (file_name.c_str());
97  this->read_implementation (in_stream);
98  }
99 }
100 
101 
102 
103 void UCDIO::write (const std::string & file_name)
104 {
105  if (file_name.rfind(".gz") < file_name.size())
106  {
107 #ifdef LIBMESH_HAVE_GZSTREAM
108  ogzstream out_stream (file_name.c_str());
109  this->write_implementation (out_stream);
110 #else
111  libmesh_error_msg("ERROR: You must have the zlib.h header files and libraries to read and write compressed streams.");
112 #endif
113  }
114 
115  else
116  {
117  std::ofstream out_stream (file_name.c_str());
118  this->write_implementation (out_stream);
119  }
120 }
121 
122 
123 
124 void UCDIO::read_implementation (std::istream & in)
125 {
126  // This is a serial-only process for now;
127  // the Mesh should be read on processor 0 and
128  // broadcast later
129  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
130 
131  // Check input buffer
132  libmesh_assert (in.good());
133 
135 
136  // Keep track of what kinds of elements this file contains
137  elems_of_dimension.clear();
138  elems_of_dimension.resize(4, false);
139 
140  this->skip_comment_lines (in, '#');
141 
142  unsigned int nNodes=0, nElem=0, dummy=0;
143 
144  in >> nNodes // Read the number of nodes from the stream
145  >> nElem // Read the number of elements from the stream
146  >> dummy
147  >> dummy
148  >> dummy;
149 
150 
151  // Read the nodal coordinates. Note that UCD format always
152  // stores (x,y,z), and in 2D z=0. We don't need to store this,
153  // however. So, we read in x,y,z for each node and make a point
154  // in the proper way based on what dimension we're in
155  {
156  for (unsigned int i=0; i<nNodes; i++)
157  {
158  libmesh_assert (in.good());
159 
160  std::array<Real, 3> xyz;
161 
162  in >> dummy // Point number
163  >> xyz[0] // x-coordinate value
164  >> xyz[1] // y-coordinate value
165  >> xyz[2]; // z-coordinate value
166 
167  Point p(xyz[0]);
168 #if LIBMESH_DIM > 1
169  p(1) = xyz[1];
170 #else
171  libmesh_assert_equal_to(xyz[1], 0);
172 #endif
173 #if LIBMESH_DIM > 2
174  p(2) = xyz[2];
175 #else
176  libmesh_assert_equal_to(xyz[2], 0);
177 #endif
178 
179  // Build the node
180  mesh.add_point (p, i);
181  }
182  }
183 
184  // Read the elements from the stream. Notice that the UCD node-numbering
185  // scheme is 1-based, and we just created a 0-based scheme above
186  // (which is of course what we want). So, when we read in the nodal
187  // connectivity for each element we need to take 1 off the value of
188  // each node so that we get the right thing.
189  {
190  unsigned int material_id=0, node=0;
191  std::string type;
192 
193  for (unsigned int i=0; i<nElem; i++)
194  {
195  libmesh_assert (in.good());
196 
197  // The cell type can be either tri, quad, tet, hex, or prism.
198  in >> dummy // Cell number, means nothing to us
199  >> material_id // We'll use this for the element subdomain id.
200  >> type; // string describing cell type
201 
202  // Convert the UCD type string to a libmesh ElemType
203  ElemType et = libmesh_map_find(_reading_element_map, type);
204 
205  // Build the required type and release it into our custody.
206  auto elem = Elem::build(et);
207 
208  for (auto n : elem->node_index_range())
209  {
210  libmesh_assert (in.good());
211 
212  in >> node; // read the current node
213  node -= 1; // UCD is 1-based, so subtract
214 
215  libmesh_assert_less (node, mesh.n_nodes());
216 
217  // assign the node
218  elem->set_node(n, mesh.node_ptr(node));
219  }
220 
221  elems_of_dimension[elem->dim()] = true;
222 
223  // Set the element's subdomain ID based on the material_id.
224  elem->subdomain_id() = cast_int<subdomain_id_type>(material_id);
225 
226  // Add the element to the mesh
227  elem->set_id(i);
228  mesh.add_elem (std::move(elem));
229  }
230 
231  // Set the mesh dimension to the largest encountered for an element
232  for (unsigned char i=0; i!=4; ++i)
233  if (elems_of_dimension[i])
235 
236 #if LIBMESH_DIM < 3
237  libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
238  "Cannot open dimension "
239  << mesh.mesh_dimension()
240  << " mesh file when configured without "
241  << mesh.mesh_dimension()
242  << "D support.");
243 #endif
244  }
245 }
246 
247 
248 
249 void UCDIO::write_implementation (std::ostream & out_stream)
250 {
251  libmesh_assert (out_stream.good());
252 
254 
255  // UCD doesn't work any dimension except 3?
256  libmesh_error_msg_if(mesh.mesh_dimension() != 3,
257  "Error: Can't write boundary elements for meshes of dimension less than 3. "
258  "Mesh dimension = " << mesh.mesh_dimension());
259 
260  // Write header
261  this->write_header(out_stream, mesh, mesh.n_elem(), 0);
262 
263  // Write the node coordinates
264  this->write_nodes(out_stream, mesh);
265 
266  // Write the elements
267  this->write_interior_elems(out_stream, mesh);
268 }
269 
270 
271 
272 void UCDIO::write_header(std::ostream & out_stream,
273  const MeshBase & mesh,
274  dof_id_type n_elems,
275  unsigned int n_vars)
276 {
277  libmesh_assert (out_stream.good());
278  // TODO: We could print out the libmesh revision used to write this file here.
279  out_stream << "# For a description of the UCD format see the AVS Developer's guide.\n"
280  << "#\n";
281 
282  // Write the mesh info
283  out_stream << mesh.n_nodes() << " "
284  << n_elems << " "
285  << n_vars << " "
286  << " 0 0\n";
287 }
288 
289 
290 
291 void UCDIO::write_nodes(std::ostream & out_stream,
292  const MeshBase & mesh)
293 {
294  // 1-based node number for UCD
295  unsigned int n=1;
296 
297  // Write the node coordinates
298  for (auto & node : mesh.node_ptr_range())
299  {
300  libmesh_assert (out_stream.good());
301 
302  out_stream << n++ << "\t";
303  node->write_unformatted(out_stream);
304  }
305 }
306 
307 
308 
309 void UCDIO::write_interior_elems(std::ostream & out_stream,
310  const MeshBase & mesh)
311 {
312  // 1-based element number for UCD
313  unsigned int e=1;
314 
315  // Write element information
316  for (const auto & elem : mesh.element_ptr_range())
317  {
318  libmesh_assert (out_stream.good());
319 
320  // Look up the corresponding UCD element type in the static map.
321  std::string elem_string = libmesh_map_find(_writing_element_map, elem->type());
322 
323  // Write the element's subdomain ID as the UCD "material_id".
324  out_stream << e++ << " " << elem->subdomain_id() << " " << elem_string << "\t";
325  elem->write_connectivity(out_stream, UCD);
326  }
327 }
328 
329 
330 
331 void UCDIO::write_nodal_data(const std::string & fname,
332  const std::vector<Number> & soln,
333  const std::vector<std::string> & names)
334 {
336 
337  const dof_id_type n_elem = mesh.n_elem();
338 
339  // Only processor 0 does the writing
340  if (mesh.processor_id())
341  return;
342 
343  std::ofstream out_stream(fname.c_str());
344 
345  // UCD doesn't work in 1D
347 
348  // Write header
349  this->write_header(out_stream,mesh,n_elem,
350  cast_int<unsigned int>(names.size()));
351 
352  // Write the node coordinates
353  this->write_nodes(out_stream, mesh);
354 
355  // Write the elements
356  this->write_interior_elems(out_stream, mesh);
357 
358  // Write the solution
359  this->write_soln(out_stream, mesh, names, soln);
360 }
361 
362 
363 
364 void UCDIO::write_soln(std::ostream & out_stream,
365  const MeshBase & mesh,
366  const std::vector<std::string> & names,
367  const std::vector<Number> & soln)
368 {
369  libmesh_assert (out_stream.good());
370 
371  // First write out how many variables and how many components per variable
372  out_stream << names.size();
373  for (std::size_t i = 0, ns = names.size(); i < ns; i++)
374  {
375  libmesh_assert (out_stream.good());
376  // Each named variable has only 1 component
377  out_stream << " 1";
378  }
379  out_stream << std::endl;
380 
381  // Now write out variable names and units. Since we don't store units
382  // We just write out dummy.
383  for (const auto & name : names)
384  {
385  libmesh_assert (out_stream.good());
386  out_stream << name << ", dummy" << std::endl;
387  }
388 
389  // Now, for each node, write out the solution variables.
390  // We use a 1-based node numbering for UCD.
391  std::size_t nv = names.size();
392  for (auto n : IntRange<std::size_t>(1, mesh.n_nodes()))
393  {
394  libmesh_assert (out_stream.good());
395  out_stream << n;
396 
397  for (std::size_t var = 0; var != nv; var++)
398  {
399  std::size_t idx = nv*(n-1) + var;
400 
401  out_stream << " " << soln[idx];
402  }
403  out_stream << std::endl;
404  }
405 }
406 
407 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
const MT & mesh() const
Definition: mesh_output.h:259
virtual void read(const std::string &) override
This method implements reading a mesh from a specified file in UCD format.
Definition: ucd_io.C:82
ElemType
Defines an enum for geometric element types.
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
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
std::vector< bool > elems_of_dimension
A vector of bools describing what dimension elements have been encountered when reading a mesh...
Definition: mesh_input.h:107
MeshBase & mesh
void write_soln(std::ostream &out, const MeshBase &mesh, const std::vector< std::string > &names, const std::vector< Number > &soln)
Writes all nodal solution variables.
Definition: ucd_io.C:364
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
static std::map< ElemType, std::string > _writing_element_map
Definition: ucd_io.h:141
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
void read_implementation(std::istream &in_stream)
The actual implementation of the read function.
Definition: ucd_io.C:124
unsigned int n_vars
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
static std::map< ElemType, std::string > build_writing_element_map()
Definition: ucd_io.C:52
libmesh_assert(ctx)
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
void write_nodes(std::ostream &out, const MeshBase &mesh)
Write node information.
Definition: ucd_io.C:291
static std::map< std::string, ElemType > build_reading_element_map()
Definition: ucd_io.C:68
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void write_interior_elems(std::ostream &out, const MeshBase &mesh)
Write element information.
Definition: ucd_io.C:309
void write_implementation(std::ostream &out_stream)
The actual implementation of the write function.
Definition: ucd_io.C:249
static std::map< std::string, ElemType > _reading_element_map
Definition: ucd_io.h:145
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
virtual void write_nodal_data(const std::string &fname, const std::vector< Number > &soln, const std::vector< std::string > &names) override
This method implements writing a mesh and solution to a specified file in UCD format.
Definition: ucd_io.C:331
virtual void write(const std::string &) override
This method implements writing a mesh to a specified file in UCD format.
Definition: ucd_io.C:103
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual dof_id_type n_nodes() const =0
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
uint8_t dof_id_type
Definition: id_types.h:67
void write_header(std::ostream &out, const MeshBase &mesh, dof_id_type n_elems, unsigned int n_vars)
Write UCD format header.
Definition: ucd_io.C:272