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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo