Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local includes
19 : #include "libmesh/libmesh_config.h"
20 : #include "libmesh/libmesh_logging.h"
21 : #include "libmesh/boundary_info.h"
22 : #include "libmesh/bounding_box.h"
23 : #include "libmesh/elem.h"
24 : #include "libmesh/gmsh_io.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/utility.h" // map_find
28 : #include "libmesh/enum_to_string.h"
29 :
30 : // C++ includes
31 : #include <fstream>
32 : #include <set>
33 : #include <cstring> // std::memcpy
34 : #include <numeric>
35 : #include <unordered_map>
36 : #include <cstddef>
37 :
38 : namespace libMesh
39 : {
40 :
41 : // Initialize the static data member
42 : GmshIO::ElementMaps GmshIO::_element_maps = GmshIO::build_element_maps();
43 :
44 :
45 :
46 : // Definition of the static function which constructs the ElementMaps object.
47 16389 : GmshIO::ElementMaps GmshIO::build_element_maps()
48 : {
49 : // Object to be filled up
50 462 : ElementMaps em;
51 :
52 : // POINT (import only)
53 16389 : em.in.emplace(15, ElementDefinition(NODEELEM, 15, 0, 1));
54 :
55 : // Add elements with trivial node mappings
56 32316 : em.add_def(ElementDefinition(EDGE2, 1, 1, 2));
57 32316 : em.add_def(ElementDefinition(EDGE3, 8, 1, 3));
58 32316 : em.add_def(ElementDefinition(TRI3, 2, 2, 3));
59 32316 : em.add_def(ElementDefinition(TRI6, 9, 2, 6));
60 32316 : em.add_def(ElementDefinition(QUAD4, 3, 2, 4));
61 32316 : em.add_def(ElementDefinition(QUAD8, 16, 2, 8));
62 32316 : em.add_def(ElementDefinition(QUAD9, 10, 2, 9));
63 32316 : em.add_def(ElementDefinition(HEX8, 5, 3, 8));
64 32316 : em.add_def(ElementDefinition(TET4, 4, 3, 4));
65 32316 : em.add_def(ElementDefinition(PRISM6, 6, 3, 6));
66 32316 : em.add_def(ElementDefinition(PYRAMID5, 7, 3, 5));
67 :
68 : // Add elements with non-trivial node mappings
69 :
70 : // HEX20
71 : {
72 924 : ElementDefinition eledef(HEX20, 17, 3, 20);
73 16389 : 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 16389 : std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
75 16389 : em.add_def(eledef);
76 : }
77 :
78 : // HEX27
79 : {
80 924 : ElementDefinition eledef(HEX27, 12, 3, 27);
81 16389 : 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 16389 : std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
84 16389 : em.add_def(eledef);
85 : }
86 :
87 : // TET10
88 : {
89 924 : ElementDefinition eledef(TET10, 11, 3, 10);
90 16389 : const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
91 16389 : std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
92 16389 : em.add_def(eledef);
93 : }
94 :
95 : // PRISM15
96 : {
97 924 : ElementDefinition eledef(PRISM15, 18, 3, 15);
98 16389 : const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13};
99 16389 : std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
100 16389 : em.add_def(eledef);
101 : }
102 :
103 : // PRISM18
104 : {
105 924 : ElementDefinition eledef(PRISM18, 13, 3, 18);
106 16389 : const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,12,14,13,15,17,16};
107 16389 : std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes); // swap trick
108 16389 : em.add_def(eledef);
109 : }
110 :
111 16389 : return em;
112 0 : }
113 :
114 :
115 :
116 0 : GmshIO::GmshIO (const MeshBase & mesh) :
117 : MeshOutput<MeshBase>(mesh),
118 0 : _binary(false),
119 0 : _write_lower_dimensional_elements(true)
120 : {
121 0 : }
122 :
123 :
124 :
125 213 : GmshIO::GmshIO (MeshBase & mesh) :
126 : MeshInput<MeshBase> (mesh),
127 : MeshOutput<MeshBase> (mesh),
128 201 : _binary (false),
129 213 : _write_lower_dimensional_elements(true)
130 : {
131 213 : }
132 :
133 :
134 :
135 0 : bool & GmshIO::binary ()
136 : {
137 0 : return _binary;
138 : }
139 :
140 :
141 :
142 0 : bool & GmshIO::write_lower_dimensional_elements ()
143 : {
144 0 : return _write_lower_dimensional_elements;
145 : }
146 :
147 :
148 :
149 36 : void GmshIO::read (const std::string & name)
150 : {
151 42 : std::ifstream in (name.c_str());
152 36 : this->read_mesh (in);
153 34 : }
154 :
155 :
156 :
157 36 : void GmshIO::read_mesh(std::istream & in)
158 : {
159 : // This is a serial-only process for now;
160 : // the Mesh should be read on processor 0 and
161 : // broadcast later
162 3 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
163 :
164 3 : libmesh_assert(in.good());
165 :
166 : // clear any data in the mesh
167 36 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
168 36 : mesh.clear();
169 :
170 : // some variables
171 36 : int format=0, size=0;
172 36 : 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 6 : 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 6 : 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 6 : 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 6 : 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 6 : std::map<std::pair<unsigned, int>, BoundingBox> entity_to_bounding_box;
197 :
198 : // For reading the file line by line
199 6 : std::string s;
200 :
201 : while (true)
202 : {
203 : // Try to read something. This may set EOF!
204 348 : std::getline(in, s);
205 :
206 377 : if (in)
207 : {
208 : // Process s...
209 :
210 324 : if (s.find("$MeshFormat") == static_cast<std::string::size_type>(0))
211 : {
212 36 : 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 36 : libmesh_error_msg_if(version < 2.0, "Error: Unknown msh file version " << version);
230 36 : libmesh_error_msg_if(format, "Error: Unknown data format for mesh in Gmsh reader.");
231 : }
232 :
233 : // Read and process the "PhysicalNames" section.
234 288 : 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 36 : unsigned int num_physical_groups = 0;
244 3 : in >> num_physical_groups;
245 :
246 : // Read rest of line including newline character.
247 36 : std::getline(in, s);
248 :
249 348 : for (unsigned int i=0; i<num_physical_groups; ++i)
250 : {
251 : // Read an entire line of the PhysicalNames section.
252 312 : 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 364 : std::istringstream s_stream(s);
258 : unsigned phys_dim;
259 : int phys_id;
260 52 : std::string phys_name;
261 312 : 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 312 : 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 390 : 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 312 : if (s.find("lower_dimensional_block") != std::string::npos)
276 : {
277 39 : 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 36 : mesh.subdomain_name(cast_int<subdomain_id_type>(phys_id)) = phys_name;
283 : }
284 260 : }
285 : }
286 :
287 252 : else if (s.find("$Entities") == static_cast<std::string::size_type>(0))
288 : {
289 36 : if (version >= 4.0)
290 : {
291 : std::size_t num_point_entities, num_curve_entities, num_surface_entities, num_volume_entities;
292 3 : in >> num_point_entities >> num_curve_entities >> num_surface_entities >> num_volume_entities;
293 :
294 240 : 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 204 : in >> point_tag >> x >> y >> z >> num_physical_tags;
300 :
301 204 : 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 307 : entity_to_bounding_box[std::make_pair(0, point_tag)] =
306 221 : BoundingBox(Point(x,y,z),Point(x,y,z));
307 :
308 204 : if (num_physical_tags)
309 : {
310 36 : in >> physical_tag;
311 39 : entity_to_physical_id[std::make_pair(0, point_tag)] = physical_tag;
312 : }
313 : }
314 360 : 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 324 : in >> curve_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
320 :
321 324 : 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 324 : entity_to_bounding_box[std::make_pair(1, curve_tag)] =
326 351 : BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
327 :
328 324 : if (num_physical_tags)
329 : {
330 60 : in >> physical_tag;
331 65 : 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 324 : std::getline(in, s);
336 : }
337 372 : 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 336 : in >> surface_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
343 :
344 336 : 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 336 : entity_to_bounding_box[std::make_pair(2, surface_tag)] =
349 364 : BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
350 :
351 336 : if (num_physical_tags)
352 : {
353 204 : in >> physical_tag;
354 221 : 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 336 : std::getline(in, s);
359 : }
360 72 : 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 36 : in >> volume_tag >> minx >> miny >> minz >> maxx >> maxy >> maxz >> num_physical_tags;
366 :
367 36 : 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 36 : entity_to_bounding_box[std::make_pair(3, volume_tag)] =
372 39 : BoundingBox(Point(minx,miny,minz),Point(maxx,maxy,maxz));
373 :
374 36 : if (num_physical_tags)
375 : {
376 36 : in >> physical_tag;
377 39 : 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 36 : std::getline(in, s);
382 : }
383 : // Read the $EndEntities
384 36 : std::getline(in, s);
385 : } // end if (version >= 4.0)
386 :
387 : else
388 0 : libmesh_error_msg("The Entities block was introduced in mesh version 4.0");
389 : } // end if $Entities
390 :
391 : // read the node block
392 414 : else if (s.find("$NOD") == static_cast<std::string::size_type>(0) ||
393 432 : s.find("$NOE") == static_cast<std::string::size_type>(0) ||
394 198 : s.find("$Nodes") == static_cast<std::string::size_type>(0))
395 : {
396 36 : if (version < 4.0)
397 : {
398 0 : unsigned int num_nodes = 0;
399 0 : in >> num_nodes;
400 0 : 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 0 : for (unsigned int i=0; i<num_nodes; ++i)
408 : {
409 0 : in >> id >> x >> y >> z;
410 0 : mesh.add_point (Point(x, y, z), i);
411 0 : nodetrans[id] = i;
412 : }
413 : }
414 : else
415 : {
416 : // Read numEntityBlocks line
417 36 : std::size_t num_entities = 0, num_nodes = 0, min_node_tag, max_node_tag;
418 3 : in >> num_entities >> num_nodes >> min_node_tag >> max_node_tag;
419 :
420 36 : mesh.reserve_nodes(num_nodes);
421 :
422 3 : std::size_t node_counter = 0;
423 :
424 : // Now loop over entities
425 588 : for (std::size_t i = 0; i < num_entities; ++i)
426 : {
427 : int entity_dim, entity_tag, parametric;
428 552 : std::size_t num_nodes_in_block = 0;
429 552 : in >> entity_dim >> entity_tag >> parametric >> num_nodes_in_block;
430 552 : 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 1908 : for (std::size_t n = 0; n < num_nodes_in_block; ++n)
435 : {
436 :
437 113 : in >> gmsh_id;
438 1356 : 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 552 : for (std::size_t libmesh_id = node_counter - num_nodes_in_block;
444 1908 : libmesh_id < node_counter;
445 : ++libmesh_id)
446 : {
447 113 : in >> x >> y >> z;
448 1469 : mesh.add_point(Point(x, y, z), libmesh_id);
449 : }
450 : }
451 : }
452 : // read the $ENDNOD delimiter
453 36 : std::getline(in, s);
454 : }
455 :
456 : // Read the element block
457 360 : else if (s.find("$ELM") == static_cast<std::string::size_type>(0) ||
458 165 : s.find("$Elements") == static_cast<std::string::size_type>(0))
459 : {
460 : // Keep track of element dimensions seen
461 50 : std::vector<unsigned> elem_dimensions_seen(3);
462 :
463 36 : if (version < 4.0)
464 : {
465 : // For reading the number of elements and the node ids from the stream
466 : unsigned int
467 0 : num_elem = 0,
468 0 : node_id = 0;
469 :
470 : // read how many elements are there, and reserve space in the mesh
471 0 : in >> num_elem;
472 0 : 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 0 : for (unsigned int iel=0; iel<num_elem; ++iel)
491 : {
492 : unsigned int
493 : id, type,
494 0 : physical=1, elementary=1,
495 0 : 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 0 : if (version <= 1.0)
502 0 : in >> id >> type >> physical >> elementary >> nnodes;
503 :
504 : else
505 : {
506 0 : in >> id >> type >> ntags;
507 :
508 0 : if (ntags > 2)
509 0 : libmesh_do_once(libMesh::err << "Warning, ntags=" << ntags << ", but we currently only support reading 2 flags." << std::endl;);
510 :
511 0 : for (unsigned int j = 0; j < ntags; j++)
512 : {
513 0 : in >> tag;
514 0 : if (j == 0)
515 0 : physical = tag;
516 0 : else if (j == 1)
517 0 : elementary = tag;
518 : }
519 : }
520 :
521 : // Get a reference to the ElementDefinition, throw an error if not found.
522 : const GmshIO::ElementDefinition & eletype =
523 0 : libmesh_map_find(_element_maps.in, type);
524 :
525 : // If we read nnodes, make sure it matches the number in eletype.nnodes
526 0 : 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 0 : 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 0 : 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 0 : elem_dimensions_seen[eletype.dim-1] = 1;
543 :
544 : // Add the element to the mesh
545 : {
546 : Elem * elem =
547 0 : mesh.add_elem(Elem::build_with_id(eletype.type, iel));
548 :
549 : // Make sure that the libmesh element we added has nnodes nodes.
550 0 : 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 0 : if (eletype.nodes.size() > 0)
562 0 : for (unsigned int i=0; i<nnodes; i++)
563 : {
564 0 : in >> node_id;
565 0 : elem->set_node(eletype.nodes[i], mesh.node_ptr(nodetrans[node_id]));
566 : }
567 : else
568 : {
569 0 : for (unsigned int i=0; i<nnodes; i++)
570 : {
571 0 : in >> node_id;
572 0 : 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 0 : 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 0 : in >> node_id;
592 0 : mesh.get_boundary_info().add_node
593 0 : (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 36 : std::size_t num_entity_blocks = 0, num_elem = 0, min_element_tag, max_element_tag;
602 :
603 : // Read entity information
604 3 : in >> num_entity_blocks >> num_elem >> min_element_tag >> max_element_tag;
605 :
606 36 : mesh.reserve_elem(num_elem);
607 :
608 3 : std::size_t iel = 0;
609 :
610 : // Loop over entity blocks
611 156 : for (std::size_t i = 0; i < num_entity_blocks; ++i)
612 : {
613 : int entity_dim, entity_tag;
614 : unsigned int element_type;
615 132 : std::size_t num_elems_in_block = 0;
616 132 : 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 132 : 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 132 : 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 130 : elem_dimensions_seen[eletype.dim-1] = 1;
633 :
634 : // Loop over elements with dim > 0
635 456 : for (std::size_t n = 0; n < num_elems_in_block; ++n)
636 : {
637 : Elem * elem =
638 364 : mesh.add_elem(Elem::build_with_id(eletype.type, iel++));
639 :
640 : std::size_t gmsh_element_id;
641 28 : 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 308 : BoundingBox expected_bounding_box;
649 28 : std::swap(expected_bounding_box.min(),
650 : expected_bounding_box.max());
651 :
652 336 : if (auto it = entity_to_bounding_box.find
653 344 : (std::make_pair(entity_dim, entity_tag));
654 28 : it != entity_to_bounding_box.end())
655 28 : expected_bounding_box = it->second;
656 :
657 : // Get the remainder of the line that includes the nodes ids
658 336 : std::getline(in, s);
659 336 : std::istringstream is(s);
660 28 : std::size_t local_node_counter = 0, gmsh_node_id;
661 2704 : while (is >> gmsh_node_id)
662 : {
663 2160 : 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 2160 : 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 2160 : if (eletype.nodes.size() > 0)
683 720 : elem->set_node(eletype.nodes[local_node_counter++],
684 120 : node);
685 : else
686 1440 : elem->set_node(local_node_counter++, node);
687 : }
688 :
689 : // Make sure that the libmesh element we added has nnodes nodes.
690 336 : 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 364 : elem->subdomain_id() = static_cast<subdomain_id_type>(
702 364 : entity_to_physical_id[std::make_pair(entity_dim, entity_tag)]);
703 :
704 280 : } // end for (loop over elements in entity block)
705 : } // end if (eletype.dim > 0)
706 :
707 : else
708 : {
709 12 : for (std::size_t n = 0; n < num_elems_in_block; ++n)
710 : {
711 : std::size_t gmsh_element_id, gmsh_node_id;
712 1 : in >> gmsh_element_id;
713 1 : 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 11 : BoundingBox expected_bounding_box;
721 1 : std::swap(expected_bounding_box.min(),
722 : expected_bounding_box.max());
723 :
724 12 : if (auto it = entity_to_bounding_box.find
725 12 : (std::make_pair(entity_dim, entity_tag));
726 1 : it != entity_to_bounding_box.end())
727 1 : expected_bounding_box = it->second;
728 :
729 13 : 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 47 : 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 :
743 0 : mesh.get_boundary_info().add_node(
744 : node,
745 : static_cast<boundary_id_type>(entity_to_physical_id[
746 10 : 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 24 : std::getline(in, s);
754 :
755 : // Record the max and min element dimension seen while reading the file.
756 : unsigned char
757 24 : max_elem_dimension_seen=1,
758 24 : min_elem_dimension_seen=3;
759 :
760 96 : for (auto i : index_range(elem_dimensions_seen))
761 72 : if (elem_dimensions_seen[i])
762 : {
763 : // Debugging
764 : // libMesh::out << "Seen elements of dimension " << i+1 << std::endl;
765 60 : max_elem_dimension_seen =
766 65 : std::max(max_elem_dimension_seen, cast_int<unsigned char>(i+1));
767 60 : min_elem_dimension_seen =
768 82 : 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 2 : 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 44 : 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 120 : for (const auto & pr : gmsh_physicals)
790 : {
791 : // Extract data
792 96 : auto [phys_id, phys_dim] = pr.first;
793 96 : 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 96 : if (phys_dim == max_elem_dimension_seen)
798 36 : 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 60 : else if (phys_dim == 0)
802 0 : 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 65 : else if (phys_dim < max_elem_dimension_seen &&
807 65 : !lower_dimensional_blocks.count(cast_int<boundary_id_type>(phys_id)))
808 60 : mesh.get_boundary_info().sideset_name(cast_int<boundary_id_type>(phys_id)) = phys_name;
809 : }
810 :
811 24 : 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 4 : provide_container_t provide_bcs;
827 :
828 : // 1st loop over active elements - get info about lower-dimensional elements.
829 662 : for (auto & elem : mesh.active_element_ptr_range())
830 344 : if (elem->dim() < max_elem_dimension_seen &&
831 96 : !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 528 : for (auto n : elem->node_index_range())
840 468 : mesh.get_boundary_info().add_node(elem->node_id(n),
841 432 : elem->subdomain_id());
842 :
843 : // Store this elem in a quickly-searchable
844 : // container to use it to assign boundary
845 : // conditions later.
846 96 : provide_bcs.emplace(elem->key(), elem);
847 20 : }
848 :
849 : // 2nd loop over active elements - use lower dimensional element data to set BCs for higher dimensional elements
850 662 : for (auto & elem : mesh.active_element_ptr_range())
851 336 : 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 1052 : for (auto sn : elem->side_index_range())
862 864 : 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 78 : std::unique_ptr<Elem> side (elem->build_side_ptr(sn));
867 :
868 : // Construct the lower-dimensional element to compare to the side.
869 72 : Elem * lower_dim_elem = pr.second;
870 :
871 : // This was a hash, so it might not be perfect. Let's verify...
872 72 : 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 72 : boundary_id_type bid = cast_int<boundary_id_type>(lower_dim_elem->subdomain_id());
879 72 : mesh.get_boundary_info().add_side(elem, sn, bid);
880 : }
881 60 : }
882 20 : }
883 :
884 : // 3rd loop over active elements - Remove the lower-dimensional elements
885 662 : for (auto & elem : mesh.active_element_ptr_range())
886 344 : if (elem->dim() < max_elem_dimension_seen &&
887 96 : !lower_dimensional_blocks.count(elem->subdomain_id()))
888 116 : mesh.delete_elem(elem);
889 : } // end if (n_dims_seen > 1)
890 : } // if $ELM
891 :
892 26 : continue;
893 260 : } // if (in)
894 :
895 :
896 : // If !in, check to see if EOF was set. If so, break out
897 : // of while loop.
898 24 : if (in.eof())
899 2 : break;
900 :
901 : // If !in and !in.eof(), stream is in a bad state!
902 0 : libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
903 :
904 312 : } // while true
905 24 : }
906 :
907 :
908 :
909 0 : void GmshIO::write (const std::string & name)
910 : {
911 0 : if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
912 : {
913 : // Open the output file stream
914 0 : std::ofstream out_stream (name.c_str());
915 :
916 : // Make sure it opened correctly
917 0 : if (!out_stream.good())
918 0 : libmesh_file_error(name.c_str());
919 :
920 0 : this->write_mesh (out_stream);
921 0 : }
922 0 : }
923 :
924 :
925 :
926 0 : void GmshIO::write_nodal_data (const std::string & fname,
927 : const std::vector<Number> & soln,
928 : const std::vector<std::string> & names)
929 : {
930 0 : LOG_SCOPE("write_nodal_data()", "GmshIO");
931 :
932 0 : if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
933 0 : this->write_post (fname, &soln, &names);
934 0 : }
935 :
936 :
937 :
938 0 : void GmshIO::write_mesh (std::ostream & out_stream)
939 : {
940 : // Be sure that the stream is valid.
941 0 : libmesh_assert (out_stream.good());
942 :
943 : // Get a const reference to the mesh
944 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
945 :
946 : // If requested, write out lower-dimensional elements for
947 : // element-side-based boundary conditions.
948 0 : std::size_t n_boundary_faces = 0;
949 0 : if (this->write_lower_dimensional_elements())
950 0 : 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 0 : out_stream << "$MeshFormat\n";
956 0 : out_stream << "2.0 0 " << sizeof(Real) << '\n';
957 0 : out_stream << "$EndMeshFormat\n";
958 :
959 : // write the nodes in (n x y z) format
960 0 : out_stream << "$Nodes\n";
961 0 : out_stream << mesh.n_nodes() << '\n';
962 :
963 0 : for (auto v : make_range(mesh.n_nodes()))
964 0 : out_stream << mesh.node_ref(v).id()+1 << " "
965 0 : << mesh.node_ref(v)(0) << " "
966 0 : << mesh.node_ref(v)(1) << " "
967 0 : << mesh.node_ref(v)(2) << '\n';
968 0 : out_stream << "$EndNodes\n";
969 :
970 : {
971 : // write the connectivity
972 0 : out_stream << "$Elements\n";
973 0 : out_stream << mesh.n_active_elem() + n_boundary_faces << '\n';
974 :
975 : // loop over the elements
976 0 : for (const auto & elem : mesh.active_element_ptr_range())
977 : {
978 : // Get a reference to the ElementDefinition object
979 : const ElementDefinition & eletype =
980 0 : 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 0 : libmesh_assert_less_equal (eletype.nodes.size(), elem->n_nodes());
985 :
986 : // elements ids are 1 based in Gmsh
987 0 : out_stream << elem->id()+1 << " ";
988 : // element type
989 0 : 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 0 : out_stream << " 3 "
996 0 : << static_cast<unsigned int>(elem->subdomain_id())
997 0 : << " 0 "
998 0 : << elem->processor_id()+1
999 0 : << " ";
1000 :
1001 : // if there is a node translation table, use it
1002 0 : if (eletype.nodes.size() > 0)
1003 0 : for (auto i : elem->node_index_range())
1004 0 : out_stream << elem->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
1005 : // otherwise keep the same node order
1006 : else
1007 0 : for (auto i : elem->node_index_range())
1008 0 : out_stream << elem->node_id(i)+1 << " "; // gmsh is 1-based
1009 0 : out_stream << "\n";
1010 0 : } // 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 0 : unsigned int e_id = mesh.max_elem_id();
1022 :
1023 : // loop over the elements, writing out boundary faces
1024 0 : if (n_boundary_faces)
1025 : {
1026 : // Build a list of (elem, side, bc) tuples.
1027 0 : auto bc_triples = mesh.get_boundary_info().build_side_list();
1028 :
1029 : // Loop over these lists, writing data to the file.
1030 0 : for (const auto & t : bc_triples)
1031 : {
1032 0 : const Elem & elem = mesh.elem_ref(std::get<0>(t));
1033 :
1034 0 : 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 0 : 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 0 : libmesh_assert_less_equal (eletype.nodes.size(), side->n_nodes());
1043 :
1044 : // elements ids are 1-based in Gmsh
1045 0 : out_stream << e_id+1 << " ";
1046 :
1047 : // element type
1048 0 : 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 0 : out_stream << " 3 "
1055 0 : << std::get<2>(t)
1056 0 : << " 0 "
1057 0 : << elem.processor_id()+1
1058 0 : << " ";
1059 :
1060 : // if there is a node translation table, use it
1061 0 : if (eletype.nodes.size() > 0)
1062 0 : for (auto i : side->node_index_range())
1063 0 : out_stream << side->node_id(eletype.nodes[i])+1 << " "; // gmsh is 1-based
1064 :
1065 : // otherwise keep the same node order
1066 : else
1067 0 : for (auto i : side->node_index_range())
1068 0 : out_stream << side->node_id(i)+1 << " "; // gmsh is 1-based
1069 :
1070 : // Go to the next line
1071 0 : out_stream << "\n";
1072 :
1073 : // increment this index too...
1074 0 : ++e_id;
1075 0 : }
1076 : }
1077 :
1078 0 : out_stream << "$EndElements\n";
1079 : }
1080 0 : }
1081 :
1082 :
1083 :
1084 0 : void GmshIO::write_post (const std::string & fname,
1085 : const std::vector<Number> * v,
1086 : const std::vector<std::string> * solution_names)
1087 : {
1088 :
1089 : // Should only do this on processor 0!
1090 0 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
1091 :
1092 : // Create an output stream
1093 0 : std::ofstream out_stream(fname.c_str());
1094 :
1095 : // Make sure it opened correctly
1096 0 : if (!out_stream.good())
1097 0 : 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 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1104 :
1105 : // write the data
1106 0 : if ((solution_names != nullptr) && (v != nullptr))
1107 : {
1108 : const unsigned int n_vars =
1109 0 : cast_int<unsigned int>(solution_names->size());
1110 :
1111 0 : if (!(v->size() == mesh.n_nodes()*n_vars))
1112 0 : libMesh::err << "ERROR: v->size()=" << v->size()
1113 0 : << ", mesh.n_nodes()=" << mesh.n_nodes()
1114 0 : << ", n_vars=" << n_vars
1115 0 : << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
1116 0 : << "\n";
1117 :
1118 0 : libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
1119 :
1120 : // write the header
1121 0 : out_stream << "$PostFormat\n";
1122 0 : if (this->binary())
1123 0 : out_stream << "1.2 1 " << sizeof(double) << "\n";
1124 : else
1125 0 : out_stream << "1.2 0 " << sizeof(double) << "\n";
1126 0 : out_stream << "$EndPostFormat\n";
1127 :
1128 : // Loop over the elements to see how much of each type there are
1129 0 : unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
1130 0 : n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
1131 0 : unsigned int n_scalar=0, n_vector=0, n_tensor=0;
1132 0 : unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;
1133 :
1134 : {
1135 0 : for (const auto & elem : mesh.active_element_ptr_range())
1136 : {
1137 0 : const ElemType elemtype = elem->type();
1138 :
1139 0 : switch (elemtype)
1140 : {
1141 0 : case EDGE2:
1142 : case EDGE3:
1143 : case EDGE4:
1144 : {
1145 0 : n_lines += 1;
1146 0 : break;
1147 : }
1148 0 : case TRI3:
1149 : case TRI6:
1150 : {
1151 0 : n_triangles += 1;
1152 0 : break;
1153 : }
1154 0 : case QUAD4:
1155 : case QUAD8:
1156 : case QUAD9:
1157 : {
1158 0 : n_quadrangles += 1;
1159 0 : break;
1160 : }
1161 0 : case TET4:
1162 : case TET10:
1163 : {
1164 0 : n_tetrahedra += 1;
1165 0 : break;
1166 : }
1167 0 : case HEX8:
1168 : case HEX20:
1169 : case HEX27:
1170 : {
1171 0 : n_hexahedra += 1;
1172 0 : break;
1173 : }
1174 0 : case PRISM6:
1175 : case PRISM15:
1176 : case PRISM18:
1177 : {
1178 0 : n_prisms += 1;
1179 0 : break;
1180 : }
1181 0 : case PYRAMID5:
1182 : {
1183 0 : n_pyramids += 1;
1184 0 : break;
1185 : }
1186 0 : default:
1187 0 : libmesh_error_msg("ERROR: Nonexistent element type " << Utility::enum_to_string(elem->type()));
1188 : }
1189 0 : }
1190 : }
1191 :
1192 : // create a view for each variable
1193 0 : for (unsigned int ivar=0; ivar < n_vars; ivar++)
1194 : {
1195 0 : 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 0 : n_scalar = 1;
1201 :
1202 : // write the variable as a view, and the number of time steps
1203 0 : out_stream << "$View\n" << varname << " " << 1 << "\n";
1204 :
1205 : // write how many of each geometry type are written
1206 0 : out_stream << n_points * n_scalar << " "
1207 0 : << n_points * n_vector << " "
1208 0 : << n_points * n_tensor << " "
1209 0 : << n_lines * n_scalar << " "
1210 0 : << n_lines * n_vector << " "
1211 0 : << n_lines * n_tensor << " "
1212 0 : << n_triangles * n_scalar << " "
1213 0 : << n_triangles * n_vector << " "
1214 0 : << n_triangles * n_tensor << " "
1215 0 : << n_quadrangles * n_scalar << " "
1216 0 : << n_quadrangles * n_vector << " "
1217 0 : << n_quadrangles * n_tensor << " "
1218 0 : << n_tetrahedra * n_scalar << " "
1219 0 : << n_tetrahedra * n_vector << " "
1220 0 : << n_tetrahedra * n_tensor << " "
1221 0 : << n_hexahedra * n_scalar << " "
1222 0 : << n_hexahedra * n_vector << " "
1223 0 : << n_hexahedra * n_tensor << " "
1224 0 : << n_prisms * n_scalar << " "
1225 0 : << n_prisms * n_vector << " "
1226 0 : << n_prisms * n_tensor << " "
1227 0 : << n_pyramids * n_scalar << " "
1228 0 : << n_pyramids * n_vector << " "
1229 0 : << n_pyramids * n_tensor << " "
1230 0 : << nb_text2d << " "
1231 0 : << nb_text2d_chars << " "
1232 0 : << nb_text3d << " "
1233 0 : << nb_text3d_chars << "\n";
1234 :
1235 : // if binary, write a marker to identify the endianness of the file
1236 0 : if (this->binary())
1237 : {
1238 0 : const int one = 1;
1239 0 : std::memcpy(buf, &one, sizeof(int));
1240 0 : out_stream.write(buf, sizeof(int));
1241 : }
1242 :
1243 : // the time steps (there is just 1 at the moment)
1244 0 : if (this->binary())
1245 : {
1246 0 : double one = 1;
1247 0 : std::memcpy(buf, &one, sizeof(double));
1248 0 : out_stream.write(buf, sizeof(double));
1249 : }
1250 : else
1251 0 : out_stream << "1\n";
1252 :
1253 : // Loop over the elements and write out the data
1254 0 : for (const auto & elem : mesh.active_element_ptr_range())
1255 : {
1256 0 : const unsigned int nv = elem->n_vertices();
1257 : // this is quite crappy, but I did not invent that file format!
1258 0 : for (unsigned int d=0; d<3; d++) // loop over the dimensions
1259 : {
1260 0 : for (unsigned int n=0; n < nv; n++) // loop over vertices
1261 : {
1262 0 : const Point & vertex = elem->point(n);
1263 0 : 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 0 : double tmp = double(vertex(d));
1269 0 : std::memcpy(buf, &tmp, sizeof(double));
1270 0 : out_stream.write(reinterpret_cast<char *>(buf), sizeof(double));
1271 : }
1272 : else
1273 0 : out_stream << vertex(d) << " ";
1274 : }
1275 0 : if (!this->binary())
1276 0 : out_stream << "\n";
1277 : }
1278 :
1279 : // now finally write out the data
1280 0 : for (unsigned int i=0; i < nv; i++) // loop over vertices
1281 0 : 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 0 : double tmp = double(libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]));
1289 0 : std::memcpy(buf, &tmp, sizeof(double));
1290 0 : 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 0 : out_stream << libmesh_real((*v)[elem->node_id(i)*n_vars + ivar]) << "\n";
1300 : }
1301 0 : }
1302 0 : if (this->binary())
1303 0 : out_stream << "\n";
1304 0 : out_stream << "$EndView\n";
1305 :
1306 : } // end variable loop (writing the views)
1307 : }
1308 0 : }
1309 :
1310 : } // namespace libMesh
|
