OGS
GmshReader.cpp
Go to the documentation of this file.
1 
10 #include "GmshReader.h"
11 
12 #include <algorithm>
13 #include <array>
14 #include <fstream>
15 #include <map>
16 #include <type_traits>
17 #include <vector>
18 
19 #include "BaseLib/FileTools.h"
21 #include "MeshLib/Elements/Hex.h"
22 #include "MeshLib/Elements/Line.h"
23 #include "MeshLib/Elements/Prism.h"
25 #include "MeshLib/Elements/Quad.h"
26 #include "MeshLib/Elements/Tet.h"
27 #include "MeshLib/Elements/Tri.h"
28 #include "MeshLib/Mesh.h"
30 #include "MeshLib/Node.h"
31 
32 namespace FileIO
33 {
34 namespace GMSH
35 {
36 bool isGMSHMeshFile(const std::string& fname)
37 {
38  std::ifstream input(fname.c_str());
39 
40  if (!input)
41  {
42  ERR("isGMSHMeshFile(): Could not open file {:s}.", fname);
43  return false;
44  }
45 
46  std::string header_first_line;
47  input >> header_first_line;
48  if (header_first_line.find("$MeshFormat") != std::string::npos)
49  {
50  // read version
51  std::string version;
52  getline(input, version);
53  getline(input, version);
54  INFO("isGMSHMeshFile(): Found GMSH mesh file version: {:s}.", version);
55  input.close();
56  return true;
57  }
58 
59  return false;
60 }
61 
62 void readNodeIDs(std::ifstream& in, unsigned n_nodes,
63  std::vector<unsigned>& node_ids,
64  std::map<unsigned, unsigned> const& id_map)
65 {
66  unsigned idx;
67  for (unsigned i = 0; i < n_nodes; i++)
68  {
69  in >> idx;
70  node_ids.push_back(id_map.at(idx));
71  }
72 }
73 
74 template <typename ElementType>
75 std::pair<MeshLib::Element*, int> createElement(
76  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
77  int const mat_id, std::map<unsigned, unsigned> const& id_map)
78 {
79  std::vector<unsigned> node_ids;
80  readNodeIDs(in, ElementType::n_all_nodes, node_ids, id_map);
81 
82  std::array<MeshLib::Node*, ElementType::n_all_nodes> element_nodes;
83 
84  std::transform(begin(node_ids), end(node_ids), begin(element_nodes),
85  [&nodes](auto const id) { return nodes[id]; });
86 
87  return std::make_pair(new ElementType(element_nodes), mat_id);
88 }
89 
90 template <>
91 std::pair<MeshLib::Element*, int> createElement<MeshLib::Tri>(
92  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
93  int const mat_id, std::map<unsigned, unsigned> const& id_map)
94 {
95  std::vector<unsigned> node_ids;
96  readNodeIDs(in, 3, node_ids, id_map);
97 
98  std::array<MeshLib::Node*, 3> element_nodes;
99 
100  std::transform(std::rbegin(node_ids), std::rend(node_ids),
101  begin(element_nodes),
102  [&nodes](auto const id) { return nodes[id]; });
103 
104  return std::make_pair(new MeshLib::Tri(element_nodes), mat_id);
105 }
106 
107 template <>
108 std::pair<MeshLib::Element*, int> createElement<MeshLib::Tet10>(
109  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
110  int const mat_id, std::map<unsigned, unsigned> const& id_map)
111 {
112  std::vector<unsigned> node_ids;
113  readNodeIDs(in, MeshLib::Tet10::n_all_nodes, node_ids, id_map);
114 
115  std::swap(node_ids[8], node_ids[9]);
116 
117  std::array<MeshLib::Node*, MeshLib::Tet10::n_all_nodes> element_nodes;
118 
119  std::transform(begin(node_ids), end(node_ids), begin(element_nodes),
120  [&nodes](auto const id) { return nodes[id]; });
121 
122  return std::make_pair(new MeshLib::Tet10(element_nodes), mat_id);
123 }
124 
125 template <>
126 std::pair<MeshLib::Element*, int> createElement<MeshLib::Hex20>(
127  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
128  int const mat_id, std::map<unsigned, unsigned> const& id_map)
129 {
130  std::vector<unsigned> node_ids;
131  readNodeIDs(in, MeshLib::Hex20::n_all_nodes, node_ids, id_map);
132 
133  std::array<MeshLib::Node*, MeshLib::Hex20::n_all_nodes> element_nodes;
134 
135  constexpr std::array node_order = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11,
136  13, 9, 16, 18, 19, 17, 10, 12, 14, 15};
137 
138  std::transform(begin(node_order), end(node_order), begin(element_nodes),
139  [&node_ids, &nodes](auto const id)
140  { return nodes[node_ids[id]]; });
141 
142  return std::make_pair(new MeshLib::Hex20(element_nodes), mat_id);
143 }
144 
145 template <>
146 std::pair<MeshLib::Element*, int> createElement<MeshLib::Prism15>(
147  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
148  int const mat_id, std::map<unsigned, unsigned> const& id_map)
149 {
150  std::vector<unsigned> node_ids;
151  readNodeIDs(in, MeshLib::Prism15::n_all_nodes, node_ids, id_map);
152 
153  std::array<MeshLib::Node*, MeshLib::Prism15::n_all_nodes> element_nodes;
154 
155  constexpr std::array node_order = {0, 1, 2, 3, 4, 5, 6, 9,
156  7, 12, 14, 13, 8, 10, 11};
157 
158  std::transform(begin(node_order), end(node_order), begin(element_nodes),
159  [&node_ids, &nodes](auto const id)
160  { return nodes[node_ids[id]]; });
161 
162  return std::make_pair(new MeshLib::Prism15(element_nodes), mat_id);
163 }
164 
165 template <>
166 std::pair<MeshLib::Element*, int> createElement<MeshLib::Pyramid13>(
167  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
168  int const mat_id, std::map<unsigned, unsigned> const& id_map)
169 {
170  std::vector<unsigned> node_ids;
171  readNodeIDs(in, MeshLib::Pyramid13::n_all_nodes, node_ids, id_map);
172  std::array<MeshLib::Node*, MeshLib::Pyramid13::n_all_nodes> element_nodes;
173 
174  constexpr std::array node_order = {0, 1, 2, 3, 4, 5, 8,
175  10, 6, 7, 9, 11, 12};
176 
177  std::transform(begin(node_order), end(node_order), begin(element_nodes),
178  [&node_ids, &nodes](auto const id)
179  { return nodes[node_ids[id]]; });
180 
181  return std::make_pair(new MeshLib::Pyramid13(element_nodes), mat_id);
182 }
183 
184 std::pair<MeshLib::Element*, int> readElement(
185  std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
186  std::map<unsigned, unsigned> const& id_map)
187 {
188  unsigned idx;
189  unsigned type;
190  unsigned n_tags;
191  unsigned dummy;
192  int mat_id;
193 
194  // element format is structured like this:
195  // element-id element-type n-tags physical-entity elementary entity node-ids
196  in >> idx >> type >> n_tags >> dummy >> mat_id;
197 
198  switch (type)
199  {
200  case 1:
201  {
202  return createElement<MeshLib::Line>(in, nodes, mat_id, id_map);
203  }
204  case 2:
205  {
206  return createElement<MeshLib::Tri>(in, nodes, mat_id, id_map);
207  }
208  case 3:
209  {
210  return createElement<MeshLib::Quad>(in, nodes, mat_id, id_map);
211  }
212  case 4:
213  {
214  return createElement<MeshLib::Tet>(in, nodes, mat_id, id_map);
215  }
216  case 5:
217  {
218  return createElement<MeshLib::Hex>(in, nodes, mat_id, id_map);
219  }
220  case 6:
221  {
222  return createElement<MeshLib::Prism>(in, nodes, mat_id, id_map);
223  }
224  case 7:
225  {
226  return createElement<MeshLib::Pyramid>(in, nodes, mat_id, id_map);
227  }
228  case 8: // 3-node second order line.
229  {
230  return createElement<MeshLib::Line3>(in, nodes, mat_id, id_map);
231  }
232  case 9: // 6-node second order triangle.
233  {
234  return createElement<MeshLib::Tri6>(in, nodes, mat_id, id_map);
235  }
236  case 10: // 9-node second order quadrangle.
237  {
238  return createElement<MeshLib::Quad9>(in, nodes, mat_id, id_map);
239  }
240  case 11: // 10-node second order tetrahedron.
241  {
242  return createElement<MeshLib::Tet10>(in, nodes, mat_id, id_map);
243  }
244  case 16: // 8-node second order quadrangle.
245  {
246  return createElement<MeshLib::Quad8>(in, nodes, mat_id, id_map);
247  }
248  case 17: // 20-node second order hexahedron.
249  {
250  return createElement<MeshLib::Hex20>(in, nodes, mat_id, id_map);
251  }
252  case 18: // 15-node second order prism.
253  {
254  return createElement<MeshLib::Prism15>(in, nodes, mat_id, id_map);
255  }
256  case 19: // 13-node second order pyramid.
257  {
258  return createElement<MeshLib::Pyramid13>(in, nodes, mat_id, id_map);
259  }
260  case 15:
261  in >> dummy; // skip rest of line
262  break;
263  default:
264  WARN("readGMSHMesh(): Unknown element type {:d}.", type);
265  break;
266  }
267  return std::make_pair(nullptr, -1);
268 }
269 
270 MeshLib::Mesh* readGMSHMesh(std::string const& fname)
271 {
272  std::string line;
273  std::ifstream in(fname.c_str(), std::ios::in);
274  if (!in.is_open())
275  {
276  WARN("readGMSHMesh() - Could not open file {:s}.", fname);
277  return nullptr;
278  }
279 
280  getline(in, line); // $MeshFormat keyword
281  if (line.find("$MeshFormat") == std::string::npos)
282  {
283  in.close();
284  WARN("No GMSH file format recognized.");
285  return nullptr;
286  }
287 
288  getline(in, line); // version-number file-type data-size
289  if (line.substr(0, 3) != "2.2")
290  {
291  WARN("Wrong gmsh file format version '{:s}'.", line.substr(0, 3));
292  return nullptr;
293  }
294 
295  if (line[4] != '0')
296  {
297  WARN("Currently reading gmsh binary file type is not supported.");
298  return nullptr;
299  }
300  getline(in, line); //$EndMeshFormat
301 
302  std::vector<MeshLib::Node*> nodes;
303  std::vector<MeshLib::Element*> elements;
304  std::vector<int> materials;
305  std::map<unsigned, unsigned> id_map;
306  while (line.find("$EndElements") == std::string::npos)
307  {
308  // Node data
309  getline(in, line); //$Nodes Keywords
310  if (line.find("$Nodes") != std::string::npos)
311  {
312  std::size_t n_nodes(0);
313  long id;
314  double x;
315  double y;
316  double z;
317  in >> n_nodes >> std::ws;
318  nodes.resize(n_nodes);
319  for (std::size_t i = 0; i < n_nodes; i++)
320  {
321  in >> id >> x >> y >> z >> std::ws;
322  id_map.insert(std::map<unsigned, unsigned>::value_type(id, i));
323  nodes[i] = new MeshLib::Node(x, y, z, id);
324  }
325  getline(in, line); // End Node keyword $EndNodes
326  }
327 
328  // Element data
329  if (line.find("$Elements") != std::string::npos)
330  {
331  std::size_t n_elements(0);
332  if (!(in >> n_elements >> std::ws))
333  { // number-of-elements
334  ERR("Read GMSH mesh does not contain any elements");
335  }
336  elements.reserve(n_elements);
337  materials.reserve(n_elements);
338  for (std::size_t i = 0; i < n_elements; i++)
339  {
340  MeshLib::Element* elem(nullptr);
341  int mat_id(0);
342  std::tie(elem, mat_id) = readElement(in, nodes, id_map);
343 
344  if (elem)
345  {
346  elements.push_back(elem);
347  materials.push_back(mat_id);
348  }
349  }
350  getline(in, line); // END keyword
351  }
352 
353  if (line.find("PhysicalNames") != std::string::npos)
354  {
355  std::size_t n_lines(0);
356  in >> n_lines >> std::ws; // number-of-lines
357  for (std::size_t i = 0; i < n_lines; i++)
358  {
359  getline(in, line);
360  }
361  getline(in, line); // END keyword
362  }
363  }
364  in.close();
365  if (elements.empty())
366  {
367  for (auto& node : nodes)
368  {
369  delete node;
370  }
371  return nullptr;
372  }
373 
374  MeshLib::Mesh* mesh(new MeshLib::Mesh(
375  BaseLib::extractBaseNameWithoutExtension(fname), nodes, elements));
376 
377  auto* const material_ids =
379  "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
380  if (!material_ids)
381  {
382  WARN("Could not create PropertyVector for MaterialIDs in Mesh.");
383  }
384  else
385  {
386  material_ids->insert(material_ids->end(), materials.cbegin(),
387  materials.cend());
388  }
389 
391 
392  INFO("\t... finished.");
393  INFO("Nr. Nodes: {:d}.", nodes.size());
394  INFO("Nr. Elements: {:d}.", elements.size());
395 
396  return mesh;
397 }
398 
399 } // end namespace GMSH
400 } // end namespace FileIO
Definition of the ElementValueModification class.
Definition of the Element class.
Filename manipulation routines.
Definition of the Hex class.
Definition of the Line class.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of the Mesh class.
Definition of the Node class.
Definition of the Prism class.
Definition of the Pyramid class.
Definition of the Quad class.
Definition of the Tet class.
Definition of the Tri class.
static std::size_t condense(MeshLib::Mesh &mesh)
Properties & getProperties()
Definition: Mesh.h:123
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
static const unsigned n_all_nodes
Constant: The number of all nodes for this element.
std::string extractBaseNameWithoutExtension(std::string const &pathname)
Definition: FileTools.cpp:180
std::pair< MeshLib::Element *, int > createElement(std::ifstream &in, std::vector< MeshLib::Node * > const &nodes, int const mat_id, std::map< unsigned, unsigned > const &id_map)
Definition: GmshReader.cpp:75
std::pair< MeshLib::Element *, int > readElement(std::ifstream &in, std::vector< MeshLib::Node * > const &nodes, std::map< unsigned, unsigned > const &id_map)
Definition: GmshReader.cpp:184
bool isGMSHMeshFile(const std::string &fname)
Definition: GmshReader.cpp:36
MeshLib::Mesh * readGMSHMesh(std::string const &fname)
Definition: GmshReader.cpp:270
void readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector< unsigned > &node_ids, std::map< unsigned, unsigned > const &id_map)
Definition: GmshReader.cpp:62