diva_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 // 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&>
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 (the_mesh.comm(),
100  the_mesh.mesh_dimension()-1);
101  the_mesh.boundary_info->sync(boundary_mesh);
102 
103 
107  out_file << the_mesh.n_nodes()
108  << ' '
109  << (MeshTools::n_active_elem_of_type(boundary_mesh,TRI3) +
110  MeshTools::n_active_elem_of_type(boundary_mesh,TRI6)*4)
111  << ' '
112  << (MeshTools::n_active_elem_of_type(boundary_mesh, QUAD4) +
113  MeshTools::n_active_elem_of_type(boundary_mesh, QUAD8) +
114  MeshTools::n_active_elem_of_type(boundary_mesh, QUAD9)*4)
115  << ' '
116  << (MeshTools::n_active_elem_of_type(the_mesh, TET4) +
118  << ' '
120  << ' '
123  << ' '
124  << (MeshTools::n_active_elem_of_type(the_mesh, HEX8) +
127  << ' '
128  << '\n';
129 
130  boundary_mesh.clear();
131 
132 
136  for (unsigned int v=0; v<the_mesh.n_nodes(); v++)
137  the_mesh.point(v).write_unformatted(out_file);
138 
139 
143  {
147  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
148  if (the_mesh.elem(e)->active())
149  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
150  if (the_mesh.elem(e)->neighbor(s) == NULL)
151  {
152  const AutoPtr<Elem> side(the_mesh.elem(e)->build_side(s));
153 
154  if (side->type() == TRI3)
155  {
156  out_file << side->node(0)+1 << " "
157  << side->node(1)+1 << " "
158  << side->node(2)+1 << '\n';
159  }
160  else if (side->type() == TRI6)
161  {
162  out_file << side->node(0)+1 << " "
163  << side->node(3)+1 << " "
164  << side->node(5)+1 << '\n'
165 
166  << side->node(3)+1 << " "
167  << side->node(1)+1 << " "
168  << side->node(4)+1 << '\n'
169 
170  << side->node(5)+1 << " "
171  << side->node(4)+1 << " "
172  << side->node(2)+1 << '\n'
173 
174  << side->node(3)+1 << " "
175  << side->node(4)+1 << " "
176  << side->node(5)+1 << '\n';
177  }
178  }
179 
180 
184  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
185  if (the_mesh.elem(e)->active())
186  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
187  if (the_mesh.elem(e)->neighbor(s) == NULL)
188  {
189  const AutoPtr<Elem> side(the_mesh.elem(e)->build_side(s));
190 
191  if ((side->type() == QUAD4) ||
192  (side->type() == QUAD8) )
193  {
194  out_file << side->node(0)+1 << " "
195  << side->node(1)+1 << " "
196  << side->node(2)+1 << " "
197  << side->node(3)+1 << '\n';
198  }
199  else if (side->type() == QUAD9)
200  {
201  out_file << side->node(0)+1 << " "
202  << side->node(4)+1 << " "
203  << side->node(8)+1 << " "
204  << side->node(7)+1 << '\n'
205 
206  << side->node(4)+1 << " "
207  << side->node(1)+1 << " "
208  << side->node(5)+1 << " "
209  << side->node(8)+1 << '\n'
210 
211  << side->node(7)+1 << " "
212  << side->node(8)+1 << " "
213  << side->node(6)+1 << " "
214  << side->node(3)+1 << '\n'
215 
216  << side->node(8)+1 << " "
217  << side->node(5)+1 << " "
218  << side->node(2)+1 << " "
219  << side->node(6)+1 << '\n';
220  }
221  }
222  }
223 
224 
225 
229  {
233  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
234  if (the_mesh.elem(e)->active())
235  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
236  if (the_mesh.elem(e)->neighbor(s) == NULL)
237  {
238  const AutoPtr<Elem> side(the_mesh.elem(e)->build_side(s));
239 
240  if ((side->type() == TRI3) ||
241  (side->type() == TRI6) )
242 
243  out_file << the_mesh.boundary_info->boundary_id(the_mesh.elem(e), s)
244  << '\n';
245  }
246 
247 
251  for(unsigned int e=0; e<the_mesh.n_elem(); e++)
252  if (the_mesh.elem(e)->active())
253  for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
254  if (the_mesh.elem(e)->neighbor(s) == NULL)
255  {
256  const AutoPtr<Elem> side(the_mesh.elem(e)->build_side(s));
257 
258  if ((side->type() == QUAD4) ||
259  (side->type() == QUAD8) ||
260  (side->type() == QUAD9))
261 
262  out_file << the_mesh.boundary_info->boundary_id(the_mesh.elem(e), s);
263  }
264  }
265 
266 
267 
271  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
272  if (the_mesh.elem(e)->active())
273  {
274  if (the_mesh.elem(e)->type() == TET4)
275  {
276  out_file << the_mesh.elem(e)->node(0)+1 << " "
277  << the_mesh.elem(e)->node(1)+1 << " "
278  << the_mesh.elem(e)->node(2)+1 << " "
279  << the_mesh.elem(e)->node(3)+1 << '\n';
280  }
281  else if (the_mesh.elem(e)->type() == TET10)
282  {
283  out_file << the_mesh.elem(e)->node(0)+1 << " "
284  << the_mesh.elem(e)->node(4)+1 << " "
285  << the_mesh.elem(e)->node(6)+1 << " "
286  << the_mesh.elem(e)->node(7)+1 << '\n';
287 
288  out_file << the_mesh.elem(e)->node(4)+1 << " "
289  << the_mesh.elem(e)->node(1)+1 << " "
290  << the_mesh.elem(e)->node(5)+1 << " "
291  << the_mesh.elem(e)->node(8)+1 << '\n';
292 
293  out_file << the_mesh.elem(e)->node(6)+1 << " "
294  << the_mesh.elem(e)->node(5)+1 << " "
295  << the_mesh.elem(e)->node(2)+1 << " "
296  << the_mesh.elem(e)->node(9)+1 << '\n';
297 
298  out_file << the_mesh.elem(e)->node(7)+1 << " "
299  << the_mesh.elem(e)->node(8)+1 << " "
300  << the_mesh.elem(e)->node(9)+1 << " "
301  << the_mesh.elem(e)->node(3)+1 << '\n';
302 
303  out_file << the_mesh.elem(e)->node(4)+1 << " "
304  << the_mesh.elem(e)->node(8)+1 << " "
305  << the_mesh.elem(e)->node(6)+1 << " "
306  << the_mesh.elem(e)->node(7)+1 << '\n';
307 
308  out_file << the_mesh.elem(e)->node(4)+1 << " "
309  << the_mesh.elem(e)->node(5)+1 << " "
310  << the_mesh.elem(e)->node(6)+1 << " "
311  << the_mesh.elem(e)->node(8)+1 << '\n';
312 
313  out_file << the_mesh.elem(e)->node(6)+1 << " "
314  << the_mesh.elem(e)->node(5)+1 << " "
315  << the_mesh.elem(e)->node(9)+1 << " "
316  << the_mesh.elem(e)->node(8)+1 << '\n';
317 
318  out_file << the_mesh.elem(e)->node(6)+1 << " "
319  << the_mesh.elem(e)->node(8)+1 << " "
320  << the_mesh.elem(e)->node(9)+1 << " "
321  << the_mesh.elem(e)->node(7)+1 << '\n';
322  }
323  }
324 
325 
329  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
330  if (the_mesh.elem(e)->active())
331  if (the_mesh.elem(e)->type() == PYRAMID5)
332  {
333  out_file << the_mesh.elem(e)->node(0)+1 << " "
334  << the_mesh.elem(e)->node(1)+1 << " "
335  << the_mesh.elem(e)->node(2)+1 << " "
336  << the_mesh.elem(e)->node(3)+1 << " "
337  << the_mesh.elem(e)->node(4)+1 << '\n';
338  }
339 
340 
341 
345  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
346  if (the_mesh.elem(e)->active())
347  {
348  if (the_mesh.elem(e)->type() == PRISM6)
349  {
350  out_file << the_mesh.elem(e)->node(0)+1 << " "
351  << the_mesh.elem(e)->node(1)+1 << " "
352  << the_mesh.elem(e)->node(2)+1 << " "
353  << the_mesh.elem(e)->node(3)+1 << " "
354  << the_mesh.elem(e)->node(4)+1 << " "
355  << the_mesh.elem(e)->node(5)+1 << '\n';
356  }
357  else if (the_mesh.elem(e)->type() == PRISM18)
358  {
359  libmesh_error();
360  }
361  }
362 
363 
367  for (unsigned int e=0; e<the_mesh.n_elem(); e++)
368  if (the_mesh.elem(e)->active())
369  if ((the_mesh.elem(e)->type() == HEX8) ||
370  (the_mesh.elem(e)->type() == HEX20) ||
371  (the_mesh.elem(e)->type() == HEX27) )
372  {
373  std::vector<dof_id_type> conn;
374  for (unsigned int se=0; se<the_mesh.elem(e)->n_sub_elem(); se++)
375  {
376  the_mesh.elem(e)->connectivity(se, TECPLOT, conn);
377 
378  out_file << conn[0] << ' '
379  << conn[1] << ' '
380  << conn[2] << ' '
381  << conn[3] << ' '
382  << conn[4] << ' '
383  << conn[5] << ' '
384  << conn[6] << ' '
385  << conn[7] << '\n';
386  }
387  }
388 }
389 
390 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo