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