libMesh
gmsh_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
19 #include <fstream>
20 #include <set>
21 #include <cstring> // std::memcpy
22 #include <numeric>
23 #include <unordered_map>
24 
25 // Local includes
26 #include "libmesh/libmesh_config.h"
27 #include "libmesh/libmesh_logging.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 namespace libMesh
34 {
35 
36 // Initialize the static data member
38 
39 
40 
41 // Definition of the static function which constructs the ElementMaps object.
43 {
44  // Object to be filled up
45  ElementMaps em;
46 
47  // POINT (import only)
48  em.in.insert(std::make_pair(15, ElementDefinition(NODEELEM, 15, 0, 1)));
49 
50  // Add elements with trivial node mappings
51  em.add_def(ElementDefinition(EDGE2, 1, 1, 2));
52  em.add_def(ElementDefinition(EDGE3, 8, 1, 3));
53  em.add_def(ElementDefinition(TRI3, 2, 2, 3));
54  em.add_def(ElementDefinition(TRI6, 9, 2, 6));
55  em.add_def(ElementDefinition(QUAD4, 3, 2, 4));
56  em.add_def(ElementDefinition(QUAD8, 16, 2, 8));
57  em.add_def(ElementDefinition(QUAD9, 10, 2, 9));
58  em.add_def(ElementDefinition(HEX8, 5, 3, 8));
59  em.add_def(ElementDefinition(TET4, 4, 3, 4));
60  em.add_def(ElementDefinition(PRISM6, 6, 3, 6));
61  em.add_def(ElementDefinition(PYRAMID5, 7, 3, 5));
62 
63  // Add elements with non-trivial node mappings
64 
65  // HEX20
66  {
67  ElementDefinition eledef(HEX20, 17, 3, 20);
68  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,15,16,19,17,18};
69  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
70  em.add_def(eledef);
71  }
72 
73  // HEX27
74  {
75  ElementDefinition eledef(HEX27, 12, 3, 27);
76  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
77  15,16,19,17,18,20,21,24,22,23,25,26};
78  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
79  em.add_def(eledef);
80  }
81 
82  // TET10
83  {
84  ElementDefinition eledef(TET10, 11, 3, 10);
85  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
86  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
87  em.add_def(eledef);
88  }
89 
90  // PRISM15
91  {
92  ElementDefinition eledef(PRISM15, 18, 3, 15);
93  const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13};
94  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
95  em.add_def(eledef);
96  }
97 
98  // PRISM18
99  {
100  ElementDefinition eledef(PRISM18, 13, 3, 18);
101  const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13,15,17,16};
102  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
103  em.add_def(eledef);
104  }
105 
106  return em;
107 }
108 
109 
110 
112  MeshOutput<MeshBase>(mesh),
113  _binary(false),
115 {
116 }
117 
118 
119 
121  MeshInput<MeshBase> (mesh),
122  MeshOutput<MeshBase> (mesh),
123  _binary (false),
125 {
126 }
127 
128 
129 
130 bool & GmshIO::binary ()
131 {
132  return _binary;
133 }
134 
135 
136 
138 {
140 }
141 
142 
143 
144 void GmshIO::read (const std::string & name)
145 {
146  std::ifstream in (name.c_str());
147  this->read_mesh (in);
148 }
149 
150 
151 
152 void GmshIO::read_mesh(std::istream & in)
153 {
154  // This is a serial-only process for now;
155  // the Mesh should be read on processor 0 and
156  // broadcast later
157  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
158 
159  libmesh_assert(in.good());
160 
161  // clear any data in the mesh
163  mesh.clear();
164 
165  // some variables
166  int format=0, size=0;
167  Real version = 1.0;
168 
169  // Keep track of lower-dimensional blocks which are not BCs, but
170  // actually blocks of lower-dimensional elements.
171  std::set<subdomain_id_type> lower_dimensional_blocks;
172 
173  // Mapping from physical id -> (physical dim, physical name) pairs.
174  // These can refer to either "sidesets" or "subdomains"; we need to
175  // wait until the Mesh has been read to know which is which. Note
176  // that we are using 'int' as the key here rather than
177  // subdomain_id_type or boundary_id_type, since at this point, it
178  // could be either.
179  typedef std::pair<unsigned, std::string> GmshPhysical;
180  std::map<int, GmshPhysical> gmsh_physicals;
181 
182  // map to hold the node numbers for translation
183  // note the the nodes can be non-consecutive
184  std::map<unsigned int, unsigned int> nodetrans;
185 
186  // For reading the file line by line
187  std::string s;
188 
189  while (true)
190  {
191  // Try to read something. This may set EOF!
192  std::getline(in, s);
193 
194  if (in)
195  {
196  // Process s...
197 
198  if (s.find("$MeshFormat") == static_cast<std::string::size_type>(0))
199  {
200  in >> version >> format >> size;
201  if ((version != 2.0) && (version != 2.1) && (version != 2.2))
202  {
203  // Some notes on gmsh mesh versions:
204  //
205  // Mesh version 2.0 goes back as far as I know. It's not explicitly
206  // mentioned here: http://www.geuz.org/gmsh/doc/VERSIONS.txt
207  //
208  // As of gmsh-2.4.0:
209  // bumped mesh version format to 2.1 (small change in the $PhysicalNames
210  // section, where the group dimension is now required);
211  // [Since we don't even parse the PhysicalNames section at the time
212  // of this writing, I don't think this change affects us.]
213  //
214  // Mesh version 2.2 tested by Manav Bhatia; no other
215  // libMesh code changes were required for support
216  libmesh_error_msg("Error: Unknown msh file version " << version);
217  }
218 
219  if (format)
220  libmesh_error_msg("Error: Unknown data format for mesh in Gmsh reader.");
221  }
222 
223  // Read and process the "PhysicalNames" section.
224  else if (s.find("$PhysicalNames") == static_cast<std::string::size_type>(0))
225  {
226  // The lines in the PhysicalNames section should look like the following:
227  // 2 1 "frac" lower_dimensional_block
228  // 2 3 "top"
229  // 2 4 "bottom"
230  // 3 2 "volume"
231 
232  // Read in the number of physical groups to expect in the file.
233  unsigned int num_physical_groups = 0;
234  in >> num_physical_groups;
235 
236  // Read rest of line including newline character.
237  std::getline(in, s);
238 
239  for (unsigned int i=0; i<num_physical_groups; ++i)
240  {
241  // Read an entire line of the PhysicalNames section.
242  std::getline(in, s);
243 
244  // Use an istringstream to extract the physical
245  // dimension, physical id, and physical name from
246  // this line.
247  std::istringstream s_stream(s);
248  unsigned phys_dim;
249  int phys_id;
250  std::string phys_name;
251  s_stream >> phys_dim >> phys_id >> phys_name;
252 
253  // Not sure if this is true for all Gmsh files, but
254  // my test file has quotes around the phys_name
255  // string. So let's erase any quotes now...
256  phys_name.erase(std::remove(phys_name.begin(), phys_name.end(), '"'), phys_name.end());
257 
258  // Record this ID for later assignment of subdomain/sideset names.
259  gmsh_physicals[phys_id] = std::make_pair(phys_dim, phys_name);
260 
261  // If 's' also contains the libmesh-specific string
262  // "lower_dimensional_block", add this block ID to
263  // the list of blocks which are not boundary
264  // conditions.
265  if (s.find("lower_dimensional_block") != std::string::npos)
266  {
267  lower_dimensional_blocks.insert(cast_int<subdomain_id_type>(phys_id));
268 
269  // The user has explicitly told us that this
270  // block is a subdomain, so set that association
271  // in the Mesh.
272  mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
273  }
274  }
275  }
276 
277  // read the node block
278  else if (s.find("$NOD") == static_cast<std::string::size_type>(0) ||
279  s.find("$NOE") == static_cast<std::string::size_type>(0) ||
280  s.find("$Nodes") == static_cast<std::string::size_type>(0))
281  {
282  unsigned int num_nodes = 0;
283  in >> num_nodes;
284  mesh.reserve_nodes (num_nodes);
285 
286  // read in the nodal coordinates and form points.
287  Real x, y, z;
288  unsigned int id;
289 
290  // add the nodal coordinates to the mesh
291  for (unsigned int i=0; i<num_nodes; ++i)
292  {
293  in >> id >> x >> y >> z;
294  mesh.add_point (Point(x, y, z), i);
295  nodetrans[id] = i;
296  }
297 
298  // read the $ENDNOD delimiter
299  std::getline(in, s);
300  }
301 
302 
303  // Read the element block
304  else if (s.find("$ELM") == static_cast<std::string::size_type>(0) ||
305  s.find("$Elements") == static_cast<std::string::size_type>(0))
306  {
307  // For reading the number of elements and the node ids from the stream
308  unsigned int
309  num_elem = 0,
310  node_id = 0;
311 
312  // read how many elements are there, and reserve space in the mesh
313  in >> num_elem;
314  mesh.reserve_elem (num_elem);
315 
316  // As of version 2.2, the format for each element line is:
317  // elm-number elm-type number-of-tags < tag > ... node-number-list
318  // From the Gmsh docs:
319  // * the first tag is the number of the
320  // physical entity to which the element belongs
321  // * the second is the number of the elementary geometrical
322  // entity to which the element belongs
323  // * the third is the number of mesh partitions to which the element
324  // belongs
325  // * The rest of the tags are the partition ids (negative
326  // partition ids indicate ghost cells). A zero tag is
327  // equivalent to no tag. Gmsh and most codes using the
328  // MSH 2 format require at least the first two tags
329  // (physical and elementary tags).
330 
331  // Keep track of all the element dimensions seen
332  std::vector<unsigned> elem_dimensions_seen(3);
333 
334  // read the elements
335  for (unsigned int iel=0; iel<num_elem; ++iel)
336  {
337  unsigned int
338  id, type,
339  physical=1, elementary=1,
340  nnodes=0, ntags;
341 
342  // Note: tag has to be an int because it could be negative,
343  // see above.
344  int tag;
345 
346  if (version <= 1.0)
347  in >> id >> type >> physical >> elementary >> nnodes;
348 
349  else
350  {
351  in >> id >> type >> ntags;
352 
353  if (ntags > 2)
354  libmesh_do_once(libMesh::err << "Warning, ntags=" << ntags << ", but we currently only support reading 2 flags." << std::endl;);
355 
356  for (unsigned int j = 0; j < ntags; j++)
357  {
358  in >> tag;
359  if (j == 0)
360  physical = tag;
361  else if (j == 1)
362  elementary = tag;
363  }
364  }
365 
366  // Consult the import element table to determine which element to build
367  std::map<unsigned int, GmshIO::ElementDefinition>::iterator eletypes_it = _element_maps.in.find(type);
368 
369  // Make sure we actually found something
370  if (eletypes_it == _element_maps.in.end())
371  libmesh_error_msg("Element type " << type << " not found!");
372 
373  // Get a reference to the ElementDefinition
374  const GmshIO::ElementDefinition & eletype = eletypes_it->second;
375 
376  // If we read nnodes, make sure it matches the number in eletype.nnodes
377  if (nnodes != 0 && nnodes != eletype.nnodes)
378  libmesh_error_msg("nnodes = " << nnodes << " and eletype.nnodes = " << eletype.nnodes << " do not match.");
379 
380  // Assign the value from the eletype object.
381  nnodes = eletype.nnodes;
382 
383  // Don't add 0-dimensional "point" elements to the
384  // Mesh. They should *always* be treated as boundary
385  // "nodeset" data.
386  if (eletype.dim > 0)
387  {
388  // Record this element dimension as being "seen".
389  // We will treat all elements with dimension <
390  // max(dimension) as specifying boundary conditions,
391  // but we won't know what max_elem_dimension_seen is
392  // until we read the entire file.
393  elem_dimensions_seen[eletype.dim-1] = 1;
394 
395  // Add the element to the mesh
396  {
397  Elem * elem = Elem::build(eletype.type).release();
398  elem->set_id(iel);
399  elem = mesh.add_elem(elem);
400 
401  // Make sure that the libmesh element we added has nnodes nodes.
402  if (elem->n_nodes() != nnodes)
403  libmesh_error_msg("Number of nodes for element " \
404  << id \
405  << " of type " << eletype.type \
406  << " (Gmsh type " << type \
407  << ") does not match Libmesh definition. " \
408  << "I expected " << elem->n_nodes() \
409  << " nodes, but got " << nnodes);
410 
411  // Add node pointers to the elements.
412  // If there is a node translation table, use it.
413  if (eletype.nodes.size() > 0)
414  for (unsigned int i=0; i<nnodes; i++)
415  {
416  in >> node_id;
417  elem->set_node(eletype.nodes[i]) = mesh.node_ptr(nodetrans[node_id]);
418  }
419  else
420  {
421  for (unsigned int i=0; i<nnodes; i++)
422  {
423  in >> node_id;
424  elem->set_node(i) = mesh.node_ptr(nodetrans[node_id]);
425  }
426  }
427 
428  // Finally, set the subdomain ID to physical. If this is a lower-dimension element, this ID will
429  // eventually go into the Mesh's BoundaryInfo object.
430  elem->subdomain_id() = static_cast<subdomain_id_type>(physical);
431  }
432  }
433 
434  // Handle 0-dimensional elements (points) by adding
435  // them to the BoundaryInfo object with
436  // boundary_id == physical.
437  else
438  {
439  // This seems like it should always be the same
440  // number as the 'id' we already read in on this
441  // line. At least it was in the example gmsh
442  // file I had...
443  in >> node_id;
444  mesh.get_boundary_info().add_node
445  (nodetrans[node_id],
446  static_cast<boundary_id_type>(physical));
447  }
448  } // element loop
449 
450  // read the $ENDELM delimiter
451  std::getline(in, s);
452 
453  // Record the max and min element dimension seen while reading the file.
454  unsigned char
455  max_elem_dimension_seen=1,
456  min_elem_dimension_seen=3;
457 
458  for (std::size_t i=0; i<elem_dimensions_seen.size(); ++i)
459  if (elem_dimensions_seen[i])
460  {
461  // Debugging
462  // libMesh::out << "Seen elements of dimension " << i+1 << std::endl;
463  max_elem_dimension_seen =
464  std::max(max_elem_dimension_seen, cast_int<unsigned char>(i+1));
465  min_elem_dimension_seen =
466  std::min(min_elem_dimension_seen, cast_int<unsigned char>(i+1));
467  }
468 
469  // Debugging:
470  // libMesh::out << "max_elem_dimension_seen=" << max_elem_dimension_seen << std::endl;
471  // libMesh::out << "min_elem_dimension_seen=" << min_elem_dimension_seen << std::endl;
472 
473  // If the difference between the max and min element dimension seen is larger than
474  // 1, (e.g. the file has 1D and 3D elements only) we don't handle this case.
475  if (max_elem_dimension_seen - min_elem_dimension_seen > 1)
476  libmesh_error_msg("Cannot handle meshes with dimension mismatch greater than 1.");
477 
478  // How many different element dimensions did we see while reading from file?
479  unsigned n_dims_seen = std::accumulate(elem_dimensions_seen.begin(),
480  elem_dimensions_seen.end(),
481  static_cast<unsigned>(0),
482  std::plus<unsigned>());
483 
484  // Have not yet tested a case where 1, 2, and 3D elements are all in the same Mesh,
485  // though it should theoretically be possible to handle.
486  if (n_dims_seen == 3)
487  libmesh_error_msg("Reading meshes with 1, 2, and 3D elements not currently supported.");
488 
489  // Set mesh_dimension based on the largest element dimension seen.
490  mesh.set_mesh_dimension(max_elem_dimension_seen);
491 
492  // Now that we know the maximum element dimension seen,
493  // we know whether the physical names are subdomain
494  // names or sideset names.
495  {
496  std::map<int, GmshPhysical>::iterator it = gmsh_physicals.begin();
497  for (; it != gmsh_physicals.end(); ++it)
498  {
499  // Extract data
500  int phys_id = it->first;
501  unsigned phys_dim = it->second.first;
502  std::string phys_name = it->second.second;
503 
504  // If the physical's dimension matches the largest
505  // dimension we've seen, it's a subdomain name.
506  if (phys_dim == max_elem_dimension_seen)
507  mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
508 
509  // Otherwise, if it's not a lower-dimensional
510  // block, it's a sideset name.
511  else if (phys_dim < max_elem_dimension_seen &&
512  !lower_dimensional_blocks.count(cast_int<boundary_id_type>(phys_id)))
513  mesh.get_boundary_info().sideset_name(cast_int<boundary_id_type>(phys_id)) = phys_name;
514  }
515  }
516 
517  if (n_dims_seen > 1)
518  {
519  // Store lower-dimensional elements in a map sorted
520  // by Elem::key(). We use a multimap for two reasons:
521  // 1.) The hash function is not guaranteed to be
522  // unique, so different lower-dimensional elements
523  // could theoretically hash to the same value,
524  // although this is pretty unlikely.
525  // 2.) The Gmsh file may contain multiple
526  // lower-dimensional elements for a single side in
527  // order to implement multiple boundary ids for a
528  // single side. These lower-dimensional elements
529  // will all hash to the same value, and we need to
530  // be able to store all of them.
531  typedef std::unordered_multimap<dof_id_type, Elem *> provide_container_t;
532  provide_container_t provide_bcs;
533 
534  // 1st loop over active elements - get info about lower-dimensional elements.
535  for (auto & elem : mesh.active_element_ptr_range())
536  if (elem->dim() < max_elem_dimension_seen &&
537  !lower_dimensional_blocks.count(elem->subdomain_id()))
538  {
539  // To be consistent with the previous
540  // GmshIO behavior, add all the
541  // lower-dimensional elements' nodes to
542  // the Mesh's BoundaryInfo object with the
543  // lower-dimensional element's subdomain
544  // ID.
545  for (auto n : elem->node_index_range())
546  mesh.get_boundary_info().add_node(elem->node_id(n),
547  elem->subdomain_id());
548 
549  // Store this elem in a quickly-searchable
550  // container to use it to assign boundary
551  // conditions later.
552  provide_bcs.insert(std::make_pair(elem->key(), elem));
553  }
554 
555  // 2nd loop over active elements - use lower dimensional element data to set BCs for higher dimensional elements
556  for (auto & elem : mesh.active_element_ptr_range())
557  if (elem->dim() == max_elem_dimension_seen)
558  {
559  // This is a max-dimension element that
560  // may require BCs. For each of its
561  // sides, including internal sides, we'll
562  // see if one more more lower-dimensional elements
563  // provides boundary information for it.
564  // Note that we have not yet called
565  // find_neighbors(), so we can't use
566  // elem->neighbor(sn) in this algorithm...
567  for (auto sn : elem->side_index_range())
568  for (const auto & pr : as_range(provide_bcs.equal_range(elem->key(sn))))
569  {
570  // For each side side in the provide_bcs multimap...
571  // Construct the side for hash verification.
572  UniquePtr<Elem> side (elem->build_side_ptr(sn));
573 
574  // Construct the lower-dimensional element to compare to the side.
575  Elem * lower_dim_elem = pr.second;
576 
577  // This was a hash, so it might not be perfect. Let's verify...
578  if (*lower_dim_elem == *side)
579  {
580  // Add the lower-dimensional
581  // element's subdomain_id as a
582  // boundary_id for the
583  // higher-dimensional element.
584  boundary_id_type bid = cast_int<boundary_id_type>(lower_dim_elem->subdomain_id());
585  mesh.get_boundary_info().add_side(elem, sn, bid);
586  }
587  }
588  }
589 
590  // 3rd loop over active elements - Remove the lower-dimensional elements
591  for (auto & elem : mesh.active_element_ptr_range())
592  if (elem->dim() < max_elem_dimension_seen &&
593  !lower_dimensional_blocks.count(elem->subdomain_id()))
594  mesh.delete_elem(elem);
595  } // end if (n_dims_seen > 1)
596  } // if $ELM
597 
598  continue;
599  } // if (in)
600 
601 
602  // If !in, check to see if EOF was set. If so, break out
603  // of while loop.
604  if (in.eof())
605  break;
606 
607  // If !in and !in.eof(), stream is in a bad state!
608  libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
609 
610  } // while true
611 }
612 
613 
614 
615 void GmshIO::write (const std::string & name)
616 {
618  {
619  // Open the output file stream
620  std::ofstream out_stream (name.c_str());
621 
622  // Make sure it opened correctly
623  if (!out_stream.good())
624  libmesh_file_error(name.c_str());
625 
626  this->write_mesh (out_stream);
627  }
628 }
629 
630 
631 
632 void GmshIO::write_nodal_data (const std::string & fname,
633  const std::vector<Number> & soln,
634  const std::vector<std::string> & names)
635 {
636  LOG_SCOPE("write_nodal_data()", "GmshIO");
637 
639  this->write_post (fname, &soln, &names);
640 }
641 
642 
643 
644 void GmshIO::write_mesh (std::ostream & out_stream)
645 {
646  // Be sure that the stream is valid.
647  libmesh_assert (out_stream.good());
648 
649  // Get a const reference to the mesh
651 
652  // If requested, write out lower-dimensional elements for
653  // element-side-based boundary conditions.
654  unsigned int n_boundary_faces = 0;
656  n_boundary_faces = mesh.get_boundary_info().n_boundary_conds();
657 
658  // Note: we are using version 2.0 of the gmsh output format.
659 
660  // Write the file header.
661  out_stream << "$MeshFormat\n";
662  out_stream << "2.0 0 " << sizeof(Real) << '\n';
663  out_stream << "$EndMeshFormat\n";
664 
665  // write the nodes in (n x y z) format
666  out_stream << "$Nodes\n";
667  out_stream << mesh.n_nodes() << '\n';
668 
669  for (unsigned int v=0; v<mesh.n_nodes(); v++)
670  out_stream << mesh.node_ref(v).id()+1 << " "
671  << mesh.node_ref(v)(0) << " "
672  << mesh.node_ref(v)(1) << " "
673  << mesh.node_ref(v)(2) << '\n';
674  out_stream << "$EndNodes\n";
675 
676  {
677  // write the connectivity
678  out_stream << "$Elements\n";
679  out_stream << mesh.n_active_elem() + n_boundary_faces << '\n';
680 
681  // loop over the elements
682  for (const auto & elem : mesh.active_element_ptr_range())
683  {
684  // Make sure we have a valid entry for
685  // the current element type.
686  libmesh_assert (_element_maps.out.count(elem->type()));
687 
688  // consult the export element table
689  std::map<ElemType, ElementDefinition>::iterator def_it =
690  _element_maps.out.find(elem->type());
691 
692  // Assert that we found it
693  if (def_it == _element_maps.out.end())
694  libmesh_error_msg("Element type " << elem->type() << " not found in _element_maps.");
695 
696  // Get a reference to the ElementDefinition object
697  const ElementDefinition & eletype = def_it->second;
698 
699  // The element mapper better not require any more nodes
700  // than are present in the current element!
701  libmesh_assert_less_equal (eletype.nodes.size(), elem->n_nodes());
702 
703  // elements ids are 1 based in Gmsh
704  out_stream << elem->id()+1 << " ";
705  // element type
706  out_stream << eletype.gmsh_type;
707 
708  // write the number of tags (3) and their values:
709  // 1 (physical entity)
710  // 2 (geometric entity)
711  // 3 (partition entity)
712  out_stream << " 3 "
713  << static_cast<unsigned int>(elem->subdomain_id())
714  << " 0 "
715  << elem->processor_id()+1
716  << " ";
717 
718  // if there is a node translation table, use it
719  if (eletype.nodes.size() > 0)
720  for (auto i : elem->node_index_range())
721  out_stream << elem->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
722  // otherwise keep the same node order
723  else
724  for (auto i : elem->node_index_range())
725  out_stream << elem->node_id(i)+1 << " "; // gmsh is 1-based
726  out_stream << "\n";
727  } // element loop
728  }
729 
730  {
731  // A counter for writing surface elements to the Gmsh file
732  // sequentially. We start numbering them with a number strictly
733  // larger than the largest element ID in the mesh. Note: the
734  // MeshBase docs say "greater than or equal to" the maximum
735  // element id in the mesh, so technically we might need a +1 here,
736  // but all of the implementations return an ID strictly greater
737  // than the largest element ID in the Mesh.
738  unsigned int e_id = mesh.max_elem_id();
739 
740  // loop over the elements, writing out boundary faces
741  if (n_boundary_faces)
742  {
743  // Construct the list of boundary sides
744  std::vector<dof_id_type> element_id_list;
745  std::vector<unsigned short int> side_list;
746  std::vector<boundary_id_type> bc_id_list;
747 
748  mesh.get_boundary_info().build_side_list(element_id_list, side_list, bc_id_list);
749 
750  // Loop over these lists, writing data to the file.
751  for (std::size_t idx=0; idx<element_id_list.size(); ++idx)
752  {
753  const Elem & elem = mesh.elem_ref(element_id_list[idx]);
754 
755  UniquePtr<const Elem> side = elem.build_side_ptr(side_list[idx]);
756 
757  // Map from libmesh elem type to gmsh elem type.
758  std::map<ElemType, ElementDefinition>::iterator def_it =
759  _element_maps.out.find(side->type());
760 
761  // If we didn't find it, that's an error
762  if (def_it == _element_maps.out.end())
763  libmesh_error_msg("Element type " << side->type() << " not found in _element_maps.");
764 
765  // consult the export element table
766  const GmshIO::ElementDefinition & eletype = def_it->second;
767 
768  // The element mapper better not require any more nodes
769  // than are present in the current element!
770  libmesh_assert_less_equal (eletype.nodes.size(), side->n_nodes());
771 
772  // elements ids are 1-based in Gmsh
773  out_stream << e_id+1 << " ";
774 
775  // element type
776  out_stream << eletype.gmsh_type;
777 
778  // write the number of tags:
779  // 1 (physical entity)
780  // 2 (geometric entity)
781  // 3 (partition entity)
782  out_stream << " 3 "
783  << bc_id_list[idx]
784  << " 0 "
785  << elem.processor_id()+1
786  << " ";
787 
788  // if there is a node translation table, use it
789  if (eletype.nodes.size() > 0)
790  for (auto i : side->node_index_range())
791  out_stream << side->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
792 
793  // otherwise keep the same node order
794  else
795  for (auto i : side->node_index_range())
796  out_stream << side->node_id(i)+1 << " "; // gmsh is 1-based
797 
798  // Go to the next line
799  out_stream << "\n";
800 
801  // increment this index too...
802  ++e_id;
803  }
804  }
805 
806  out_stream << "$EndElements\n";
807  }
808 }
809 
810 
811 
812 void GmshIO::write_post (const std::string & fname,
813  const std::vector<Number> * v,
814  const std::vector<std::string> * solution_names)
815 {
816 
817  // Should only do this on processor 0!
818  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
819 
820  // Create an output stream
821  std::ofstream out_stream(fname.c_str());
822 
823  // Make sure it opened correctly
824  if (!out_stream.good())
825  libmesh_file_error(fname.c_str());
826 
827  // create a character buffer
828  char buf[80];
829 
830  // Get a constant reference to the mesh.
832 
833  // write the data
834  if ((solution_names != libmesh_nullptr) && (v != libmesh_nullptr))
835  {
836  const unsigned int n_vars =
837  cast_int<unsigned int>(solution_names->size());
838 
839  if (!(v->size() == mesh.n_nodes()*n_vars))
840  libMesh::err << "ERROR: v->size()=" << v->size()
841  << ", mesh.n_nodes()=" << mesh.n_nodes()
842  << ", n_vars=" << n_vars
843  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
844  << "\n";
845 
846  libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
847 
848  // write the header
849  out_stream << "$PostFormat\n";
850  if (this->binary())
851  out_stream << "1.2 1 " << sizeof(double) << "\n";
852  else
853  out_stream << "1.2 0 " << sizeof(double) << "\n";
854  out_stream << "$EndPostFormat\n";
855 
856  // Loop over the elements to see how much of each type there are
857  unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
858  n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
859  unsigned int n_scalar=0, n_vector=0, n_tensor=0;
860  unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;
861 
862  {
863  for (const auto & elem : mesh.active_element_ptr_range())
864  {
865  const ElemType elemtype = elem->type();
866 
867  switch (elemtype)
868  {
869  case EDGE2:
870  case EDGE3:
871  case EDGE4:
872  {
873  n_lines += 1;
874  break;
875  }
876  case TRI3:
877  case TRI6:
878  {
879  n_triangles += 1;
880  break;
881  }
882  case QUAD4:
883  case QUAD8:
884  case QUAD9:
885  {
886  n_quadrangles += 1;
887  break;
888  }
889  case TET4:
890  case TET10:
891  {
892  n_tetrahedra += 1;
893  break;
894  }
895  case HEX8:
896  case HEX20:
897  case HEX27:
898  {
899  n_hexahedra += 1;
900  break;
901  }
902  case PRISM6:
903  case PRISM15:
904  case PRISM18:
905  {
906  n_prisms += 1;
907  break;
908  }
909  case PYRAMID5:
910  {
911  n_pyramids += 1;
912  break;
913  }
914  default:
915  libmesh_error_msg("ERROR: Nonexistent element type " << elem->type());
916  }
917  }
918  }
919 
920  // create a view for each variable
921  for (unsigned int ivar=0; ivar < n_vars; ivar++)
922  {
923  std::string varname = (*solution_names)[ivar];
924 
925  // at the moment, we just write out scalar quantities
926  // later this should be made configurable through
927  // options to the writer class
928  n_scalar = 1;
929 
930  // write the variable as a view, and the number of time steps
931  out_stream << "$View\n" << varname << " " << 1 << "\n";
932 
933  // write how many of each geometry type are written
934  out_stream << n_points * n_scalar << " "
935  << n_points * n_vector << " "
936  << n_points * n_tensor << " "
937  << n_lines * n_scalar << " "
938  << n_lines * n_vector << " "
939  << n_lines * n_tensor << " "
940  << n_triangles * n_scalar << " "
941  << n_triangles * n_vector << " "
942  << n_triangles * n_tensor << " "
943  << n_quadrangles * n_scalar << " "
944  << n_quadrangles * n_vector << " "
945  << n_quadrangles * n_tensor << " "
946  << n_tetrahedra * n_scalar << " "
947  << n_tetrahedra * n_vector << " "
948  << n_tetrahedra * n_tensor << " "
949  << n_hexahedra * n_scalar << " "
950  << n_hexahedra * n_vector << " "
951  << n_hexahedra * n_tensor << " "
952  << n_prisms * n_scalar << " "
953  << n_prisms * n_vector << " "
954  << n_prisms * n_tensor << " "
955  << n_pyramids * n_scalar << " "
956  << n_pyramids * n_vector << " "
957  << n_pyramids * n_tensor << " "
958  << nb_text2d << " "
959  << nb_text2d_chars << " "
960  << nb_text3d << " "
961  << nb_text3d_chars << "\n";
962 
963  // if binary, write a marker to identify the endianness of the file
964  if (this->binary())
965  {
966  const int one = 1;
967  std::memcpy(buf, &one, sizeof(int));
968  out_stream.write(buf, sizeof(int));
969  }
970 
971  // the time steps (there is just 1 at the moment)
972  if (this->binary())
973  {
974  double one = 1;
975  std::memcpy(buf, &one, sizeof(double));
976  out_stream.write(buf, sizeof(double));
977  }
978  else
979  out_stream << "1\n";
980 
981  // Loop over the elements and write out the data
982  for (const auto & elem : mesh.active_element_ptr_range())
983  {
984  // this is quite crappy, but I did not invent that file format!
985  for (unsigned int d=0; d<3; d++) // loop over the dimensions
986  {
987  for (unsigned int n=0; n < elem->n_vertices(); n++) // loop over vertices
988  {
989  const Point & vertex = elem->point(n);
990  if (this->binary())
991  {
992  double tmp = vertex(d);
993  std::memcpy(buf, &tmp, sizeof(double));
994  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
995  }
996  else
997  out_stream << vertex(d) << " ";
998  }
999  if (!this->binary())
1000  out_stream << "\n";
1001  }
1002 
1003  // now finally write out the data
1004  for (unsigned int i=0; i < elem->n_vertices(); i++) // loop over vertices
1005  if (this->binary())
1006  {
1007 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1008  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1009  << "complex numbers. Will only write the real part of "
1010  << "variable " << varname << std::endl;
1011 #endif
1012  double tmp = libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]);
1013  std::memcpy(buf, &tmp, sizeof(double));
1014  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1015  }
1016  else
1017  {
1018 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1019  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1020  << "complex numbers. Will only write the real part of "
1021  << "variable " << varname << std::endl;
1022 #endif
1023  out_stream << libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]) << "\n";
1024  }
1025  }
1026  if (this->binary())
1027  out_stream << "\n";
1028  out_stream << "$EndView\n";
1029 
1030  } // end variable loop (writing the views)
1031  }
1032 }
1033 
1034 } // namespace libMesh
T libmesh_real(T a)
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
OStreamProxy err
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
bool _binary
Flag to write binary data.
Definition: gmsh_io.h:142
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:238
ElemType
Defines an enum for geometric element types.
bool & write_lower_dimensional_elements()
Access to the flag which controls whether boundary elements are written to the Mesh file...
Definition: gmsh_io.C:137
void write_post(const std::string &, const std::vector< Number > *=libmesh_nullptr, const std::vector< std::string > *=libmesh_nullptr)
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: gmsh_io.C:812
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
const MT & mesh() const
Definition: mesh_output.h:216
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
bool & binary()
Flag indicating whether or not to write a binary file.
Definition: gmsh_io.C:130
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
std::map< unsigned int, ElementDefinition > in
Definition: gmsh_io.h:186
long double max(long double a, double b)
This is the MeshBase class.
Definition: mesh_base.h:68
dof_id_type & set_id()
Definition: dof_object.h:641
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
std::map< ElemType, ElementDefinition > out
Definition: gmsh_io.h:185
PetscErrorCode Vec x
int8_t boundary_id_type
Definition: id_types.h:51
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:50
SimpleRange< I > as_range(const std::pair< I, I > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
void write_mesh(std::ostream &out)
This method implements writing a mesh to a specified file.
Definition: gmsh_io.C:644
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
bool _write_lower_dimensional_elements
If true, lower-dimensional elements based on the boundary conditions get written to the output file...
Definition: gmsh_io.h:148
Defines mapping from libMesh element types to Gmsh element types or vice-versa.
Definition: gmsh_io.h:153
static ElementMaps build_element_maps()
A static function used to construct the _element_maps struct, statically.
Definition: gmsh_io.C:42
virtual void write(const std::string &name) libmesh_override
This method implements writing a mesh to a specified file in the Gmsh *.msh format.
Definition: gmsh_io.C:615
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
X_input swap(X_system)
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) libmesh_override
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: gmsh_io.C:632
OStreamProxy out
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
Definition: gmsh_io.h:193
virtual void read(const std::string &name) libmesh_override
Reads in a mesh in the Gmsh *.msh format from the ASCII file given by name.
Definition: gmsh_io.C:144
GmshIO(MeshBase &mesh)
Constructor.
Definition: gmsh_io.C:120
void add_def(const ElementDefinition &eledef)
Definition: gmsh_io.h:179
virtual dof_id_type n_nodes() const =0
processor_id_type processor_id()
Definition: libmesh_base.h:96
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
struct which holds a map from Gmsh to libMesh element numberings and vice-versa.
Definition: gmsh_io.h:176
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
void read_mesh(std::istream &in)
Implementation of the read() function.
Definition: gmsh_io.C:152
processor_id_type processor_id() const
Definition: dof_object.h:694
std::vector< unsigned int > nodes
Definition: gmsh_io.h:169