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