libMesh
Classes | Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
libMesh::GmshIO Class Reference

Reading and writing meshes in the Gmsh format. More...

#include <gmsh_io.h>

Inheritance diagram for libMesh::GmshIO:
[legend]

Classes

struct  ElementDefinition
 Defines mapping from libMesh element types to Gmsh element types or vice-versa. More...
 
struct  ElementMaps
 struct which holds a map from Gmsh to libMesh element numberings and vice-versa. More...
 

Public Member Functions

 GmshIO (MeshBase &mesh)
 Constructor. More...
 
 GmshIO (const MeshBase &mesh)
 Constructor. More...
 
virtual void read (const std::string &name) override
 Reads in a mesh in the Gmsh *.msh format from the ASCII file given by name. More...
 
virtual void write (const std::string &name) override
 This method implements writing a mesh to a specified file in the Gmsh *.msh format. More...
 
virtual void write_nodal_data (const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
 This method implements writing a mesh with nodal data to a specified file where the nodal data and variable names are provided. More...
 
bool & binary ()
 Flag indicating whether or not to write a binary file. More...
 
bool & write_lower_dimensional_elements ()
 Access to the flag which controls whether boundary elements are written to the Mesh file. More...
 
bool is_parallel_format () const
 Returns true iff this mesh file format and input class are parallelized, so that all processors can read their share of the data at once. More...
 
virtual void write_equation_systems (const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
 This method implements writing a mesh with data to a specified file where the data is taken from the EquationSystems object. More...
 
virtual void write_discontinuous_equation_systems (const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
 This method implements writing a mesh with discontinuous data to a specified file where the data is taken from the EquationSystems object. More...
 
virtual void write_nodal_data (const std::string &, const NumericVector< Number > &, const std::vector< std::string > &)
 This method may be overridden by "parallel" output formats for writing nodal data. More...
 
virtual void write_nodal_data (const std::string &, const EquationSystems &, const std::set< std::string > *)
 This method should be overridden by "parallel" output formats for writing nodal data. More...
 
virtual void write_nodal_data_discontinuous (const std::string &, const std::vector< Number > &, const std::vector< std::string > &)
 This method implements writing a mesh with discontinuous data to a specified file where the nodal data and variables names are provided. More...
 
unsigned intascii_precision ()
 Return/set the precision to use when writing ASCII files. More...
 

Protected Member Functions

MeshBasemesh ()
 
void set_n_partitions (unsigned int n_parts)
 Sets the number of partitions in the mesh. More...
 
void skip_comment_lines (std::istream &in, const char comment_start)
 Reads input from in, skipping all the lines that start with the character comment_start. More...
 
const MeshBasemesh () const
 
virtual bool get_add_sides ()
 

Protected Attributes

std::vector< bool > elems_of_dimension
 A vector of bools describing what dimension elements have been encountered when reading a mesh. More...
 
const bool _is_parallel_format
 Flag specifying whether this format is parallel-capable. More...
 
const bool _serial_only_needed_on_proc_0
 Flag specifying whether this format can be written by only serializing the mesh to processor zero. More...
 

Private Member Functions

void read_mesh (std::istream &in)
 Implementation of the read() function. More...
 
void write_mesh (std::ostream &out)
 This method implements writing a mesh to a specified file. More...
 
void write_post (const std::string &, const std::vector< Number > *=nullptr, const std::vector< std::string > *=nullptr)
 This method implements writing a mesh with nodal data to a specified file where the nodal data and variable names are optionally provided. More...
 

Static Private Member Functions

static ElementMaps build_element_maps ()
 A static function used to construct the _element_maps struct, statically. More...
 

Private Attributes

bool _binary
 Flag to write binary data. More...
 
bool _write_lower_dimensional_elements
 If true, lower-dimensional elements based on the boundary conditions get written to the output file. More...
 

Static Private Attributes

static ElementMaps _element_maps = GmshIO::build_element_maps()
 A static ElementMaps object that is built statically and used by all instances of this class. More...
 

Detailed Description

Reading and writing meshes in the Gmsh format.

For a full description of the Gmsh format and to obtain the GMSH software see the Gmsh home page

Author
John W. Peterson
Date
2004, 2014
Author
Martin Luthi
Date
2005

Definition at line 51 of file gmsh_io.h.

Constructor & Destructor Documentation

◆ GmshIO() [1/2]

libMesh::GmshIO::GmshIO ( MeshBase mesh)
explicit

Constructor.

Takes a non-const Mesh reference which it will fill up with elements via the read() command.

Definition at line 125 of file gmsh_io.C.

125  :
126  MeshInput<MeshBase> (mesh),
128  _binary (false),
130 {
131 }
bool _binary
Flag to write binary data.
Definition: gmsh_io.h:149
template class LIBMESH_EXPORT MeshOutput< MeshBase >
Definition: mesh_output.C:180
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:155

◆ GmshIO() [2/2]

libMesh::GmshIO::GmshIO ( const MeshBase mesh)
explicit

Constructor.

Takes a reference to a constant mesh object. This constructor will only allow us to write the mesh.

Definition at line 116 of file gmsh_io.C.

116  :
118  _binary(false),
120 {
121 }
bool _binary
Flag to write binary data.
Definition: gmsh_io.h:149
template class LIBMESH_EXPORT MeshOutput< MeshBase >
Definition: mesh_output.C:180
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:155

Member Function Documentation

◆ ascii_precision()

unsigned int & libMesh::MeshOutput< MeshBase >::ascii_precision ( )
inlineinherited

Return/set the precision to use when writing ASCII files.

By default we use numeric_limits<Real>::max_digits10, which should be enough to write out to ASCII and get the exact same Real back when reading in.

Definition at line 269 of file mesh_output.h.

Referenced by libMesh::UNVIO::nodes_out(), libMesh::FroIO::write(), libMesh::STLIO::write(), libMesh::MEDITIO::write_ascii(), libMesh::TecplotIO::write_ascii(), libMesh::GMVIO::write_ascii_new_impl(), and libMesh::GMVIO::write_ascii_old_impl().

270 {
271  return _ascii_precision;
272 }
unsigned int _ascii_precision
Precision to use when writing ASCII files.
Definition: mesh_output.h:207

◆ binary()

bool & libMesh::GmshIO::binary ( )

Flag indicating whether or not to write a binary file.

While binary files may end up being smaller than equivalent ASCII files, they will almost certainly take longer to write. The reason for this is that the ostream::write() function which is used to write "binary" data to streams, only takes a pointer to char as its first argument. This means if you want to write anything other than a buffer of chars, you first have to use a strange memcpy hack to get the data into the desired format. See the templated to_binary_stream() function below.

Definition at line 135 of file gmsh_io.C.

References _binary.

Referenced by write_post().

136 {
137  return _binary;
138 }
bool _binary
Flag to write binary data.
Definition: gmsh_io.h:149

◆ build_element_maps()

GmshIO::ElementMaps libMesh::GmshIO::build_element_maps ( )
staticprivate

A static function used to construct the _element_maps struct, statically.

Definition at line 47 of file gmsh_io.C.

References libMesh::GmshIO::ElementMaps::add_def(), libMesh::EDGE2, libMesh::EDGE3, libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::GmshIO::ElementMaps::in, libMesh::GmshIO::ElementDefinition::nnodes, libMesh::NODEELEM, libMesh::GmshIO::ElementDefinition::nodes, libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::TET10, libMesh::TET4, libMesh::TRI3, and libMesh::TRI6.

48 {
49  // Object to be filled up
50  ElementMaps em;
51 
52  // POINT (import only)
53  em.in.emplace(15, ElementDefinition(NODEELEM, 15, 0, 1));
54 
55  // Add elements with trivial node mappings
56  em.add_def(ElementDefinition(EDGE2, 1, 1, 2));
57  em.add_def(ElementDefinition(EDGE3, 8, 1, 3));
58  em.add_def(ElementDefinition(TRI3, 2, 2, 3));
59  em.add_def(ElementDefinition(TRI6, 9, 2, 6));
60  em.add_def(ElementDefinition(QUAD4, 3, 2, 4));
61  em.add_def(ElementDefinition(QUAD8, 16, 2, 8));
62  em.add_def(ElementDefinition(QUAD9, 10, 2, 9));
63  em.add_def(ElementDefinition(HEX8, 5, 3, 8));
64  em.add_def(ElementDefinition(TET4, 4, 3, 4));
65  em.add_def(ElementDefinition(PRISM6, 6, 3, 6));
66  em.add_def(ElementDefinition(PYRAMID5, 7, 3, 5));
67 
68  // Add elements with non-trivial node mappings
69 
70  // HEX20
71  {
72  ElementDefinition eledef(HEX20, 17, 3, 20);
73  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,15,16,19,17,18};
74  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
75  em.add_def(eledef);
76  }
77 
78  // HEX27
79  {
80  ElementDefinition eledef(HEX27, 12, 3, 27);
81  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
82  15,16,19,17,18,20,21,24,22,23,25,26};
83  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
84  em.add_def(eledef);
85  }
86 
87  // TET10
88  {
89  ElementDefinition eledef(TET10, 11, 3, 10);
90  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
91  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
92  em.add_def(eledef);
93  }
94 
95  // PRISM15
96  {
97  ElementDefinition eledef(PRISM15, 18, 3, 15);
98  const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13};
99  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
100  em.add_def(eledef);
101  }
102 
103  // PRISM18
104  {
105  ElementDefinition eledef(PRISM18, 13, 3, 18);
106  const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13,15,17,16};
107  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
108  em.add_def(eledef);
109  }
110 
111  return em;
112 }

