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