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