ucd_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
33 #include "libmesh/enum_elem_type.h"
34 
35 #ifdef LIBMESH_HAVE_GZSTREAM
36 # include "libmesh/ignore_warnings.h" // shadowing in gzstream.h
37 # include "gzstream.h" // For reading/writing compressed streams
39 #endif
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  Point xyz;
157 
158  for (unsigned int i=0; i<nNodes; i++)
159  {
160  libmesh_assert (in.good());
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  // Build the node
168  mesh.add_point (xyz, i);
169  }
170  }
171 
172  // Read the elements from the stream. Notice that the UCD node-numbering
173  // scheme is 1-based, and we just created a 0-based scheme above
174  // (which is of course what we want). So, when we read in the nodal
175  // connectivity for each element we need to take 1 off the value of
176  // each node so that we get the right thing.
177  {
178  unsigned int material_id=0, node=0;
179  std::string type;
180 
181  for (unsigned int i=0; i<nElem; i++)
182  {
183  libmesh_assert (in.good());
184 
185  // The cell type can be either tri, quad, tet, hex, or prism.
186  in >> dummy // Cell number, means nothing to us
187  >> material_id // We'll use this for the element subdomain id.
188  >> type; // string describing cell type
189 
190  // Convert the UCD type string to a libmesh ElementType
191  auto it = _reading_element_map.find(type);
192  if (it == _reading_element_map.end())
193  libmesh_error_msg("Unsupported element type = " << type);
194 
195  // Build the required type and release it into our custody.
196  Elem * elem = Elem::build(it->second).release();
197 
198  for (unsigned int n=0; n<elem->n_nodes(); n++)
199  {
200  libmesh_assert (in.good());
201 
202  in >> node; // read the current node
203  node -= 1; // UCD is 1-based, so subtract
204 
205  libmesh_assert_less (node, mesh.n_nodes());
206 
207  // assign the node
208  elem->set_node(n) = mesh.node_ptr(node);
209  }
210 
211  elems_of_dimension[elem->dim()] = true;
212 
213  // Set the element's subdomain ID based on the material_id.
214  elem->subdomain_id() = cast_int<subdomain_id_type>(material_id);
215 
216  // Add the element to the mesh
217  elem->set_id(i);
218  mesh.add_elem (elem);
219  }
220 
221  // Set the mesh dimension to the largest encountered for an element
222  for (unsigned char i=0; i!=4; ++i)
223  if (elems_of_dimension[i])
225 
226 #if LIBMESH_DIM < 3
227  if (mesh.mesh_dimension() > LIBMESH_DIM)
228  libmesh_error_msg("Cannot open dimension " \
229  << mesh.mesh_dimension() \
230  << " mesh file when configured without " \
231  << mesh.mesh_dimension() \
232  << "D support.");
233 #endif
234  }
235 }
236 
237 
238 
239 void UCDIO::write_implementation (std::ostream & out_stream)
240 {
241  libmesh_assert (out_stream.good());
242 
244 
245  // UCD doesn't work any dimension except 3?
246  if (mesh.mesh_dimension() != 3)
247  libmesh_error_msg("Error: Can't write boundary elements for meshes of dimension less than 3. " \
248  << "Mesh dimension = " << mesh.mesh_dimension());
249 
250  // Write header
251  this->write_header(out_stream, mesh, mesh.n_elem(), 0);
252 
253  // Write the node coordinates
254  this->write_nodes(out_stream, mesh);
255 
256  // Write the elements
257  this->write_interior_elems(out_stream, mesh);
258 }
259 
260 
261 
262 void UCDIO::write_header(std::ostream & out_stream,
263  const MeshBase & mesh,
264  dof_id_type n_elems,
265  unsigned int n_vars)
266 {
267  libmesh_assert (out_stream.good());
268  // TODO: We could print out the libmesh revision used to write this file here.
269  out_stream << "# For a description of the UCD format see the AVS Developer's guide.\n"
270  << "#\n";
271 
272  // Write the mesh info
273  out_stream << mesh.n_nodes() << " "
274  << n_elems << " "
275  << n_vars << " "
276  << " 0 0\n";
277 }
278 
279 
280 
281 void UCDIO::write_nodes(std::ostream & out_stream,
282  const MeshBase & mesh)
283 {
284  // 1-based node number for UCD
285  unsigned int n=1;
286 
287  // Write the node coordinates
288  for (auto & node : mesh.node_ptr_range())
289  {
290  libmesh_assert (out_stream.good());
291 
292  out_stream << n++ << "\t";
293  node->write_unformatted(out_stream);
294  }
295 }
296 
297 
298 
299 void UCDIO::write_interior_elems(std::ostream & out_stream,
300  const MeshBase & mesh)
301 {
302  // 1-based element number for UCD
303  unsigned int e=1;
304 
305  // Write element information
306  for (const auto & elem : mesh.element_ptr_range())
307  {
308  libmesh_assert (out_stream.good());
309 
310  // Look up the corresponding UCD element type in the static map.
311  const ElemType etype = elem->type();
312  auto it = _writing_element_map.find(etype);
313  if (it == _writing_element_map.end())
314  libmesh_error_msg("Error: Unsupported ElemType " << etype << " for UCDIO.");
315 
316  // Write the element's subdomain ID as the UCD "material_id".
317  out_stream << e++ << " " << elem->subdomain_id() << " " << it->second << "\t";
318  elem->write_connectivity(out_stream, UCD);
319  }
320 }
321 
322 
323 
324 void UCDIO::write_nodal_data(const std::string & fname,
325  const std::vector<Number> & soln,
326  const std::vector<std::string> & names)
327 {
329 
330  const dof_id_type n_elem = mesh.n_elem();
331 
332  // Only processor 0 does the writing
333  if (mesh.processor_id())
334  return;
335 
336  std::ofstream out_stream(fname.c_str());
337 
338  // UCD doesn't work in 1D
339  libmesh_assert (mesh.mesh_dimension() != 1);
340 
341  // Write header
342  this->write_header(out_stream,mesh,n_elem,
343  cast_int<unsigned int>(names.size()));
344 
345  // Write the node coordinates
346  this->write_nodes(out_stream, mesh);
347 
348  // Write the elements
349  this->write_interior_elems(out_stream, mesh);
350 
351  // Write the solution
352  this->write_soln(out_stream, mesh, names, soln);
353 }
354 
355 
356 
357 void UCDIO::write_soln(std::ostream & out_stream,
358  const MeshBase & mesh,
359  const std::vector<std::string> & names,
360  const std::vector<Number> & soln)
361 {
362  libmesh_assert (out_stream.good());
363 
364  // First write out how many variables and how many components per variable
365  out_stream << names.size();
366  for (std::size_t i = 0; i < names.size(); i++)
367  {
368  libmesh_assert (out_stream.good());
369  // Each named variable has only 1 component
370  out_stream << " 1";
371  }
372  out_stream << std::endl;
373 
374  // Now write out variable names and units. Since we don't store units
375  // We just write out dummy.
376  for (const auto & name : names)
377  {
378  libmesh_assert (out_stream.good());
379  out_stream << name << ", dummy" << std::endl;
380  }
381 
382  // Now, for each node, write out the solution variables.
383  // We use a 1-based node numbering for UCD.
384  std::size_t nv = names.size();
385  for (std::size_t n = 1; n <= mesh.n_nodes(); n++)
386  {
387  libmesh_assert (out_stream.good());
388  out_stream << n;
389 
390  for (std::size_t var = 0; var != nv; var++)
391  {
392  std::size_t idx = nv*(n-1) + var;
393 
394  out_stream << " " << soln[idx];
395  }
396  out_stream << std::endl;
397  }
398 }
399 
400 } // namespace libMesh
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
const MT & mesh() const
Definition: mesh_output.h:234
virtual void read(const std::string &) override
Definition: ucd_io.C:82
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2024
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Definition: mesh_tools.C:702
std::vector< bool > elems_of_dimension
Definition: mesh_input.h:97
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
void write_soln(std::ostream &out, const MeshBase &mesh, const std::vector< std::string > &names, const std::vector< Number > &soln)
Definition: ucd_io.C:357
static std::map< ElemType, std::string > _writing_element_map
Definition: ucd_io.h:149
const unsigned int n_vars
Definition: tecplot_io.C:69
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
void skip_comment_lines(std::istream &in, const char comment_start)
Definition: mesh_input.h:179
Base class for Mesh.
Definition: mesh_base.h:77
dof_id_type & set_id()
Definition: dof_object.h:664
void read_implementation(std::istream &in_stream)
Definition: ucd_io.C:124
virtual SimpleRange< element_iterator > element_ptr_range()=0
virtual unsigned int n_nodes() const =0
virtual Elem * add_elem(Elem *e)=0
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:245
static std::map< ElemType, std::string > build_writing_element_map()
Definition: ucd_io.C:52
virtual SimpleRange< node_iterator > node_ptr_range()=0
void set_mesh_dimension(unsigned char d)
Definition: mesh_base.h:213
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
virtual unsigned short dim() const =0
void write_nodes(std::ostream &out, const MeshBase &mesh)
Definition: ucd_io.C:281
static std::map< std::string, ElemType > build_reading_element_map()
Definition: ucd_io.C:68
unsigned int mesh_dimension() const
Definition: mesh_base.C:126
void write_interior_elems(std::ostream &out, const MeshBase &mesh)
Definition: ucd_io.C:299
void write_implementation(std::ostream &out_stream)
Definition: ucd_io.C:239
static std::map< std::string, ElemType > _reading_element_map
Definition: ucd_io.h:153
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
Definition: ucd_io.C:324
virtual void write(const std::string &) override
Definition: ucd_io.C:103
A geometric point in (x,y,z) space.
Definition: point.h:38
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)
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)
Definition: ucd_io.C:262