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