ucd_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 
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 
42 
43 
44 // ------------------------------------------------------------
45 // UCDIO class members
46 void UCDIO::read (const std::string& file_name)
47 {
48  if (file_name.rfind(".gz") < file_name.size())
49  {
50 #ifdef LIBMESH_HAVE_GZSTREAM
51 
52  igzstream in_stream (file_name.c_str());
53  this->read_implementation (in_stream);
54 
55 #else
56 
57  libMesh::err << "ERROR: You must have the zlib.h header "
58  << "files and libraries to read and write "
59  << "compressed streams."
60  << std::endl;
61  libmesh_error();
62 
63 #endif
64  return;
65  }
66 
67  else
68  {
69  std::ifstream in_stream (file_name.c_str());
70  this->read_implementation (in_stream);
71  return;
72  }
73 }
74 
75 
76 
77 void UCDIO::write (const std::string& file_name)
78 {
79  if (file_name.rfind(".gz") < file_name.size())
80  {
81 #ifdef LIBMESH_HAVE_GZSTREAM
82 
83  ogzstream out_stream (file_name.c_str());
84  this->write_implementation (out_stream);
85 
86 #else
87 
88  libMesh::err << "ERROR: You must have the zlib.h header "
89  << "files and libraries to read and write "
90  << "compressed streams."
91  << std::endl;
92  libmesh_error();
93 
94 #endif
95  return;
96  }
97 
98  else
99  {
100  std::ofstream out_stream (file_name.c_str());
101  this->write_implementation (out_stream);
102  return;
103  }
104 }
105 
106 
107 
108 void UCDIO::read_implementation (std::istream& in)
109 {
110  // This is a serial-only process for now;
111  // the Mesh should be read on processor 0 and
112  // broadcast later
113  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
114 
115  // Check input buffer
116  libmesh_assert (in.good());
117 
119 
120  // Keep track of what kinds of elements this file contains
121  elems_of_dimension.clear();
122  elems_of_dimension.resize(4, false);
123 
124  this->skip_comment_lines (in, '#');
125 
126  unsigned int nNodes=0, nElem=0, dummy=0;
127 
128  in >> nNodes // Read the number of nodes from the stream
129  >> nElem // Read the number of elements from the stream
130  >> dummy
131  >> dummy
132  >> dummy;
133 
134 
135  // Read the nodal coordinates. Note that UCD format always
136  // stores (x,y,z), and in 2D z=0. We don't need to store this,
137  // however. So, we read in x,y,z for each node and make a point
138  // in the proper way based on what dimension we're in
139  {
140  Point xyz;
141 
142  for (unsigned int i=0; i<nNodes; i++)
143  {
144  libmesh_assert (in.good());
145 
146  in >> dummy // Point number
147  >> xyz(0) // x-coordinate value
148  >> xyz(1) // y-coordinate value
149  >> xyz(2); // z-coordinate value
150 
151  // Build the node
152  mesh.add_point (xyz, i);
153  }
154  }
155 
156 
157 
158  // Read the elements from the stream. Notice that the UCD node-numbering
159  // scheme is 1-based, and we just created a 0-based scheme above
160  // (which is of course what we want). So, when we read in the nodal
161  // connectivity for each element we need to take 1 off the value of
162  // each node so that we get the right thing.
163  {
164  unsigned int material_id=0, node=0;
165  std::string type;
166 
167  for (unsigned int i=0; i<nElem; i++)
168  {
169  Elem* elem = NULL;
170 
171  libmesh_assert (in.good());
172 
173  in >> dummy // Cell number, means nothing to us
174  >> material_id // doesn't mean anything at present, might later
175  >> type; // string describing cell type:
176  // either tri, quad, tet, hex, or prism for the
177  // obvious cases
178 
179 
180  // Now read the connectivity.
181  if (type == "quad")
182  elem = new Quad4;
183  else if (type == "tri")
184  elem = new Tri3;
185  else if (type == "hex")
186  elem = new Hex8;
187  else if (type == "tet")
188  elem = new Tet4;
189  else if (type == "prism")
190  elem = new Prism6;
191  else
192  libmesh_error();
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  elem->set_node(n) =
204  mesh.node_ptr(node); // assign the node
205  }
206 
207  elems_of_dimension[elem->dim()] = true;
208 
209  // Add the element to the mesh
210  elem->set_id(i);
211  mesh.add_elem (elem);
212  }
213 
214  // Set the mesh dimension to the largest encountered for an element
215  for (unsigned int i=0; i!=4; ++i)
216  if (elems_of_dimension[i])
217  mesh.set_mesh_dimension(i);
218 
219 #if LIBMESH_DIM < 3
220  if (mesh.mesh_dimension() > LIBMESH_DIM)
221  {
222  libMesh::err << "Cannot open dimension " <<
223  mesh.mesh_dimension() <<
224  " mesh file when configured without " <<
225  mesh.mesh_dimension() << "D support." <<
226  std::endl;
227  libmesh_error();
228  }
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 in 1D
242  libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
243  if(mesh.mesh_dimension() != 3)
244  {
245  libMesh::err << "Error: Can't write boundary elements for meshes of dimension less than 3"
246  << "Mesh dimension = " << mesh.mesh_dimension()
247  << std::endl;
248  libmesh_error();
249  }
250 
251  // Write header
252  this->write_header(out_stream,mesh,mesh.n_elem(),0);
253 
254  // Write the node coordinates
255  this->write_nodes(out_stream,mesh);
256 
257  // Write the elements
258  this->write_interior_elems(out_stream,mesh);
259 
260  return;
261 }
262 
263 void UCDIO::write_header(std::ostream& out_stream, const MeshBase& mesh,
264  dof_id_type n_elems, unsigned int n_vars )
265 {
266  libmesh_assert (out_stream.good());
267  // TODO: We used to print out the SVN revision here when we did keyword expansions...
268  out_stream << "# For a description of the UCD format see the AVS Developer's guide.\n"
269  << "#\n";
270 
271  // Write the mesh info
272  out_stream << mesh.n_nodes() << " "
273  << n_elems << " "
274  << n_vars << " "
275  << " 0 0\n";
276  return;
277 }
278 
279 void UCDIO::write_nodes(std::ostream& out_stream, const MeshBase& mesh)
280 {
283 
284  unsigned int n=1; // 1-based node number for UCD
285 
286  // Write the node coordinates
287  for (; it != end; ++it)
288  {
289  libmesh_assert (out_stream.good());
290 
291  out_stream << n++ << "\t";
292  (*it)->write_unformatted(out_stream);
293  }
294 
295  return;
296 }
297 
298 void UCDIO::write_interior_elems(std::ostream& out_stream, const MeshBase& mesh)
299 {
300  std::string type[] =
301  { "edge", "edge", "edge",
302  "tri", "tri",
303  "quad", "quad", "quad",
304  "tet", "tet",
305  "hex", "hex", "hex",
306  "prism", "prism", "prism",
307  "pyramid" };
308 
311 
312  unsigned int e=1; // 1-based element number for UCD
313 
314  // Write element information
315  for (; it != end; ++it)
316  {
317  libmesh_assert (out_stream.good());
318 
319  // PB: I believe these are the only supported ElemTypes.
320  const ElemType etype = (*it)->type();
321  if( (etype != TRI3) && (etype != QUAD4) &&
322  (etype != TET4) && (etype != HEX8) &&
323  (etype != PRISM6) && (etype != PYRAMID5) )
324  {
325  libMesh::err << "Error: Unsupported ElemType for UCDIO."
326  << std::endl;
327  libmesh_error();
328  }
329 
330  out_stream << e++ << " 0 " << type[etype] << "\t";
331  // (*it)->write_ucd_connectivity(out_stream);
332  (*it)->write_connectivity(out_stream, UCD);
333  }
334 
335  return;
336 }
337 
338 void UCDIO::write_nodal_data(const std::string& fname,
339  const std::vector<Number>&soln,
340  const std::vector<std::string>& names)
341 {
343 
344  const dof_id_type n_elem = mesh.n_elem();
345 
346  // Only processor 0 does the writing
347  if (mesh.processor_id())
348  return;
349 
350  std::ofstream out_stream(fname.c_str());
351 
352  // UCD doesn't work in 1D
353  libmesh_assert (mesh.mesh_dimension() != 1);
354 
355  // Write header
356  this->write_header(out_stream,mesh,n_elem,names.size());
357 
358  // Write the node coordinates
359  this->write_nodes(out_stream,mesh);
360 
361  // Write the elements
362  this->write_interior_elems(out_stream,mesh);
363 
364  // Write the solution
365  this->write_soln(out_stream,mesh,names,soln);
366 
367  return;
368 }
369 
370 void UCDIO::write_soln(std::ostream& out_stream, const MeshBase& mesh,
371  const std::vector<std::string>& names,
372  const std::vector<Number>&soln)
373 {
374  libmesh_assert (out_stream.good());
375 
376  // First write out how many variables and how many components per variable
377  out_stream << names.size();
378  for( unsigned int i = 0; i < names.size(); i++ )
379  {
380  libmesh_assert (out_stream.good());
381  // Each named variable has only 1 component
382  out_stream << " 1";
383  }
384  out_stream << std::endl;
385 
386  // Now write out variable names and units. Since we don't store units
387  // We just write out dummy.
388  for( std::vector<std::string>::const_iterator var = names.begin();
389  var != names.end();
390  var++)
391  {
392  libmesh_assert (out_stream.good());
393  out_stream << (*var) << ", dummy" << std::endl;
394  }
395 
396  // Now, for each node, write out the solution variables
397  std::size_t nv = names.size();
398  for( std::size_t n = 1; // 1-based node number for UCD
399  n <= mesh.n_nodes(); n++)
400  {
401  libmesh_assert (out_stream.good());
402  out_stream << n;
403 
404  for( std::size_t var = 0; var != nv; var++ )
405  {
406  std::size_t idx = nv*(n-1) + var;
407 
408  out_stream << " " << soln[idx];
409  }
410  out_stream << std::endl;
411  }
412 
413  return;
414 }
415 
416 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo