diva_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2016 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 // Local includes
23 #include "libmesh/diva_io.h"
24 #include "libmesh/boundary_mesh.h"
25 #include "libmesh/mesh_tools.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/boundary_info.h"
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // DivaIO class members
34 void DivaIO::write (const std::string & fname)
35 {
36  // We may need to gather a ParallelMesh to output it, making that
37  // const qualifier in our constructor a dirty lie
38  MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
39 
40  // Open the output file stream
41  std::ofstream out_file(fname.c_str());
42 
43  // Make sure it opened correctly
44  if (!out_file.good())
45  libmesh_file_error(fname.c_str());
46 
47  this->write_stream (out_file);
48 }
49 
50 
51 
52 
53 void DivaIO::write_stream (std::ostream & out_file)
54 {
55  /*
56  From Kelly: (kelly@tacc.utexas.edu)
57 
58  Ok, the following is the format:
59 
60  #points #triangles #quads #tets #prisms #pyramids #hexs
61  loop over all points (written out x y z x y z ...)
62  loop over all triangles (written out i1 i2 i3) (These are indices into
63  the points going from
64  1 to #points)
65  loop over all quads (written out i1 i2 i3 i4) (Same numbering scheme)
66  loop over all triangles and quads (write out b1) (This is a boundary
67  condition for each
68  triangle and each
69  hex. You can put
70  anything you want
71  here)
72  loop over all tets (written out i1 i2 i3 i4) (Same)
73  loop over all pyramids (written out i1 i2 i3 i4 i5) (Same)
74  loop over all prisms (written out i1 i2 i3 i4 i5 i6) (Same)
75  loop over all hexs (written out i1 i2 i3 i4 i5 i6 i7 i8) (Same)
76 
77  */
78 
79  // Be sure the stream has been created successfully.
80  libmesh_assert (out_file.good());
81 
82  // Can't use a constant mesh reference since we have to
83  // sync the boundary info.
84  libmesh_here();
85  libMesh::err << "WARNING... Sure you want to do this?"
86  << std::endl;
87  MeshBase & the_mesh = const_cast<MeshBase &>
88  (MeshOutput<MeshBase>::mesh());
89 
90  if (the_mesh.mesh_dimension() < 3)
91  {
92  libMesh::err << "WARNING: DIVA only supports 3D meshes.\n\n"
93  << "Exiting without producing output.\n";
94  return;
95  }
96 
97 
98 
99  BoundaryMesh boundary_mesh
100  (the_mesh.comm(),
101  cast_int<unsigned char>(the_mesh.mesh_dimension()-1));
102  the_mesh.get_boundary_info().sync(boundary_mesh);
103 
104 
108  out_file << the_mesh.n_nodes()
109  << ' '
110  << (MeshTools::n_active_elem_of_type(boundary_mesh,TRI3) +
111  MeshTools::n_active_elem_of_type(boundary_mesh,TRI6)*4)
112  << ' '
113  << (MeshTools::n_active_elem_of_type(boundary_mesh, QUAD4) +
114  MeshTools::n_active_elem_of_type(boundary_mesh, QUAD8) +
115  MeshTools::n_active_elem_of_type(boundary_mesh, QUAD9)*4)
116  << ' '
117  << (MeshTools::n_active_elem_of_type(the_mesh, TET4) +
118  MeshTools::n_active_elem_of_type(the_mesh, TET10)*8)
119  << ' '
120  << MeshTools::n_active_elem_of_type(the_mesh, PYRAMID5)
121  << ' '
122  << (MeshTools::n_active_elem_of_type(the_mesh, PRISM6) +
123  MeshTools::n_active_elem_of_type(the_mesh, PRISM18)*8)
124  << ' '
125  << (MeshTools::n_active_elem_of_type(the_mesh, HEX8) +
126  MeshTools::n_active_elem_of_type(the_mesh, HEX20) +
127  MeshTools::n_active_elem_of_type(the_mesh, HEX27)*8)
128  << ' '
129  << '\n';
130 
131  boundary_mesh.clear();
132 
133 
137  for (unsigned int v=0; v<the_mesh.n_nodes(); v++)
138  the_mesh.point(v).write_unformatted(out_file);
139 
140 
144  {
148  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
149  if (the_mesh.elem(e)->active())
150  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
151  if (the_mesh.elem(e)->neighbor(s) == libmesh_nullptr)
152  {
153  const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
154 
155  if (side->type() == TRI3)
156  {
157  out_file << side->node(0)+1 << " "
158  << side->node(1)+1 << " "
159  << side->node(2)+1 << '\n';
160  }
161  else if (side->type() == TRI6)
162  {
163  out_file << side->node(0)+1 << " "
164  << side->node(3)+1 << " "
165  << side->node(5)+1 << '\n'
166 
167  << side->node(3)+1 << " "
168  << side->node(1)+1 << " "
169  << side->node(4)+1 << '\n'
170 
171  << side->node(5)+1 << " "
172  << side->node(4)+1 << " "
173  << side->node(2)+1 << '\n'
174 
175  << side->node(3)+1 << " "
176  << side->node(4)+1 << " "
177  << side->node(5)+1 << '\n';
178  }
179  }
180 
181 
185  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
186  if (the_mesh.elem(e)->active())
187  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
188  if (the_mesh.elem(e)->neighbor(s) == libmesh_nullptr)
189  {
190  const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
191 
192  if ((side->type() == QUAD4) ||
193  (side->type() == QUAD8) )
194  {
195  out_file << side->node(0)+1 << " "
196  << side->node(1)+1 << " "
197  << side->node(2)+1 << " "
198  << side->node(3)+1 << '\n';
199  }
200  else if (side->type() == QUAD9)
201  {
202  out_file << side->node(0)+1 << " "
203  << side->node(4)+1 << " "
204  << side->node(8)+1 << " "
205  << side->node(7)+1 << '\n'
206 
207  << side->node(4)+1 << " "
208  << side->node(1)+1 << " "
209  << side->node(5)+1 << " "
210  << side->node(8)+1 << '\n'
211 
212  << side->node(7)+1 << " "
213  << side->node(8)+1 << " "
214  << side->node(6)+1 << " "
215  << side->node(3)+1 << '\n'
216 
217  << side->node(8)+1 << " "
218  << side->node(5)+1 << " "
219  << side->node(2)+1 << " "
220  << side->node(6)+1 << '\n';
221  }
222  }
223  }
224 
225 
226 
230  {
234  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
235  if (the_mesh.elem(e)->active())
236  for (unsigned short s=0; s<the_mesh.elem(e)->n_sides(); s++)
237  if (the_mesh.elem(e)->neighbor(s) == libmesh_nullptr)
238  {
239  const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
240 
241  if ((side->type() == TRI3) ||
242  (side->type() == TRI6) )
243 
244  out_file << the_mesh.get_boundary_info().boundary_id(the_mesh.elem(e), s)
245  << '\n';
246  }
247 
248 
252  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
253  if (the_mesh.elem(e)->active())
254  for (unsigned short s=0; s<the_mesh.elem(e)->n_sides(); s++)
255  if (the_mesh.elem(e)->neighbor(s) == libmesh_nullptr)
256  {
257  const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
258 
259  if ((side->type() == QUAD4) ||
260  (side->type() == QUAD8) ||
261  (side->type() == QUAD9))
262 
263  out_file << the_mesh.get_boundary_info().boundary_id(the_mesh.elem(e), s);
264  }
265  }
266 
267 
268 
272  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
273  if (the_mesh.elem(e)->active())
274  {
275  if (the_mesh.elem(e)->type() == TET4)
276  {
277  out_file << the_mesh.elem(e)->node(0)+1 << " "
278  << the_mesh.elem(e)->node(1)+1 << " "
279  << the_mesh.elem(e)->node(2)+1 << " "
280  << the_mesh.elem(e)->node(3)+1 << '\n';
281  }
282  else if (the_mesh.elem(e)->type() == TET10)
283  {
284  out_file << the_mesh.elem(e)->node(0)+1 << " "
285  << the_mesh.elem(e)->node(4)+1 << " "
286  << the_mesh.elem(e)->node(6)+1 << " "
287  << the_mesh.elem(e)->node(7)+1 << '\n';
288 
289  out_file << the_mesh.elem(e)->node(4)+1 << " "
290  << the_mesh.elem(e)->node(1)+1 << " "
291  << the_mesh.elem(e)->node(5)+1 << " "
292  << the_mesh.elem(e)->node(8)+1 << '\n';
293 
294  out_file << the_mesh.elem(e)->node(6)+1 << " "
295  << the_mesh.elem(e)->node(5)+1 << " "
296  << the_mesh.elem(e)->node(2)+1 << " "
297  << the_mesh.elem(e)->node(9)+1 << '\n';
298 
299  out_file << the_mesh.elem(e)->node(7)+1 << " "
300  << the_mesh.elem(e)->node(8)+1 << " "
301  << the_mesh.elem(e)->node(9)+1 << " "
302  << the_mesh.elem(e)->node(3)+1 << '\n';
303 
304  out_file << the_mesh.elem(e)->node(4)+1 << " "
305  << the_mesh.elem(e)->node(8)+1 << " "
306  << the_mesh.elem(e)->node(6)+1 << " "
307  << the_mesh.elem(e)->node(7)+1 << '\n';
308 
309  out_file << the_mesh.elem(e)->node(4)+1 << " "
310  << the_mesh.elem(e)->node(5)+1 << " "
311  << the_mesh.elem(e)->node(6)+1 << " "
312  << the_mesh.elem(e)->node(8)+1 << '\n';
313 
314  out_file << the_mesh.elem(e)->node(6)+1 << " "
315  << the_mesh.elem(e)->node(5)+1 << " "
316  << the_mesh.elem(e)->node(9)+1 << " "
317  << the_mesh.elem(e)->node(8)+1 << '\n';
318 
319  out_file << the_mesh.elem(e)->node(6)+1 << " "
320  << the_mesh.elem(e)->node(8)+1 << " "
321  << the_mesh.elem(e)->node(9)+1 << " "
322  << the_mesh.elem(e)->node(7)+1 << '\n';
323  }
324  }
325 
326 
330  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
331  if (the_mesh.elem(e)->active())
332  if (the_mesh.elem(e)->type() == PYRAMID5)
333  {
334  out_file << the_mesh.elem(e)->node(0)+1 << " "
335  << the_mesh.elem(e)->node(1)+1 << " "
336  << the_mesh.elem(e)->node(2)+1 << " "
337  << the_mesh.elem(e)->node(3)+1 << " "
338  << the_mesh.elem(e)->node(4)+1 << '\n';
339  }
340 
341 
342 
346  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
347  if (the_mesh.elem(e)->active())
348  {
349  if (the_mesh.elem(e)->type() == PRISM6)
350  {
351  out_file << the_mesh.elem(e)->node(0)+1 << " "
352  << the_mesh.elem(e)->node(1)+1 << " "
353  << the_mesh.elem(e)->node(2)+1 << " "
354  << the_mesh.elem(e)->node(3)+1 << " "
355  << the_mesh.elem(e)->node(4)+1 << " "
356  << the_mesh.elem(e)->node(5)+1 << '\n';
357  }
358  else if (the_mesh.elem(e)->type() == PRISM18)
359  libmesh_error_msg("PRISM18 element type not supported.");
360  }
361 
362 
366  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
367  if (the_mesh.elem(e)->active())
368  if ((the_mesh.elem(e)->type() == HEX8) ||
369  (the_mesh.elem(e)->type() == HEX20) ||
370  (the_mesh.elem(e)->type() == HEX27) )
371  {
372  std::vector<dof_id_type> conn;
373  for (unsigned int se=0; se<the_mesh.elem(e)->n_sub_elem(); se++)
374  {
375  the_mesh.elem(e)->connectivity(se, TECPLOT, conn);
376 
377  out_file << conn[0] << ' '
378  << conn[1] << ' '
379  << conn[2] << ' '
380  << conn[3] << ' '
381  << conn[4] << ' '
382  << conn[5] << ' '
383  << conn[6] << ' '
384  << conn[7] << '\n';
385  }
386  }
387 }
388 
389 } // namespace libMesh
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:111
libMesh::Elem::active
bool active() const
Definition: elem.h:1721
libMesh::libmesh_assert
libmesh_assert(remote_elem)
boundary_mesh.h
side
unsigned short int side
Definition: xdr_io.C:50
libMesh::BoundaryInfo::sync
void sync(UnstructuredMesh &boundary_mesh, MeshData *boundary_mesh_data=libmesh_nullptr, MeshData *this_mesh_data=libmesh_nullptr)
Definition: boundary_info.C:129
libMesh::TypeVector::write_unformatted
void write_unformatted(std::ostream &out, const bool newline=true) const
Definition: type_vector.C:92
libmesh_nullptr
const class libmesh_nullptr_t libmesh_nullptr
libMesh::DivaIO::write
virtual void write(const std::string &) libmesh_override
Definition: diva_io.C:34
libMesh::MeshOutput< MeshBase >::mesh
const MeshBase & mesh() const
libMesh::ParallelMesh::clear
virtual void clear() libmesh_override
Definition: parallel_mesh.C:800
libMesh::UniquePtr
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
libMesh::MeshBase::elem
virtual const Elem * elem(const dof_id_type i) const =0
libMesh::DivaIO::write_stream
void write_stream(std::ostream &out)
Definition: diva_io.C:53
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:235
libMesh::BoundaryInfo::boundary_id
boundary_id_type boundary_id(const Elem *const elem, const unsigned short int side) const
Definition: boundary_info.C:985
libMesh::Elem::node
dof_id_type node(const unsigned int i) const
Definition: elem.h:1499
libMesh::Elem::build_side
virtual UniquePtr< Elem > build_side(const unsigned int i, bool proxy=true) const =0
libMesh::Elem::n_sides
virtual unsigned int n_sides() const =0
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:87
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:127
libMesh::MeshTools::n_active_elem_of_type
dof_id_type n_active_elem_of_type(const MeshBase &mesh, const ElemType type)
Definition: mesh_tools.C:551
libMesh::Elem::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const =0
mesh_tools.h
diva_io.h
boundary_info.h
libMesh::MeshBase::point
virtual const Point & point(const dof_id_type i) const =0
libMesh::Elem::n_sub_elem
virtual unsigned int n_sub_elem() const =0
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::Elem::neighbor
Elem * neighbor(const unsigned int i) const
Definition: elem.h:1580
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
generated by doxygen