1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18 // C++ includes
19 #include <fstream>
20 #include <set>
21 #include <cstring> // std::memcpy
22 #include <numeric>
23 #include <unordered_map>
24 #include <cstddef>
25
26 // Local includes
27 #include "libmesh/libmesh_config.h"
28 #include "libmesh/libmesh_logging.h"
29 #include "libmesh/boundary_info.h"
30 #include "libmesh/elem.h"
31 #include "libmesh/gmsh_io.h"
32 #include "libmesh/mesh_base.h"
33 #include "libmesh/int_range.h"
34 #include "libmesh/utility.h" // map_find
35
36 namespace libMesh
37 {
38
39 // Initialize the static data member
40 GmshIO::ElementMaps GmshIO::_element_maps = GmshIO::build_element_maps();
41
42
43
44 // Definition of the static function which constructs the ElementMaps object.
build_element_maps()45 GmshIO::ElementMaps GmshIO::build_element_maps()
46 {
47 // Object to be filled up
48 ElementMaps em;
49
50 // POINT (import only)
51 em.in.emplace(15, ElementDefinition(NODEELEM, 15, 0, 1));
52
53 // Add elements with trivial node mappings
54 em.add_def(ElementDefinition(EDGE2, 1, 1, 2));
55 em.add_def(ElementDefinition(EDGE3, 8, 1, 3));
56 em.add_def(ElementDefinition(TRI3, 2, 2, 3));
57 em.add_def(ElementDefinition(TRI6, 9, 2, 6));
58 em.add_def(ElementDefinition(QUAD4, 3, 2, 4));
59 em.add_def(ElementDefinition(QUAD8, 16, 2, 8));
60 em.add_def(ElementDefinition(QUAD9, 10, 2, 9));
61 em.add_def(ElementDefinition(HEX8, 5, 3, 8));
62 em.add_def(ElementDefinition(TET4, 4, 3, 4));
63 em.add_def(ElementDefinition(PRISM6, 6, 3, 6));
64 em.add_def(ElementDefinition(PYRAMID5, 7, 3, 5));
65
66 // Add elements with non-trivial node mappings
67
68 // HEX20
69 {
70 ElementDefinition eledef(HEX20, 17, 3, 20);
71 const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,15,16,19,17,18};
72 std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
73 em.add_def(eledef);
74 }
75
76 // HEX27
77 {
78 ElementDefinition eledef(HEX27, 12, 3, 27);
79 const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
80 15,16,19,17,18,20,21,24,22,23,25,26};
81 std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
82 em.add_def(eledef);
83 }
84
85 // TET10
86 {
87 ElementDefinition eledef(TET10, 11, 3, 10);
88 const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
89 std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
90 em.add_def(eledef);
91 }
92
93 // PRISM15
94 {
95 ElementDefinition eledef(PRISM15, 18, 3, 15);
96 const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13};
97 std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
98 em.add_def(eledef);
99 }
100
101 // PRISM18
102 {
103 ElementDefinition eledef(PRISM18, 13, 3, 18);
104 const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13,15,17,16};
105 std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
106 em.add_def(eledef);
107 }
108
109 return em;
110 }
111
112
113
GmshIO(const MeshBase & mesh)114 GmshIO::GmshIO (const MeshBase & mesh) :
115 MeshOutput<MeshBase>(mesh),
116 _binary(false),
117 _write_lower_dimensional_elements(true)
118 {
119 }
120
121
122
GmshIO(MeshBase & mesh)123 GmshIO::GmshIO (MeshBase & mesh) :
124 MeshInput<MeshBase> (mesh),
125 MeshOutput<MeshBase> (mesh),
126 _binary (false),
127 _write_lower_dimensional_elements(true)
128 {
129 }
130
131
132
binary()133 bool & GmshIO::binary ()
134 {
135 return _binary;
136 }
137
138
139
write_lower_dimensional_elements()140 bool & GmshIO::write_lower_dimensional_elements ()
141 {
142 return _write_lower_dimensional_elements;
143 }
144
145
146
read(const std::string & name)147 void GmshIO::read (const std::string & name)
148 {
149 std::ifstream in (name.c_str());
150 this->read_mesh (in);
151 }
152
153
154
read_mesh(std::istream & in)155 void GmshIO::read_mesh(std::istream & in)
156 {
157 // This is a serial-only process for now;
158 // the Mesh should be read on processor 0 and
159 // broadcast later
160 libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
161
162 libmesh_assert(in.good());
163
164 // clear any data in the mesh
165 MeshBase & mesh = MeshInput<MeshBase>::mesh();
166 mesh.clear();
167
168 // some variables
169 int format=0, size=0;
170 Real version = 1.0;
171
172 // Keep track of lower-dimensional blocks which are not BCs, but
173 // actually blocks of lower-dimensional elements.
174 std::set<subdomain_id_type> lower_dimensional_blocks;
175
176 // Mapping from physical id -> (physical dim, physical name) pairs.
177 // These can refer to either "sidesets" or "subdomains"; we need to
178 // wait until the Mesh has been read to know which is which. Note
179 // that we are using 'int' as the key here rather than
180 // subdomain_id_type or boundary_id_type, since at this point, it
181 // could be either.
182 typedef std::pair<unsigned, std::string> GmshPhysical;
183 std::map<int, GmshPhysical> gmsh_physicals;
184
185 // map to hold the node numbers for translation
186 // note the the nodes can be non-consecutive
187 std::map<unsigned int, unsigned int> nodetrans;
188
189 // Map from entity tag to physical id. The key is a pair with the first
190 // item being the dimension of the entity and the second item being
191 // the entity tag/id
192 std::map<std::pair<unsigned, int>, int> entity_to_physical_id;
193
194 // For reading the file line by line
195 std::string s;
196
197 while (true)
198 {
199 // Try to read something. This may set EOF!
200 std::getline(in, s);
201
202 if (in)
203 {
204 // Process s...
205
206 if (s.find("$MeshFormat") == static_cast<std::string::size_type>(0))
207 {
208 in >> version >> format >> size;
209
210 // Some notes on gmsh mesh versions:
211 //
212 // Mesh version 2.0 goes back as far as I know. It's not explicitly
213 // mentioned here: http://www.geuz.org/gmsh/doc/VERSIONS.txt
214 //
215 // As of gmsh-2.4.0:
216 // bumped mesh version format to 2.1 (small change in the $PhysicalNames
217 // section, where the group dimension is now required);
218 // [Since we don't even parse the PhysicalNames section at the time
219 // of this writing, I don't think this change affects us.]
220 //
221 // Mesh version 2.2 tested by Manav Bhatia; no other
222 // libMesh code changes were required for support
223 //
224 // Mesh version 4.0 is a near complete rewrite of the previous mesh version
225 libmesh_error_msg_if(version < 2.0, "Error: Unknown msh file version " << version);
226 libmesh_error_msg_if(format, "Error: Unknown data format for mesh in Gmsh reader.");
227 }
228
229 // Read and process the "PhysicalNames" section.
230 else if (s.find("$PhysicalNames") == static_cast<std::string::size_type>(0))
231 {
232 // The lines in the PhysicalNames section should look like the following:
233 // 2 1 "frac" lower_dimensional_block
234 // 2 3 "top"
235 // 2 4 "bottom"
236 // 3 2 "volume"
237
238 // Read in the number of physical groups to expect in the file.
239 unsigned int num_physical_groups = 0;
240 in >> num_physical_groups;
241
242 // Read rest of line including newline character.
243 std::getline(in, s);
244
245 for (unsigned int i=0; i<num_physical_groups; ++i)
246 {
247 // Read an entire line of the PhysicalNames section.
248 std::getline(in, s);
249
250 // Use an istringstream to extract the physical
251 // dimension, physical id, and physical name from
252 // this line.
253 std::istringstream s_stream(s);
254 unsigned phys_dim;
255 int phys_id;
256 std::string phys_name;
257 s_stream >> phys_dim >> phys_id >> phys_name;
258
259 // Not sure if this is true for all Gmsh files, but
260 // my test file has quotes around the phys_name
261 // string. So let's erase any quotes now...
262 phys_name.erase(std::remove(phys_name.begin(), phys_name.end(), '"'), phys_name.end());
263
264 // Record this ID for later assignment of subdomain/sideset names.
265 gmsh_physicals[phys_id] = std::make_pair(phys_dim, phys_name);
266
267 // If 's' also contains the libmesh-specific string
268 // "lower_dimensional_block", add this block ID to
269 // the list of blocks which are not boundary
270 // conditions.
271 if (s.find("lower_dimensional_block") != std::string::npos)
272 {
273 lower_dimensional_blocks.insert(cast_int<subdomain_id_type>(phys_id));
274
275 // The user has explicitly told us that this
276 // block is a subdomain, so set that association
277 // in the Mesh.
278 mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
279 }
280 }
281 }
282
283 else if (s.find("$Entities") == static_cast<std::string::size_type>(0))
284 {
285 if (version >= 4.0)
286 {
287 std::size_t num_point_entities, num_curve_entities, num_surface_entities, num_volume_entities;
288 in >> num_point_entities >> num_curve_entities >> num_surface_entities >> num_volume_entities;
289
290 for (std::size_t n = 0; n < num_point_entities; ++n)
291 {
292 int point_tag, physical_tag;
293 Real x, y, z;
294 std::size_t num_physical_tags;
295 in >> point_tag >> x >> y >> z >> num_physical_tags;
296
297 libmesh_error_msg_if(num_physical_tags > 1,
298 "Sorry, you cannot currently specify multiple subdomain or "
299 "boundary ids for a given geometric entity");
300
301 if (num_physical_tags)
302 {
303 in >> physical_tag;
304 entity_to_physical_id[std::make_pair(0, point_tag)] = physical_tag;
305 }
306 }
307 for (std::size_t n = 0; n < num_curve_entities; ++n)
308 {
309 int curve_tag, physical_tag;
310 Real minx, miny, minz, maxx, maxy, maxz;
311 std::size_t num_physical_tags;
312 in >> curve_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
313
314 libmesh_error_msg_if(num_physical_tags > 1,
315 "I don't believe that we can specify multiple subdomain or "
316 "boundary ids for a given geometric entity");
317
318 if (num_physical_tags)
319 {
320 in >> physical_tag;
321 entity_to_physical_id[std::make_pair(1, curve_tag)] = physical_tag;
322 }
323
324 // Read to end of line; this captures bounding information that we don't care about
325 std::getline(in, s);
326 }
327 for (std::size_t n = 0; n < num_surface_entities; ++n)
328 {
329 int surface_tag, physical_tag;
330 Real minx, miny, minz, maxx, maxy, maxz;
331 std::size_t num_physical_tags;
332 in >> surface_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
333
334 libmesh_error_msg_if(num_physical_tags > 1,
335 "I don't believe that we can specify multiple subdomain or "
336 "boundary ids for a given geometric entity");
337
338 if (num_physical_tags)
339 {
340 in >> physical_tag;
341 entity_to_physical_id[std::make_pair(2, surface_tag)] = physical_tag;
342 }
343
344 // Read to end of line; this captures bounding information that we don't care about
345 std::getline(in, s);
346 }
347 for (std::size_t n = 0; n < num_volume_entities; ++n)
348 {
349 int volume_tag, physical_tag;
350 Real minx, miny, minz, maxx, maxy, maxz;
351 std::size_t num_physical_tags;
352 in >> volume_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
353
354 libmesh_error_msg_if(num_physical_tags > 1,
355 "I don't believe that we can specify multiple subdomain or "
356 "boundary ids for a given geometric entity");
357
358 if (num_physical_tags)
359 {
360 in >> physical_tag;
361 entity_to_physical_id[std::make_pair(3, volume_tag)] = physical_tag;
362 }
363
364 // Read to end of line; this captures bounding information that we don't care about
365 std::getline(in, s);
366 }
367 // Read the $EndEntities
368 std::getline(in, s);
369 } // end if (version >= 4.0)
370
371 else
372 libmesh_error_msg("The Entities block was introduced in mesh version 4.0");
373 } // end if $Entities
374
375 // read the node block
376 else if (s.find("$NOD") == static_cast<std::string::size_type>(0) ||
377 s.find("$NOE") == static_cast<std::string::size_type>(0) ||
378 s.find("$Nodes") == static_cast<std::string::size_type>(0))
379 {
380 if (version < 4.0)
381 {
382 unsigned int num_nodes = 0;
383 in >> num_nodes;
384 mesh.reserve_nodes (num_nodes);
385
386 // read in the nodal coordinates and form points.
387 Real x, y, z;
388 unsigned int id;
389
390 // add the nodal coordinates to the mesh
391 for (unsigned int i=0; i<num_nodes; ++i)
392 {
393 in >> id >> x >> y >> z;
394 mesh.add_point (Point(x, y, z), i);
395 nodetrans[id] = i;
396 }
397 }
398 else
399 {
400 // Read numEntityBlocks line
401 std::size_t num_entities = 0, num_nodes = 0, min_node_tag, max_node_tag;
402 in >> num_entities >> num_nodes >> min_node_tag >> max_node_tag;
403
404 mesh.reserve_nodes(num_nodes);
405
406 std::size_t node_counter = 0;
407
408 // Now loop over entities
409 for (std::size_t i = 0; i < num_entities; ++i)
410 {
411 int entity_dim, entity_tag, parametric;
412 std::size_t num_nodes_in_block = 0;
413 in >> entity_dim >> entity_tag >> parametric >> num_nodes_in_block;
414 libmesh_error_msg_if(parametric, "We don't currently support reading parametric gmsh entities");
415
416 // Read the node tags/ids
417 std::size_t gmsh_id;
418 for (std::size_t n = 0; n < num_nodes_in_block; ++n)
419 {
420
421 in >> gmsh_id;
422 nodetrans[gmsh_id] = node_counter++;
423 }
424
425 // Read the node coordinates and add the nodes to the mesh
426 Real x, y, z;
427 for (std::size_t libmesh_id = node_counter - num_nodes_in_block;
428 libmesh_id < node_counter;
429 ++libmesh_id)
430 {
431 in >> x >> y >> z;
432 mesh.add_point(Point(x, y, z), libmesh_id);
433 }
434 }
435 }
436 // read the $ENDNOD delimiter
437 std::getline(in, s);
438 }
439
440 // Read the element block
441 else if (s.find("$ELM") == static_cast<std::string::size_type>(0) ||
442 s.find("$Elements") == static_cast<std::string::size_type>(0))
443 {
444 // Keep track of element dimensions seen
445 std::vector<unsigned> elem_dimensions_seen(3);
446
447 if (version < 4.0)
448 {
449 // For reading the number of elements and the node ids from the stream
450 unsigned int
451 num_elem = 0,
452 node_id = 0;
453
454 // read how many elements are there, and reserve space in the mesh
455 in >> num_elem;
456 mesh.reserve_elem (num_elem);
457
458 // As of version 2.2, the format for each element line is:
459 // elm-number elm-type number-of-tags < tag > ... node-number-list
460 // From the Gmsh docs:
461 // * the first tag is the number of the
462 // physical entity to which the element belongs
463 // * the second is the number of the elementary geometrical
464 // entity to which the element belongs
465 // * the third is the number of mesh partitions to which the element
466 // belongs
467 // * The rest of the tags are the partition ids (negative
468 // partition ids indicate ghost cells). A zero tag is
469 // equivalent to no tag. Gmsh and most codes using the
470 // MSH 2 format require at least the first two tags
471 // (physical and elementary tags).
472
473 // read the elements
474 for (unsigned int iel=0; iel<num_elem; ++iel)
475 {
476 unsigned int
477 id, type,
478 physical=1, elementary=1,
479 nnodes=0, ntags;
480
481 // Note: tag has to be an int because it could be negative,
482 // see above.
483 int tag;
484
485 if (version <= 1.0)
486 in >> id >> type >> physical >> elementary >> nnodes;
487
488 else
489 {
490 in >> id >> type >> ntags;
491
492 if (ntags > 2)
493 libmesh_do_once(libMesh::err << "Warning, ntags=" << ntags << ", but we currently only support reading 2 flags." << std::endl;);
494
495 for (unsigned int j = 0; j < ntags; j++)
496 {
497 in >> tag;
498 if (j == 0)
499 physical = tag;
500 else if (j == 1)
501 elementary = tag;
502 }
503 }
504
505 // Get a reference to the ElementDefinition, throw an error if not found.
506 const GmshIO::ElementDefinition & eletype =
507 libmesh_map_find(_element_maps.in, type);
508
509 // If we read nnodes, make sure it matches the number in eletype.nnodes
510 libmesh_error_msg_if(nnodes != 0 && nnodes != eletype.nnodes,
511 "nnodes = " << nnodes << " and eletype.nnodes = " << eletype.nnodes << " do not match.");
512
513 // Assign the value from the eletype object.
514 nnodes = eletype.nnodes;
515
516 // Don't add 0-dimensional "point" elements to the
517 // Mesh. They should *always* be treated as boundary
518 // "nodeset" data.
519 if (eletype.dim > 0)
520 {
521 // Record this element dimension as being "seen".
522 // We will treat all elements with dimension <
523 // max(dimension) as specifying boundary conditions,
524 // but we won't know what max_elem_dimension_seen is
525 // until we read the entire file.
526 elem_dimensions_seen[eletype.dim-1] = 1;
527
528 // Add the element to the mesh
529 {
530 Elem * elem =
531 mesh.add_elem(Elem::build_with_id(eletype.type, iel));
532
533 // Make sure that the libmesh element we added has nnodes nodes.
534 libmesh_error_msg_if(elem->n_nodes() != nnodes,
535 "Number of nodes for element "
536 << id
537 << " of type " << eletype.type
538 << " (Gmsh type " << type
539 << ") does not match Libmesh definition. "
540 << "I expected " << elem->n_nodes()
541 << " nodes, but got " << nnodes);
542
543 // Add node pointers to the elements.
544 // If there is a node translation table, use it.
545 if (eletype.nodes.size() > 0)
546 for (unsigned int i=0; i<nnodes; i++)
547 {
548 in >> node_id;
549 elem->set_node(eletype.nodes[i]) = mesh.node_ptr(nodetrans[node_id]);
550 }
551 else
552 {
553 for (unsigned int i=0; i<nnodes; i++)
554 {
555 in >> node_id;
556 elem->set_node(i) = mesh.node_ptr(nodetrans[node_id]);
557 }
558 }
559
560 // Finally, set the subdomain ID to physical. If this is a lower-dimension element, this ID will
561 // eventually go into the Mesh's BoundaryInfo object.
562 elem->subdomain_id() = static_cast<subdomain_id_type>(physical);
563 }
564 }
565
566 // Handle 0-dimensional elements (points) by adding
567 // them to the BoundaryInfo object with
568 // boundary_id == physical.
569 else
570 {
571 // This seems like it should always be the same
572 // number as the 'id' we already read in on this
573 // line. At least it was in the example gmsh
574 // file I had...
575 in >> node_id;
576 mesh.get_boundary_info().add_node
577 (nodetrans[node_id],
578 static_cast<boundary_id_type>(physical));
579 }
580 } // element loop
581 } // end if (version < 4.0)
582
583 else
584 {
585 std::size_t num_entity_blocks = 0, num_elem = 0, min_element_tag, max_element_tag;
586
587 // Read entity information
588 in >> num_entity_blocks >> num_elem >> min_element_tag >> max_element_tag;
589
590 mesh.reserve_elem(num_elem);
591
592 std::size_t iel = 0;
593
594 // Loop over entity blocks
595 for (std::size_t i = 0; i < num_entity_blocks; ++i)
596 {
597 int entity_dim, entity_tag;
598 unsigned int element_type;
599 std::size_t num_elems_in_block = 0;
600 in >> entity_dim >> entity_tag >> element_type >> num_elems_in_block;
601
602 // Get a reference to the ElementDefinition
603 const GmshIO::ElementDefinition & eletype =
604 libmesh_map_find(_element_maps.in, element_type);
605
606 // Don't add 0-dimensional "point" elements to the
607 // Mesh. They should *always* be treated as boundary
608 // "nodeset" data.
609 if (eletype.dim > 0)
610 {
611 // Record this element dimension as being "seen".
612 // We will treat all elements with dimension <
613 // max(dimension) as specifying boundary conditions,
614 // but we won't know what max_elem_dimension_seen is
615 // until we read the entire file.
616 elem_dimensions_seen[eletype.dim-1] = 1;
617
618 // Loop over elements with dim > 0
619 for (std::size_t n = 0; n < num_elems_in_block; ++n)
620 {
621 Elem * elem =
622 mesh.add_elem(Elem::build_with_id(eletype.type, iel++));
623
624 std::size_t gmsh_element_id;
625 in >> gmsh_element_id;
626
627 // Get the remainer of the line that includes the nodes ids
628 std::getline(in, s);
629 std::istringstream is(s);
630 std::size_t local_node_counter = 0, gmsh_node_id;
631 while (is >> gmsh_node_id)
632 {
633 // Add node pointers to the elements.
634 // If there is a node translation table, use it.
635 if (eletype.nodes.size() > 0)
636 elem->set_node(eletype.nodes[local_node_counter++]) =
637 mesh.node_ptr(nodetrans[gmsh_node_id]);
638 else
639 elem->set_node(local_node_counter++) = mesh.node_ptr(nodetrans[gmsh_node_id]);
640 }
641
642 // Make sure that the libmesh element we added has nnodes nodes.
643 libmesh_error_msg_if(elem->n_nodes() != local_node_counter,
644 "Number of nodes for element "
645 << gmsh_element_id
646 << " of type " << eletype.type
647 << " (Gmsh type " << element_type
648 << ") does not match Libmesh definition. "
649 << "I expected " << elem->n_nodes()
650 << " nodes, but got " << local_node_counter);
651
652 // Finally, set the subdomain ID to physical. If this is a lower-dimension element, this ID will
653 // eventually go into the Mesh's BoundaryInfo object.
654 elem->subdomain_id() = static_cast<subdomain_id_type>(
655 entity_to_physical_id[std::make_pair(entity_dim, entity_tag)]);
656
657 } // end for (loop over elements in entity block)
658 } // end if (eletype.dim > 0)
659
660 else
661 {
662 for (std::size_t n = 0; n < num_elems_in_block; ++n)
663 {
664 std::size_t gmsh_element_id, gmsh_node_id;
665 in >> gmsh_element_id;
666 in >> gmsh_node_id;
667 mesh.get_boundary_info().add_node(
668 nodetrans[gmsh_node_id],
669 static_cast<boundary_id_type>(entity_to_physical_id[
670 std::make_pair(entity_dim, entity_tag)]));
671 } // end for (loop over elements in entity block)
672 } // end if (eletype.dim == 0)
673 } // end for (loop over entity blocks)
674 } // end if (version >= 4.0)
675
676 // read the $ENDELM delimiter
677 std::getline(in, s);
678
679 // Record the max and min element dimension seen while reading the file.
680 unsigned char
681 max_elem_dimension_seen=1,
682 min_elem_dimension_seen=3;
683
684 for (auto i : index_range(elem_dimensions_seen))
685 if (elem_dimensions_seen[i])
686 {
687 // Debugging
688 // libMesh::out << "Seen elements of dimension " << i+1 << std::endl;
689 max_elem_dimension_seen =
690 std::max(max_elem_dimension_seen, cast_int<unsigned char>(i+1));
691 min_elem_dimension_seen =
692 std::min(min_elem_dimension_seen, cast_int<unsigned char>(i+1));
693 }
694
695 // Debugging:
696 // libMesh::out << "max_elem_dimension_seen=" << max_elem_dimension_seen << std::endl;
697 // libMesh::out << "min_elem_dimension_seen=" << min_elem_dimension_seen << std::endl;
698
699
700 // How many different element dimensions did we see while reading from file?
701 unsigned n_dims_seen = std::accumulate(elem_dimensions_seen.begin(),
702 elem_dimensions_seen.end(),
703 static_cast<unsigned>(0),
704 std::plus<unsigned>());
705
706
707 // Set mesh_dimension based on the largest element dimension seen.
708 mesh.set_mesh_dimension(max_elem_dimension_seen);
709
710 // Now that we know the maximum element dimension seen,
711 // we know whether the physical names are subdomain
712 // names or sideset names.
713 for (const auto & pr : gmsh_physicals)
714 {
715 // Extract data
716 int phys_id = pr.first;
717 unsigned phys_dim = pr.second.first;
718 const std::string & phys_name = pr.second.second;
719
720 // If the physical's dimension matches the largest
721 // dimension we've seen, it's a subdomain name.
722 if (phys_dim == max_elem_dimension_seen)
723 mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
724
725 // Otherwise, if it's not a lower-dimensional
726 // block, it's a sideset name.
727 else if (phys_dim < max_elem_dimension_seen &&
728 !lower_dimensional_blocks.count(cast_int<boundary_id_type>(phys_id)))
729 mesh.get_boundary_info().sideset_name(cast_int<boundary_id_type>(phys_id)) = phys_name;
730 }
731
732 if (n_dims_seen > 1)
733 {
734 // Store lower-dimensional elements in a map sorted
735 // by Elem::key(). We use a multimap for two reasons:
736 // 1.) The hash function is not guaranteed to be
737 // unique, so different lower-dimensional elements
738 // could theoretically hash to the same value,
739 // although this is pretty unlikely.
740 // 2.) The Gmsh file may contain multiple
741 // lower-dimensional elements for a single side in
742 // order to implement multiple boundary ids for a
743 // single side. These lower-dimensional elements
744 // will all hash to the same value, and we need to
745 // be able to store all of them.
746 typedef std::unordered_multimap<dof_id_type, Elem *> provide_container_t;
747 provide_container_t provide_bcs;
748
749 // 1st loop over active elements - get info about lower-dimensional elements.
750 for (auto & elem : mesh.active_element_ptr_range())
751 if (elem->dim() < max_elem_dimension_seen &&
752 !lower_dimensional_blocks.count(elem->subdomain_id()))
753 {
754 // To be consistent with the previous
755 // GmshIO behavior, add all the
756 // lower-dimensional elements' nodes to
757 // the Mesh's BoundaryInfo object with the
758 // lower-dimensional element's subdomain
759 // ID.
760 for (auto n : elem->node_index_range())
761 mesh.get_boundary_info().add_node(elem->node_id(n),
762 elem->subdomain_id());
763
764 // Store this elem in a quickly-searchable
765 // container to use it to assign boundary
766 // conditions later.
767 provide_bcs.emplace(elem->key(), elem);
768 }
769
770 // 2nd loop over active elements - use lower dimensional element data to set BCs for higher dimensional elements
771 for (auto & elem : mesh.active_element_ptr_range())
772 if (elem->dim() == max_elem_dimension_seen)
773 {
774 // This is a max-dimension element that
775 // may require BCs. For each of its
776 // sides, including internal sides, we'll
777 // see if one more more lower-dimensional elements
778 // provides boundary information for it.
779 // Note that we have not yet called
780 // find_neighbors(), so we can't use
781 // elem->neighbor(sn) in this algorithm...
782 for (auto sn : elem->side_index_range())
783 for (const auto & pr : as_range(provide_bcs.equal_range(elem->key(sn))))
784 {
785 // For each side side in the provide_bcs multimap...
786 // Construct the side for hash verification.
787 std::unique_ptr<Elem> side (elem->build_side_ptr(sn));
788
789 // Construct the lower-dimensional element to compare to the side.
790 Elem * lower_dim_elem = pr.second;
791
792 // This was a hash, so it might not be perfect. Let's verify...
793 if (*lower_dim_elem == *side)
794 {
795 // Add the lower-dimensional
796 // element's subdomain_id as a
797 // boundary_id for the
798 // higher-dimensional element.
799 boundary_id_type bid = cast_int<boundary_id_type>(lower_dim_elem->subdomain_id());
800 mesh.get_boundary_info().add_side(elem, sn, bid);
801 }
802 }
803 }
804
805 // 3rd loop over active elements - Remove the lower-dimensional elements
806 for (auto & elem : mesh.active_element_ptr_range())
807 if (elem->dim() < max_elem_dimension_seen &&
808 !lower_dimensional_blocks.count(elem->subdomain_id()))
809 mesh.delete_elem(elem);
810 } // end if (n_dims_seen > 1)
811 } // if $ELM
812
813 continue;
814 } // if (in)
815
816
817 // If !in, check to see if EOF was set. If so, break out
818 // of while loop.
819 if (in.eof())
820 break;
821
822 // If !in and !in.eof(), stream is in a bad state!
823 libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
824
825 } // while true
826 }
827
828
829
write(const std::string & name)830 void GmshIO::write (const std::string & name)
831 {
832 if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
833 {
834 // Open the output file stream
835 std::ofstream out_stream (name.c_str());
836
837 // Make sure it opened correctly
838 if (!out_stream.good())
839 libmesh_file_error(name.c_str());
840
841 this->write_mesh (out_stream);
842 }
843 }
844
845
846
write_nodal_data(const std::string & fname,const std::vector<Number> & soln,const std::vector<std::string> & names)847 void GmshIO::write_nodal_data (const std::string & fname,
848 const std::vector<Number> & soln,
849 const std::vector<std::string> & names)
850 {
851 LOG_SCOPE("write_nodal_data()", "GmshIO");
852
853 if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
854 this->write_post (fname, &soln, &names);
855 }
856
857
858
write_mesh(std::ostream & out_stream)859 void GmshIO::write_mesh (std::ostream & out_stream)
860 {
861 // Be sure that the stream is valid.
862 libmesh_assert (out_stream.good());
863
864 // Get a const reference to the mesh
865 const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
866
867 // If requested, write out lower-dimensional elements for
868 // element-side-based boundary conditions.
869 std::size_t n_boundary_faces = 0;
870 if (this->write_lower_dimensional_elements())
871 n_boundary_faces = mesh.get_boundary_info().n_boundary_conds();
872
873 // Note: we are using version 2.0 of the gmsh output format.
874
875 // Write the file header.
876 out_stream << "$MeshFormat\n";
877 out_stream << "2.0 0 " << sizeof(Real) << '\n';
878 out_stream << "$EndMeshFormat\n";
879
880 // write the nodes in (n x y z) format
881 out_stream << "$Nodes\n";
882 out_stream << mesh.n_nodes() << '\n';
883
884 for (auto v : make_range(mesh.n_nodes()))
885 out_stream << mesh.node_ref(v).id()+1 << " "
886 << mesh.node_ref(v)(0) << " "
887 << mesh.node_ref(v)(1) << " "
888 << mesh.node_ref(v)(2) << '\n';
889 out_stream << "$EndNodes\n";
890
891 {
892 // write the connectivity
893 out_stream << "$Elements\n";
894 out_stream << mesh.n_active_elem() + n_boundary_faces << '\n';
895
896 // loop over the elements
897 for (const auto & elem : mesh.active_element_ptr_range())
898 {
899 // Get a reference to the ElementDefinition object
900 const ElementDefinition & eletype =
901 libmesh_map_find(_element_maps.out, elem->type());
902
903 // The element mapper better not require any more nodes
904 // than are present in the current element!
905 libmesh_assert_less_equal (eletype.nodes.size(), elem->n_nodes());
906
907 // elements ids are 1 based in Gmsh
908 out_stream << elem->id()+1 << " ";
909 // element type
910 out_stream << eletype.gmsh_type;
911
912 // write the number of tags (3) and their values:
913 // 1 (physical entity)
914 // 2 (geometric entity)
915 // 3 (partition entity)
916 out_stream << " 3 "
917 << static_cast<unsigned int>(elem->subdomain_id())
918 << " 0 "
919 << elem->processor_id()+1
920 << " ";
921
922 // if there is a node translation table, use it
923 if (eletype.nodes.size() > 0)
924 for (auto i : elem->node_index_range())
925 out_stream << elem->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
926 // otherwise keep the same node order
927 else
928 for (auto i : elem->node_index_range())
929 out_stream << elem->node_id(i)+1 << " "; // gmsh is 1-based
930 out_stream << "\n";
931 } // element loop
932 }
933
934 {
935 // A counter for writing surface elements to the Gmsh file
936 // sequentially. We start numbering them with a number strictly
937 // larger than the largest element ID in the mesh. Note: the
938 // MeshBase docs say "greater than or equal to" the maximum
939 // element id in the mesh, so technically we might need a +1 here,
940 // but all of the implementations return an ID strictly greater
941 // than the largest element ID in the Mesh.
942 unsigned int e_id = mesh.max_elem_id();
943
944 // loop over the elements, writing out boundary faces
945 if (n_boundary_faces)
946 {
947 // Build a list of (elem, side, bc) tuples.
948 auto bc_triples = mesh.get_boundary_info().build_side_list();
949
950 // Loop over these lists, writing data to the file.
951 for (const auto & t : bc_triples)
952 {
953 const Elem & elem = mesh.elem_ref(std::get<0>(t));
954
955 std::unique_ptr<const Elem> side = elem.build_side_ptr(std::get<1>(t));
956
957 // consult the export element table
958 const GmshIO::ElementDefinition & eletype =
959 libmesh_map_find(_element_maps.out, side->type());
960
961 // The element mapper better not require any more nodes
962 // than are present in the current element!
963 libmesh_assert_less_equal (eletype.nodes.size(), side->n_nodes());
964
965 // elements ids are 1-based in Gmsh
966 out_stream << e_id+1 << " ";
967
968 // element type
969 out_stream << eletype.gmsh_type;
970
971 // write the number of tags:
972 // 1 (physical entity)
973 // 2 (geometric entity)
974 // 3 (partition entity)
975 out_stream << " 3 "
976 << std::get<2>(t)
977 << " 0 "
978 << elem.processor_id()+1
979 << " ";
980
981 // if there is a node translation table, use it
982 if (eletype.nodes.size() > 0)
983 for (auto i : side->node_index_range())
984 out_stream << side->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
985
986 // otherwise keep the same node order
987 else
988 for (auto i : side->node_index_range())
989 out_stream << side->node_id(i)+1 << " "; // gmsh is 1-based
990
991 // Go to the next line
992 out_stream << "\n";
993
994 // increment this index too...
995 ++e_id;
996 }
997 }
998
999 out_stream << "$EndElements\n";
1000 }
1001 }
1002
1003
1004
write_post(const std::string & fname,const std::vector<Number> * v,const std::vector<std::string> * solution_names)1005 void GmshIO::write_post (const std::string & fname,
1006 const std::vector<Number> * v,
1007 const std::vector<std::string> * solution_names)
1008 {
1009
1010 // Should only do this on processor 0!
1011 libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
1012
1013 // Create an output stream
1014 std::ofstream out_stream(fname.c_str());
1015
1016 // Make sure it opened correctly
1017 if (!out_stream.good())
1018 libmesh_file_error(fname.c_str());
1019
1020 // create a character buffer
1021 char buf[80];
1022
1023 // Get a constant reference to the mesh.
1024 const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1025
1026 // write the data
1027 if ((solution_names != nullptr) && (v != nullptr))
1028 {
1029 const unsigned int n_vars =
1030 cast_int<unsigned int>(solution_names->size());
1031
1032 if (!(v->size() == mesh.n_nodes()*n_vars))
1033 libMesh::err << "ERROR: v->size()=" << v->size()
1034 << ", mesh.n_nodes()=" << mesh.n_nodes()
1035 << ", n_vars=" << n_vars
1036 << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
1037 << "\n";
1038
1039 libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
1040
1041 // write the header
1042 out_stream << "$PostFormat\n";
1043 if (this->binary())
1044 out_stream << "1.2 1 " << sizeof(double) << "\n";
1045 else
1046 out_stream << "1.2 0 " << sizeof(double) << "\n";
1047 out_stream << "$EndPostFormat\n";
1048
1049 // Loop over the elements to see how much of each type there are
1050 unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
1051 n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
1052 unsigned int n_scalar=0, n_vector=0, n_tensor=0;
1053 unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;
1054
1055 {
1056 for (const auto & elem : mesh.active_element_ptr_range())
1057 {
1058 const ElemType elemtype = elem->type();
1059
1060 switch (elemtype)
1061 {
1062 case EDGE2:
1063 case EDGE3:
1064 case EDGE4:
1065 {
1066 n_lines += 1;
1067 break;
1068 }
1069 case TRI3:
1070 case TRI6:
1071 {
1072 n_triangles += 1;
1073 break;
1074 }
1075 case QUAD4:
1076 case QUAD8:
1077 case QUAD9:
1078 {
1079 n_quadrangles += 1;
1080 break;
1081 }
1082 case TET4:
1083 case TET10:
1084 {
1085 n_tetrahedra += 1;
1086 break;
1087 }
1088 case HEX8:
1089 case HEX20:
1090 case HEX27:
1091 {
1092 n_hexahedra += 1;
1093 break;
1094 }
1095 case PRISM6:
1096 case PRISM15:
1097 case PRISM18:
1098 {
1099 n_prisms += 1;
1100 break;
1101 }
1102 case PYRAMID5:
1103 {
1104 n_pyramids += 1;
1105 break;
1106 }
1107 default:
1108 libmesh_error_msg("ERROR: Nonexistent element type " << elem->type());
1109 }
1110 }
1111 }
1112
1113 // create a view for each variable
1114 for (unsigned int ivar=0; ivar < n_vars; ivar++)
1115 {
1116 std::string varname = (*solution_names)[ivar];
1117
1118 // at the moment, we just write out scalar quantities
1119 // later this should be made configurable through
1120 // options to the writer class
1121 n_scalar = 1;
1122
1123 // write the variable as a view, and the number of time steps
1124 out_stream << "$View\n" << varname << " " << 1 << "\n";
1125
1126 // write how many of each geometry type are written
1127 out_stream << n_points * n_scalar << " "
1128 << n_points * n_vector << " "
1129 << n_points * n_tensor << " "
1130 << n_lines * n_scalar << " "
1131 << n_lines * n_vector << " "
1132 << n_lines * n_tensor << " "
1133 << n_triangles * n_scalar << " "
1134 << n_triangles * n_vector << " "
1135 << n_triangles * n_tensor << " "
1136 << n_quadrangles * n_scalar << " "
1137 << n_quadrangles * n_vector << " "
1138 << n_quadrangles * n_tensor << " "
1139 << n_tetrahedra * n_scalar << " "
1140 << n_tetrahedra * n_vector << " "
1141 << n_tetrahedra * n_tensor << " "
1142 << n_hexahedra * n_scalar << " "
1143 << n_hexahedra * n_vector << " "
1144 << n_hexahedra * n_tensor << " "
1145 << n_prisms * n_scalar << " "
1146 << n_prisms * n_vector << " "
1147 << n_prisms * n_tensor << " "
1148 << n_pyramids * n_scalar << " "
1149 << n_pyramids * n_vector << " "
1150 << n_pyramids * n_tensor << " "
1151 << nb_text2d << " "
1152 << nb_text2d_chars << " "
1153 << nb_text3d << " "
1154 << nb_text3d_chars << "\n";
1155
1156 // if binary, write a marker to identify the endianness of the file
1157 if (this->binary())
1158 {
1159 const int one = 1;
1160 std::memcpy(buf, &one, sizeof(int));
1161 out_stream.write(buf, sizeof(int));
1162 }
1163
1164 // the time steps (there is just 1 at the moment)
1165 if (this->binary())
1166 {
1167 double one = 1;
1168 std::memcpy(buf, &one, sizeof(double));
1169 out_stream.write(buf, sizeof(double));
1170 }
1171 else
1172 out_stream << "1\n";
1173
1174 // Loop over the elements and write out the data
1175 for (const auto & elem : mesh.active_element_ptr_range())
1176 {
1177 const unsigned int nv = elem->n_vertices();
1178 // this is quite crappy, but I did not invent that file format!
1179 for (unsigned int d=0; d<3; d++) // loop over the dimensions
1180 {
1181 for (unsigned int n=0; n < nv; n++) // loop over vertices
1182 {
1183 const Point & vertex = elem->point(n);
1184 if (this->binary())
1185 {
1186 #if defined(LIBMESH_DEFAULT_TRIPLE_PRECISION) || defined(LIBMESH_DEFAULT_QUADRUPLE_PRECISION)
1187 libmesh_warning("Gmsh binary writes use only double precision!");
1188 #endif
1189 double tmp = double(vertex(d));
1190 std::memcpy(buf, &tmp, sizeof(double));
1191 out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1192 }
1193 else
1194 out_stream << vertex(d) << " ";
1195 }
1196 if (!this->binary())
1197 out_stream << "\n";
1198 }
1199
1200 // now finally write out the data
1201 for (unsigned int i=0; i < nv; i++) // loop over vertices
1202 if (this->binary())
1203 {
1204 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1205 libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1206 << "complex numbers. Will only write the real part of "
1207 << "variable " << varname << std::endl;
1208 #endif
1209 double tmp = double(libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]));
1210 std::memcpy(buf, &tmp, sizeof(double));
1211 out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1212 }
1213 else
1214 {
1215 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1216 libMesh::out << "WARNING: Gmsh::write_post does not fully support "
1217 << "complex numbers. Will only write the real part of "
1218 << "variable " << varname << std::endl;
1219 #endif
1220 out_stream << libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]) << "\n";
1221 }
1222 }
1223 if (this->binary())
1224 out_stream << "\n";
1225 out_stream << "$EndView\n";
1226
1227 } // end variable loop (writing the views)
1228 }
1229 }
1230
1231 } // namespace libMesh
1232