libMesh
ensight_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 
19 
20 #include <sstream>
21 #include <fstream>
22 #include <string>
23 #include <iomanip>
24 
25 #include "libmesh/dof_map.h"
26 #include "libmesh/ensight_io.h"
27 #include "libmesh/equation_systems.h"
28 #include "libmesh/fe_interface.h"
29 #include "libmesh/libmesh.h"
30 #include "libmesh/system.h"
31 #include "libmesh/elem.h"
32 
33 namespace libMesh
34 {
35 
36 // Initialize the static data members by calling the static build functions.
37 std::map<ElemType, std::string> EnsightIO::_element_map = EnsightIO::build_element_map();
38 
39 // Static function used to build the _element_map.
40 std::map<ElemType, std::string> EnsightIO::build_element_map()
41 {
42  std::map<ElemType, std::string> ret;
43  ret[EDGE2] = "bar2";
44  ret[EDGE3] = "bar3";
45  ret[QUAD4] = "quad4";
46  ret[QUAD8] = "quad8";
47  // ret[QUAD9] = "quad9"; // not supported
48  ret[TRI3] = "tria3";
49  ret[TRI6] = "tria6";
50  ret[TET4] = "tetra4";
51  ret[TET10] = "tetra10";
52  ret[HEX8] = "hexa8";
53  ret[HEX20] = "hexa20";
54  // ret[HEX27] = "HEX27"; // not supported
55  ret[PYRAMID5] = "pyramid5";
56  return ret;
57 }
58 
59 
60 EnsightIO::EnsightIO (const std::string & filename,
61  const EquationSystems & eq) :
62  MeshOutput<MeshBase> (eq.get_mesh()),
64 {
66  _ensight_file_name = filename;
67  else
68  {
69  std::ostringstream tmp_file;
70  tmp_file << filename << "_rank" << _equation_systems.processor_id();
71  _ensight_file_name = tmp_file.str();
72  }
73 }
74 
75 
76 
77 void EnsightIO::add_vector (const std::string & system_name,
78  const std::string & vec_description,
79  const std::string & u,
80  const std::string & v)
81 {
83  libmesh_assert (_equation_systems.get_system(system_name).has_variable(u));
84  libmesh_assert (_equation_systems.get_system(system_name).has_variable(v));
85 
86  Vectors vec;
87  vec.description = vec_description;
88  vec.components.push_back(u);
89  vec.components.push_back(v);
90 
91  _system_vars_map[system_name].EnsightVectors.push_back(vec);
92 }
93 
94 
95 
96 void EnsightIO::add_vector (const std::string & system_name,
97  const std::string & vec_name,
98  const std::string & u,
99  const std::string & v,
100  const std::string & w)
101 {
103  libmesh_assert(_equation_systems.get_system(system_name).has_variable(u));
104  libmesh_assert(_equation_systems.get_system(system_name).has_variable(v));
105  libmesh_assert(_equation_systems.get_system(system_name).has_variable(w));
106 
107  Vectors vec;
108  vec.description = vec_name;
109  vec.components.push_back(u);
110  vec.components.push_back(v);
111  vec.components.push_back(w);
112  _system_vars_map[system_name].EnsightVectors.push_back(vec);
113 }
114 
115 
116 
117 void EnsightIO::add_scalar(const std::string & system_name,
118  const std::string & scl_description,
119  const std::string & s)
120 {
122  libmesh_assert(_equation_systems.get_system(system_name).has_variable(s));
123 
124  Scalars scl;
125  scl.description = scl_description;
126  scl.scalar_name = s;
127 
128  _system_vars_map[system_name].EnsightScalars.push_back(scl);
129 }
130 
131 
132 
133 // This method must be implemented as it is pure virtual in
134 // the MeshOutput base class.
135 void EnsightIO::write (const std::string & name)
136 {
137  // We may need to gather a DistributedMesh to output it, making that
138  // const qualifier in our constructor a dirty lie
139  MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
140 
142  this->write();
143 }
144 
145 
146 
148 {
149  this->write_ascii(time);
150  this->write_case();
151 }
152 
153 
154 
156 {
157  _time_steps.push_back(time);
158 
159  this->write_geometry_ascii();
160  this->write_solution_ascii();
161 }
162 
163 
164 
166 {
167  std::ostringstream file;
168  file << _ensight_file_name
169  << ".geo"
170  << std::setw(3)
171  << std::setprecision(0)
172  << std::setfill('0')
173  << std::right
174  << _time_steps.size()-1;
175 
176  // Open a stream to write the mesh
177  std::ofstream mesh_stream(file.str().c_str());
178 
179  mesh_stream << "EnSight Gold Geometry File Format\n";
180  mesh_stream << "Generated by \n";
181  mesh_stream << "node id off\n";
182  mesh_stream << "element id given\n";
183  mesh_stream << "part\n";
184  mesh_stream << std::setw(10) << 1 << "\n";
185  mesh_stream << "uns-elements\n";
186  mesh_stream << "coordinates\n";
187 
188  // mapping between nodal index and your coordinates
189  typedef std::map<int, Point> mesh_nodes_map_t;
190  typedef mesh_nodes_map_t::iterator mesh_nodes_iterator;
191  mesh_nodes_map_t mesh_nodes_map;
192 
193  // Map for grouping elements of the same type
194  typedef std::map<ElemType, std::vector<const Elem *>> ensight_parts_map_t;
195  typedef ensight_parts_map_t::iterator ensight_parts_iterator;
196  ensight_parts_map_t ensight_parts_map;
197 
198  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
199 
200  // Construct the various required maps
201  for (const auto & elem : the_mesh.active_local_element_ptr_range())
202  {
203  ensight_parts_map[elem->type()].push_back(elem);
204 
205  for (unsigned int i = 0; i < elem->n_nodes(); i++)
206  mesh_nodes_map[elem->node_id(i)] = elem->point(i);
207  }
208 
209  // Write number of local points
210  mesh_stream << std::setw(10) << mesh_nodes_map.size() << "\n";
211 
212  // write x, y, and z node positions, build mapping between
213  // ensight and libmesh node numbers.
214  std::map <int, int> ensight_node_index;
215  {
216  mesh_nodes_iterator no_it = mesh_nodes_map.begin();
217  const mesh_nodes_iterator no_end_it = mesh_nodes_map.end();
218 
219  for (unsigned direction=0; direction<3; ++direction)
220  {
221  for (int i = 1; no_it != no_end_it; ++no_it, i++)
222  {
223  mesh_stream << std::setw(12)
224  << std::setprecision(5)
225  << std::scientific
226  << no_it->second(direction)
227  << "\n";
228  ensight_node_index[no_it->first] = i;
229  }
230 
231  // Reset iterator to the beginning of the map
232  no_it = mesh_nodes_map.begin();
233  }
234  }
235 
236  // Write parts
237  {
238  ensight_parts_iterator parts_it = ensight_parts_map.begin();
239  const ensight_parts_iterator end_parts_it = ensight_parts_map.end();
240 
241  for (; parts_it != end_parts_it; ++parts_it)
242  {
243  // Look up this ElemType in the map, error if not present.
244  std::map<ElemType, std::string>::iterator name_it = _element_map.find(parts_it->first);
245  if (name_it == _element_map.end())
246  libmesh_error_msg("Error: Unsupported ElemType " << parts_it->first << " for EnsightIO.");
247 
248  // Write element type
249  mesh_stream << "\n" << name_it->second << "\n";
250 
251  std::vector<const Elem *> elem_ref = parts_it->second;
252 
253  // Write number of element
254  mesh_stream << std::setw(10) << elem_ref.size() << "\n";
255 
256  // Write element id
257  for (std::size_t i = 0; i < elem_ref.size(); i++)
258  mesh_stream << std::setw(10) << elem_ref[i]->id() << "\n";
259 
260  // Write connectivity
261  for (std::size_t i = 0; i < elem_ref.size(); i++)
262  {
263  for (unsigned int j = 0; j < elem_ref[i]->n_nodes(); j++)
264  {
265  // tests!
266  if (parts_it->first == QUAD9 && i==4)
267  continue;
268 
269  // tests!
270  if (parts_it->first == HEX27 &&
271  (i==4 || i ==10 || i == 12 ||
272  i == 13 || i ==14 || i == 16 || i == 22))
273  continue;
274 
275  mesh_stream << std::setw(10) << ensight_node_index[elem_ref[i]->node_id(j)];
276  }
277  mesh_stream << "\n";
278  }
279  }
280  }
281 }
282 
283 
284 
285 
286 
288 {
289  std::ostringstream case_file;
290  case_file << _ensight_file_name << ".case";
291 
292  // Open a stream for writing the case file.
293  std::ofstream case_stream(case_file.str().c_str());
294 
295  case_stream << "FORMAT\n";
296  case_stream << "type: ensight gold\n\n";
297  case_stream << "GEOMETRY\n";
298  case_stream << "model: 1 " << _ensight_file_name << ".geo" << "***\n";
299 
300  system_vars_map_t::iterator sys_it = _system_vars_map.begin();
301  const system_vars_map_t::iterator sys_end = _system_vars_map.end();
302 
303  // Write Variable per node section
304  if (sys_it != sys_end)
305  case_stream << "\n\nVARIABLE\n";
306 
307  for (; sys_it != sys_end; ++sys_it)
308  {
309  for (std::size_t i=0; i < sys_it->second.EnsightScalars.size(); i++)
310  {
311  Scalars scalar = sys_it->second.EnsightScalars[i];
312  case_stream << "scalar per node: 1 "
313  << scalar.description << " "
314  << _ensight_file_name << "_" << scalar.scalar_name << ".scl***\n";
315  }
316 
317  for (std::size_t i=0; i < sys_it->second.EnsightVectors.size(); i++)
318  {
319  Vectors vec = sys_it->second.EnsightVectors[i];
320  case_stream << "vector per node: 1 "
321  << vec.description << " "
322  << _ensight_file_name << "_" << vec.description << ".vec***\n";
323  }
324 
325  // Write time step section
326  if (_time_steps.size() != 0)
327  {
328  case_stream << "\n\nTIME\n";
329  case_stream << "time set: 1\n";
330  case_stream << "number of steps: " << std::setw(10) << _time_steps.size() << "\n";
331  case_stream << "filename start number: " << std::setw(10) << 0 << "\n";
332  case_stream << "filename increment: " << std::setw(10) << 1 << "\n";
333  case_stream << "time values:\n";
334  for (std::size_t i = 0; i < _time_steps.size(); i++)
335  case_stream << std::setw(12) << std::setprecision(5) << std::scientific << _time_steps[i] << "\n";
336  }
337  }
338 }
339 
340 
341 // Write scalar and vector solution
343 {
344  system_vars_map_t::iterator sys_it = _system_vars_map.begin();
345  const system_vars_map_t::iterator sys_end = _system_vars_map.end();
346 
347  for (; sys_it != sys_end; ++sys_it)
348  {
349  for (std::size_t i = 0; i < sys_it->second.EnsightScalars.size(); i++)
350  this->write_scalar_ascii(sys_it->first,
351  sys_it->second.EnsightScalars[i].scalar_name);
352 
353  for (std::size_t i = 0; i < sys_it->second.EnsightVectors.size(); i++)
354  this->write_vector_ascii(sys_it->first,
355  sys_it->second.EnsightVectors[i].components,
356  sys_it->second.EnsightVectors[i].description);
357  }
358 }
359 
360 
361 void EnsightIO::write_scalar_ascii(const std::string & sys,
362  const std::string & var_name)
363 {
364  // Construct scalar variable filename
365  std::ostringstream scl_file;
366  scl_file << _ensight_file_name
367  << "_"
368  << var_name
369  << ".scl"
370  << std::setw(3)
371  << std::setprecision(0)
372  << std::setfill('0')
373  << std::right
374  << _time_steps.size()-1;
375 
376  // Open a stream and start writing scalar variable info.
377  std::ofstream scl_stream(scl_file.str().c_str());
378  scl_stream << "Per node scalar value\n";
379  scl_stream << "part\n";
380  scl_stream << std::setw(10) << 1 << "\n";
381  scl_stream << "coordinates\n";
382 
383  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
384  const unsigned int dim = the_mesh.mesh_dimension();
385  const System & system = _equation_systems.get_system(sys);
386  const DofMap & dof_map = system.get_dof_map();
387  int var = system.variable_number(var_name);
388 
389  std::vector<dof_id_type> dof_indices_scl;
390 
391  // Map from node id -> solution value. We end up just writing this
392  // map out in order, not sure what would happen if there were holes
393  // in the numbering...
394  typedef std::map<int, Real> map_local_soln;
395  typedef map_local_soln::iterator local_soln_iterator;
396  map_local_soln local_soln;
397 
398  std::vector<Number> elem_soln;
399  std::vector<Number> nodal_soln;
400 
401  // Loop over active local elements, construct the nodal solution, and write it to file.
402  for (const auto & elem : the_mesh.active_local_element_ptr_range())
403  {
404  const FEType & fe_type = system.variable_type(var);
405 
406  dof_map.dof_indices (elem, dof_indices_scl, var);
407 
408  elem_soln.resize(dof_indices_scl.size());
409 
410  for (std::size_t i = 0; i < dof_indices_scl.size(); i++)
411  elem_soln[i] = system.current_solution(dof_indices_scl[i]);
412 
413  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln, nodal_soln);
414 
415  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
416 
417 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
418  libmesh_error_msg("Complex-valued Ensight output not yet supported");
419 #endif
420 
421  for (unsigned int n=0; n<elem->n_nodes(); n++)
422  local_soln[elem->node_id(n)] = libmesh_real(nodal_soln[n]);
423  }
424 
425  {
426  local_soln_iterator it = local_soln.begin();
427  const local_soln_iterator it_end = local_soln.end();
428  for ( ; it != it_end; ++it)
429  scl_stream << std::setw(12)
430  << std::setprecision(5)
431  << std::scientific
432  << it->second
433  << "\n";
434  }
435 }
436 
437 
438 void EnsightIO::write_vector_ascii(const std::string & sys,
439  const std::vector<std::string> & vec,
440  const std::string & var_name)
441 {
442  // Construct vector variable filename
443  std::ostringstream vec_file;
444  vec_file << _ensight_file_name
445  << "_"
446  << var_name
447  << ".vec"
448  << std::setw(3)
449  << std::setprecision(0)
450  << std::setfill('0')
451  << std::right
452  << _time_steps.size()-1;
453 
454  // Open a stream and start writing vector variable info.
455  std::ofstream vec_stream(vec_file.str().c_str());
456  vec_stream << "Per vector per value\n";
457  vec_stream << "part\n";
458  vec_stream << std::setw(10) << 1 << "\n";
459  vec_stream << "coordinates\n";
460 
461  // Get a constant reference to the mesh object.
462  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
463 
464  // The dimension that we are running
465  const unsigned int dim = the_mesh.mesh_dimension();
466 
467  const System & system = _equation_systems.get_system(sys);
468 
469  const DofMap & dof_map = system.get_dof_map();
470 
471  const unsigned int u_var = system.variable_number(vec[0]);
472  const unsigned int v_var = system.variable_number(vec[1]);
473  const unsigned int w_var = (dim==3) ? system.variable_number(vec[2]) : 0;
474 
475  std::vector<dof_id_type> dof_indices_u;
476  std::vector<dof_id_type> dof_indices_v;
477  std::vector<dof_id_type> dof_indices_w;
478 
479  // Map from node id -> solution value. We end up just writing this
480  // map out in order, not sure what would happen if there were holes
481  // in the numbering...
482  typedef std::map<int,std::vector<Real>> map_local_soln;
483  typedef map_local_soln::iterator local_soln_iterator;
484  map_local_soln local_soln;
485 
486  // Now we will loop over all the elements in the mesh.
487  for (const auto & elem : the_mesh.active_local_element_ptr_range())
488  {
489  const FEType & fe_type = system.variable_type(u_var);
490 
491  dof_map.dof_indices (elem, dof_indices_u, u_var);
492  dof_map.dof_indices (elem, dof_indices_v, v_var);
493  if (dim==3)
494  dof_map.dof_indices (elem, dof_indices_w, w_var);
495 
496  std::vector<Number> elem_soln_u;
497  std::vector<Number> elem_soln_v;
498  std::vector<Number> elem_soln_w;
499 
500  std::vector<Number> nodal_soln_u;
501  std::vector<Number> nodal_soln_v;
502  std::vector<Number> nodal_soln_w;
503 
504  elem_soln_u.resize(dof_indices_u.size());
505  elem_soln_v.resize(dof_indices_v.size());
506  if (dim == 3)
507  elem_soln_w.resize(dof_indices_w.size());
508 
509  for (std::size_t i = 0; i < dof_indices_u.size(); i++)
510  {
511  elem_soln_u[i] = system.current_solution(dof_indices_u[i]);
512  elem_soln_v[i] = system.current_solution(dof_indices_v[i]);
513  if (dim==3)
514  elem_soln_w[i] = system.current_solution(dof_indices_w[i]);
515  }
516 
517  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_u, nodal_soln_u);
518  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_v, nodal_soln_v);
519  if (dim == 3)
520  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_w, nodal_soln_w);
521 
522  libmesh_assert_equal_to (nodal_soln_u.size(), elem->n_nodes());
523  libmesh_assert_equal_to (nodal_soln_v.size(), elem->n_nodes());
524 
525 #ifdef LIBMESH_ENABLE_COMPLEX
526  libmesh_error_msg("Complex-valued Ensight output not yet supported");
527 #endif
528 
529  for (unsigned int n=0; n<elem->n_nodes(); n++)
530  {
531  std::vector<Real> node_vec(3);
532  node_vec[0] = libmesh_real(nodal_soln_u[n]);
533  node_vec[1] = libmesh_real(nodal_soln_v[n]);
534  node_vec[2] = 0.0;
535  if (dim==3)
536  node_vec[2] = libmesh_real(nodal_soln_w[n]);
537  local_soln[elem->node_id(n)] = node_vec;
538  }
539  }
540 
541  {
542  local_soln_iterator it = local_soln.begin();
543  const local_soln_iterator it_end = local_soln.end();
544 
545  for (unsigned dir=0; dir<3; ++dir)
546  {
547  for (; it != it_end; ++it)
548  vec_stream << std::setw(12)
549  << std::scientific
550  << std::setprecision(5)
551  << it->second[dir]
552  << "\n";
553 
554  // Reset the iterator to the beginning of the map
555  it = local_soln.begin();
556  }
557  }
558 }
559 
560 } // namespace libMesh
T libmesh_real(T a)
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
void write_solution_ascii()
Definition: ensight_io.C:342
void write_geometry_ascii()
Definition: ensight_io.C:165
std::string _ensight_file_name
Definition: ensight_io.h:144
This is the EquationSystems class.
bool has_system(const std::string &name) const
EnsightIO(const std::string &filename, const EquationSystems &eq)
Constructor.
Definition: ensight_io.C:60
static std::map< ElemType, std::string > build_element_map()
Definition: ensight_io.C:40
unsigned int dim
ImplicitSystem & sys
const bool _is_parallel_format
Flag specifying whether this format is parallel-capable.
Definition: mesh_output.h:141
processor_id_type n_processors() const
PetscBool eq
void write_scalar_ascii(const std::string &sys, const std::string &var)
Definition: ensight_io.C:361
const MeshBase & mesh() const
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
The libMesh namespace provides an interface to certain functionality in the library.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
void add_vector(const std::string &system, const std::string &vec_description, const std::string &u, const std::string &v)
Tell the EnsightIO interface that the variables (u,v) constitute a vector.
Definition: ensight_io.C:77
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
const DofMap & get_dof_map() const
Definition: system.h:2030
static std::map< ElemType, std::string > _element_map
Definition: ensight_io.h:155
system_vars_map_t _system_vars_map
Definition: ensight_io.h:149
void write_vector_ascii(const std::string &sys, const std::vector< std::string > &vec, const std::string &var_name)
Definition: ensight_io.C:438
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
void write(Real time=0)
Calls write_ascii() and write_case().
Definition: ensight_io.C:147
const EquationSystems & _equation_systems
Definition: ensight_io.h:152
void add_scalar(const std::string &system, const std::string &scalar_description, const std::string &s)
Tell the EnsightIO interface to output the finite element (not SCALAR) variable named "s"...
Definition: ensight_io.C:117
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
Temporarily serialize a DistributedMesh for output; a distributed mesh is allgathered by the MeshSeri...
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
const T_sys & get_system(const std::string &name) const
static void nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
Definition: fe_interface.C:526
std::vector< Real > _time_steps
Definition: ensight_io.h:145
void write_ascii(Real time=0)
Definition: ensight_io.C:155
processor_id_type processor_id() const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917