◆ get_add_sides()

virtual bool libMesh::MeshOutput< MeshBase >::get_add_sides ( )
inlineprotectedvirtualinherited
Returns
Whether or not added sides are expected to be output, to plot SIDE_DISCONTINUOUS data. Subclasses should override this if they are capable of plotting such data.

Reimplemented in libMesh::ExodusII_IO.

Definition at line 176 of file mesh_output.h.

176 { return false; }

◆ is_parallel_format()

bool libMesh::MeshInput< MeshBase >::is_parallel_format ( ) const
inlineinherited

Returns true iff this mesh file format and input class are parallelized, so that all processors can read their share of the data at once.

Definition at line 87 of file mesh_input.h.

References libMesh::MeshInput< MT >::_is_parallel_format.

87 { return this->_is_parallel_format; }
const bool _is_parallel_format
Flag specifying whether this format is parallel-capable.
Definition: mesh_input.h:130

◆ mesh() [1/2]

MeshBase & libMesh::MeshInput< MeshBase >::mesh ( )
inlineprotectedinherited
Returns
The object as a writable reference.

Definition at line 178 of file mesh_input.h.

Referenced by libMesh::GMVIO::_read_one_cell(), libMesh::VTKIO::cells_to_vtk(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::Nemesis_IO::copy_elemental_solution(), libMesh::ExodusII_IO::copy_nodal_solution(), libMesh::TetGenIO::element_in(), libMesh::UNVIO::elements_in(), libMesh::UNVIO::elements_out(), libMesh::VTKIO::get_local_node_values(), libMesh::ExodusII_IO::get_sideset_data_indices(), libMesh::UNVIO::groups_in(), libMesh::TetGenIO::node_in(), libMesh::UNVIO::nodes_in(), libMesh::UNVIO::nodes_out(), libMesh::VTKIO::nodes_to_vtk(), libMesh::Nemesis_IO::prepare_to_write_nodal_data(), libMesh::GMVIO::read(), libMesh::STLIO::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::ExodusII_IO::read(), libMesh::CheckpointIO::read(), libMesh::VTKIO::read(), libMesh::STLIO::read_ascii(), libMesh::CheckpointIO::read_bcs(), libMesh::STLIO::read_binary(), libMesh::CheckpointIO::read_connectivity(), libMesh::ExodusII_IO::read_header(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::UCDIO::read_implementation(), libMesh::UNVIO::read_implementation(), read_mesh(), libMesh::DynaIO::read_mesh(), libMesh::CheckpointIO::read_nodes(), libMesh::CheckpointIO::read_nodesets(), libMesh::CheckpointIO::read_remote_elem(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::ExodusII_IO::read_sideset_data(), libMesh::OFFIO::read_stream(), libMesh::MatlabIO::read_stream(), libMesh::CheckpointIO::read_subdomain_names(), libMesh::STLIO::write(), libMesh::TetGenIO::write(), libMesh::Nemesis_IO::write(), libMesh::XdrIO::write(), libMesh::CheckpointIO::write(), libMesh::ExodusII_IO::write(), libMesh::GMVIO::write_ascii_new_impl(), libMesh::GMVIO::write_ascii_old_impl(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::Nemesis_IO::write_element_data(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO::write_elemsets(), libMesh::UCDIO::write_header(), libMesh::UCDIO::write_implementation(), libMesh::UCDIO::write_interior_elems(), write_mesh(), libMesh::UCDIO::write_nodal_data(), libMesh::VTKIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_common(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::UCDIO::write_nodes(), libMesh::CheckpointIO::write_nodesets(), libMesh::XdrIO::write_parallel(), write_post(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::ExodusII_IO::write_sideset_data(), libMesh::UCDIO::write_soln(), and libMesh::CheckpointIO::write_subdomain_names().

179 {
180  libmesh_error_msg_if(_obj == nullptr, "ERROR: _obj should not be nullptr!");
181  return *_obj;
182 }
MeshBase * _obj
A pointer to a non-const object object.
Definition: mesh_input.h:123

◆ mesh() [2/2]

const MeshBase & libMesh::MeshOutput< MeshBase >::mesh ( ) const
inlineprotectedinherited

◆ read()

void libMesh::GmshIO::read ( const std::string &  name)
overridevirtual

Reads in a mesh in the Gmsh *.msh format from the ASCII file given by name.

Note
The user is responsible for calling Mesh::prepare_for_use() after reading the mesh and before using it.

The physical group names defined in the Gmsh-file are stored, depending on the element dimension, as either subdomain name or as side name. The IDs of the former can be retrieved by using the MeshBase::get_id_by_name method; the IDs of the latter can be retrieved by using the MeshBase::get_boundary_info and the BoundaryInfo::get_id_by_name methods.

Implements libMesh::MeshInput< MeshBase >.

Definition at line 149 of file gmsh_io.C.

References libMesh::Quality::name(), and read_mesh().

Referenced by libMesh::NameBasedIO::read(), MeshInputTest::testBadGmsh(), MeshInputTest::testGmshBCIDOverlap(), and MeshInputTest::testGoodGmsh().

150 {
151  std::ifstream in (name.c_str());
152  this->read_mesh (in);
153 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
void read_mesh(std::istream &in)
Implementation of the read() function.
Definition: gmsh_io.C:157

◆ read_mesh()

void libMesh::GmshIO::read_mesh ( std::istream &  in)
private

Implementation of the read() function.

This function is called by the public interface function and implements reading the file.

Definition at line 157 of file gmsh_io.C.

References _element_maps, libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_node(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::as_range(), libMesh::Elem::build_with_id(), libMesh::MeshBase::clear(), libMesh::BoundingBox::contains_point(), libMesh::MeshBase::delete_elem(), libMesh::GmshIO::ElementDefinition::dim, libMesh::Utility::enum_to_string(), libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::GmshIO::ElementMaps::in, libMesh::index_range(), libMesh::is, libMesh::libmesh_assert(), libMesh::BoundingBox::max(), libMesh::MeshInput< MeshBase >::mesh(), libMesh::MeshInput< MT >::mesh(), libMesh::BoundingBox::min(), libMesh::Elem::n_nodes(), libMesh::GmshIO::ElementDefinition::nnodes, libMesh::MeshBase::node_ptr(), libMesh::GmshIO::ElementDefinition::nodes, libMesh::BoundaryInfo::nodeset_name(), libMesh::Real, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), libMesh::BoundaryInfo::sideset_name(), libMesh::Elem::subdomain_id(), libMesh::MeshBase::subdomain_name(), libMesh::TOLERANCE, and libMesh::GmshIO::ElementDefinition::type.

Referenced by read().

158 {
159  // This is a serial-only process for now;
160  // the Mesh should be read on processor 0 and
161  // broadcast later
162  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
163 
164  libmesh_assert(in.good());
165 
166  // clear any data in the mesh
167  MeshBase & mesh = MeshInput<MeshBase>::mesh();
168  mesh.clear();
169 
170  // some variables
171  int format=0, size=0;
172  Real version = 1.0;
173 
174  // Keep track of lower-dimensional blocks which are not BCs, but
175  // actually blocks of lower-dimensional elements.
176  std::set<subdomain_id_type> lower_dimensional_blocks;
177 
178  // Mapping from (physical id, physical dim) -> physical name.
179  // These can refer to either "sidesets" or "subdomains"; we need to
180  // wait until the Mesh has been read to know which is which. Note
181  // that we are using ptrdiff_t as the key here rather than
182  // subdomain_id_type or boundary_id_type, since at this point, it
183  // could be either.
184  std::map<std::pair<std::ptrdiff_t, unsigned>, std::string> gmsh_physicals;
185 
186  // map to hold the node numbers for translation
187  // note the the nodes can be non-consecutive
188  std::map<unsigned int, unsigned int> nodetrans;
189 
190  // Map from entity tag to physical id. The key is a pair with the first
191  // item being the dimension of the entity and the second item being
192  // the entity tag/id
193  std::map<std::pair<unsigned, int>, int> entity_to_physical_id;
194 
195  // Map from entities to bounding boxes. The key is the same.
196  std::map<std::pair<unsigned, int>, BoundingBox> entity_to_bounding_box;
197 
198  // For reading the file line by line
199  std::string s;
200 
201  while (true)
202  {
203  // Try to read something. This may set EOF!
204  std::getline(in, s);
205 
206  if (in)
207  {
208  // Process s...
209 
210  if (s.find("$MeshFormat") == static_cast<std::string::size_type>(0))
211  {
212  in >> version >> format >> size;
213 
214  // Some notes on gmsh mesh versions:
215  //
216  // Mesh version 2.0 goes back as far as I know. It's not explicitly
217  // mentioned here: http://www.geuz.org/gmsh/doc/VERSIONS.txt
218  //
219  // As of gmsh-2.4.0:
220  // bumped mesh version format to 2.1 (small change in the $PhysicalNames
221  // section, where the group dimension is now required);
222  // [Since we don't even parse the PhysicalNames section at the time
223  // of this writing, I don't think this change affects us.]
224  //
225  // Mesh version 2.2 tested by Manav Bhatia; no other
226  // libMesh code changes were required for support
227  //
228  // Mesh version 4.0 is a near complete rewrite of the previous mesh version
229  libmesh_error_msg_if(version < 2.0, "Error: Unknown msh file version " << version);
230  libmesh_error_msg_if(format, "Error: Unknown data format for mesh in Gmsh reader.");
231  }
232 
233  // Read and process the "PhysicalNames" section.
234  else if (s.find("$PhysicalNames") == static_cast<std::string::size_type>(0))
235  {
236  // The lines in the PhysicalNames section should look like the following:
237  // 2 1 "frac" lower_dimensional_block
238  // 2 3 "top"
239  // 2 4 "bottom"
240  // 3 2 "volume"
241 
242  // Read in the number of physical groups to expect in the file.
243  unsigned int num_physical_groups = 0;
244  in >> num_physical_groups;
245 
246  // Read rest of line including newline character.
247  std::getline(in, s);
248 
249  for (unsigned int i=0; i<num_physical_groups; ++i)
250  {
251  // Read an entire line of the PhysicalNames section.
252  std::getline(in, s);
253 
254  // Use an istringstream to extract the physical
255  // dimension, physical id, and physical name from
256  // this line.
257  std::istringstream s_stream(s);
258  unsigned phys_dim;
259  int phys_id;
260  std::string phys_name;
261  s_stream >> phys_dim >> phys_id >> phys_name;
262 
263  // Not sure if this is true for all Gmsh files, but
264  // my test file has quotes around the phys_name
265  // string. So let's erase any quotes now...
266  phys_name.erase(std::remove(phys_name.begin(), phys_name.end(), '"'), phys_name.end());
267 
268  // Record this ID for later assignment of subdomain/sideset names.
269  gmsh_physicals[std::make_pair(phys_id, phys_dim)] = phys_name;
270 
271  // If 's' also contains the libmesh-specific string
272  // "lower_dimensional_block", add this block ID to
273  // the list of blocks which are not boundary
274  // conditions.
275  if (s.find("lower_dimensional_block") != std::string::npos)
276  {
277  lower_dimensional_blocks.insert(cast_int<subdomain_id_type>(phys_id));
278 
279  // The user has explicitly told us that this
280  // block is a subdomain, so set that association
281  // in the Mesh.
282  mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
283  }
284  }
285  }
286 
287  else if (s.find("$Entities") == static_cast<std::string::size_type>(0))
288  {
289  if (version >= 4.0)
290  {
291  std::size_t num_point_entities, num_curve_entities, num_surface_entities, num_volume_entities;
292  in >> num_point_entities >> num_curve_entities >> num_surface_entities >> num_volume_entities;
293 
294  for (std::size_t n = 0; n < num_point_entities; ++n)
295  {
296  int point_tag, physical_tag;
297  Real x, y, z;
298  std::size_t num_physical_tags;
299  in >> point_tag >> x >> y >> z >> num_physical_tags;
300 
301  libmesh_error_msg_if(num_physical_tags > 1,
302  "Sorry, you cannot currently specify multiple subdomain or "
303  "boundary ids for a given geometric entity");
304 
305  entity_to_bounding_box[std::make_pair(0, point_tag)] =
306  BoundingBox(Point(x,y,z),Point(x,y,z));
307 
308  if (num_physical_tags)
309  {
310  in >> physical_tag;
311  entity_to_physical_id[std::make_pair(0, point_tag)] = physical_tag;
312  }
313  }
314  for (std::size_t n = 0; n < num_curve_entities; ++n)
315  {
316  int curve_tag, physical_tag;
317  Real minx, miny, minz, maxx, maxy, maxz;
318  std::size_t num_physical_tags;
319  in >> curve_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
320 
321  libmesh_error_msg_if(num_physical_tags > 1,
322  "I don't believe that we can specify multiple subdomain or "
323  "boundary ids for a given geometric entity");
324 
325  entity_to_bounding_box[std::make_pair(1, curve_tag)] =
326  BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
327 
328  if (num_physical_tags)
329  {
330  in >> physical_tag;
331  entity_to_physical_id[std::make_pair(1, curve_tag)] = physical_tag;
332  }
333 
334  // Read to end of line; this captures bounding information that we don't care about
335  std::getline(in, s);
336  }
337  for (std::size_t n = 0; n < num_surface_entities; ++n)
338  {
339  int surface_tag, physical_tag;
340  Real minx, miny, minz, maxx, maxy, maxz;
341  std::size_t num_physical_tags;
342  in >> surface_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
343 
344  libmesh_error_msg_if(num_physical_tags > 1,
345  "I don't believe that we can specify multiple subdomain or "
346  "boundary ids for a given geometric entity");
347 
348  entity_to_bounding_box[std::make_pair(2, surface_tag)] =
349  BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
350 
351  if (num_physical_tags)
352  {
353  in >> physical_tag;
354  entity_to_physical_id[std::make_pair(2, surface_tag)] = physical_tag;
355  }
356 
357  // Read to end of line; this captures bounding information that we don't care about
358  std::getline(in, s);
359  }
360  for (std::size_t n = 0; n < num_volume_entities; ++n)
361  {
362  int volume_tag, physical_tag;
363  Real minx, miny, minz, maxx, maxy, maxz;
364  std::size_t num_physical_tags;
365  in >> volume_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
366 
367  libmesh_error_msg_if(num_physical_tags > 1,
368  "I don't believe that we can specify multiple subdomain or "
369  "boundary ids for a given geometric entity");
370 
371  entity_to_bounding_box[std::make_pair(3, volume_tag)] =
372  BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
373 
374  if (num_physical_tags)
375  {
376  in >> physical_tag;
377  entity_to_physical_id[std::make_pair(3, volume_tag)] = physical_tag;
378  }
379 
380  // Read to end of line; this captures bounding information that we don't care about
381  std::getline(in, s);
382  }
383  // Read the $EndEntities
384  std::getline(in, s);
385  } // end if (version >= 4.0)
386 
387  else
388  libmesh_error_msg("The Entities block was introduced in mesh version 4.0");
389  } // end if $Entities
390 
391  // read the node block
392  else if (s.find("$NOD") == static_cast<std::string::size_type>(0) ||
393  s.find("$NOE") == static_cast<std::string::size_type>(0) ||
394  s.find("$Nodes") == static_cast<std::string::size_type>(0))
395  {
396  if (version < 4.0)
397  {
398  unsigned int num_nodes = 0;
399  in >> num_nodes;
400  mesh.reserve_nodes (num_nodes);
401 
402  // read in the nodal coordinates and form points.
403  Real x, y, z;
404  unsigned int id;
405 
406  // add the nodal coordinates to the mesh
407  for (unsigned int i=0; i<num_nodes; ++i)
408  {
409  in >> id >> x >> y >> z;
410  mesh.add_point (Point(x, y, z), i);
411  nodetrans[id] = i;
412  }
413  }
414  else
415  {
416  // Read numEntityBlocks line
417  std::size_t num_entities = 0, num_nodes = 0, min_node_tag, max_node_tag;
418  in >> num_entities >> num_nodes >> min_node_tag >> max_node_tag;
419 
420  mesh.reserve_nodes(num_nodes);
421 
422  std::size_t node_counter = 0;
423 
424  // Now loop over entities
425  for (std::size_t i = 0; i < num_entities; ++i)
426  {
427  int entity_dim, entity_tag, parametric;
428  std::size_t num_nodes_in_block = 0;
429  in >> entity_dim >> entity_tag >> parametric >> num_nodes_in_block;
430  libmesh_error_msg_if(parametric, "We don't currently support reading parametric gmsh entities");
431 
432  // Read the node tags/ids
433  std::size_t gmsh_id;
434  for (std::size_t n = 0; n < num_nodes_in_block; ++n)
435  {
436 
437  in >> gmsh_id;
438  nodetrans[gmsh_id] = node_counter++;
439  }
440 
441  // Read the node coordinates and add the nodes to the mesh
442  Real x, y, z;
443  for (std::size_t libmesh_id = node_counter - num_nodes_in_block;
444  libmesh_id < node_counter;
445  ++libmesh_id)
446  {
447  in >> x >> y >> z;
448  mesh.add_point(Point(x, y, z), libmesh_id);
449  }
450  }
451  }
452  // read the $ENDNOD delimiter
453  std::getline(in, s);
454  }
455 
456  // Read the element block
457  else if (s.find("$ELM") == static_cast<std::string::size_type>(0) ||
458  s.find("$Elements") == static_cast<std::string::size_type>(0))
459  {
460  // Keep track of element dimensions seen
461  std::vector<unsigned> elem_dimensions_seen(3);
462 
463  if (version < 4.0)
464  {
465  // For reading the number of elements and the node ids from the stream
466  unsigned int
467  num_elem = 0,
468  node_id = 0;
469 
470  // read how many elements are there, and reserve space in the mesh
471  in >> num_elem;
472  mesh.reserve_elem (num_elem);
473 
474  // As of version 2.2, the format for each element line is:
475  // elm-number elm-type number-of-tags < tag > ... node-number-list
476  // From the Gmsh docs:
477  // * the first tag is the number of the
478  // physical entity to which the element belongs
479  // * the second is the number of the elementary geometrical
480  // entity to which the element belongs
481  // * the third is the number of mesh partitions to which the element
482  // belongs
483  // * The rest of the tags are the partition ids (negative
484  // partition ids indicate ghost cells). A zero tag is
485  // equivalent to no tag. Gmsh and most codes using the
486  // MSH 2 format require at least the first two tags
487  // (physical and elementary tags).
488 
489  // read the elements
490  for (unsigned int iel=0; iel<num_elem; ++iel)
491  {
492  unsigned int
493  id, type,
494  physical=1, elementary=1,
495  nnodes=0, ntags;
496 
497  // Note: tag has to be an int because it could be negative,
498  // see above.
499  int tag;
500 
501  if (version <= 1.0)
502  in >> id >> type >> physical >> elementary >> nnodes;
503 
504  else
505  {
506  in >> id >> type >> ntags;
507 
508  if (ntags > 2)
509  libmesh_do_once(libMesh::err << "Warning, ntags=" << ntags << ", but we currently only support reading 2 flags." << std::endl;);
510 
511  for (unsigned int j = 0; j < ntags; j++)
512  {
513  in >> tag;
514  if (j == 0)
515  physical = tag;
516  else if (j == 1)
517  elementary = tag;
518  }
519  }
520 
521  // Get a reference to the ElementDefinition, throw an error if not found.
522  const GmshIO::ElementDefinition & eletype =
523  libmesh_map_find(_element_maps.in, type);
524 
525  // If we read nnodes, make sure it matches the number in eletype.nnodes
526  libmesh_error_msg_if(nnodes != 0 && nnodes != eletype.nnodes,
527  "nnodes = " << nnodes << " and eletype.nnodes = " << eletype.nnodes << " do not match.");
528 
529  // Assign the value from the eletype object.
530  nnodes = eletype.nnodes;
531 
532  // Don't add 0-dimensional "point" elements to the
533  // Mesh. They should *always* be treated as boundary
534  // "nodeset" data.
535  if (eletype.dim > 0)
536  {
537  // Record this element dimension as being "seen".
538  // We will treat all elements with dimension <
539  // max(dimension) as specifying boundary conditions,
540  // but we won't know what max_elem_dimension_seen is
541  // until we read the entire file.
542  elem_dimensions_seen[eletype.dim-1] = 1;
543 
544  // Add the element to the mesh
545  {
546  Elem * elem =
547  mesh.add_elem(Elem::build_with_id(eletype.type, iel));
548 
549  // Make sure that the libmesh element we added has nnodes nodes.
550  libmesh_error_msg_if(elem->n_nodes() != nnodes,
551  "Number of nodes for element "
552  << id
553  << " of type " << eletype.type
554  << " (Gmsh type " << type
555  << ") does not match Libmesh definition. "
556  << "I expected " << elem->n_nodes()
557  << " nodes, but got " << nnodes);
558 
559  // Add node pointers to the elements.
560  // If there is a node translation table, use it.
561  if (eletype.nodes.size() > 0)
562  for (unsigned int i=0; i<nnodes; i++)
563  {
564  in >> node_id;
565  elem->set_node(eletype.nodes[i], mesh.node_ptr(nodetrans[node_id]));
566  }
567  else
568  {
569  for (unsigned int i=0; i<nnodes; i++)
570  {
571  in >> node_id;
572  elem->set_node(i, mesh.node_ptr(nodetrans[node_id]));
573  }
574  }
575 
576  // Finally, set the subdomain ID to physical. If this is a lower-dimension element, this ID will
577  // eventually go into the Mesh's BoundaryInfo object.
578  elem->subdomain_id() = static_cast<subdomain_id_type>(physical);
579  }
580  }
581 
582  // Handle 0-dimensional elements (points) by adding
583  // them to the BoundaryInfo object with
584  // boundary_id == physical.
585  else
586  {
587  // This seems like it should always be the same
588  // number as the 'id' we already read in on this
589  // line. At least it was in the example gmsh
590  // file I had...
591  in >> node_id;
593  (nodetrans[node_id],
594  static_cast<boundary_id_type>(physical));
595  }
596  } // element loop
597  } // end if (version < 4.0)
598 
599  else
600  {
601  std::size_t num_entity_blocks = 0, num_elem = 0, min_element_tag, max_element_tag;
602 
603  // Read entity information
604  in >> num_entity_blocks >> num_elem >> min_element_tag >> max_element_tag;
605 
606  mesh.reserve_elem(num_elem);
607 
608  std::size_t iel = 0;
609 
610  // Loop over entity blocks
611  for (std::size_t i = 0; i < num_entity_blocks; ++i)
612  {
613  int entity_dim, entity_tag;
614  unsigned int element_type;
615  std::size_t num_elems_in_block = 0;
616  in >> entity_dim >> entity_tag >> element_type >> num_elems_in_block;
617 
618  // Get a reference to the ElementDefinition
619  const GmshIO::ElementDefinition & eletype =
620  libmesh_map_find(_element_maps.in, element_type);
621 
622  // Don't add 0-dimensional "point" elements to the
623  // Mesh. They should *always* be treated as boundary
624  // "nodeset" data.
625  if (eletype.dim > 0)
626  {
627  // Record this element dimension as being "seen".
628  // We will treat all elements with dimension <
629  // max(dimension) as specifying boundary conditions,
630  // but we won't know what max_elem_dimension_seen is
631  // until we read the entire file.
632  elem_dimensions_seen[eletype.dim-1] = 1;
633 
634  // Loop over elements with dim > 0
635  for (std::size_t n = 0; n < num_elems_in_block; ++n)
636  {
637  Elem * elem =
638  mesh.add_elem(Elem::build_with_id(eletype.type, iel++));
639 
640  std::size_t gmsh_element_id;
641  in >> gmsh_element_id;
642 
643  // Make sure this element isn't somewhere
644  // unexpected
645 
646  // A default bounding box is [inf,-inf] (empty);
647  // swap that and we get [-inf,inf] (everything)
648  BoundingBox expected_bounding_box;
649  std::swap(expected_bounding_box.min(),
650  expected_bounding_box.max());
651 
652  if (auto it = entity_to_bounding_box.find
653  (std::make_pair(entity_dim, entity_tag));
654  it != entity_to_bounding_box.end())
655  expected_bounding_box = it->second;
656 
657  // Get the remainder of the line that includes the nodes ids
658  std::getline(in, s);
659  std::istringstream is(s);
660  std::size_t local_node_counter = 0, gmsh_node_id;
661  while (is >> gmsh_node_id)
662  {
663  Node * node = mesh.node_ptr(nodetrans[gmsh_node_id]);
664 
665  // Make sure the file is consistent about entity
666  // placement. Well, mostly consistent. We have
667  // files that claim a bounding box but still
668  // have points epsilon outside it.
669  libmesh_error_msg_if
670  (!expected_bounding_box.contains_point
671  (*node, /* abs */ 0, /* relative */ TOLERANCE),
672  "$Elements dim " << entity_dim << " element "
673  << gmsh_element_id << " (entity " << entity_tag
674  << ", " <<
675  Utility::enum_to_string(eletype.type) <<
676  ") has node at " << *static_cast<Node*>(node)
677  << "\n outside entity physical bounding box " <<
678  expected_bounding_box);
679 
680  // Add node pointers to the elements.
681  // If there is a node translation table, use it.
682  if (eletype.nodes.size() > 0)
683  elem->set_node(eletype.nodes[local_node_counter++],
684  node);
685  else
686  elem->set_node(local_node_counter++, node);
687  }
688 
689  // Make sure that the libmesh element we added has nnodes nodes.
690  libmesh_error_msg_if(elem->n_nodes() != local_node_counter,
691  "Number of nodes for element "
692  << gmsh_element_id
693  << " of type " << eletype.type
694  << " (Gmsh type " << element_type
695  << ") does not match Libmesh definition. "
696  << "I expected " << elem->n_nodes()
697  << " nodes, but got " << local_node_counter);
698 
699  // Finally, set the subdomain ID to physical. If this is a lower-dimension element, this ID will
700  // eventually go into the Mesh's BoundaryInfo object.
701  elem->subdomain_id() = static_cast<subdomain_id_type>(
702  entity_to_physical_id[std::make_pair(entity_dim, entity_tag)]);
703 
704  } // end for (loop over elements in entity block)
705  } // end if (eletype.dim > 0)
706 
707  else
708  {
709  for (std::size_t n = 0; n < num_elems_in_block; ++n)
710  {
711  std::size_t gmsh_element_id, gmsh_node_id;
712  in >> gmsh_element_id;
713  in >> gmsh_node_id;
714 
715  // Make sure the file is consistent about entity
716  // placement.
717 
718  // A default bounding box is [inf,-inf] (empty);
719  // swap that and we get [-inf,inf] (everything)
720  BoundingBox expected_bounding_box;
721  std::swap(expected_bounding_box.min(),
722  expected_bounding_box.max());
723 
724  if (auto it = entity_to_bounding_box.find
725  (std::make_pair(entity_dim, entity_tag));
726  it != entity_to_bounding_box.end())
727  expected_bounding_box = it->second;
728 
729  Node * node = mesh.node_ptr(nodetrans[gmsh_node_id]);
730 
731  // We'll accept *mostly* consistent. We have
732  // files that claim a bounding box but still have
733  // points epsilon outside it.
734  libmesh_error_msg_if
735  (!expected_bounding_box.contains_point
736  (*node, /* abs */ 0, /* relative */ TOLERANCE),
737  "$Elements dim " << entity_dim << " element "
738  << gmsh_element_id << " (entity " << entity_tag <<
739  ") has node at " << *static_cast<Node*>(node)
740  << "\n outside entity physical bounding box " <<
741  expected_bounding_box);
742 
744  node,
745  static_cast<boundary_id_type>(entity_to_physical_id[
746  std::make_pair(entity_dim, entity_tag)]));
747  } // end for (loop over elements in entity block)
748  } // end if (eletype.dim == 0)
749  } // end for (loop over entity blocks)
750  } // end if (version >= 4.0)
751 
752  // read the $ENDELM delimiter
753  std::getline(in, s);
754 
755  // Record the max and min element dimension seen while reading the file.
756  unsigned char
757  max_elem_dimension_seen=1,
758  min_elem_dimension_seen=3;
759 
760  for (auto i : index_range(elem_dimensions_seen))
761  if (elem_dimensions_seen[i])
762  {
763  // Debugging
764  // libMesh::out << "Seen elements of dimension " << i+1 << std::endl;
765  max_elem_dimension_seen =
766  std::max(max_elem_dimension_seen, cast_int<unsigned char>(i+1));
767  min_elem_dimension_seen =
768  std::min(min_elem_dimension_seen, cast_int<unsigned char>(i+1));
769  }
770 
771  // Debugging:
772  // libMesh::out << "max_elem_dimension_seen=" << max_elem_dimension_seen << std::endl;
773  // libMesh::out << "min_elem_dimension_seen=" << min_elem_dimension_seen << std::endl;
774 
775 
776  // How many different element dimensions did we see while reading from file?
777  unsigned n_dims_seen = std::accumulate(elem_dimensions_seen.begin(),
778  elem_dimensions_seen.end(),
779  static_cast<unsigned>(0),
780  std::plus<unsigned>());
781 
782 
783  // Set mesh_dimension based on the largest element dimension seen.
784  mesh.set_mesh_dimension(max_elem_dimension_seen);
785 
786  // Now that we know the maximum element dimension seen,
787  // we know whether the physical names are subdomain
788  // names or sideset names.
789  for (const auto & pr : gmsh_physicals)
790  {
791  // Extract data
792  auto [phys_id, phys_dim] = pr.first;
793  const std::string & phys_name = pr.second;
794 
795  // If the physical's dimension matches the largest
796  // dimension we've seen, it's a subdomain name.
797  if (phys_dim == max_elem_dimension_seen)
798  mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
799 
800  // If it's zero-dimensional then it's a nodeset
801  else if (phys_dim == 0)
802  mesh.get_boundary_info().nodeset_name(cast_int<boundary_id_type>(phys_id)) = phys_name;
803 
804  // Otherwise, if it's not a lower-dimensional
805  // block, it's a sideset name.
806  else if (phys_dim < max_elem_dimension_seen &&
807  !lower_dimensional_blocks.count(cast_int<boundary_id_type>(phys_id)))
808  mesh.get_boundary_info().sideset_name(cast_int<boundary_id_type>(phys_id)) = phys_name;
809  }
810 
811  if (n_dims_seen > 1)
812  {
813  // Store lower-dimensional elements in a map sorted
814  // by Elem::key(). We use a multimap for two reasons:
815  // 1.) The hash function is not guaranteed to be
816  // unique, so different lower-dimensional elements
817  // could theoretically hash to the same value,
818  // although this is pretty unlikely.
819  // 2.) The Gmsh file may contain multiple
820  // lower-dimensional elements for a single side in
821  // order to implement multiple boundary ids for a
822  // single side. These lower-dimensional elements
823  // will all hash to the same value, and we need to
824  // be able to store all of them.
825  typedef std::unordered_multimap<dof_id_type, Elem *> provide_container_t;
826  provide_container_t provide_bcs;
827 
828  // 1st loop over active elements - get info about lower-dimensional elements.
829  for (auto & elem : mesh.active_element_ptr_range())
830  if (elem->dim() < max_elem_dimension_seen &&
831  !lower_dimensional_blocks.count(elem->subdomain_id()))
832  {
833  // To be consistent with the previous
834  // GmshIO behavior, add all the
835  // lower-dimensional elements' nodes to
836  // the Mesh's BoundaryInfo object with the
837  // lower-dimensional element's subdomain
838  // ID.
839  for (auto n : elem->node_index_range())
840  mesh.get_boundary_info().add_node(elem->node_id(n),
841  elem->subdomain_id());
842 
843  // Store this elem in a quickly-searchable
844  // container to use it to assign boundary
845  // conditions later.
846  provide_bcs.emplace(elem->key(), elem);
847  }
848 
849  // 2nd loop over active elements - use lower dimensional element data to set BCs for higher dimensional elements
850  for (auto & elem : mesh.active_element_ptr_range())
851  if (elem->dim() == max_elem_dimension_seen)
852  {
853  // This is a max-dimension element that
854  // may require BCs. For each of its
855  // sides, including internal sides, we'll
856  // see if one more more lower-dimensional elements
857  // provides boundary information for it.
858  // Note that we have not yet called
859  // find_neighbors(), so we can't use
860  // elem->neighbor(sn) in this algorithm...
861  for (auto sn : elem->side_index_range())
862  for (const auto & pr : as_range(provide_bcs.equal_range(elem->key(sn))))
863  {
864  // For each side side in the provide_bcs multimap...
865  // Construct the side for hash verification.
866  std::unique_ptr<Elem> side (elem->build_side_ptr(sn));
867 
868  // Construct the lower-dimensional element to compare to the side.
869  Elem * lower_dim_elem = pr.second;
870 
871  // This was a hash, so it might not be perfect. Let's verify...
872  if (*lower_dim_elem == *side)
873  {
874  // Add the lower-dimensional
875  // element's subdomain_id as a
876  // boundary_id for the
877  // higher-dimensional element.
878  boundary_id_type bid = cast_int<boundary_id_type>(lower_dim_elem->subdomain_id());
879  mesh.get_boundary_info().add_side(elem, sn, bid);
880  }
881  }
882  }
883 
884  // 3rd loop over active elements - Remove the lower-dimensional elements
885  for (auto & elem : mesh.active_element_ptr_range())
886  if (elem->dim() < max_elem_dimension_seen &&
887  !lower_dimensional_blocks.count(elem->subdomain_id()))
888  mesh.delete_elem(elem);
889  } // end if (n_dims_seen > 1)
890  } // if $ELM
891 
892  continue;
893  } // if (in)
894 
895 
896  // If !in, check to see if EOF was set. If so, break out
897  // of while loop.
898  if (in.eof())
899  break;
900 
901  // If !in and !in.eof(), stream is in a bad state!
902  libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
903 
904  } // while true
905 }
OStreamProxy err
const MeshBase & mesh() const
Definition: mesh_output.h:259
virtual void reserve_nodes(const dof_id_type nn)=0
Reserves space for a known number of nodes.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
std::string & nodeset_name(boundary_id_type id)
static constexpr Real TOLERANCE
TestClass subdomain_id_type
Based on the 4-byte comment warning above, this probably doesn&#39;t work with exodusII at all...
Definition: id_types.h:43
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
std::map< unsigned int, ElementDefinition > in
Definition: gmsh_io.h:193
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
int8_t boundary_id_type
Definition: id_types.h:51
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
libmesh_assert(ctx)
PetscErrorCode PetscInt const PetscInt IS * is
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:1692
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
std::string & sideset_name(boundary_id_type id)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem&#39;s ...
Definition: elem.C:558
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
Definition: gmsh_io.h:200
virtual const Node * node_ptr(const dof_id_type i) const =0
virtual void reserve_elem(const dof_id_type ne)=0
Reserves space for a known number of elements.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ set_n_partitions()

void libMesh::MeshInput< MeshBase >::set_n_partitions ( unsigned int  n_parts)
inlineprotectedinherited

Sets the number of partitions in the mesh.

Typically this gets done by the partitioner, but some parallel file formats begin "pre-partitioned".

Definition at line 101 of file mesh_input.h.

References libMesh::MeshInput< MT >::mesh().

Referenced by libMesh::Nemesis_IO::read(), and libMesh::XdrIO::read_header().

101 { this->mesh().set_n_partitions() = n_parts; }
unsigned int & set_n_partitions()
Definition: mesh_base.h:1865

◆ skip_comment_lines()

void libMesh::MeshInput< MeshBase >::skip_comment_lines ( std::istream &  in,
const char  comment_start 
)
protectedinherited

Reads input from in, skipping all the lines that start with the character comment_start.

Definition at line 187 of file mesh_input.h.

Referenced by libMesh::TetGenIO::read(), and libMesh::UCDIO::read_implementation().

189 {
190  char c, line[256];
191 
192  while (in.get(c), c==comment_start)
193  in.getline (line, 255);
194 
195  // put back first character of
196  // first non-comment line
197  in.putback (c);
198 }

◆ write()

void libMesh::GmshIO::write ( const std::string &  name)
overridevirtual

This method implements writing a mesh to a specified file in the Gmsh *.msh format.

Implements libMesh::MeshOutput< MeshBase >.

Definition at line 909 of file gmsh_io.C.

References libMesh::Quality::name(), and write_mesh().

Referenced by libMesh::NameBasedIO::write().

910 {
911  if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
912  {
913  // Open the output file stream
914  std::ofstream out_stream (name.c_str());
915 
916  // Make sure it opened correctly
917  if (!out_stream.good())
918  libmesh_file_error(name.c_str());
919 
920  this->write_mesh (out_stream);
921  }
922 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
const MeshBase & mesh() const
Definition: mesh_output.h:259
void write_mesh(std::ostream &out)
This method implements writing a mesh to a specified file.
Definition: gmsh_io.C:938

◆ write_discontinuous_equation_systems()

void libMesh::MeshOutput< MeshBase >::write_discontinuous_equation_systems ( const std::string &  fname,
const EquationSystems es,
const std::set< std::string > *  system_names = nullptr 
)
virtualinherited

This method implements writing a mesh with discontinuous data to a specified file where the data is taken from the EquationSystems object.

Definition at line 89 of file mesh_output.C.

References libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_variable_names(), libMesh::EquationSystems::get_mesh(), libMesh::libmesh_assert(), and libMesh::out.

Referenced by libMesh::ExodusII_IO::write_timestep_discontinuous().

92 {
93  LOG_SCOPE("write_discontinuous_equation_systems()", "MeshOutput");
94 
95  // We may need to gather and/or renumber a DistributedMesh to output
96  // it, making that const qualifier in our constructor a dirty lie
97  MT & my_mesh = const_cast<MT &>(*_obj);
98 
99  // If we're asked to write data that's associated with a different
100  // mesh, output files full of garbage are the result.
101  libmesh_assert_equal_to(&es.get_mesh(), _obj);
102 
103  // A non-renumbered mesh may not have a contiguous numbering, and
104  // that needs to be fixed before we can build a solution vector.
105  if (my_mesh.max_elem_id() != my_mesh.n_elem() ||
106  my_mesh.max_node_id() != my_mesh.n_nodes())
107  {
108  // If we were allowed to renumber then we should have already
109  // been properly renumbered...
110  libmesh_assert(!my_mesh.allow_renumbering());
111 
112  libmesh_do_once(libMesh::out <<
113  "Warning: This MeshOutput subclass only supports meshes which are contiguously renumbered!"
114  << std::endl;);
115 
116  my_mesh.allow_renumbering(true);
117 
118  my_mesh.renumber_nodes_and_elements();
119 
120  // Not sure what good going back to false will do here, the
121  // renumbering horses have already left the barn...
122  my_mesh.allow_renumbering(false);
123  }
124 
125  MeshSerializer serialize(const_cast<MT &>(*_obj), !_is_parallel_format, _serial_only_needed_on_proc_0);
126 
127  // Build the list of variable names that will be written.
128  std::vector<std::string> names;
129  es.build_variable_names (names, nullptr, system_names);
130 
131  if (!_is_parallel_format)
132  {
133  // Build the nodal solution values & get the variable
134  // names from the EquationSystems object
135  std::vector<Number> soln;
136  es.build_discontinuous_solution_vector (soln, system_names,
137  nullptr, false, /* defaults */
138  this->get_add_sides());
139 
140  this->write_nodal_data_discontinuous (fname, soln, names);
141  }
142  else // _is_parallel_format
143  {
144  libmesh_not_implemented();
145  }
146 }
virtual void write_nodal_data_discontinuous(const std::string &, const std::vector< Number > &, const std::vector< std::string > &)
This method implements writing a mesh with discontinuous data to a specified file where the nodal dat...
Definition: mesh_output.h:118
const MeshBase *const _obj
A pointer to a constant object.
Definition: mesh_output.h:202
const bool _is_parallel_format
Flag specifying whether this format is parallel-capable.
Definition: mesh_output.h:184
libmesh_assert(ctx)
OStreamProxy out
const bool _serial_only_needed_on_proc_0
Flag specifying whether this format can be written by only serializing the mesh to processor zero...
Definition: mesh_output.h:193

◆ write_equation_systems()

void libMesh::MeshOutput< MeshBase >::write_equation_systems ( const std::string &  fname,
const EquationSystems es,
const std::set< std::string > *  system_names = nullptr 
)
virtualinherited

This method implements writing a mesh with data to a specified file where the data is taken from the EquationSystems object.

Reimplemented in libMesh::ExodusII_IO, and libMesh::NameBasedIO.

Definition at line 31 of file mesh_output.C.

References libMesh::EquationSystems::build_solution_vector(), libMesh::EquationSystems::build_variable_names(), libMesh::EquationSystems::get_mesh(), libMesh::libmesh_assert(), and libMesh::out.

Referenced by libMesh::Nemesis_IO::write_timestep().

34 {
35  LOG_SCOPE("write_equation_systems()", "MeshOutput");
36 
37  // We may need to gather and/or renumber a DistributedMesh to output
38  // it, making that const qualifier in our constructor a dirty lie
39  MT & my_mesh = const_cast<MT &>(*_obj);
40 
41  // If we're asked to write data that's associated with a different
42  // mesh, output files full of garbage are the result.
43  libmesh_assert_equal_to(&es.get_mesh(), _obj);
44 
45  // A non-parallel format, non-renumbered mesh may not have a contiguous
46  // numbering, and that needs to be fixed before we can build a solution vector.
47  if (!_is_parallel_format &&
48  (my_mesh.max_elem_id() != my_mesh.n_elem() ||
49  my_mesh.max_node_id() != my_mesh.n_nodes()))
50  {
51  // If we were allowed to renumber then we should have already
52  // been properly renumbered...
53  libmesh_assert(!my_mesh.allow_renumbering());
54 
55  libmesh_do_once(libMesh::out <<
56  "Warning: This MeshOutput subclass only supports meshes which are contiguously renumbered!"
57  << std::endl;);
58 
59  my_mesh.allow_renumbering(true);
60 
61  my_mesh.renumber_nodes_and_elements();
62 
63  // Not sure what good going back to false will do here, the
64  // renumbering horses have already left the barn...
65  my_mesh.allow_renumbering(false);
66  }
67 
69  {
70  MeshSerializer serialize(const_cast<MT &>(*_obj), !_is_parallel_format, _serial_only_needed_on_proc_0);
71 
72  // Build the list of variable names that will be written.
73  std::vector<std::string> names;
74  es.build_variable_names (names, nullptr, system_names);
75 
76  // Build the nodal solution values & get the variable
77  // names from the EquationSystems object
78  std::vector<Number> soln;
79  es.build_solution_vector (soln, system_names,
80  this->get_add_sides());
81 
82  this->write_nodal_data (fname, soln, names);
83  }
84  else // _is_parallel_format
85  this->write_nodal_data (fname, es, system_names);
86 }
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &)
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: mesh_output.h:109
const MeshBase *const _obj
A pointer to a constant object.
Definition: mesh_output.h:202
const bool _is_parallel_format
Flag specifying whether this format is parallel-capable.
Definition: mesh_output.h:184
libmesh_assert(ctx)
OStreamProxy out
const bool _serial_only_needed_on_proc_0
Flag specifying whether this format can be written by only serializing the mesh to processor zero...
Definition: mesh_output.h:193

◆ write_lower_dimensional_elements()

bool & libMesh::GmshIO::write_lower_dimensional_elements ( )

Access to the flag which controls whether boundary elements are written to the Mesh file.

Definition at line 142 of file gmsh_io.C.

References _write_lower_dimensional_elements.

Referenced by write_mesh().

143 {
145 }
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:155

◆ write_mesh()

void libMesh::GmshIO::write_mesh ( std::ostream &  out)
private

This method implements writing a mesh to a specified file.

This will write an ASCII *.msh file.

Definition at line 938 of file gmsh_io.C.

References _element_maps, libMesh::BoundaryInfo::build_side_list(), libMesh::Elem::build_side_ptr(), libMesh::MeshBase::elem_ref(), libMesh::MeshBase::get_boundary_info(), libMesh::GmshIO::ElementDefinition::gmsh_type, libMesh::DofObject::id(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_elem_id(), libMesh::MeshInput< MeshBase >::mesh(), libMesh::MeshOutput< MT >::mesh(), libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ref(), libMesh::GmshIO::ElementDefinition::nodes, libMesh::GmshIO::ElementMaps::out, libMesh::DofObject::processor_id(), libMesh::Real, and write_lower_dimensional_elements().

Referenced by write().

939 {
940  // Be sure that the stream is valid.
941  libmesh_assert (out_stream.good());
942 
943  // Get a const reference to the mesh
944  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
945 
946  // If requested, write out lower-dimensional elements for
947  // element-side-based boundary conditions.
948  std::size_t n_boundary_faces = 0;
950  n_boundary_faces = mesh.get_boundary_info().n_boundary_conds();
951 
952  // Note: we are using version 2.0 of the gmsh output format.
953 
954  // Write the file header.
955  out_stream << "$MeshFormat\n";
956  out_stream << "2.0 0 " << sizeof(Real) << '\n';
957  out_stream << "$EndMeshFormat\n";
958 
959  // write the nodes in (n x y z) format
960  out_stream << "$Nodes\n";
961  out_stream << mesh.n_nodes() << '\n';
962 
963  for (auto v : make_range(mesh.n_nodes()))
964  out_stream << mesh.node_ref(v).id()+1 << " "
965  << mesh.node_ref(v)(0) << " "
966  << mesh.node_ref(v)(1) << " "
967  << mesh.node_ref(v)(2) << '\n';
968  out_stream << "$EndNodes\n";
969 
970  {
971  // write the connectivity
972  out_stream << "$Elements\n";
973  out_stream << mesh.n_active_elem() + n_boundary_faces << '\n';
974 
975  // loop over the elements
976  for (const auto & elem : mesh.active_element_ptr_range())
977  {
978  // Get a reference to the ElementDefinition object
979  const ElementDefinition & eletype =
980  libmesh_map_find(_element_maps.out, elem->type());
981 
982  // The element mapper better not require any more nodes
983  // than are present in the current element!
984  libmesh_assert_less_equal (eletype.nodes.size(), elem->n_nodes());
985 
986  // elements ids are 1 based in Gmsh
987  out_stream << elem->id()+1 << " ";
988  // element type
989  out_stream << eletype.gmsh_type;
990 
991  // write the number of tags (3) and their values:
992  // 1 (physical entity)
993  // 2 (geometric entity)
994  // 3 (partition entity)
995  out_stream << " 3 "
996  << static_cast<unsigned int>(elem->subdomain_id())
997  << " 0 "
998  << elem->processor_id()+1
999  << " ";
1000 
1001  // if there is a node translation table, use it
1002  if (eletype.nodes.size() > 0)
1003  for (auto i : elem->node_index_range())
1004  out_stream << elem->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
1005  // otherwise keep the same node order
1006  else
1007  for (auto i : elem->node_index_range())
1008  out_stream << elem->node_id(i)+1 << " "; // gmsh is 1-based
1009  out_stream << "\n";
1010  } // element loop
1011  }
1012 
1013  {
1014  // A counter for writing surface elements to the Gmsh file
1015  // sequentially. We start numbering them with a number strictly
1016  // larger than the largest element ID in the mesh. Note: the
1017  // MeshBase docs say "greater than or equal to" the maximum
1018  // element id in the mesh, so technically we might need a +1 here,
1019  // but all of the implementations return an ID strictly greater
1020  // than the largest element ID in the Mesh.
1021  unsigned int e_id = mesh.max_elem_id();
1022 
1023  // loop over the elements, writing out boundary faces
1024  if (n_boundary_faces)
1025  {
1026  // Build a list of (elem, side, bc) tuples.
1027  auto bc_triples = mesh.get_boundary_info().build_side_list();
1028 
1029  // Loop over these lists, writing data to the file.
1030  for (const auto & t : bc_triples)
1031  {
1032  const Elem & elem = mesh.elem_ref(std::get<0>(t));
1033 
1034  std::unique_ptr<const Elem> side = elem.build_side_ptr(std::get<1>(t));
1035 
1036  // consult the export element table
1037  const GmshIO::ElementDefinition & eletype =
1038  libmesh_map_find(_element_maps.out, side->type());
1039 
1040  // The element mapper better not require any more nodes
1041  // than are present in the current element!
1042  libmesh_assert_less_equal (eletype.nodes.size(), side->n_nodes());
1043 
1044  // elements ids are 1-based in Gmsh
1045  out_stream << e_id+1 << " ";
1046 
1047  // element type
1048  out_stream << eletype.gmsh_type;
1049 
1050  // write the number of tags:
1051  // 1 (physical entity)
1052  // 2 (geometric entity)
1053  // 3 (partition entity)
1054  out_stream << " 3 "
1055  << std::get<2>(t)
1056  << " 0 "
1057  << elem.processor_id()+1
1058  << " ";
1059 
1060  // if there is a node translation table, use it
1061  if (eletype.nodes.size() > 0)
1062  for (auto i : side->node_index_range())
1063  out_stream << side->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
1064 
1065  // otherwise keep the same node order
1066  else
1067  for (auto i : side->node_index_range())
1068  out_stream << side->node_id(i)+1 << " "; // gmsh is 1-based
1069 
1070  // Go to the next line
1071  out_stream << "\n";
1072 
1073  // increment this index too...
1074  ++e_id;
1075  }
1076  }
1077 
1078  out_stream << "$EndElements\n";
1079  }
1080 }
const MT & mesh() const
Definition: mesh_output.h:259
std::size_t n_boundary_conds() const
virtual dof_id_type n_active_elem() const =0
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
bool & write_lower_dimensional_elements()
Access to the flag which controls whether boundary elements are written to the Mesh file...
Definition: gmsh_io.C:142
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
std::map< ElemType, ElementDefinition > out
Definition: gmsh_io.h:192
dof_id_type id() const
Definition: dof_object.h:828
virtual dof_id_type max_elem_id() const =0
libmesh_assert(ctx)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
static ElementMaps _element_maps
A static ElementMaps object that is built statically and used by all instances of this class...
Definition: gmsh_io.h:200
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
virtual dof_id_type n_nodes() const =0

◆ write_nodal_data() [1/3]

void libMesh::GmshIO::write_nodal_data ( const std::string &  fname,
const std::vector< Number > &  soln,
const std::vector< std::string > &  names 
)
overridevirtual

This method implements writing a mesh with nodal data to a specified file where the nodal data and variable names are provided.

Reimplemented from libMesh::MeshOutput< MeshBase >.

Definition at line 926 of file gmsh_io.C.

References write_post().

Referenced by libMesh::NameBasedIO::write_nodal_data().

929 {
930  LOG_SCOPE("write_nodal_data()", "GmshIO");
931 
932  if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
933  this->write_post (fname, &soln, &names);
934 }
const MeshBase & mesh() const
Definition: mesh_output.h:259
void write_post(const std::string &, const std::vector< Number > *=nullptr, const std::vector< std::string > *=nullptr)
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: gmsh_io.C:1084

◆ write_nodal_data() [2/3]

void libMesh::MeshOutput< MeshBase >::write_nodal_data ( const std::string &  fname,
const NumericVector< Number > &  parallel_soln,
const std::vector< std::string > &  names 
)
virtualinherited

This method may be overridden by "parallel" output formats for writing nodal data.

Instead of getting a localized copy of the nodal solution vector, it is passed a NumericVector of type=PARALLEL which is in node-major order i.e. (u0,v0,w0, u1,v1,w1, u2,v2,w2, u3,v3,w3, ...) and contains n_nodes*n_vars total entries. Then, it is up to the individual I/O class to extract the required solution values from this vector and write them in parallel.

If not implemented, localizes the parallel vector into a std::vector and calls the other version of this function.

Reimplemented in libMesh::Nemesis_IO.

Definition at line 149 of file mesh_output.C.

References libMesh::NumericVector< T >::localize().

152 {
153  // This is the fallback implementation for parallel I/O formats that
154  // do not yet implement proper writing in parallel, and instead rely
155  // on the full solution vector being available on all processors.
156  std::vector<Number> soln;
157  parallel_soln.localize(soln);
158  this->write_nodal_data(fname, soln, names);
159 }
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &)
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: mesh_output.h:109
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.

◆ write_nodal_data() [3/3]

void libMesh::MeshOutput< MeshBase >::write_nodal_data ( const std::string &  fname,
const EquationSystems es,
const std::set< std::string > *  system_names 
)
virtualinherited

This method should be overridden by "parallel" output formats for writing nodal data.

Instead of getting a localized copy of the nodal solution vector, it directly uses EquationSystems current_local_solution vectors to look up nodal values.

If not implemented, reorders the solutions into a nodal-only NumericVector and calls the above version of this function.

Reimplemented in libMesh::Nemesis_IO.

Definition at line 162 of file mesh_output.C.

References libMesh::EquationSystems::build_parallel_solution_vector(), and libMesh::EquationSystems::build_variable_names().

165 {
166  std::vector<std::string> names;
167  es.build_variable_names (names, nullptr, system_names);
168 
169  std::unique_ptr<NumericVector<Number>> parallel_soln =
170  es.build_parallel_solution_vector(system_names);
171 
172  this->write_nodal_data (fname, *parallel_soln, names);
173 }
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &)
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: mesh_output.h:109

◆ write_nodal_data_discontinuous()

virtual void libMesh::MeshOutput< MeshBase >::write_nodal_data_discontinuous ( const std::string &  ,
const std::vector< Number > &  ,
const std::vector< std::string > &   
)
inlinevirtualinherited

This method implements writing a mesh with discontinuous data to a specified file where the nodal data and variables names are provided.

Reimplemented in libMesh::ExodusII_IO.

Definition at line 118 of file mesh_output.h.

121  { libmesh_not_implemented(); }

◆ write_post()

void libMesh::GmshIO::write_post ( const std::string &  fname,
const std::vector< Number > *  v = nullptr,
const std::vector< std::string > *  solution_names = nullptr 
)
private

This method implements writing a mesh with nodal data to a specified file where the nodal data and variable names are optionally provided.

This will write an ASCII or binary *.pos file, depending on the binary flag.

Definition at line 1084 of file gmsh_io.C.

References binary(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::Utility::enum_to_string(), libMesh::err, libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::libmesh_real(), libMesh::MeshInput< MeshBase >::mesh(), libMesh::MeshOutput< MT >::mesh(), libMesh::MeshBase::n_nodes(), n_vars, libMesh::out, libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::TET10, libMesh::TET4, libMesh::TRI3, and libMesh::TRI6.

Referenced by write_nodal_data().

1087 {
1088 
1089  // Should only do this on processor 0!
1090  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
1091 
1092  // Create an output stream
1093  std::ofstream out_stream(fname.c_str());
1094 
1095  // Make sure it opened correctly
1096  if (!out_stream.good())
1097  libmesh_file_error(fname.c_str());
1098 
1099  // create a character buffer
1100  char buf[80];
1101 
1102  // Get a constant reference to the mesh.
1103  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1104 
1105  // write the data
1106  if ((solution_names != nullptr) && (v != nullptr))
1107  {
1108  const unsigned int n_vars =
1109  cast_int<unsigned int>(solution_names->size());
1110 
1111  if (!(v->size() == mesh.n_nodes()*n_vars))
1112  libMesh::err << "ERROR: v->size()=" << v->size()
1113  << ", mesh.n_nodes()=" << mesh.n_nodes()
1114  << ", n_vars=" << n_vars
1115  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
1116  << "\n";
1117 
1118  libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
1119 
1120  // write the header
1121  out_stream << "$PostFormat\n";
1122  if (this->binary())
1123  out_stream << "1.2 1 " << sizeof(double) << "\n";
1124  else
1125  out_stream << "1.2 0 " << sizeof(double) << "\n";
1126  out_stream << "$EndPostFormat\n";
1127 
1128  // Loop over the elements to see how much of each type there are
1129  unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
1130  n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
1131  unsigned int n_scalar=0, n_vector=0, n_tensor=0;
1132  unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;
1133 
1134  {
1135  for (const auto & elem : mesh.active_element_ptr_range())
1136  {
1137  const ElemType elemtype = elem->type();
1138 
1139  switch (elemtype)
1140  {
1141  case EDGE2:
1142  case EDGE3:
1143  case EDGE4:
1144  {
1145  n_lines += 1;
1146  break;
1147  }
1148  case TRI3:
1149  case TRI6:
1150  {
1151  n_triangles += 1;
1152  break;
1153  }
1154  case QUAD4:
1155  case QUAD8:
1156  case QUAD9:
1157  {
1158  n_quadrangles += 1;
1159  break;
1160  }
1161  case TET4:
1162  case TET10:
1163  {
1164  n_tetrahedra += 1;
1165  break;
1166  }
1167  case HEX8:
1168  case HEX20:
1169  case HEX27:
1170  {
1171  n_hexahedra += 1;
1172  break;
1173  }
1174  case PRISM6:
1175  case PRISM15:
1176  case PRISM18:
1177  {
1178  n_prisms += 1;
1179  break;
1180  }
1181  case PYRAMID5:
1182  {
1183  n_pyramids += 1;
1184  break;
1185  }
1186  default:
1187  libmesh_error_msg("ERROR: Nonexistent element type " << Utility::enum_to_string(elem->type()));
1188  }
1189  }
1190  }
1191 
1192  // create a view for each variable
1193  for (unsigned int ivar=0; ivar < n_vars; ivar++)
1194  {
1195  std::string varname = (*solution_names)[ivar];
1196 
1197  // at the moment, we just write out scalar quantities
1198  // later this should be made configurable through
1199  // options to the writer class
1200  n_scalar = 1;
1201 
1202  // write the variable as a view, and the number of time steps
1203  out_stream << "$View\n" << varname << " " << 1 << "\n";
1204 
1205  // write how many of each geometry type are written
1206  out_stream << n_points * n_scalar << " "
1207  << n_points * n_vector << " "
1208  << n_points * n_tensor << " "
1209  << n_lines * n_scalar << " "
1210  << n_lines * n_vector << " "
1211  << n_lines * n_tensor << " "
1212  << n_triangles * n_scalar << " "
1213  << n_triangles * n_vector << " "
1214  << n_triangles * n_tensor << " "
1215  << n_quadrangles * n_scalar << " "
1216  << n_quadrangles * n_vector << " "
1217  << n_quadrangles * n_tensor << " "
1218  << n_tetrahedra * n_scalar << " "
1219  << n_tetrahedra * n_vector << " "
1220  << n_tetrahedra * n_tensor << " "
1221  << n_hexahedra * n_scalar << " "
1222  << n_hexahedra * n_vector << " "
1223  << n_hexahedra * n_tensor << " "
1224  << n_prisms * n_scalar << " "
1225  << n_prisms * n_vector << " "
1226  << n_prisms * n_tensor << " "
1227  << n_pyramids * n_scalar << " "
1228  << n_pyramids * n_vector << " "
1229  << n_pyramids * n_tensor << " "
1230  << nb_text2d << " "
1231  << nb_text2d_chars << " "
1232  << nb_text3d << " "
1233  << nb_text3d_chars << "\n";
1234 
1235  // if binary, write a marker to identify the endianness of the file
1236  if (this->binary())
1237  {
1238  const int one = 1;
1239  std::memcpy(buf, &one, sizeof(int));
1240  out_stream.write(buf, sizeof(int));
1241  }
1242 
1243  // the time steps (there is just 1 at the moment)
1244  if (this->binary())
1245  {
1246  double one = 1;
1247  std::memcpy(buf, &one, sizeof(double));
1248  out_stream.write(buf, sizeof(double));
1249  }
1250  else
1251  out_stream << "1\n";
1252 
1253  // Loop over the elements and write out the data
1254  for (const auto & elem : mesh.active_element_ptr_range())
1255  {
1256  const unsigned int nv = elem->n_vertices();
1257  // this is quite crappy, but I did not invent that file format!
1258  for (unsigned int d=0; d<3; d++) // loop over the dimensions
1259  {
1260  for (unsigned int n=0; n < nv; n++) // loop over vertices
1261  {
1262  const Point & vertex = elem->point(n);
1263  if (this->binary())
1264  {
1265 #if defined(LIBMESH_DEFAULT_TRIPLE_PRECISION) || defined(LIBMESH_DEFAULT_QUADRUPLE_PRECISION)
1266  libmesh_warning("Gmsh binary writes use only double precision!");
1267 #endif
1268  double tmp = double(vertex(d));
1269  std::memcpy(buf, &tmp, sizeof(double));
1270  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1271  }
1272  else
1273  out_stream << vertex(d) << " ";
1274  }
1275  if (!this->binary())
1276  out_stream << "\n";
1277  }
1278 
1279  // now finally write out the data
1280  for (unsigned int i=0; i < nv; i++) // loop over vertices
1281  if (this->binary())
1282  {
1283 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1284  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1285  << "complex numbers. Will only write the real part of "
1286  << "variable " << varname << std::endl;
1287 #endif
1288  double tmp = double(libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]));
1289  std::memcpy(buf, &tmp, sizeof(double));
1290  out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1291  }
1292  else
1293  {
1294 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1295  libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1296  << "complex numbers. Will only write the real part of "
1297  << "variable " << varname << std::endl;
1298 #endif
1299  out_stream << libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]) << "\n";
1300  }
1301  }
1302  if (this->binary())
1303  out_stream << "\n";
1304  out_stream << "$EndView\n";
1305 
1306  } // end variable loop (writing the views)
1307  }
1308 }
T libmesh_real(T a)
OStreamProxy err
const MeshBase & mesh() const
Definition: mesh_output.h:259
ElemType
Defines an enum for geometric element types.
bool & binary()
Flag indicating whether or not to write a binary file.
Definition: gmsh_io.C:135
unsigned int n_vars
std::string enum_to_string(const T e)
OStreamProxy out
virtual dof_id_type n_nodes() const =0

Member Data Documentation

◆ _binary

bool libMesh::GmshIO::_binary
private

Flag to write binary data.

Definition at line 149 of file gmsh_io.h.

Referenced by binary().

◆ _element_maps

GmshIO::ElementMaps libMesh::GmshIO::_element_maps = GmshIO::build_element_maps()
staticprivate

A static ElementMaps object that is built statically and used by all instances of this class.

Definition at line 200 of file gmsh_io.h.

Referenced by read_mesh(), and write_mesh().

◆ _is_parallel_format

const bool libMesh::MeshOutput< MeshBase >::_is_parallel_format
protectedinherited

Flag specifying whether this format is parallel-capable.

If this is false (default) I/O is only permitted when the mesh has been serialized.

Definition at line 184 of file mesh_output.h.

Referenced by libMesh::FroIO::write(), libMesh::PostscriptIO::write(), and libMesh::EnsightIO::write().

◆ _serial_only_needed_on_proc_0

const bool libMesh::MeshOutput< MeshBase >::_serial_only_needed_on_proc_0
protectedinherited

Flag specifying whether this format can be written by only serializing the mesh to processor zero.

If this is false (default) the mesh will be serialized to all processors

Definition at line 193 of file mesh_output.h.

◆ _write_lower_dimensional_elements

bool libMesh::GmshIO::_write_lower_dimensional_elements
private

If true, lower-dimensional elements based on the boundary conditions get written to the output file.

Definition at line 155 of file gmsh_io.h.

Referenced by write_lower_dimensional_elements().

◆ elems_of_dimension

std::vector<bool> libMesh::MeshInput< MeshBase >::elems_of_dimension
protectedinherited

The documentation for this class was generated from the following files: