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 : // C++ includes
19 : #include <iomanip>
20 : #include <fstream>
21 : #include <vector>
22 : #include <cctype> // isspace
23 :
24 : // Local includes
25 : #include "libmesh/libmesh_config.h"
26 : #include "libmesh/libmesh_logging.h"
27 : #include "libmesh/gmv_io.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/elem.h"
30 : #include "libmesh/equation_systems.h"
31 : #include "libmesh/numeric_vector.h"
32 : #include "libmesh/enum_to_string.h"
33 : #include "libmesh/enum_io_package.h"
34 : #include "libmesh/enum_elem_type.h"
35 : #include "libmesh/int_range.h"
36 : #include "libmesh/utility.h"
37 :
38 : // Wrap everything in a GMVLib namespace and
39 : // use extern "C" to avoid name mangling.
40 : #ifdef LIBMESH_HAVE_GMV
41 : namespace GMVLib
42 : {
43 : extern "C"
44 : {
45 : #include "gmvread.h"
46 : }
47 : }
48 : #endif
49 :
50 : // anonymous namespace to hold local data
51 : namespace
52 : {
53 : using namespace libMesh;
54 :
55 : /**
56 : * Defines mapping from libMesh element types to GMV element types.
57 : * Note: Not all of the GMV element types have an identity mapping
58 : * to libmesh node numbering, but the node mappings do all happen to
59 : * be their own inverse, that is, pairs of nodes are simply swapped
60 : * between the two definitions. Therefore we need only one node map
61 : * for both reading and writing.
62 : */
63 0 : struct ElementDefinition {
64 : // GMV element name
65 : std::string label;
66 :
67 : // Used to map libmesh nodes to GMV for writing
68 : std::vector<unsigned> node_map;
69 : };
70 :
71 :
72 : // maps from a libMesh element type to the proper GMV
73 : // ElementDefinition. Placing the data structure here in this
74 : // anonymous namespace gives us the benefits of a global variable
75 : // without the nasty side-effects.
76 : std::map<ElemType, ElementDefinition> eletypes;
77 :
78 : // Helper function to fill up eletypes map
79 0 : void add_eletype_entry(ElemType libmesh_elem_type,
80 : const unsigned * node_map,
81 : const std::string & gmv_label,
82 : unsigned nodes_size )
83 : {
84 : // If map entry does not exist, this will create it
85 0 : ElementDefinition & map_entry = eletypes[libmesh_elem_type];
86 :
87 : // Set the label
88 0 : map_entry.label = gmv_label;
89 :
90 : // Use the "swap trick" from Scott Meyer's "Effective STL" to swap
91 : // an unnamed temporary vector into the map_entry's vector. Note:
92 : // the vector(iter, iter) constructor is used.
93 0 : std::vector<unsigned int>(node_map,
94 0 : node_map+nodes_size).swap(map_entry.node_map);
95 0 : }
96 :
97 :
98 : // ------------------------------------------------------------
99 : // helper function to initialize the eletypes map
100 0 : void init_eletypes ()
101 : {
102 0 : if (eletypes.empty())
103 : {
104 : // This should happen only once. The first time this method
105 : // is called the eletypes data structure will be empty, and
106 : // we will fill it. Any subsequent calls will find an initialized
107 : // eletypes map and will do nothing.
108 :
109 : // EDGE2
110 : {
111 0 : const unsigned int node_map[] = {0,1};
112 0 : add_eletype_entry(EDGE2, node_map, "line 2", 2);
113 : }
114 :
115 : // LINE3
116 : {
117 0 : const unsigned int node_map[] = {0,1,2};
118 0 : add_eletype_entry(EDGE3, node_map, "3line 3", 3);
119 : }
120 :
121 : // TRI3
122 : {
123 0 : const unsigned int node_map[] = {0,1,2};
124 0 : add_eletype_entry(TRI3, node_map, "tri3 3", 3);
125 : }
126 :
127 : // TRI6
128 : {
129 0 : const unsigned int node_map[] = {0,1,2,3,4,5};
130 0 : add_eletype_entry(TRI6, node_map, "6tri 6", 6);
131 : }
132 :
133 : // QUAD4
134 : {
135 0 : const unsigned int node_map[] = {0,1,2,3};
136 0 : add_eletype_entry(QUAD4, node_map, "quad 4", 4);
137 : }
138 :
139 : // QUAD8, QUAD9
140 : {
141 0 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7};
142 0 : add_eletype_entry(QUAD8, node_map, "8quad 8", 8);
143 :
144 : // QUAD9 was not supported by GMV but it gets the same entry, even the label (is that correct?)
145 0 : eletypes[QUAD9] = eletypes[QUAD8];
146 : }
147 :
148 : // HEX8
149 : {
150 0 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7};
151 0 : add_eletype_entry(HEX8, node_map, "phex8 8", 8);
152 : }
153 :
154 : // HEX20, HEX27
155 : {
156 : // Note: This map is its own inverse
157 0 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15};
158 0 : add_eletype_entry(HEX20, node_map, "phex20 20", 20);
159 :
160 : // HEX27 was not supported by GMV but it gets the same entry, even the label (is that correct?)
161 0 : eletypes[HEX27] = eletypes[HEX20];
162 : }
163 :
164 : // TET4
165 : {
166 : // This is correct, see write_ascii_old_impl() to confirm.
167 : // This map is also its own inverse.
168 0 : const unsigned node_map[] = {0,2,1,3};
169 0 : add_eletype_entry(TET4, node_map, "tet 4", 4);
170 : }
171 :
172 : // TET10
173 : {
174 0 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,9};
175 0 : add_eletype_entry(TET10, node_map, "ptet10 10", 10);
176 : }
177 :
178 : // PRISM6
179 : {
180 0 : const unsigned int node_map[] = {0,1,2,3,4,5};
181 0 : add_eletype_entry(PRISM6, node_map, "pprism6 6", 6);
182 : }
183 :
184 : // PRISM15, PRISM18
185 : {
186 : // Note: This map is its own inverse
187 0 : const unsigned int node_map[] = {0,1,2,3,4,5,6,7,8,12,13,14, 9,10,11};
188 0 : add_eletype_entry(PRISM15, node_map, "pprism15 15", 15);
189 :
190 : // PRISM18 was not supported by GMV but it gets the same entry, even the label (is that correct?)
191 0 : eletypes[PRISM18] = eletypes[PRISM15];
192 : }
193 : //==============================
194 : }
195 0 : }
196 :
197 : } // end anonymous namespace
198 :
199 :
200 : namespace libMesh
201 : {
202 :
203 : // Initialize the static data members by calling the static build functions.
204 : std::map<std::string, ElemType> GMVIO::_reading_element_map = GMVIO::build_reading_element_map();
205 :
206 :
207 :
208 : // Static function used to build the _reading_element_map.
209 16389 : std::map<std::string, ElemType> GMVIO::build_reading_element_map()
210 : {
211 462 : std::map<std::string, ElemType> ret;
212 :
213 16389 : ret["line"] = EDGE2;
214 16389 : ret["tri"] = TRI3;
215 16389 : ret["tri3"] = TRI3;
216 16389 : ret["quad"] = QUAD4;
217 16389 : ret["tet"] = TET4;
218 16389 : ret["ptet4"] = TET4;
219 16389 : ret["hex"] = HEX8;
220 16389 : ret["phex8"] = HEX8;
221 16389 : ret["prism"] = PRISM6;
222 16389 : ret["pprism6"] = PRISM6;
223 16389 : ret["phex20"] = HEX20;
224 16389 : ret["phex27"] = HEX27;
225 16389 : ret["pprism15"] = PRISM15;
226 16389 : ret["ptet10"] = TET10;
227 16389 : ret["6tri"] = TRI6;
228 16389 : ret["8quad"] = QUAD8;
229 16389 : ret["3line"] = EDGE3;
230 :
231 : // Unsupported/Unused types
232 : // ret["vface2d"]
233 : // ret["vface3d"]
234 : // ret["pyramid"]
235 : // ret["ppyrmd5"]
236 : // ret["ppyrmd13"]
237 :
238 16389 : return ret;
239 : }
240 :
241 :
242 0 : GMVIO::GMVIO (const MeshBase & mesh) :
243 : MeshOutput<MeshBase> (mesh),
244 0 : _binary (false),
245 0 : _discontinuous (false),
246 0 : _partitioning (true),
247 0 : _write_subdomain_id_as_material (false),
248 0 : _subdivide_second_order (true),
249 0 : _p_levels (true),
250 0 : _next_elem_id (0)
251 : {
252 0 : }
253 :
254 :
255 :
256 35664 : GMVIO::GMVIO (MeshBase & mesh) :
257 : MeshInput<MeshBase> (mesh),
258 : MeshOutput<MeshBase>(mesh),
259 33196 : _binary (false),
260 33196 : _discontinuous (false),
261 33196 : _partitioning (true),
262 33196 : _write_subdomain_id_as_material (false),
263 33196 : _subdivide_second_order (true),
264 33196 : _p_levels (true),
265 35664 : _next_elem_id (0)
266 : {
267 35664 : }
268 :
269 :
270 :
271 0 : void GMVIO::write (const std::string & fname)
272 : {
273 0 : if (this->binary())
274 0 : this->write_binary (fname);
275 : else
276 0 : this->write_ascii_old_impl (fname);
277 0 : }
278 :
279 :
280 :
281 31718 : void GMVIO::write_nodal_data (const std::string & fname,
282 : const std::vector<Number> & soln,
283 : const std::vector<std::string> & names)
284 : {
285 2096 : LOG_SCOPE("write_nodal_data()", "GMVIO");
286 :
287 31718 : if (this->binary())
288 0 : this->write_binary (fname, &soln, &names);
289 : else
290 31718 : this->write_ascii_old_impl (fname, &soln, &names);
291 31718 : }
292 :
293 :
294 :
295 0 : void GMVIO::write_ascii_new_impl (const std::string & fname,
296 : const std::vector<Number> * v,
297 : const std::vector<std::string> * solution_names)
298 : {
299 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
300 :
301 0 : libMesh::err << "WARNING: GMVIO::write_ascii_new_impl() not infinite-element aware!"
302 0 : << std::endl;
303 0 : libmesh_here();
304 :
305 : // Set it to our current precision
306 0 : this->write_ascii_old_impl (fname, v, solution_names);
307 :
308 : #else
309 :
310 : // Get a reference to the mesh
311 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
312 :
313 : // This is a parallel_only function
314 0 : const unsigned int n_active_elem = mesh.n_active_elem();
315 :
316 0 : if (MeshOutput<MeshBase>::mesh().processor_id() != 0)
317 0 : return;
318 :
319 : // Open the output file stream
320 0 : std::ofstream out_stream (fname.c_str());
321 :
322 0 : out_stream << std::setprecision(this->ascii_precision());
323 :
324 : // Make sure it opened correctly
325 0 : if (!out_stream.good())
326 0 : libmesh_file_error(fname.c_str());
327 :
328 0 : unsigned int mesh_max_p_level = 0;
329 :
330 : // Begin interfacing with the GMV data file
331 : {
332 0 : out_stream << "gmvinput ascii\n\n";
333 :
334 : // write the nodes
335 0 : out_stream << "nodes " << mesh.n_nodes() << "\n";
336 0 : for (auto n : make_range(mesh.n_nodes()))
337 0 : out_stream << mesh.point(n)(0) << " ";
338 0 : out_stream << "\n";
339 :
340 0 : for (auto n : make_range(mesh.n_nodes()))
341 : #if LIBMESH_DIM > 1
342 0 : out_stream << mesh.point(n)(1) << " ";
343 : #else
344 : out_stream << 0. << " ";
345 : #endif
346 0 : out_stream << "\n";
347 :
348 0 : for (auto n : make_range(mesh.n_nodes()))
349 : #if LIBMESH_DIM > 2
350 0 : out_stream << mesh.point(n)(2) << " ";
351 : #else
352 : out_stream << 0. << " ";
353 : #endif
354 0 : out_stream << "\n\n";
355 : }
356 :
357 : {
358 : // write the connectivity
359 0 : out_stream << "cells " << n_active_elem << "\n";
360 :
361 : // initialize the eletypes map (eletypes is a file-global variable)
362 0 : init_eletypes();
363 :
364 0 : for (const auto & elem : mesh.active_element_ptr_range())
365 : {
366 0 : mesh_max_p_level = std::max(mesh_max_p_level,
367 0 : elem->p_level());
368 :
369 : // Make sure we have a valid entry for
370 : // the current element type.
371 : libmesh_assert (eletypes.count(elem->type()));
372 :
373 0 : const ElementDefinition & ele = eletypes[elem->type()];
374 :
375 : // The element mapper better not require any more nodes
376 : // than are present in the current element!
377 : libmesh_assert_less_equal (ele.node_map.size(), elem->n_nodes());
378 :
379 0 : out_stream << ele.label << "\n";
380 0 : for (auto i : index_range(ele.node_map))
381 0 : out_stream << elem->node_id(ele.node_map[i])+1 << " ";
382 0 : out_stream << "\n";
383 0 : }
384 0 : out_stream << "\n";
385 : }
386 :
387 : // optionally write the partition information
388 0 : if (this->partitioning())
389 : {
390 0 : libmesh_error_msg_if(this->write_subdomain_id_as_material(),
391 : "Not yet supported in GMVIO::write_ascii_new_impl");
392 :
393 0 : out_stream << "material "
394 : << mesh.n_partitions()
395 : // Note: GMV may give you errors like
396 : // Error, material for cell 1 is greater than 1
397 : // Error, material for cell 2 is greater than 1
398 : // Error, material for cell 3 is greater than 1
399 : // ... because you put the wrong number of partitions here.
400 : // To ensure you write the correct number of materials, call
401 : // mesh.recalculate_n_partitions() before writing out the
402 : // mesh.
403 : // Note: we can't call it now because the Mesh is const here and
404 : // it is a non-const function.
405 0 : << " 0\n";
406 :
407 0 : for (auto proc : make_range(mesh.n_partitions()))
408 0 : out_stream << "proc_" << proc << "\n";
409 :
410 : // FIXME - don't we need to use an ElementDefinition here? - RHS
411 0 : for (const auto & elem : mesh.active_element_ptr_range())
412 0 : out_stream << elem->processor_id()+1 << "\n";
413 0 : out_stream << "\n";
414 : }
415 :
416 : // If there are *any* variables at all in the system (including
417 : // p level, or arbitrary cell-based data)
418 : // to write, the gmv file needs to contain the word "variable"
419 : // on a line by itself.
420 : bool write_variable = false;
421 :
422 : // 1.) p-levels
423 0 : if (this->p_levels() && mesh_max_p_level)
424 : write_variable = true;
425 :
426 : // 2.) solution data
427 0 : if ((solution_names != nullptr) && (v != nullptr))
428 : write_variable = true;
429 :
430 : // 3.) cell-centered data
431 0 : if (!(this->_cell_centered_data.empty()))
432 : write_variable = true;
433 :
434 0 : if (write_variable)
435 0 : out_stream << "variable\n";
436 :
437 : // if ((this->p_levels() && mesh_max_p_level) ||
438 : // ((solution_names != nullptr) && (v != nullptr)))
439 : // out_stream << "variable\n";
440 :
441 : // optionally write the polynomial degree information
442 0 : if (this->p_levels() && mesh_max_p_level)
443 : {
444 0 : out_stream << "p_level 0\n";
445 :
446 0 : for (const auto & elem : mesh.active_element_ptr_range())
447 : {
448 0 : const ElementDefinition & ele = eletypes[elem->type()];
449 :
450 : // The element mapper better not require any more nodes
451 : // than are present in the current element!
452 : libmesh_assert_less_equal (ele.node_map.size(), elem->n_nodes());
453 :
454 0 : for (std::size_t i=0, enms=ele.node_map.size(); i < enms; i++)
455 0 : out_stream << elem->p_level() << " ";
456 0 : }
457 0 : out_stream << "\n\n";
458 : }
459 :
460 :
461 : // optionally write cell-centered data
462 0 : if (!(this->_cell_centered_data.empty()))
463 : {
464 0 : for (const auto & [var_name, the_array] : this->_cell_centered_data)
465 : {
466 : // write out the variable name, followed by a zero.
467 0 : out_stream << var_name << " 0\n";
468 :
469 : // Loop over active elements, write out cell data. If second-order cells
470 : // are split into sub-elements, the sub-elements inherit their parent's
471 : // cell-centered data.
472 0 : for (const auto & elem : mesh.active_element_ptr_range())
473 : {
474 : // Use the element's ID to find the value.
475 : libmesh_assert_less (elem->id(), the_array->size());
476 0 : const Real the_value = (*the_array)[elem->id()];
477 :
478 0 : if (this->subdivide_second_order())
479 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
480 0 : out_stream << the_value << " ";
481 : else
482 0 : out_stream << the_value << " ";
483 0 : }
484 :
485 0 : out_stream << "\n\n";
486 : }
487 : }
488 :
489 :
490 : // optionally write the data
491 0 : if ((solution_names != nullptr) && (v != nullptr))
492 : {
493 0 : const unsigned int n_vars = solution_names->size();
494 :
495 0 : if (!(v->size() == mesh.n_nodes()*n_vars))
496 : libMesh::err << "ERROR: v->size()=" << v->size()
497 0 : << ", mesh.n_nodes()=" << mesh.n_nodes()
498 : << ", n_vars=" << n_vars
499 0 : << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
500 : << "\n";
501 :
502 : libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
503 :
504 0 : for (unsigned int c=0; c<n_vars; c++)
505 : {
506 :
507 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
508 :
509 : // in case of complex data, write _three_ data sets
510 : // for each component
511 :
512 : // this is the real part
513 : out_stream << "r_" << (*solution_names)[c] << " 1\n";
514 :
515 : for (auto n : make_range(mesh.n_nodes()))
516 : out_stream << (*v)[n*n_vars + c].real() << " ";
517 :
518 : out_stream << "\n\n";
519 :
520 : // this is the imaginary part
521 : out_stream << "i_" << (*solution_names)[c] << " 1\n";
522 :
523 : for (auto n : make_range(mesh.n_nodes()))
524 : out_stream << (*v)[n*n_vars + c].imag() << " ";
525 :
526 : out_stream << "\n\n";
527 :
528 : // this is the magnitude
529 : out_stream << "a_" << (*solution_names)[c] << " 1\n";
530 : for (auto n : make_range(mesh.n_nodes()))
531 : out_stream << std::abs((*v)[n*n_vars + c]) << " ";
532 :
533 : out_stream << "\n\n";
534 :
535 : #else
536 :
537 0 : out_stream << (*solution_names)[c] << " 1\n";
538 :
539 0 : for (auto n : make_range(mesh.n_nodes()))
540 0 : out_stream << (*v)[n*n_vars + c] << " ";
541 :
542 0 : out_stream << "\n\n";
543 :
544 : #endif
545 : }
546 :
547 : }
548 :
549 : // If we wrote any variables, we have to close the variable section now
550 0 : if (write_variable)
551 0 : out_stream << "endvars\n";
552 :
553 :
554 : // end of the file
555 0 : out_stream << "\nendgmv\n";
556 :
557 : #endif
558 0 : }
559 :
560 :
561 :
562 :
563 :
564 :
565 31718 : void GMVIO::write_ascii_old_impl (const std::string & fname,
566 : const std::vector<Number> * v,
567 : const std::vector<std::string> * solution_names)
568 : {
569 : // Get a reference to the mesh
570 2096 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
571 :
572 : // Use a MeshSerializer object to gather a parallel mesh before outputting it.
573 : // Note that we cast away constness here (which is bad), but the destructor of
574 : // the MeshSerializer object reparallelizes the Mesh, hopefully keeping it
575 : // "logically const" outside the context of this function...
576 : MeshSerializer serialize(const_cast<MeshBase &>(mesh),
577 32242 : !MeshOutput<MeshBase>::_is_parallel_format);
578 :
579 : // These are parallel_only functions
580 31718 : const dof_id_type n_active_elem = mesh.n_active_elem(),
581 31718 : n_active_sub_elem = mesh.n_active_sub_elem();
582 :
583 32766 : if (MeshOutput<MeshBase>::mesh().processor_id() != 0)
584 1048 : return;
585 :
586 : // Open the output file stream
587 6266 : std::ofstream out_stream (fname.c_str());
588 :
589 : // Set it to our current precision
590 5218 : out_stream << std::setprecision(this->ascii_precision());
591 :
592 : // Make sure it opened correctly
593 5218 : if (!out_stream.good())
594 0 : libmesh_file_error(fname.c_str());
595 :
596 : // Make sure our nodes are contiguous and serialized
597 524 : libmesh_assert_equal_to (mesh.n_nodes(), mesh.max_node_id());
598 :
599 : // libmesh_assert (mesh.is_serial());
600 : // if (!mesh.is_serial())
601 : // {
602 : // if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
603 : // libMesh::err << "Error: GMVIO cannot yet write a DistributedMesh solution"
604 : // << std::endl;
605 : // return;
606 : // }
607 :
608 5218 : unsigned int mesh_max_p_level = 0;
609 :
610 : // Begin interfacing with the GMV data file
611 :
612 : // FIXME - if subdivide_second_order() is off,
613 : // we probably should only be writing the
614 : // vertex nodes - RHS
615 : {
616 : // write the nodes
617 :
618 5218 : out_stream << "gmvinput ascii\n\n";
619 9388 : out_stream << "nodes " << mesh.n_nodes() << '\n';
620 5223851 : for (auto n : make_range(mesh.n_nodes()))
621 5218633 : out_stream << mesh.point(n)(0) << " ";
622 :
623 5218 : out_stream << '\n';
624 :
625 5223851 : for (auto n : make_range(mesh.n_nodes()))
626 : #if LIBMESH_DIM > 1
627 5218633 : out_stream << mesh.point(n)(1) << " ";
628 : #else
629 : out_stream << 0. << " ";
630 : #endif
631 :
632 5218 : out_stream << '\n';
633 :
634 5223851 : for (auto n : make_range(mesh.n_nodes()))
635 : #if LIBMESH_DIM > 2
636 5218633 : out_stream << mesh.point(n)(2) << " ";
637 : #else
638 : out_stream << 0. << " ";
639 : #endif
640 :
641 5218 : out_stream << '\n' << '\n';
642 : }
643 :
644 :
645 :
646 : {
647 : // write the connectivity
648 :
649 5218 : out_stream << "cells ";
650 5218 : if (this->subdivide_second_order())
651 524 : out_stream << n_active_sub_elem;
652 : else
653 0 : out_stream << n_active_elem;
654 5218 : out_stream << '\n';
655 :
656 : // The same temporary storage will be used for each element
657 1048 : std::vector<dof_id_type> conn;
658 :
659 7654918 : for (const auto & elem : mesh.active_element_ptr_range())
660 : {
661 4207984 : mesh_max_p_level = std::max(mesh_max_p_level,
662 4593465 : elem->p_level());
663 :
664 4207984 : switch (elem->dim())
665 : {
666 72 : case 1:
667 : {
668 864 : if (this->subdivide_second_order())
669 2592 : for (auto se : make_range(elem->n_sub_elem()))
670 : {
671 1728 : out_stream << "line 2\n";
672 1728 : elem->connectivity(se, TECPLOT, conn);
673 5184 : for (const auto & idx : conn)
674 3456 : out_stream << idx << " ";
675 :
676 1728 : out_stream << '\n';
677 : }
678 : else
679 : {
680 0 : out_stream << "line 2\n";
681 0 : if (elem->default_order() == FIRST)
682 0 : elem->connectivity(0, TECPLOT, conn);
683 : else
684 : {
685 0 : std::unique_ptr<Elem> lo_elem = Elem::build(Elem::first_order_equivalent_type(elem->type()));
686 0 : for (auto i : lo_elem->node_index_range())
687 0 : lo_elem->set_node(i, const_cast<Node*>(elem->node_ptr(i)));
688 0 : lo_elem->connectivity(0, TECPLOT, conn);
689 0 : }
690 0 : for (const auto & idx : conn)
691 0 : out_stream << idx << " ";
692 :
693 0 : out_stream << '\n';
694 : }
695 :
696 72 : break;
697 : }
698 :
699 363236 : case 2:
700 : {
701 3949137 : if (this->subdivide_second_order())
702 8696775 : for (auto se : make_range(elem->n_sub_elem()))
703 : {
704 : // Quad elements
705 5846836 : if ((elem->type() == QUAD4) ||
706 5186954 : (elem->type() == QUAD8) || // Note: QUAD8 will be output as one central quad and
707 : // four surrounding triangles (though they will be written
708 : // to GMV as QUAD4s).
709 1216342 : (elem->type() == QUAD9)
710 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
711 37914 : || (elem->type() == INFQUAD4)
712 1140218 : || (elem->type() == INFQUAD6)
713 : #endif
714 : )
715 : {
716 4595964 : out_stream << "quad 4\n";
717 4595964 : elem->connectivity(se, TECPLOT, conn);
718 22979820 : for (const auto & idx : conn)
719 18383856 : out_stream << idx << " ";
720 : }
721 :
722 : // Triangle elements
723 151674 : else if ((elem->type() == TRI3) ||
724 0 : (elem->type() == TRI6))
725 : {
726 151674 : out_stream << "tri 3\n";
727 151674 : elem->connectivity(se, TECPLOT, conn);
728 606696 : for (unsigned int i=0; i<3; i++)
729 492942 : out_stream << conn[i] << " ";
730 : }
731 : else
732 0 : libmesh_error_msg("Unsupported element type: " << Utility::enum_to_string(elem->type()));
733 : }
734 : else // !this->subdivide_second_order()
735 : {
736 : // Quad elements
737 0 : if ((elem->type() == QUAD4)
738 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
739 0 : || (elem->type() == INFQUAD4)
740 : #endif
741 : )
742 : {
743 0 : elem->connectivity(0, TECPLOT, conn);
744 0 : out_stream << "quad 4\n";
745 0 : for (const auto & idx : conn)
746 0 : out_stream << idx << " ";
747 : }
748 0 : else if ((elem->type() == QUAD8) ||
749 0 : (elem->type() == QUAD9)
750 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
751 0 : || (elem->type() == INFQUAD6)
752 : #endif
753 : )
754 : {
755 0 : std::unique_ptr<Elem> lo_elem = Elem::build(Elem::first_order_equivalent_type(elem->type()));
756 0 : for (auto i : lo_elem->node_index_range())
757 0 : lo_elem->set_node(i, const_cast<Node*>(elem->node_ptr(i)));
758 0 : lo_elem->connectivity(0, TECPLOT, conn);
759 0 : out_stream << "quad 4\n";
760 0 : for (const auto & idx : conn)
761 0 : out_stream << idx << " ";
762 0 : }
763 0 : else if (elem->type() == TRI3)
764 : {
765 0 : elem->connectivity(0, TECPLOT, conn);
766 0 : out_stream << "tri 3\n";
767 0 : for (unsigned int i=0; i<3; i++)
768 0 : out_stream << conn[i] << " ";
769 : }
770 0 : else if (elem->type() == TRI6)
771 : {
772 0 : std::unique_ptr<Elem> lo_elem = Elem::build(Elem::first_order_equivalent_type(elem->type()));
773 0 : for (auto i : lo_elem->node_index_range())
774 0 : lo_elem->set_node(i, const_cast<Node*>(elem->node_ptr(i)));
775 0 : lo_elem->connectivity(0, TECPLOT, conn);
776 0 : out_stream << "tri 3\n";
777 0 : for (unsigned int i=0; i<3; i++)
778 0 : out_stream << conn[i] << " ";
779 0 : }
780 :
781 0 : out_stream << '\n';
782 : }
783 :
784 363236 : break;
785 : }
786 :
787 22173 : case 3:
788 : {
789 257983 : if (this->subdivide_second_order())
790 515966 : for (auto se : make_range(elem->n_sub_elem()))
791 : {
792 :
793 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
794 216837 : if ((elem->type() == HEX8) ||
795 17280 : (elem->type() == HEX27))
796 : {
797 182277 : out_stream << "phex8 8\n";
798 182277 : elem->connectivity(se, TECPLOT, conn);
799 1640493 : for (const auto & idx : conn)
800 1458216 : out_stream << idx << " ";
801 : }
802 :
803 17280 : else if (elem->type() == HEX20)
804 : {
805 0 : out_stream << "phex20 20\n";
806 0 : out_stream << elem->node_id(0)+1 << " "
807 0 : << elem->node_id(1)+1 << " "
808 0 : << elem->node_id(2)+1 << " "
809 0 : << elem->node_id(3)+1 << " "
810 0 : << elem->node_id(4)+1 << " "
811 0 : << elem->node_id(5)+1 << " "
812 0 : << elem->node_id(6)+1 << " "
813 0 : << elem->node_id(7)+1 << " "
814 0 : << elem->node_id(8)+1 << " "
815 0 : << elem->node_id(9)+1 << " "
816 0 : << elem->node_id(10)+1 << " "
817 0 : << elem->node_id(11)+1 << " "
818 0 : << elem->node_id(16)+1 << " "
819 0 : << elem->node_id(17)+1 << " "
820 0 : << elem->node_id(18)+1 << " "
821 0 : << elem->node_id(19)+1 << " "
822 0 : << elem->node_id(12)+1 << " "
823 0 : << elem->node_id(13)+1 << " "
824 0 : << elem->node_id(14)+1 << " "
825 0 : << elem->node_id(15)+1 << " ";
826 : }
827 :
828 : // According to our copy of gmvread.c, this is the
829 : // mapping for the Hex27 element. Unfortunately,
830 : // I tried it and Paraview does not seem to be
831 : // able to read Hex27 elements. Since this is
832 : // unlikely to change any time soon, we'll
833 : // continue to write out Hex27 elements as 8 Hex8
834 : // sub-elements.
835 : //
836 : // TODO:
837 : // 1.) If we really wanted to use this code for
838 : // something, we'd want to avoid repeating the
839 : // hard-coded node ordering from the HEX20 case.
840 : // These should both be able to use
841 : // ElementDefinitions.
842 : // 2.) You would also need to change
843 : // Hex27::n_sub_elem() to return 1, not sure how
844 : // much other code that would affect...
845 :
846 : // else if (elem->type() == HEX27)
847 : // {
848 : // out_stream << "phex27 27\n";
849 : // out_stream << elem->node_id(0)+1 << " "
850 : // << elem->node_id(1)+1 << " "
851 : // << elem->node_id(2)+1 << " "
852 : // << elem->node_id(3)+1 << " "
853 : // << elem->node_id(4)+1 << " "
854 : // << elem->node_id(5)+1 << " "
855 : // << elem->node_id(6)+1 << " "
856 : // << elem->node_id(7)+1 << " "
857 : // << elem->node_id(8)+1 << " "
858 : // << elem->node_id(9)+1 << " "
859 : // << elem->node_id(10)+1 << " "
860 : // << elem->node_id(11)+1 << " "
861 : // << elem->node_id(16)+1 << " "
862 : // << elem->node_id(17)+1 << " "
863 : // << elem->node_id(18)+1 << " "
864 : // << elem->node_id(19)+1 << " "
865 : // << elem->node_id(12)+1 << " "
866 : // << elem->node_id(13)+1 << " "
867 : // << elem->node_id(14)+1 << " "
868 : // << elem->node_id(15)+1 << " "
869 : // // mid-face nodes
870 : // << elem->node_id(21)+1 << " " // GMV front
871 : // << elem->node_id(22)+1 << " " // GMV right
872 : // << elem->node_id(23)+1 << " " // GMV back
873 : // << elem->node_id(24)+1 << " " // GMV left
874 : // << elem->node_id(20)+1 << " " // GMV bottom
875 : // << elem->node_id(25)+1 << " " // GMV top
876 : // // center node
877 : // << elem->node_id(26)+1 << " ";
878 : // }
879 :
880 : #else // LIBMESH_ENABLE_INFINITE_ELEMENTS
881 :
882 : // In case of infinite elements, HEX20
883 : // should be handled just like the
884 : // INFHEX16, since these connect to each other
885 62266 : if ((elem->type() == HEX8) ||
886 9600 : (elem->type() == HEX27) ||
887 9600 : (elem->type() == INFHEX8) ||
888 9600 : (elem->type() == INFHEX16) ||
889 69946 : (elem->type() == INFHEX18) ||
890 5760 : (elem->type() == HEX20))
891 : {
892 52666 : out_stream << "phex8 8\n";
893 52666 : elem->connectivity(se, TECPLOT, conn);
894 473994 : for (const auto & idx : conn)
895 421328 : out_stream << idx << " ";
896 : }
897 : #endif
898 :
899 46080 : else if ((elem->type() == TET4) ||
900 23040 : (elem->type() == TET10))
901 : {
902 0 : out_stream << "tet 4\n";
903 : // Tecplot connectivity returns 8 entries for
904 : // the Tet, enough to store it as a degenerate Hex.
905 : // For GMV we only pick out the four relevant node
906 : // indices.
907 0 : elem->connectivity(se, TECPLOT, conn);
908 0 : out_stream << conn[0] << " " // libmesh tet node 0
909 0 : << conn[2] << " " // libmesh tet node 2
910 0 : << conn[1] << " " // libmesh tet node 1
911 0 : << conn[4] << " "; // libmesh tet node 3
912 : }
913 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
914 17280 : else if ((elem->type() == PRISM6) ||
915 0 : (elem->type() == PRISM15) ||
916 17280 : (elem->type() == PRISM18) ||
917 0 : (elem->type() == PYRAMID5))
918 : #else
919 5760 : else if ((elem->type() == PRISM6) ||
920 0 : (elem->type() == PRISM15) ||
921 0 : (elem->type() == PRISM18) ||
922 0 : (elem->type() == PYRAMID5) ||
923 5760 : (elem->type() == INFPRISM6) ||
924 0 : (elem->type() == INFPRISM12))
925 : #endif
926 : {
927 : // Note that the prisms are treated as
928 : // degenerated phex8's.
929 23040 : out_stream << "phex8 8\n";
930 23040 : elem->connectivity(se, TECPLOT, conn);
931 207360 : for (const auto & idx : conn)
932 184320 : out_stream << idx << " ";
933 : }
934 :
935 : else
936 0 : libmesh_error_msg("Encountered an unrecognized element " \
937 : << "type: " << elem->type() \
938 : << "\nPossibly a dim-1 dimensional " \
939 : << "element? Aborting...");
940 :
941 257983 : out_stream << '\n';
942 : }
943 : else // !this->subdivide_second_order()
944 : {
945 0 : std::unique_ptr<Elem> lo_elem = Elem::build(Elem::first_order_equivalent_type(elem->type()));
946 0 : for (auto i : lo_elem->node_index_range())
947 0 : lo_elem->set_node(i, const_cast<Node*>(elem->node_ptr(i)));
948 0 : if ((lo_elem->type() == HEX8)
949 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
950 0 : || (lo_elem->type() == HEX27)
951 : #endif
952 : )
953 : {
954 0 : out_stream << "phex8 8\n";
955 0 : lo_elem->connectivity(0, TECPLOT, conn);
956 0 : for (const auto & idx : conn)
957 0 : out_stream << idx << " ";
958 : }
959 :
960 0 : else if (lo_elem->type() == TET4)
961 : {
962 0 : out_stream << "tet 4\n";
963 0 : lo_elem->connectivity(0, TECPLOT, conn);
964 0 : out_stream << conn[0] << " "
965 0 : << conn[2] << " "
966 0 : << conn[1] << " "
967 0 : << conn[4] << " ";
968 : }
969 0 : else if ((lo_elem->type() == PRISM6)
970 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
971 0 : || (lo_elem->type() == INFPRISM6)
972 : #endif
973 : )
974 : {
975 : // Note that the prisms are treated as
976 : // degenerated phex8's.
977 0 : out_stream << "phex8 8\n";
978 0 : lo_elem->connectivity(0, TECPLOT, conn);
979 0 : for (const auto & idx : conn)
980 0 : out_stream << idx << " ";
981 : }
982 :
983 : else
984 0 : libmesh_error_msg("Encountered an unrecognized element " \
985 : << "type. Possibly a dim-1 dimensional " \
986 : << "element? Aborting...");
987 :
988 0 : out_stream << '\n';
989 0 : }
990 :
991 22173 : break;
992 : }
993 :
994 0 : default:
995 0 : libmesh_error_msg("Unsupported element dimension: " <<
996 : elem->dim());
997 : }
998 4170 : }
999 :
1000 5218 : out_stream << '\n';
1001 : }
1002 :
1003 :
1004 :
1005 : // optionally write the partition information
1006 5218 : if (this->partitioning())
1007 : {
1008 5218 : if (this->write_subdomain_id_as_material())
1009 : {
1010 : // Subdomain IDs can be non-contiguous and need not
1011 : // necessarily start at 0. Furthermore, since the user is
1012 : // free to define subdomain_id_type to be a signed type, we
1013 : // can't even assume max(subdomain_id) >= # unique subdomain ids.
1014 :
1015 : // We build a map<subdomain_id, unsigned> to associate to each
1016 : // user-selected subdomain ID a unique, contiguous unsigned value
1017 : // which we can write to file.
1018 0 : std::map<subdomain_id_type, unsigned> sbdid_map;
1019 :
1020 : // Try to insert with dummy value
1021 0 : for (const auto & elem : mesh.active_element_ptr_range())
1022 0 : sbdid_map.emplace(elem->subdomain_id(), 0);
1023 :
1024 : // Map is created, iterate through it to set indices. They will be
1025 : // used repeatedly below.
1026 : {
1027 0 : unsigned ctr=0;
1028 0 : for (auto & pr : sbdid_map)
1029 0 : pr.second = ctr++;
1030 : }
1031 :
1032 0 : out_stream << "material "
1033 0 : << sbdid_map.size()
1034 0 : << " 0\n";
1035 :
1036 0 : for (auto sbdid : make_range(sbdid_map.size()))
1037 0 : out_stream << "proc_" << sbdid << "\n";
1038 :
1039 0 : for (const auto & elem : mesh.active_element_ptr_range())
1040 : {
1041 : // Find the unique index for elem->subdomain_id(), print that to file
1042 0 : unsigned gmv_mat_number = libmesh_map_find(sbdid_map, elem->subdomain_id());
1043 :
1044 0 : if (this->subdivide_second_order())
1045 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1046 0 : out_stream << gmv_mat_number+1 << '\n';
1047 : else
1048 0 : out_stream << gmv_mat_number+1 << "\n";
1049 0 : }
1050 0 : out_stream << '\n';
1051 :
1052 : }
1053 : else // write processor IDs as materials. This is the default
1054 : {
1055 4694 : out_stream << "material "
1056 1048 : << mesh.n_partitions()
1057 5218 : << " 0"<< '\n';
1058 :
1059 27420 : for (auto proc : make_range(mesh.n_partitions()))
1060 42308 : out_stream << "proc_" << proc << '\n';
1061 :
1062 7654918 : for (const auto & elem : mesh.active_element_ptr_range())
1063 4207984 : if (this->subdivide_second_order())
1064 9215333 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1065 5471280 : out_stream << elem->processor_id()+1 << '\n';
1066 : else
1067 4170 : out_stream << elem->processor_id()+1 << '\n';
1068 :
1069 5218 : out_stream << '\n';
1070 : }
1071 : }
1072 :
1073 :
1074 : // If there are *any* variables at all in the system (including
1075 : // p level, or arbitrary cell-based data)
1076 : // to write, the gmv file needs to contain the word "variable"
1077 : // on a line by itself.
1078 524 : bool write_variable = false;
1079 :
1080 : // 1.) p-levels
1081 5218 : if (this->p_levels() && mesh_max_p_level)
1082 0 : write_variable = true;
1083 :
1084 : // 2.) solution data
1085 5218 : if ((solution_names != nullptr) && (v != nullptr))
1086 524 : write_variable = true;
1087 :
1088 : // 3.) cell-centered data
1089 5218 : if (!(this->_cell_centered_data.empty()))
1090 0 : write_variable = true;
1091 :
1092 5218 : if (write_variable)
1093 5218 : out_stream << "variable\n";
1094 :
1095 :
1096 : // optionally write the p-level information
1097 5218 : if (this->p_levels() && mesh_max_p_level)
1098 : {
1099 0 : out_stream << "p_level 0\n";
1100 :
1101 0 : for (const auto & elem : mesh.active_element_ptr_range())
1102 0 : if (this->subdivide_second_order())
1103 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1104 0 : out_stream << elem->p_level() << " ";
1105 : else
1106 0 : out_stream << elem->p_level() << " ";
1107 0 : out_stream << "\n\n";
1108 : }
1109 :
1110 :
1111 :
1112 :
1113 : // optionally write cell-centered data
1114 5218 : if (!(this->_cell_centered_data.empty()))
1115 : {
1116 0 : for (const auto & [var_name, the_array] : this->_cell_centered_data)
1117 : {
1118 : // write out the variable name, followed by a zero.
1119 0 : out_stream << var_name << " 0\n";
1120 :
1121 : // Loop over active elements, write out cell data. If second-order cells
1122 : // are split into sub-elements, the sub-elements inherit their parent's
1123 : // cell-centered data.
1124 0 : for (const auto & elem : mesh.active_element_ptr_range())
1125 : {
1126 : // Use the element's ID to find the value...
1127 0 : libmesh_assert_less (elem->id(), the_array->size());
1128 0 : const Real the_value = (*the_array)[elem->id()];
1129 :
1130 0 : if (this->subdivide_second_order())
1131 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1132 0 : out_stream << the_value << " ";
1133 : else
1134 0 : out_stream << the_value << " ";
1135 0 : }
1136 :
1137 0 : out_stream << "\n\n";
1138 : }
1139 : }
1140 :
1141 :
1142 :
1143 :
1144 : // optionally write the data
1145 5218 : if ((solution_names != nullptr) &&
1146 : (v != nullptr))
1147 : {
1148 : const unsigned int n_vars =
1149 1048 : cast_int<unsigned int>(solution_names->size());
1150 :
1151 5742 : if (!(v->size() == mesh.n_nodes()*n_vars))
1152 0 : libMesh::err << "ERROR: v->size()=" << v->size()
1153 0 : << ", mesh.n_nodes()=" << mesh.n_nodes()
1154 0 : << ", n_vars=" << n_vars
1155 0 : << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
1156 0 : << std::endl;
1157 :
1158 524 : libmesh_assert_equal_to (v->size(), mesh.n_nodes()*n_vars);
1159 :
1160 10871 : for (unsigned int c=0; c<n_vars; c++)
1161 : {
1162 :
1163 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1164 :
1165 : // in case of complex data, write _tree_ data sets
1166 : // for each component
1167 :
1168 : // this is the real part
1169 132 : out_stream << "r_" << (*solution_names)[c] << " 1\n";
1170 :
1171 138716 : for (auto n : make_range(mesh.n_nodes()))
1172 138584 : out_stream << (*v)[n*n_vars + c].real() << " ";
1173 :
1174 132 : out_stream << '\n' << '\n';
1175 :
1176 :
1177 : // this is the imaginary part
1178 264 : out_stream << "i_" << (*solution_names)[c] << " 1\n";
1179 :
1180 138716 : for (auto n : make_range(mesh.n_nodes()))
1181 138584 : out_stream << (*v)[n*n_vars + c].imag() << " ";
1182 :
1183 132 : out_stream << '\n' << '\n';
1184 :
1185 : // this is the magnitude
1186 264 : out_stream << "a_" << (*solution_names)[c] << " 1\n";
1187 138716 : for (auto n : make_range(mesh.n_nodes()))
1188 138584 : out_stream << std::abs((*v)[n*n_vars + c]) << " ";
1189 :
1190 132 : out_stream << '\n' << '\n';
1191 :
1192 : #else
1193 :
1194 5521 : out_stream << (*solution_names)[c] << " 1\n";
1195 :
1196 5218743 : for (auto n : make_range(mesh.n_nodes()))
1197 5710690 : out_stream << (*v)[n*n_vars + c] << " ";
1198 :
1199 5521 : out_stream << '\n' << '\n';
1200 :
1201 : #endif
1202 : }
1203 :
1204 : }
1205 :
1206 : // If we wrote any variables, we have to close the variable section now
1207 5218 : if (write_variable)
1208 5218 : out_stream << "endvars\n";
1209 :
1210 :
1211 : // end of the file
1212 5218 : out_stream << "\nendgmv\n";
1213 29622 : }
1214 :
1215 :
1216 :
1217 :
1218 :
1219 :
1220 :
1221 0 : void GMVIO::write_binary (const std::string & fname,
1222 : const std::vector<Number> * vec,
1223 : const std::vector<std::string> * solution_names)
1224 : {
1225 : // We use the file-scope global variable eletypes for mapping nodes
1226 : // from GMV to libmesh indices, so initialize that data now.
1227 0 : init_eletypes();
1228 :
1229 : // Get a reference to the mesh
1230 0 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1231 :
1232 : // This is a parallel_only function
1233 0 : const dof_id_type n_active_elem = mesh.n_active_elem();
1234 :
1235 0 : if (MeshOutput<MeshBase>::mesh().processor_id() != 0)
1236 0 : return;
1237 :
1238 0 : std::ofstream out_stream (fname.c_str());
1239 :
1240 0 : libmesh_assert (out_stream.good());
1241 :
1242 0 : unsigned int mesh_max_p_level = 0;
1243 :
1244 0 : std::string buffer;
1245 :
1246 : // Begin interfacing with the GMV data file
1247 : {
1248 : // write the nodes
1249 0 : buffer = "gmvinput";
1250 0 : out_stream.write(buffer.c_str(), buffer.size());
1251 :
1252 0 : buffer = "ieeei4r4";
1253 0 : out_stream.write(buffer.c_str(), buffer.size());
1254 : }
1255 :
1256 :
1257 :
1258 : // write the nodes
1259 : {
1260 0 : buffer = "nodes ";
1261 0 : out_stream.write(buffer.c_str(), buffer.size());
1262 :
1263 0 : unsigned int tempint = mesh.n_nodes();
1264 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1265 :
1266 : // write the x coordinate
1267 0 : std::vector<float> temp(mesh.n_nodes());
1268 0 : for (auto v : make_range(mesh.n_nodes()))
1269 0 : temp[v] = static_cast<float>(mesh.point(v)(0));
1270 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1271 :
1272 : // write the y coordinate
1273 0 : for (auto v : make_range(mesh.n_nodes()))
1274 : {
1275 : #if LIBMESH_DIM > 1
1276 0 : temp[v] = static_cast<float>(mesh.point(v)(1));
1277 : #else
1278 : temp[v] = 0.;
1279 : #endif
1280 : }
1281 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1282 :
1283 : // write the z coordinate
1284 0 : for (auto v : make_range(mesh.n_nodes()))
1285 : {
1286 : #if LIBMESH_DIM > 2
1287 0 : temp[v] = static_cast<float>(mesh.point(v)(2));
1288 : #else
1289 : temp[v] = 0.;
1290 : #endif
1291 : }
1292 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1293 : }
1294 :
1295 :
1296 : // write the connectivity
1297 : {
1298 0 : buffer = "cells ";
1299 0 : out_stream.write(buffer.c_str(), buffer.size());
1300 :
1301 0 : unsigned int tempint = n_active_elem;
1302 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1303 :
1304 0 : for (const auto & elem : mesh.active_element_ptr_range())
1305 : {
1306 0 : mesh_max_p_level = std::max(mesh_max_p_level,
1307 0 : elem->p_level());
1308 :
1309 : // The ElementDefinition label contains the GMV name
1310 : // and the number of nodes for the respective
1311 : // element.
1312 0 : const ElementDefinition & ed = eletypes[elem->type()];
1313 :
1314 : // The ElementDefinition labels all have a name followed by a
1315 : // whitespace and then the number of nodes. To write the
1316 : // string in binary, we want to drop everything after that
1317 : // space, and then pad the string out to 8 characters.
1318 0 : buffer = ed.label;
1319 :
1320 : // Erase everything from the first whitespace character to the end of the string.
1321 0 : buffer.erase(buffer.find_first_of(' '), std::string::npos);
1322 :
1323 : // Append whitespaces until the string is exactly 8 characters long.
1324 0 : while (buffer.size() < 8)
1325 0 : buffer.insert(buffer.end(), ' ');
1326 :
1327 : // Finally, write the 8 character stream to file.
1328 0 : out_stream.write(buffer.c_str(), buffer.size());
1329 :
1330 : // Write the number of nodes that this type of element has.
1331 : // Note: don't want to use the number of nodes of the element,
1332 : // because certain elements, like QUAD9 and HEX27 only support
1333 : // being written out as lower-order elements (QUAD8 and HEX20,
1334 : // respectively).
1335 0 : tempint = cast_int<unsigned int>(ed.node_map.size());
1336 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1337 :
1338 : // Write the element connectivity
1339 0 : for (const auto & ed_id : ed.node_map)
1340 : {
1341 0 : dof_id_type id = elem->node_id(ed_id) + 1;
1342 0 : out_stream.write(reinterpret_cast<char *>(&id), sizeof(dof_id_type));
1343 : }
1344 0 : }
1345 : }
1346 :
1347 :
1348 :
1349 : // optionally write the partition information
1350 0 : if (this->partitioning())
1351 : {
1352 0 : libmesh_error_msg_if(this->write_subdomain_id_as_material(),
1353 : "Not yet supported in GMVIO::write_binary");
1354 :
1355 0 : buffer = "material";
1356 0 : out_stream.write(buffer.c_str(), buffer.size());
1357 :
1358 0 : unsigned int tmpint = mesh.n_processors();
1359 0 : out_stream.write(reinterpret_cast<char *>(&tmpint), sizeof(unsigned int));
1360 :
1361 0 : tmpint = 0; // IDs are cell based
1362 0 : out_stream.write(reinterpret_cast<char *>(&tmpint), sizeof(unsigned int));
1363 :
1364 0 : for (auto proc : make_range(mesh.n_processors()))
1365 : {
1366 : // We have to write exactly 8 bytes. This means that
1367 : // the processor id can be up to 3 characters long, and
1368 : // we pad it with blank characters on the end.
1369 0 : std::ostringstream oss;
1370 0 : oss << "proc_" << std::setw(3) << std::left << proc;
1371 0 : out_stream.write(oss.str().c_str(), oss.str().size());
1372 0 : }
1373 :
1374 0 : std::vector<unsigned int> proc_id (n_active_elem);
1375 :
1376 0 : unsigned int n=0;
1377 :
1378 : // We no longer write sub-elems in GMV, so just assign a
1379 : // processor id value to each element.
1380 0 : for (const auto & elem : mesh.active_element_ptr_range())
1381 0 : proc_id[n++] = elem->processor_id() + 1;
1382 :
1383 0 : out_stream.write(reinterpret_cast<char *>(proc_id.data()),
1384 0 : sizeof(unsigned int)*proc_id.size());
1385 : }
1386 :
1387 : // If there are *any* variables at all in the system (including
1388 : // p level, or arbitrary cell-based data)
1389 : // to write, the gmv file needs to contain the word "variable"
1390 : // on a line by itself.
1391 0 : bool write_variable = false;
1392 :
1393 : // 1.) p-levels
1394 0 : if (this->p_levels() && mesh_max_p_level)
1395 0 : write_variable = true;
1396 :
1397 : // 2.) solution data
1398 0 : if ((solution_names != nullptr) && (vec != nullptr))
1399 0 : write_variable = true;
1400 :
1401 : // // 3.) cell-centered data - unsupported
1402 : // if (!(this->_cell_centered_data.empty()))
1403 : // write_variable = true;
1404 :
1405 0 : if (write_variable)
1406 : {
1407 0 : buffer = "variable";
1408 0 : out_stream.write(buffer.c_str(), buffer.size());
1409 : }
1410 :
1411 : // optionally write the partition information
1412 0 : if (this->p_levels() && mesh_max_p_level)
1413 : {
1414 0 : unsigned int n_floats = n_active_elem;
1415 0 : for (unsigned int i=0, d=mesh.mesh_dimension(); i != d; ++i)
1416 0 : n_floats *= 2;
1417 :
1418 0 : std::vector<float> temp(n_floats);
1419 :
1420 0 : buffer = "p_level";
1421 0 : out_stream.write(buffer.c_str(), buffer.size());
1422 :
1423 0 : unsigned int tempint = 0; // p levels are cell data
1424 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1425 :
1426 0 : unsigned int n=0;
1427 0 : for (const auto & elem : mesh.active_element_ptr_range())
1428 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1429 0 : temp[n++] = static_cast<float>( elem->p_level() );
1430 :
1431 0 : out_stream.write(reinterpret_cast<char *>(temp.data()),
1432 0 : sizeof(float)*n_floats);
1433 : }
1434 :
1435 :
1436 : // optionally write cell-centered data
1437 0 : if (!(this->_cell_centered_data.empty()))
1438 0 : libMesh::err << "Cell-centered data not (yet) supported in binary I/O mode!" << std::endl;
1439 :
1440 :
1441 :
1442 :
1443 : // optionally write the data
1444 0 : if ((solution_names != nullptr) &&
1445 : (vec != nullptr))
1446 : {
1447 0 : std::vector<float> temp(mesh.n_nodes());
1448 :
1449 : const unsigned int n_vars =
1450 0 : cast_int<unsigned int>(solution_names->size());
1451 :
1452 0 : for (unsigned int c=0; c<n_vars; c++)
1453 : {
1454 :
1455 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1456 : // for complex data, write three datasets
1457 :
1458 :
1459 : // Real part
1460 : buffer = "r_";
1461 0 : out_stream.write(buffer.c_str(), buffer.size());
1462 :
1463 0 : buffer = (*solution_names)[c];
1464 0 : out_stream.write(buffer.c_str(), buffer.size());
1465 :
1466 0 : unsigned int tempint = 1; // always do nodal data
1467 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1468 :
1469 0 : for (auto n : make_range(mesh.n_nodes()))
1470 0 : temp[n] = static_cast<float>( (*vec)[n*n_vars + c].real() );
1471 :
1472 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1473 :
1474 :
1475 : // imaginary part
1476 : buffer = "i_";
1477 0 : out_stream.write(buffer.c_str(), buffer.size());
1478 :
1479 : buffer = (*solution_names)[c];
1480 0 : out_stream.write(buffer.c_str(), buffer.size());
1481 :
1482 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1483 :
1484 0 : for (auto n : make_range(mesh.n_nodes()))
1485 0 : temp[n] = static_cast<float>( (*vec)[n*n_vars + c].imag() );
1486 :
1487 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1488 :
1489 : // magnitude
1490 : buffer = "a_";
1491 0 : out_stream.write(buffer.c_str(), buffer.size());
1492 : buffer = (*solution_names)[c];
1493 0 : out_stream.write(buffer.c_str(), buffer.size());
1494 :
1495 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1496 :
1497 0 : for (auto n : make_range(mesh.n_nodes()))
1498 0 : temp[n] = static_cast<float>(std::abs((*vec)[n*n_vars + c]));
1499 :
1500 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1501 :
1502 : #else
1503 :
1504 0 : buffer = (*solution_names)[c];
1505 0 : out_stream.write(buffer.c_str(), buffer.size());
1506 :
1507 0 : unsigned int tempint = 1; // always do nodal data
1508 0 : out_stream.write(reinterpret_cast<char *>(&tempint), sizeof(unsigned int));
1509 :
1510 0 : for (auto n : make_range(mesh.n_nodes()))
1511 0 : temp[n] = static_cast<float>((*vec)[n*n_vars + c]);
1512 :
1513 0 : out_stream.write(reinterpret_cast<char *>(temp.data()), sizeof(float)*mesh.n_nodes());
1514 :
1515 : #endif
1516 : }
1517 : }
1518 :
1519 : // If we wrote any variables, we have to close the variable section now
1520 0 : if (write_variable)
1521 : {
1522 0 : buffer = "endvars ";
1523 0 : out_stream.write(buffer.c_str(), buffer.size());
1524 : }
1525 :
1526 : // end the file
1527 0 : buffer = "endgmv ";
1528 0 : out_stream.write(buffer.c_str(), buffer.size());
1529 0 : }
1530 :
1531 :
1532 :
1533 :
1534 :
1535 :
1536 :
1537 :
1538 :
1539 3946 : void GMVIO::write_discontinuous_gmv (const std::string & name,
1540 : const EquationSystems & es,
1541 : const bool write_partitioning,
1542 : const std::set<std::string> * system_names) const
1543 : {
1544 279 : std::vector<std::string> solution_names;
1545 186 : std::vector<Number> v;
1546 :
1547 : // Get a reference to the mesh
1548 372 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
1549 :
1550 3946 : es.build_variable_names (solution_names, nullptr, system_names);
1551 3946 : es.build_discontinuous_solution_vector (v, system_names);
1552 :
1553 : // These are parallel_only functions
1554 3946 : const unsigned int n_active_elem = mesh.n_active_elem();
1555 :
1556 4132 : if (mesh.processor_id() != 0)
1557 93 : return;
1558 :
1559 909 : std::ofstream out_stream(name.c_str());
1560 :
1561 93 : libmesh_assert (out_stream.good());
1562 :
1563 : // Begin interfacing with the GMV data file
1564 : {
1565 :
1566 : // write the nodes
1567 630 : out_stream << "gmvinput ascii" << std::endl << std::endl;
1568 :
1569 : // Compute the total weight
1570 : {
1571 93 : unsigned int tw=0;
1572 :
1573 44568 : for (const auto & elem : mesh.active_element_ptr_range())
1574 43752 : tw += elem->n_nodes();
1575 :
1576 630 : out_stream << "nodes " << tw << std::endl;
1577 : }
1578 :
1579 :
1580 :
1581 : // Write all the x values
1582 : {
1583 71935 : for (const auto & elem : mesh.active_element_ptr_range())
1584 223999 : for (const Node & node : elem->node_ref_range())
1585 173397 : out_stream << node(0) << " ";
1586 :
1587 93 : out_stream << std::endl;
1588 : }
1589 :
1590 :
1591 : // Write all the y values
1592 : {
1593 71935 : for (const auto & elem : mesh.active_element_ptr_range())
1594 223999 : for (const Node & node : elem->node_ref_range())
1595 : #if LIBMESH_DIM > 1
1596 173397 : out_stream << node(1) << " ";
1597 : #else
1598 : out_stream << 0. << " ";
1599 : #endif
1600 :
1601 93 : out_stream << std::endl;
1602 : }
1603 :
1604 :
1605 : // Write all the z values
1606 : {
1607 71935 : for (const auto & elem : mesh.active_element_ptr_range())
1608 223999 : for (const Node & node : elem->node_ref_range())
1609 : #if LIBMESH_DIM > 2
1610 173397 : out_stream << node(2) << " ";
1611 : #else
1612 : out_stream << 0. << " ";
1613 : #endif
1614 :
1615 93 : out_stream << std::endl << std::endl;
1616 : }
1617 : }
1618 :
1619 :
1620 :
1621 : {
1622 : // write the connectivity
1623 :
1624 630 : out_stream << "cells " << n_active_elem << std::endl;
1625 :
1626 93 : unsigned int nn=1;
1627 :
1628 723 : switch (mesh.mesh_dimension())
1629 : {
1630 0 : case 1:
1631 : {
1632 0 : for (const auto & elem : mesh.active_element_ptr_range())
1633 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1634 : {
1635 0 : if ((elem->type() == EDGE2) ||
1636 0 : (elem->type() == EDGE3) ||
1637 0 : (elem->type() == EDGE4)
1638 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1639 0 : || (elem->type() == INFEDGE2)
1640 : #endif
1641 : )
1642 : {
1643 0 : out_stream << "line 2" << std::endl;
1644 0 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1645 0 : out_stream << nn++ << " ";
1646 :
1647 : }
1648 : else
1649 0 : libmesh_error_msg("Unsupported 1D element type: " << Utility::enum_to_string(elem->type()));
1650 :
1651 0 : out_stream << std::endl;
1652 0 : }
1653 :
1654 0 : break;
1655 : }
1656 :
1657 723 : case 2:
1658 : {
1659 71935 : for (const auto & elem : mesh.active_element_ptr_range())
1660 86430 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1661 : {
1662 43215 : if ((elem->type() == QUAD4) ||
1663 21781 : (elem->type() == QUAD8) || // Note: QUAD8 will be output as one central quad and
1664 : // four surrounding triangles (though they will be written
1665 : // to GMV as QUAD4s).
1666 0 : (elem->type() == QUAD9)
1667 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1668 0 : || (elem->type() == INFQUAD4)
1669 21434 : || (elem->type() == INFQUAD6)
1670 : #endif
1671 : )
1672 : {
1673 35291 : out_stream << "quad 4" << std::endl;
1674 216075 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1675 172860 : out_stream << nn++ << " ";
1676 :
1677 : }
1678 0 : else if ((elem->type() == TRI3) ||
1679 0 : (elem->type() == TRI6))
1680 : {
1681 0 : out_stream << "tri 3" << std::endl;
1682 0 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1683 0 : out_stream << nn++ << " ";
1684 :
1685 : }
1686 : else
1687 0 : libmesh_error_msg("Unsupported 2D element type: " << Utility::enum_to_string(elem->type()));
1688 :
1689 7924 : out_stream << std::endl;
1690 537 : }
1691 :
1692 723 : break;
1693 : }
1694 :
1695 :
1696 0 : case 3:
1697 : {
1698 0 : for (const auto & elem : mesh.active_element_ptr_range())
1699 0 : for (unsigned int se=0, nse=elem->n_sub_elem(); se<nse; se++)
1700 : {
1701 0 : if ((elem->type() == HEX8) ||
1702 0 : (elem->type() == HEX20) ||
1703 0 : (elem->type() == HEX27)
1704 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1705 0 : || (elem->type() == INFHEX8)
1706 0 : || (elem->type() == INFHEX16)
1707 0 : || (elem->type() == INFHEX18)
1708 : #endif
1709 : )
1710 : {
1711 0 : out_stream << "phex8 8" << std::endl;
1712 0 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1713 0 : out_stream << nn++ << " ";
1714 : }
1715 0 : else if ((elem->type() == PRISM6) ||
1716 0 : (elem->type() == PRISM15) ||
1717 0 : (elem->type() == PRISM18)
1718 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1719 0 : || (elem->type() == INFPRISM6)
1720 0 : || (elem->type() == INFPRISM12)
1721 : #endif
1722 : )
1723 : {
1724 0 : out_stream << "pprism6 6" << std::endl;
1725 0 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1726 0 : out_stream << nn++ << " ";
1727 : }
1728 0 : else if ((elem->type() == TET4) ||
1729 0 : (elem->type() == TET10))
1730 : {
1731 0 : out_stream << "tet 4" << std::endl;
1732 0 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1733 0 : out_stream << nn++ << " ";
1734 : }
1735 : else
1736 0 : libmesh_error_msg("Unsupported 3D element type: " << Utility::enum_to_string(elem->type()));
1737 :
1738 0 : out_stream << std::endl;
1739 0 : }
1740 :
1741 0 : break;
1742 : }
1743 :
1744 0 : default:
1745 0 : libmesh_error_msg("Unsupported mesh dimension: " << mesh.mesh_dimension());
1746 : }
1747 :
1748 93 : out_stream << std::endl;
1749 : }
1750 :
1751 :
1752 :
1753 : // optionally write the partition information
1754 723 : if (write_partitioning)
1755 : {
1756 0 : libmesh_error_msg_if(_write_subdomain_id_as_material,
1757 : "Not yet supported in GMVIO::write_discontinuous_gmv");
1758 :
1759 0 : out_stream << "material "
1760 0 : << mesh.n_processors()
1761 0 : << " 0"<< std::endl;
1762 :
1763 0 : for (auto proc : make_range(mesh.n_processors()))
1764 0 : out_stream << "proc_" << proc << std::endl;
1765 :
1766 0 : for (const auto & elem : mesh.active_element_ptr_range())
1767 0 : out_stream << elem->processor_id()+1 << std::endl;
1768 :
1769 0 : out_stream << std::endl;
1770 : }
1771 :
1772 :
1773 : // Writing cell-centered data is not yet supported in discontinuous GMV files.
1774 723 : if (!(this->_cell_centered_data.empty()))
1775 : {
1776 0 : libMesh::err << "Cell-centered data not (yet) supported for discontinuous GMV files!" << std::endl;
1777 : }
1778 :
1779 :
1780 :
1781 : // write the data
1782 : {
1783 : const unsigned int n_vars =
1784 186 : cast_int<unsigned int>(solution_names.size());
1785 :
1786 : // libmesh_assert_equal_to (v.size(), tw*n_vars);
1787 :
1788 630 : out_stream << "variable" << std::endl;
1789 :
1790 :
1791 1446 : for (unsigned int c=0; c<n_vars; c++)
1792 : {
1793 :
1794 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1795 :
1796 : // in case of complex data, write _tree_ data sets
1797 : // for each component
1798 :
1799 : // this is the real part
1800 44 : out_stream << "r_" << solution_names[c] << " 1" << std::endl;
1801 11260 : for (const auto & elem : mesh.active_element_ptr_range())
1802 27930 : for (auto n : elem->node_index_range())
1803 22388 : out_stream << v[(n++)*n_vars + c].real() << " ";
1804 : out_stream << std::endl << std::endl;
1805 :
1806 :
1807 : // this is the imaginary part
1808 88 : out_stream << "i_" << solution_names[c] << " 1" << std::endl;
1809 11260 : for (const auto & elem : mesh.active_element_ptr_range())
1810 27930 : for (auto n : elem->node_index_range())
1811 22388 : out_stream << v[(n++)*n_vars + c].imag() << " ";
1812 : out_stream << std::endl << std::endl;
1813 :
1814 : // this is the magnitude
1815 88 : out_stream << "a_" << solution_names[c] << " 1" << std::endl;
1816 11260 : for (const auto & elem : mesh.active_element_ptr_range())
1817 27930 : for (auto n : elem->node_index_range())
1818 22388 : out_stream << std::abs(v[(n++)*n_vars + c]) << " ";
1819 : out_stream << std::endl << std::endl;
1820 :
1821 : #else
1822 :
1823 679 : out_stream << solution_names[c] << " 1" << std::endl;
1824 : {
1825 93 : unsigned int nn=0;
1826 :
1827 60675 : for (const auto & elem : mesh.active_element_ptr_range())
1828 188145 : for (unsigned int i=0, enn=elem->n_nodes(); i<enn; i++)
1829 182705 : out_stream << v[(nn++)*n_vars + c] << " ";
1830 : }
1831 93 : out_stream << std::endl << std::endl;
1832 :
1833 : #endif
1834 :
1835 : }
1836 :
1837 630 : out_stream << "endvars" << std::endl;
1838 : }
1839 :
1840 :
1841 : // end of the file
1842 630 : out_stream << std::endl << "endgmv" << std::endl;
1843 4111 : }
1844 :
1845 :
1846 :
1847 :
1848 :
1849 0 : void GMVIO::add_cell_centered_data (const std::string & cell_centered_data_name,
1850 : const std::vector<Real> * cell_centered_data_vals)
1851 : {
1852 0 : libmesh_assert(cell_centered_data_vals);
1853 :
1854 : // Make sure there are *at least* enough entries for all the active elements.
1855 : // There can also be entries for inactive elements, they will be ignored.
1856 : // libmesh_assert_greater_equal (cell_centered_data_vals->size(),
1857 : // MeshOutput<MeshBase>::mesh().n_active_elem());
1858 0 : this->_cell_centered_data[cell_centered_data_name] = cell_centered_data_vals;
1859 0 : }
1860 :
1861 :
1862 :
1863 :
1864 :
1865 :
1866 0 : void GMVIO::read (const std::string & name)
1867 : {
1868 : // This is a serial-only process for now;
1869 : // the Mesh should be read on processor 0 and
1870 : // broadcast later
1871 0 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
1872 :
1873 0 : _next_elem_id = 0;
1874 :
1875 : libmesh_experimental();
1876 :
1877 : #ifndef LIBMESH_HAVE_GMV
1878 :
1879 : libmesh_error_msg("Cannot read GMV file " << name << " without the GMV API.");
1880 :
1881 : #else
1882 : // We use the file-scope global variable eletypes for mapping nodes
1883 : // from GMV to libmesh indices, so initialize that data now.
1884 0 : init_eletypes();
1885 :
1886 : // Clear the mesh so we are sure to start from a pristine state.
1887 0 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
1888 0 : mesh.clear();
1889 :
1890 : // Keep track of what kinds of elements this file contains
1891 0 : elems_of_dimension.clear();
1892 0 : elems_of_dimension.resize(4, false);
1893 :
1894 : // It is apparently possible for gmv files to contain
1895 : // a "fromfile" directive (?) But we currently don't make
1896 : // any use of this feature in LibMesh. Nonzero return val
1897 : // from any function usually means an error has occurred.
1898 0 : int ierr = GMVLib::gmvread_open_fromfileskip(const_cast<char *>(name.c_str()));
1899 0 : libmesh_error_msg_if(ierr != 0, "GMVLib::gmvread_open_fromfileskip failed!");
1900 :
1901 : // Loop through file until GMVEND.
1902 0 : int iend = 0;
1903 0 : while (iend == 0)
1904 : {
1905 0 : GMVLib::gmvread_data();
1906 :
1907 : // Check for GMVEND.
1908 0 : if (GMVLib::gmv_data.keyword == GMVEND)
1909 : {
1910 0 : iend = 1;
1911 0 : GMVLib::gmvread_close();
1912 0 : break;
1913 : }
1914 :
1915 : // Check for GMVERROR.
1916 0 : libmesh_error_msg_if(GMVLib::gmv_data.keyword == GMVERROR,
1917 : "Encountered GMVERROR while reading!");
1918 :
1919 : // Process the data.
1920 0 : switch (GMVLib::gmv_data.keyword)
1921 : {
1922 0 : case NODES:
1923 : {
1924 0 : if (GMVLib::gmv_data.num2 == NODES)
1925 0 : this->_read_nodes();
1926 :
1927 : else
1928 0 : libmesh_error_msg_if(GMVLib::gmv_data.num2 == NODE_V,
1929 : "Unsupported GMV data type NODE_V!");
1930 :
1931 0 : break;
1932 : }
1933 :
1934 0 : case CELLS:
1935 : {
1936 : // Read 1 cell at a time
1937 0 : this->_read_one_cell();
1938 0 : break;
1939 : }
1940 :
1941 0 : case MATERIAL:
1942 : {
1943 : // keyword == 6
1944 : // These are the materials, which we use to specify the mesh
1945 : // partitioning.
1946 0 : this->_read_materials();
1947 0 : break;
1948 : }
1949 :
1950 0 : case VARIABLE:
1951 : {
1952 : // keyword == 8
1953 : // This is a field variable.
1954 :
1955 : // Check to see if we're done reading variables and break out.
1956 0 : if (GMVLib::gmv_data.datatype == ENDKEYWORD)
1957 : {
1958 0 : break;
1959 : }
1960 :
1961 0 : if (GMVLib::gmv_data.datatype == NODE)
1962 : {
1963 0 : this->_read_var();
1964 0 : break;
1965 : }
1966 :
1967 : else
1968 : {
1969 0 : libMesh::err << "Warning: Skipping variable: "
1970 0 : << GMVLib::gmv_data.name1
1971 0 : << " which is of unsupported GMV datatype "
1972 0 : << GMVLib::gmv_data.datatype
1973 0 : << ". Nodal field data is currently the only type currently supported."
1974 0 : << std::endl;
1975 0 : break;
1976 : }
1977 :
1978 : }
1979 :
1980 0 : default:
1981 0 : libmesh_error_msg("Encountered unknown GMV keyword " << GMVLib::gmv_data.keyword);
1982 :
1983 : } // end switch
1984 : } // end while
1985 :
1986 : // Set the mesh dimension to the largest encountered for an element
1987 0 : for (unsigned char i=0; i!=4; ++i)
1988 0 : if (elems_of_dimension[i])
1989 0 : mesh.set_mesh_dimension(i);
1990 :
1991 : #if LIBMESH_DIM < 3
1992 : libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
1993 : "Cannot open dimension "
1994 : << mesh.mesh_dimension()
1995 : << " mesh file when configured without "
1996 : << mesh.mesh_dimension()
1997 : << "D support.");
1998 : #endif
1999 :
2000 : // Done reading in the mesh, now call find_neighbors, etc.
2001 : // mesh.find_neighbors();
2002 :
2003 : // Don't allow renumbering of nodes and elements when calling
2004 : // prepare_for_use(). It is usually confusing for users when the
2005 : // Mesh gets renumbered automatically, since sometimes there are
2006 : // other parts of the simulation which rely on the original
2007 : // ordering.
2008 0 : mesh.allow_renumbering(false);
2009 0 : mesh.prepare_for_use();
2010 : #endif
2011 0 : }
2012 :
2013 :
2014 :
2015 :
2016 0 : void GMVIO::_read_var()
2017 : {
2018 : #ifdef LIBMESH_HAVE_GMV
2019 :
2020 : // Copy all the variable's values into a local storage vector.
2021 0 : _nodal_data.emplace(std::string(GMVLib::gmv_data.name1),
2022 0 : std::vector<Number>(GMVLib::gmv_data.doubledata1, GMVLib::gmv_data.doubledata1+GMVLib::gmv_data.num));
2023 : #endif
2024 0 : }
2025 :
2026 :
2027 :
2028 0 : void GMVIO::_read_materials()
2029 : {
2030 : #ifdef LIBMESH_HAVE_GMV
2031 :
2032 : // LibMesh assigns materials on a per-cell basis
2033 0 : libmesh_assert_equal_to (GMVLib::gmv_data.datatype, CELL);
2034 :
2035 : // Material names: LibMesh has no use for these currently...
2036 : // for (int i = 0; i < GMVLib::gmv_data.num; i++)
2037 : // {
2038 : // // Build a 32-char string from the appropriate entries
2039 : // std::string mat_string(&GMVLib::gmv_data.chardata1[i*33], 32);
2040 : // }
2041 :
2042 0 : for (int i = 0; i < GMVLib::gmv_data.nlongdata1; i++)
2043 0 : MeshInput<MeshBase>::mesh().elem_ref(i).processor_id() =
2044 0 : cast_int<processor_id_type>(GMVLib::gmv_data.longdata1[i]-1);
2045 :
2046 : #endif
2047 0 : }
2048 :
2049 :
2050 :
2051 :
2052 0 : void GMVIO::_read_nodes()
2053 : {
2054 : #ifdef LIBMESH_HAVE_GMV
2055 : // LibMesh writes UNSTRUCT=100 node data
2056 0 : libmesh_assert_equal_to (GMVLib::gmv_data.datatype, UNSTRUCT);
2057 :
2058 : // The nodal data is stored in gmv_data.doubledata{1,2,3}
2059 : // and is nnodes long
2060 0 : for (int i = 0; i < GMVLib::gmv_data.num; i++)
2061 : {
2062 : // Add the point to the Mesh
2063 0 : MeshInput<MeshBase>::mesh().add_point(Point(GMVLib::gmv_data.doubledata1[i],
2064 0 : GMVLib::gmv_data.doubledata2[i],
2065 0 : GMVLib::gmv_data.doubledata3[i]), i);
2066 : }
2067 : #endif
2068 0 : }
2069 :
2070 :
2071 0 : void GMVIO::_read_one_cell()
2072 : {
2073 : #ifdef LIBMESH_HAVE_GMV
2074 : // This is either a REGULAR=111 cell or
2075 : // the ENDKEYWORD=207 of the cells
2076 : #ifndef NDEBUG
2077 0 : bool recognized =
2078 0 : (GMVLib::gmv_data.datatype==REGULAR) ||
2079 0 : (GMVLib::gmv_data.datatype==ENDKEYWORD);
2080 : #endif
2081 0 : libmesh_assert (recognized);
2082 :
2083 0 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
2084 :
2085 0 : if (GMVLib::gmv_data.datatype == REGULAR)
2086 : {
2087 : // We need a mapping from GMV element types to LibMesh
2088 : // ElemTypes. Basically the reverse of the eletypes
2089 : // std::map above.
2090 : //
2091 : // FIXME: Since Quad9's apparently don't exist for GMV, and since
2092 : // In general we write linear sub-elements to GMV files, we need
2093 : // to be careful to read back in exactly what we wrote out...
2094 : //
2095 : // The gmv_data.name1 field is padded with blank characters when
2096 : // it is read in by GMV, so the function that converts it to the
2097 : // libmesh element type needs to take that into account.
2098 0 : ElemType type = this->gmv_elem_to_libmesh_elem(GMVLib::gmv_data.name1);
2099 :
2100 0 : auto elem = Elem::build(type);
2101 0 : elem->set_id(_next_elem_id++);
2102 :
2103 : // Get the ElementDefinition object for this element type
2104 0 : const ElementDefinition & eledef = eletypes[type];
2105 :
2106 : // Print out the connectivity information for
2107 : // this cell.
2108 0 : for (int i=0; i<GMVLib::gmv_data.num2; i++)
2109 : {
2110 : // Map index i to GMV's numbering scheme
2111 0 : unsigned mapped_i = eledef.node_map[i];
2112 :
2113 : // Note: Node numbers (as stored in libmesh) are 1-based
2114 0 : elem->set_node
2115 0 : (i, mesh.node_ptr
2116 0 : (cast_int<dof_id_type>(GMVLib::gmv_data.longdata1[mapped_i]-1)));
2117 : }
2118 :
2119 0 : elems_of_dimension[elem->dim()] = true;
2120 :
2121 : // Add the newly-created element to the mesh
2122 0 : mesh.add_elem(std::move(elem));
2123 0 : }
2124 :
2125 :
2126 0 : if (GMVLib::gmv_data.datatype == ENDKEYWORD)
2127 : {
2128 : // There isn't a cell to read, so we just return
2129 0 : return;
2130 : }
2131 :
2132 : #endif
2133 : }
2134 :
2135 :
2136 0 : ElemType GMVIO::gmv_elem_to_libmesh_elem(std::string elemname)
2137 : {
2138 : // Erase any whitespace padding in name coming from gmv before performing comparison.
2139 0 : elemname.erase(std::remove_if(elemname.begin(), elemname.end(),
2140 0 : [](unsigned char const c){return std::isspace(c);}),
2141 0 : elemname.end());
2142 0 : return libmesh_map_find(_reading_element_map, elemname);
2143 : }
2144 :
2145 :
2146 :
2147 :
2148 0 : void GMVIO::copy_nodal_solution(EquationSystems & es)
2149 : {
2150 : // Check for easy return if there isn't any nodal data
2151 0 : if (_nodal_data.empty())
2152 : {
2153 0 : libMesh::err << "Unable to copy nodal solution: No nodal "
2154 0 : << "solution has been read in from file." << std::endl;
2155 0 : return;
2156 : }
2157 :
2158 : // Be sure there is at least one system
2159 0 : libmesh_assert (es.n_systems());
2160 :
2161 : // Keep track of variable names which have been found and
2162 : // copied already. This could be used to prevent us from
2163 : // e.g. copying the same var into 2 different systems ...
2164 : // but this seems unlikely. Also, it is used to tell if
2165 : // any variables which were read in were not successfully
2166 : // copied to the EquationSystems.
2167 0 : std::set<std::string> vars_copied;
2168 :
2169 : // For each entry in the nodal data map, try to find a system
2170 : // that has the same variable key name.
2171 0 : for (auto sys : make_range(es.n_systems()))
2172 : {
2173 : // Get a generic reference to the current System
2174 0 : System & system = es.get_system(sys);
2175 :
2176 : // For each var entry in the _nodal_data map, try to find
2177 : // that var in the system
2178 0 : for (const auto & [var_name, vec] : _nodal_data)
2179 : {
2180 0 : if (system.has_variable(var_name))
2181 : {
2182 : // Check if there are as many nodes in the mesh as there are entries
2183 : // in the stored nodal data vector
2184 0 : libmesh_assert_equal_to (vec.size(), MeshInput<MeshBase>::mesh().n_nodes());
2185 :
2186 0 : const unsigned int var_num = system.variable_number(var_name);
2187 :
2188 : // The only type of nodal data we can read in from GMV is for
2189 : // linear LAGRANGE type elements.
2190 0 : const FEType & fe_type = system.variable_type(var_num);
2191 0 : if ((fe_type.order.get_order() != FIRST) || (fe_type.family != LAGRANGE))
2192 : {
2193 0 : libMesh::err << "Only FIRST-order LAGRANGE variables can be read from GMV files. "
2194 0 : << "Skipping variable " << var_name << std::endl;
2195 0 : break;
2196 : }
2197 :
2198 :
2199 : // Loop over the stored vector's entries, inserting them into
2200 : // the System's solution if appropriate.
2201 0 : for (dof_id_type i=0,
2202 0 : sz = cast_int<dof_id_type>(vec.size());
2203 0 : i != sz; ++i)
2204 : {
2205 : // Since this var came from a GMV file, the index i corresponds to
2206 : // the (single) DOF value of the current variable for node i.
2207 : const unsigned int dof_index =
2208 0 : MeshInput<MeshBase>::mesh().node_ptr(i)->dof_number(sys, // system #
2209 : var_num, // var #
2210 : 0); // component #, always zero for LAGRANGE
2211 :
2212 : // If the dof_index is local to this processor, set the value
2213 0 : if ((dof_index >= system.solution->first_local_index()) &&
2214 0 : (dof_index < system.solution->last_local_index()))
2215 0 : system.solution->set (dof_index, vec[i]);
2216 : } // end loop over my GMVIO's copy of the solution
2217 :
2218 : // Add the most recently copied var to the set of copied vars
2219 0 : vars_copied.insert (var_name);
2220 : } // end if (system.has_variable)
2221 : } // end for loop over _nodal_data
2222 :
2223 : // Communicate parallel values before going to the next system.
2224 0 : system.solution->close();
2225 0 : system.update();
2226 : } // end loop over all systems
2227 :
2228 :
2229 :
2230 : // Warn the user if any GMV variables were not successfully copied over to the EquationSystems object
2231 0 : for (const auto & pr : _nodal_data)
2232 0 : if (vars_copied.find(pr.first) == vars_copied.end())
2233 0 : libMesh::err << "Warning: Variable "
2234 0 : << pr.first
2235 0 : << " was not copied to the EquationSystems object."
2236 0 : << std::endl;
2237 : }
2238 :
2239 : } // namespace libMesh
|
