gmsh_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 // This file was massively overhauled and extended by Martin Lüthi, mluthi@tnoo.net
19 
20 // C++ includes
21 #include <fstream>
22 #include <set>
23 #include <cstring> // std::memcpy, std::strncmp
24 
25 // Local includes
26 #include "libmesh/libmesh_config.h"
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/gmsh_io.h"
31 #include "libmesh/mesh_base.h"
32 
33 
34 // anonymous namespace to hold local data
35 namespace
36 {
37  using namespace libMesh;
38 
45  struct boundaryElementInfo {
46  std::set<dof_id_type> nodes;
47  unsigned int id;
48  };
49 
53  struct elementDefinition {
54  std::string label;
55  std::vector<unsigned int> nodes;
56  ElemType type;
57  unsigned int exptype;
58  unsigned int dim;
59  unsigned int nnodes;
60  };
61 
62 
63  // maps from a libMesh element type to the proper
64  // Gmsh elementDefinition. Placing the data structure
65  // here in this anonymous namespace gives us the
66  // benefits of a global variable without the nasty
67  // side-effects
68  std::map<ElemType, elementDefinition> eletypes_exp;
69  std::map<unsigned int, elementDefinition> eletypes_imp;
70 
71 
72 
73  // ------------------------------------------------------------
74  // helper function to initialize the eletypes map
75  void init_eletypes ()
76  {
77  if (eletypes_exp.empty() && eletypes_imp.empty())
78  {
79  // This should happen only once. The first time this method
80  // is called the eletypes data struture will be empty, and
81  // we will fill it. Any subsequent calls will find an initialized
82  // eletypes map and will do nothing.
83 
84  //==============================
85  // setup the element definitions
86  elementDefinition eledef;
87 
88  // use "swap trick" from Scott Meyer's "Effective STL" to initialize
89  // eledef.nodes vector
90 
91  // POINT (only Gmsh)
92  {
93  eledef.exptype = 15;
94  eledef.dim = 0;
95  eledef.nnodes = 1;
96  eledef.nodes.clear();
97 
98  // import only
99  eletypes_imp[15] = eledef;
100  }
101 
102  // EDGE2
103  {
104  eledef.type = EDGE2;
105  eledef.dim = 1;
106  eledef.nnodes = 2;
107  eledef.exptype = 1;
108  eledef.nodes.clear();
109 
110  eletypes_exp[EDGE2] = eledef;
111  eletypes_imp[1] = eledef;
112  }
113 
114  // EDGE3
115  {
116  eledef.type = EDGE3;
117  eledef.dim = 1;
118  eledef.nnodes = 3;
119  eledef.exptype = 8;
120  eledef.nodes.clear();
121 
122  eletypes_exp[EDGE3] = eledef;
123  eletypes_imp[8] = eledef;
124  }
125 
126  // TRI3
127  {
128  eledef.type = TRI3;
129  eledef.dim = 2;
130  eledef.nnodes = 3;
131  eledef.exptype = 2;
132  eledef.nodes.clear();
133 
134  eletypes_exp[TRI3] = eledef;
135  eletypes_imp[2] = eledef;
136  }
137 
138  // TRI6
139  {
140  eledef.type = TRI6;
141  eledef.dim = 2;
142  eledef.nnodes = 6;
143  eledef.exptype = 9;
144  eledef.nodes.clear();
145 
146  eletypes_exp[TRI6] = eledef;
147  eletypes_imp[9] = eledef;
148  }
149 
150  // QUAD4
151  {
152  eledef.type = QUAD4;
153  eledef.dim = 2;
154  eledef.nnodes = 4;
155  eledef.exptype = 3;
156  eledef.nodes.clear();
157 
158  eletypes_exp[QUAD4] = eledef;
159  eletypes_imp[3] = eledef;
160  }
161 
162  // QUAD8
163  // TODO: what should be done with this on writing?
164  {
165  eledef.type = QUAD8;
166  eledef.dim = 2;
167  eledef.nnodes = 8;
168  eledef.exptype = 100;
169  const unsigned int nodes[] = {1,2,3,4,5,6,7,8};
170  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
171 
172  eletypes_exp[QUAD8] = eledef;
173  eletypes_imp[10] = eledef;
174  }
175 
176  // QUAD9
177  {
178  eledef.type = QUAD9;
179  eledef.dim = 2;
180  eledef.nnodes = 9;
181  eledef.exptype = 10;
182  eledef.nodes.clear();
183 
184  eletypes_exp[QUAD9] = eledef;
185  eletypes_imp[10] = eledef;
186  }
187 
188  // HEX8
189  {
190  eledef.type = HEX8;
191  eledef.dim = 3;
192  eledef.nnodes = 8;
193  eledef.exptype = 5;
194  eledef.nodes.clear();
195 
196  eletypes_exp[HEX8] = eledef;
197  eletypes_imp[5] = eledef;
198  }
199 
200  // HEX20
201  // TODO: what should be done with this on writing?
202  {
203  eledef.type = HEX20;
204  eledef.dim = 3;
205  eledef.nnodes = 20;
206  eledef.exptype = 101;
207  const unsigned int nodes[] = {1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15,16};
208  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
209 
210  eletypes_exp[HEX20] = eledef;
211  eletypes_imp[12] = eledef;
212  }
213 
214  // HEX27
215  {
216  eledef.type = HEX27;
217  eledef.dim = 3;
218  eledef.nnodes = 27;
219  eledef.exptype = 12;
220  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
221  15,16,19,17,18,20,21,24,22,23,25,26};
222  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
223 
224  eletypes_exp[HEX27] = eledef;
225  eletypes_imp[12] = eledef;
226  }
227 
228  // TET4
229  {
230  eledef.type = TET4;
231  eledef.dim = 3;
232  eledef.nnodes = 4;
233  eledef.exptype = 4;
234  eledef.nodes.clear();
235 
236  eletypes_exp[TET4] = eledef;
237  eletypes_imp[4] = eledef;
238  }
239 
240  // TET10
241  {
242  eledef.type = TET10;
243  eledef.dim = 3;
244  eledef.nnodes = 10;
245  eledef.exptype = 11;
246  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
247  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
248  eletypes_exp[TET10] = eledef;
249  eletypes_imp[11] = eledef;
250  }
251 
252  // PRISM6
253  {
254  eledef.type = PRISM6;
255  eledef.dim = 3;
256  eledef.nnodes = 6;
257  eledef.exptype = 6;
258  eledef.nodes.clear();
259 
260  eletypes_exp[PRISM6] = eledef;
261  eletypes_imp[6] = eledef;
262  }
263 
264  // PRISM15
265  // TODO: what should be done with this on writing?
266  {
267  eledef.type = PRISM15;
268  eledef.dim = 3;
269  eledef.nnodes = 15;
270  eledef.exptype = 103;
271  eledef.nodes.clear();
272 
273  eletypes_exp[PRISM15] = eledef;
274  eletypes_imp[13] = eledef;
275  }
276 
277  // PRISM18
278  {
279  eledef.type = PRISM18;
280  eledef.dim = 3;
281  eledef.nnodes = 18;
282  eledef.exptype = 13;
283  const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,
284  12,14,13,15,17,16};
285  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
286 
287  eletypes_exp[PRISM18] = eledef;
288  eletypes_imp[13] = eledef;
289  }
290 
291  // PYRAMID5
292  {
293  eledef.type = PYRAMID5;
294  eledef.dim = 3;
295  eledef.nnodes = 5;
296  eledef.exptype = 7;
297  eledef.nodes.clear();
298 
299  eletypes_exp[PYRAMID5] = eledef;
300  eletypes_imp[7] = eledef;
301  }
302 
303  //==============================
304  }
305  }
306 
307 } // end anonymous namespace
308 
309 
310 namespace libMesh
311 {
312 
313 // ------------------------------------------------------------
314 // GmshIO members
315 void GmshIO::read (const std::string& name)
316 {
317  std::ifstream in (name.c_str());
318  this->read_mesh (in);
319 }
320 
321 
322 void GmshIO::read_mesh(std::istream& in)
323 {
324  // This is a serial-only process for now;
325  // the Mesh should be read on processor 0 and
326  // broadcast later
327  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
328 
329  libmesh_assert(in.good());
330 
331  // initialize the map with element types
332  init_eletypes();
333 
334  // clear any data in the mesh
336  mesh.clear();
337 
338  const unsigned int dim = mesh.mesh_dimension();
339 
340  // some variables
341  const int bufLen = 256;
342  char buf[bufLen+1];
343  int format=0, size=0;
344  Real version = 1.0;
345 
346  // map to hold the node numbers for translation
347  // note the the nodes can be non-consecutive
348  std::map<unsigned int, unsigned int> nodetrans;
349 
350  {
351  while (!in.eof()) {
352  in >> buf;
353 
354  if (!std::strncmp(buf,"$MeshFormat",11))
355  {
356  in >> version >> format >> size;
357  if ((version != 2.0) && (version != 2.1) && (version != 2.2)) {
358  // Some notes on gmsh mesh versions:
359  //
360  // Mesh version 2.0 goes back as far as I know. It's not explicitly
361  // mentioned here: http://www.geuz.org/gmsh/doc/VERSIONS.txt
362  //
363  // As of gmsh-2.4.0:
364  // bumped mesh version format to 2.1 (small change in the $PhysicalNames
365  // section, where the group dimension is now required);
366  // [Since we don't even parse the PhysicalNames section at the time
367  // of this writing, I don't think this change affects us.]
368  //
369  // Mesh version 2.2 tested by Manav Bhatia; no other
370  // libMesh code changes were required for support
371  libMesh::err << "Error: Wrong msh file version " << version << "\n";
372  libmesh_error();
373  }
374  if(format){
375  libMesh::err << "Error: Unknown data format for mesh\n";
376  libmesh_error();
377  }
378  }
379 
380  // read the node block
381  else if (!std::strncmp(buf,"$NOD",4) ||
382  !std::strncmp(buf,"$NOE",4) ||
383  !std::strncmp(buf,"$Nodes",6)
384  )
385  {
386  unsigned int numNodes = 0;
387  in >> numNodes;
388  mesh.reserve_nodes (numNodes);
389 
390  // read in the nodal coordinates and form points.
391  Real x, y, z;
392  unsigned int id;
393 
394  // add the nodal coordinates to the mesh
395  for (unsigned int i=0; i<numNodes; ++i)
396  {
397  in >> id >> x >> y >> z;
398  mesh.add_point (Point(x, y, z), i);
399  nodetrans[id] = i;
400  }
401  // read the $ENDNOD delimiter
402  in >> buf;
403  }
404 
415  else if (!std::strncmp(buf,"$ELM",4) ||
416  !std::strncmp(buf,"$Elements",9)
417  )
418  {
419  unsigned int numElem = 0;
420  std::vector< boundaryElementInfo > boundary_elem;
421 
422  // read how many elements are there, and reserve space in the mesh
423  in >> numElem;
424  mesh.reserve_elem (numElem);
425 
426  // read the elements
427  unsigned int elem_id_counter = 0;
428  for (unsigned int iel=0; iel<numElem; ++iel)
429  {
430  unsigned int id, type, physical, elementary,
431  /* partition = 1,*/ nnodes, ntags;
432  // note - partition was assigned but never used - BSK
433  if(version <= 1.0)
434  {
435  in >> id >> type >> physical >> elementary >> nnodes;
436  }
437  else
438  {
439  in >> id >> type >> ntags;
440  elementary = physical = /* partition = */ 1;
441  for(unsigned int j = 0; j < ntags; j++)
442  {
443  int tag;
444  in >> tag;
445  if(j == 0)
446  physical = tag;
447  else if(j == 1)
448  elementary = tag;
449  // else if(j == 2)
450  // partition = tag;
451  // ignore any other tags for now
452  }
453  }
454 
455  // consult the import element table which element to build
456  const elementDefinition& eletype = eletypes_imp[type];
457  nnodes = eletype.nnodes;
458 
459  // only elements that match the mesh dimension are added
460  // if the element dimension is one less than dim, the nodes and
461  // sides are added to the mesh.boundary_info
462  if (eletype.dim == dim)
463  {
464  // add the elements to the mesh
465  Elem* elem = Elem::build(eletype.type).release();
466  elem->set_id(elem_id_counter);
467  mesh.add_elem(elem);
468 
469  // different to iel, lower dimensional elems aren't added
470  elem_id_counter++;
471 
472  // check number of nodes. We cannot do that for version 2.0
473  if (version <= 1.0)
474  {
475  if (elem->n_nodes() != nnodes)
476  {
477  libMesh::err << "Number of nodes for element " << id
478  << " of type " << eletypes_imp[type].type
479  << " (Gmsh type " << type
480  << ") does not match Libmesh definition. "
481  << "I expected " << elem->n_nodes()
482  << " nodes, but got " << nnodes << "\n";
483  libmesh_error();
484  }
485  }
486 
487  // add node pointers to the elements
488  int nod = 0;
489  // if there is a node translation table, use it
490  if (eletype.nodes.size() > 0)
491  for (unsigned int i=0; i<nnodes; i++)
492  {
493  in >> nod;
494  elem->set_node(eletype.nodes[i]) = mesh.node_ptr(nodetrans[nod]);
495  }
496  else
497  {
498  for (unsigned int i=0; i<nnodes; i++)
499  {
500  in >> nod;
501  elem->set_node(i) = mesh.node_ptr(nodetrans[nod]);
502  }
503  }
504 
505  // Finally, set the subdomain ID to physical
506  elem->subdomain_id() = static_cast<subdomain_id_type>(physical);
507  } // if element.dim == dim
508  // if this is a boundary
509  else if (eletype.dim == dim-1)
510  {
515  boundaryElementInfo binfo;
516  std::set<dof_id_type>::iterator iter = binfo.nodes.begin();
517  int nod = 0;
518  for (unsigned int i=0; i<nnodes; i++)
519  {
520  in >> nod;
521  mesh.boundary_info->add_node(nodetrans[nod], physical);
522  binfo.nodes.insert(iter, nodetrans[nod]);
523  }
524  binfo.id = physical;
525  boundary_elem.push_back(binfo);
526  }
531  else
532  {
533  static bool seen_high_dim_element = false;
534  if (!seen_high_dim_element)
535  {
536  libMesh::err << "Warning: can't load an element of dimension "
537  << eletype.dim << " into a mesh of dimension "
538  << dim << std::endl;
539  seen_high_dim_element = true;
540  }
541  int nod = 0;
542  for (unsigned int i=0; i<nnodes; i++)
543  in >> nod;
544  }
545  }//element loop
546  // read the $ENDELM delimiter
547  in >> buf;
548 
554  if (boundary_elem.size() > 0)
555  {
556  // create a index of the boundary nodes to easily locate which
557  // element might have that boundary
558  std::map<dof_id_type, std::vector<unsigned int> > node_index;
559  for (unsigned int i=0; i<boundary_elem.size(); i++)
560  {
561  boundaryElementInfo binfo = boundary_elem[i];
562  std::set<dof_id_type>::iterator iter = binfo.nodes.begin();
563  for (;iter!= binfo.nodes.end(); iter++)
564  node_index[*iter].push_back(i);
565  }
566 
569 
570  // iterate over all elements and see which boundary element has
571  // the same set of nodes as on of the boundary elements previously read
572  for ( ; it != end; ++it)
573  {
574  const Elem* elem = *it;
575  for (unsigned int s=0; s<elem->n_sides(); s++)
576  if (elem->neighbor(s) == NULL)
577  {
578  AutoPtr<Elem> side (elem->build_side(s));
579  std::set<dof_id_type> side_nodes;
580  std::set<dof_id_type>::iterator iter = side_nodes.begin();
581 
582  // make a set with all nodes from this side
583  // this allows for easy comparison
584  for (unsigned int ns=0; ns<side->n_nodes(); ns++)
585  side_nodes.insert(iter, side->node(ns));
586 
587  // See whether one of the side node occurs in the list
588  // of tagged nodes. If we would loop over all side
589  // nodes, we would just get multiple hits, so taking
590  // node 0 is enough to do the job
591  dof_id_type sn = side->node(0);
592  if (node_index.count(sn) > 0)
593  {
594  // Loop over all tagged ("physical") "sides" which
595  // contain the node sn (typically just 1 to
596  // three). For each of these the set of nodes is
597  // compared to the current element's side nodes
598  for (std::size_t n=0; n<node_index[sn].size(); n++)
599  {
600  unsigned int bidx = node_index[sn][n];
601  if (boundary_elem[bidx].nodes == side_nodes)
602  mesh.boundary_info->add_side(elem, s, boundary_elem[bidx].id);
603  }
604  }
605  } // if elem->neighbor(s) == NULL
606  } // element loop
607  } // if boundary_elem.size() > 0
608  } // if $ELM
609 
610  } // while !in.eof()
611 
612  }
613 
614 }
615 
616 
617 
618 void GmshIO::write (const std::string& name)
619 {
621  {
622  // Open the output file stream
623  std::ofstream out_stream (name.c_str());
624 
625  // Make sure it opened correctly
626  if (!out_stream.good())
627  libmesh_file_error(name.c_str());
628 
629  this->write_mesh (out_stream);
630  }
631 }
632 
633 
634 void GmshIO::write_nodal_data (const std::string& fname,
635  const std::vector<Number>& soln,
636  const std::vector<std::string>& names)
637 {
638  START_LOG("write_nodal_data()", "GmshIO");
639 
640  //this->_binary = true;
642  this->write_post (fname, &soln, &names);
643 
644  STOP_LOG("write_nodal_data()", "GmshIO");
645 }
646 
647 
648 void GmshIO::write_mesh (std::ostream& out_stream)
649 {
650  // Be sure that the stream is valid.
651  libmesh_assert (out_stream.good());
652 
653  // initialize the map with element types
654  init_eletypes();
655 
656  // Get a const reference to the mesh
658 
659  // Note: we are using version 2.0 of the gmsh output format.
660 
661  {
662  // Write the file header.
663  out_stream << "$MeshFormat\n";
664  out_stream << "2.0 0 " << sizeof(Real) << '\n';
665  out_stream << "$EndMeshFormat\n";
666  }
667 
668  {
669  // write the nodes in (n x y z) format
670  out_stream << "$Nodes\n";
671  out_stream << mesh.n_nodes() << '\n';
672 
673  for (unsigned int v=0; v<mesh.n_nodes(); v++)
674  out_stream << mesh.node(v).id()+1 << " "
675  << mesh.node(v)(0) << " "
676  << mesh.node(v)(1) << " "
677  << mesh.node(v)(2) << '\n';
678  out_stream << "$EndNodes\n";
679  }
680 
681  {
682  // write the connectivity
683  out_stream << "$Elements\n";
684  out_stream << mesh.n_active_elem() << '\n';
685 
688 
689  // loop over the elements
690  for ( ; it != end; ++it)
691  {
692  const Elem* elem = *it;
693 
694  // Make sure we have a valid entry for
695  // the current element type.
696  libmesh_assert (eletypes_exp.count(elem->type()));
697 
698  // consult the export element table
699  const elementDefinition& eletype = eletypes_exp[elem->type()];
700 
701  // The element mapper better not require any more nodes
702  // than are present in the current element!
703  libmesh_assert_less_equal (eletype.nodes.size(), elem->n_nodes());
704 
705  // elements ids are 1 based in Gmsh
706  out_stream << elem->id()+1 << " ";
707 
708  // element type
709  out_stream << eletype.exptype;
710 
711  // write the number of tags and
712  // tag1 (physical entity), tag2 (geometric entity), and tag3 (partition entity)
713  out_stream << " 3 "
714  << static_cast<unsigned int>(elem->subdomain_id())
715  << " 1 "
716  << (elem->processor_id()+1) << " ";
717 
718  // if there is a node translation table, use it
719  if (eletype.nodes.size() > 0)
720  for (unsigned int i=0; i < elem->n_nodes(); i++)
721  out_stream << elem->node(eletype.nodes[i])+1 << " "; // gmsh is 1-based
722  // otherwise keep the same node order
723  else
724  for (unsigned int i=0; i < elem->n_nodes(); i++)
725  out_stream << elem->node(i)+1 << " "; // gmsh is 1-based
726  out_stream << "\n";
727  } // element loop
728  out_stream << "$EndElements\n";
729  }
730 }
731 
732 
733 void GmshIO::write_post (const std::string& fname,
734  const std::vector<Number>* v,
735  const std::vector<std::string>* solution_names)
736 {
737 
738  // Should only do this on processor 0!
739  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
740 
741  // Create an output stream
742  std::ofstream out_stream(fname.c_str());
743 
744  // Make sure it opened correctly
745  if (!out_stream.good())
746  libmesh_file_error(fname.c_str());
747 
748  // initialize the map with element types
749  init_eletypes();
750 
751  // create a character buffer
752  char buf[80];
753 
754  // Get a constant reference to the mesh.
756 
757  // write the data
758  if ((solution_names != NULL) && (v != NULL))
759  {
760  const unsigned int n_vars =
761  libmesh_cast_int<unsigned int>(solution_names->size());
762 
763  if (!(v->size() == mesh.n_nodes()*n_vars))
764  libMesh::err << "ERROR: v->size()=" << v->size()
765  << ", mesh.n_nodes()=" << mesh.n_nodes()
766  << ", n_vars=" << n_vars
767  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
768  << "\n";
769 
770  libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
771 
772  // write the header
773  out_stream << "$PostFormat\n";
774  if (this->binary())
775  out_stream << "1.2 1 " << sizeof(double) << "\n";
776  else
777  out_stream << "1.2 0 " << sizeof(double) << "\n";
778  out_stream << "$EndPostFormat\n";
779 
780  // Loop over the elements to see how much of each type there are
781  unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
782  n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
783  unsigned int n_scalar=0, n_vector=0, n_tensor=0;
784  unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;
785 
786  {
789 
790 
791  for ( ; it != end; ++it)
792  {
793  const ElemType elemtype = (*it)->type();
794 
795  switch (elemtype)
796  {
797  case EDGE2:
798  case EDGE3:
799  case EDGE4:
800  {
801  n_lines += 1;
802  break;
803  }
804  case TRI3:
805  case TRI6:
806  {
807  n_triangles += 1;
808  break;
809  }
810  case QUAD4:
811  case QUAD8:
812  case QUAD9:
813  {
814  n_quadrangles += 1;
815  break;
816  }
817  case TET4:
818  case TET10:
819  {
820  n_tetrahedra += 1;
821  break;
822  }
823  case HEX8:
824  case HEX20:
825  case HEX27:
826  {
827  n_hexahedra += 1;
828  break;
829  }
830  case PRISM6:
831  case PRISM15:
832  case PRISM18:
833  {
834  n_prisms += 1;
835  break;
836  }
837  case PYRAMID5:
838  {
839  n_pyramids += 1;
840  break;
841  }
842  default:
843  {
844  libMesh::err << "ERROR: Not existant element type "
845  << (*it)->type() << std::endl;
846  libmesh_error();
847  }
848  }
849  }
850  }
851 
852  // create a view for each variable
853  for (unsigned int ivar=0; ivar < n_vars; ivar++)
854  {
855  std::string varname = (*solution_names)[ivar];
856 
857  // at the moment, we just write out scalar quantities
858  // later this should be made configurable through
859  // options to the writer class
860  n_scalar = 1;
861 
862  // write the variable as a view, and the number of time steps
863  out_stream << "$View\n" << varname << " " << 1 << "\n";
864 
865  // write how many of each geometry type are written
866  out_stream << n_points * n_scalar << " "
867  << n_points * n_vector << " "
868  << n_points * n_tensor << " "
869  << n_lines * n_scalar << " "
870  << n_lines * n_vector << " "
871  << n_lines * n_tensor << " "
872  << n_triangles * n_scalar << " "
873  << n_triangles * n_vector << " "
874  << n_triangles * n_tensor << " "
875  << n_quadrangles * n_scalar << " "
876  << n_quadrangles * n_vector << " "
877  << n_quadrangles * n_tensor << " "
878  << n_tetrahedra * n_scalar << " "
879  << n_tetrahedra * n_vector << " "
880  << n_tetrahedra * n_tensor << " "
881  << n_hexahedra * n_scalar << " "
882  << n_hexahedra * n_vector << " "
883  << n_hexahedra * n_tensor << " "
884  << n_prisms * n_scalar << " "
885  << n_prisms * n_vector << " "
886  << n_prisms * n_tensor << " "
887  << n_pyramids * n_scalar << " "
888  << n_pyramids * n_vector << " "
889  << n_pyramids * n_tensor << " "
890  << nb_text2d << " "
891  << nb_text2d_chars << " "
892  << nb_text3d << " "
893  << nb_text3d_chars << "\n";
894 
895  // if binary, write a marker to identify the endianness of the file
896  if (this->binary())
897  {
898  const int one = 1;
899  std::memcpy(buf, &one, sizeof(int));
900  out_stream.write(buf, sizeof(int));
901  }
902 
903  // the time steps (there is just 1 at the moment)
904  if (this->binary())
905  {
906  double one = 1;
907  std::memcpy(buf, &one, sizeof(double));
908  out_stream.write(buf, sizeof(double));
909  }
910  else
911  out_stream << "1\n";
912 
913  // Loop over the elements and write out the data
916 
917  for ( ; it != end; ++it)
918  {
919  const Elem* elem = *it;
920 
921  // this is quite crappy, but I did not invent that file format!
922  for (unsigned int d=0; d<3; d++) // loop over the dimensions
923  {
924  for (unsigned int n=0; n < elem->n_vertices(); n++) // loop over vertices
925  {
926  const Point vertex = elem->point(n);
927  if (this->binary())
928  {
929  double tmp = vertex(d);
930  std::memcpy(buf, &tmp, sizeof(double));
931  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
932  }
933  else
934  out_stream << vertex(d) << " ";
935  }
936  if (!this->binary())
937  out_stream << "\n";
938  }
939 
940  // now finally write out the data
941  for (unsigned int i=0; i < elem->n_vertices(); i++) // loop over vertices
942  if (this->binary())
943  {
944 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
945  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
946  << "complex numbers. Will only write the real part of "
947  << "variable " << varname << std::endl;
948 #endif
949  double tmp = libmesh_real((*v)[elem->node(i)*n_vars + ivar]);
950  std::memcpy(buf, &tmp, sizeof(double));
951  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
952  }
953  else
954  {
955 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
956  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
957  << "complex numbers. Will only write the real part of "
958  << "variable " << varname << std::endl;
959 #endif
960  out_stream << libmesh_real((*v)[elem->node(i)*n_vars + ivar]) << "\n";
961  }
962  }
963  if (this->binary())
964  out_stream << "\n";
965  out_stream << "$EndView\n";
966 
967  } // end variable loop (writing the views)
968  }
969 
970 }
971 
972 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo