ensight_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // libMesh includes
21 #include "libmesh/dof_map.h"
22 #include "libmesh/ensight_io.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/libmesh.h"
26 #include "libmesh/system.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/enum_elem_type.h"
29 #include "libmesh/int_range.h"
30 
31 // C++ includes
32 #include <sstream>
33 #include <fstream>
34 #include <string>
35 #include <iomanip>
36 
37 namespace libMesh
38 {
39 
40 // Initialize the static data members by calling the static build functions.
41 std::map<ElemType, std::string> EnsightIO::_element_map = EnsightIO::build_element_map();
42 
43 // Static function used to build the _element_map.
44 std::map<ElemType, std::string> EnsightIO::build_element_map()
45 {
46  std::map<ElemType, std::string> ret;
47  ret[EDGE2] = "bar2";
48  ret[EDGE3] = "bar3";
49  ret[QUAD4] = "quad4";
50  ret[QUAD8] = "quad8";
51  // ret[QUAD9] = "quad9"; // not supported
52  ret[TRI3] = "tria3";
53  ret[TRI6] = "tria6";
54  ret[TET4] = "tetra4";
55  ret[TET10] = "tetra10";
56  ret[HEX8] = "hexa8";
57  ret[HEX20] = "hexa20";
58  // ret[HEX27] = "HEX27"; // not supported
59  ret[PYRAMID5] = "pyramid5";
60  return ret;
61 }
62 
63 
64 EnsightIO::EnsightIO (const std::string & filename,
65  const EquationSystems & eq) :
66  MeshOutput<MeshBase> (eq.get_mesh()),
67  _equation_systems(eq)
68 {
70  _ensight_file_name = filename;
71  else
72  {
73  std::ostringstream tmp_file;
74  tmp_file << filename << "_rank" << _equation_systems.processor_id();
75  _ensight_file_name = tmp_file.str();
76  }
77 }
78 
79 
80 
81 void EnsightIO::add_vector (const std::string & system_name,
82  const std::string & vec_description,
83  const std::string & u,
84  const std::string & v)
85 {
86  libmesh_assert (_equation_systems.has_system(system_name));
87  libmesh_assert (_equation_systems.get_system(system_name).has_variable(u));
88  libmesh_assert (_equation_systems.get_system(system_name).has_variable(v));
89 
90  Vectors vec;
91  vec.description = vec_description;
92  vec.components.push_back(u);
93  vec.components.push_back(v);
94 
95  _system_vars_map[system_name].EnsightVectors.push_back(vec);
96 }
97 
98 
99 
100 void EnsightIO::add_vector (const std::string & system_name,
101  const std::string & vec_name,
102  const std::string & u,
103  const std::string & v,
104  const std::string & w)
105 {
106  libmesh_assert(_equation_systems.has_system(system_name));
107  libmesh_assert(_equation_systems.get_system(system_name).has_variable(u));
108  libmesh_assert(_equation_systems.get_system(system_name).has_variable(v));
109  libmesh_assert(_equation_systems.get_system(system_name).has_variable(w));
110 
111  Vectors vec;
112  vec.description = vec_name;
113  vec.components.push_back(u);
114  vec.components.push_back(v);
115  vec.components.push_back(w);
116  _system_vars_map[system_name].EnsightVectors.push_back(vec);
117 }
118 
119 
120 
121 void EnsightIO::add_scalar(const std::string & system_name,
122  const std::string & scl_description,
123  const std::string & s)
124 {
125  libmesh_assert(_equation_systems.has_system(system_name));
126  libmesh_assert(_equation_systems.get_system(system_name).has_variable(s));
127 
128  Scalars scl;
129  scl.description = scl_description;
130  scl.scalar_name = s;
131 
132  _system_vars_map[system_name].EnsightScalars.push_back(scl);
133 }
134 
135 
136 
137 // This method must be implemented as it is pure virtual in
138 // the MeshOutput base class.
139 void EnsightIO::write (const std::string & name)
140 {
141  // We may need to gather a DistributedMesh to output it, making that
142  // const qualifier in our constructor a dirty lie
143  MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
144 
146  this->write();
147 }
148 
149 
150 
152 {
153  this->write_ascii(time);
154  this->write_case();
155 }
156 
157 
158 
160 {
161  _time_steps.push_back(time);
162 
163  this->write_geometry_ascii();
164  this->write_solution_ascii();
165 }
166 
167 
168 
170 {
171  std::ostringstream file;
172  file << _ensight_file_name
173  << ".geo"
174  << std::setw(3)
175  << std::setprecision(0)
176  << std::setfill('0')
177  << std::right
178  << _time_steps.size()-1;
179 
180  // Open a stream to write the mesh
181  std::ofstream mesh_stream(file.str().c_str());
182 
183  mesh_stream << "EnSight Gold Geometry File Format\n";
184  mesh_stream << "Generated by \n";
185  mesh_stream << "node id off\n";
186  mesh_stream << "element id given\n";
187  mesh_stream << "part\n";
188  mesh_stream << std::setw(10) << 1 << "\n";
189  mesh_stream << "uns-elements\n";
190  mesh_stream << "coordinates\n";
191 
192  // mapping between nodal index and your coordinates
193  std::map<int, Point> mesh_nodes_map;
194 
195  // Map for grouping elements of the same type
196  std::map<ElemType, std::vector<const Elem *>> 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  for (unsigned direction=0; direction<3; ++direction)
216  {
217  int i = 1;
218  for (const auto & pr : mesh_nodes_map)
219  {
220  mesh_stream << std::setw(12)
221  << std::setprecision(5)
222  << std::scientific
223  << pr.second(direction)
224  << "\n";
225  ensight_node_index[pr.first] = i++;
226  }
227  }
228 
229  // Write parts
230  for (const auto & pr : ensight_parts_map)
231  {
232  // Look up this ElemType in the map, error if not present.
233  auto name_it = _element_map.find(pr.first);
234  if (name_it == _element_map.end())
235  libmesh_error_msg("Error: Unsupported ElemType " << pr.first << " for EnsightIO.");
236 
237  // Write element type
238  mesh_stream << "\n" << name_it->second << "\n";
239 
240  const std::vector<const Elem *> & elem_ref = pr.second;
241 
242  // Write number of element
243  mesh_stream << std::setw(10) << elem_ref.size() << "\n";
244 
245  // Write element id
246  for (const auto & elem : elem_ref)
247  mesh_stream << std::setw(10) << elem->id() << "\n";
248 
249  // Write connectivity
250  for (auto i : index_range(elem_ref))
251  {
252  for (const auto & node : elem_ref[i]->node_ref_range())
253  {
254  // tests!
255  if (pr.first == QUAD9 && i==4)
256  continue;
257 
258  // tests!
259  if (pr.first == HEX27 &&
260  (i==4 || i ==10 || i == 12 ||
261  i == 13 || i ==14 || i == 16 || i == 22))
262  continue;
263 
264  mesh_stream << std::setw(10) << ensight_node_index[node.id()];
265  }
266  mesh_stream << "\n";
267  }
268  }
269 }
270 
271 
272 
273 
274 
276 {
277  std::ostringstream case_file;
278  case_file << _ensight_file_name << ".case";
279 
280  // Open a stream for writing the case file.
281  std::ofstream case_stream(case_file.str().c_str());
282 
283  case_stream << "FORMAT\n";
284  case_stream << "type: ensight gold\n\n";
285  case_stream << "GEOMETRY\n";
286  case_stream << "model: 1 " << _ensight_file_name << ".geo" << "***\n";
287 
288  // Write Variable per node section
289  if (!_system_vars_map.empty())
290  case_stream << "\n\nVARIABLE\n";
291 
292  for (const auto & pr : _system_vars_map)
293  {
294  for (const auto & scalar : pr.second.EnsightScalars)
295  case_stream << "scalar per node: 1 "
296  << scalar.description << " "
297  << _ensight_file_name << "_" << scalar.scalar_name << ".scl***\n";
298 
299  for (const auto & vec : pr.second.EnsightVectors)
300  case_stream << "vector per node: 1 "
301  << vec.description << " "
302  << _ensight_file_name << "_" << vec.description << ".vec***\n";
303 
304  // Write time step section
305  if (_time_steps.size() != 0)
306  {
307  case_stream << "\n\nTIME\n";
308  case_stream << "time set: 1\n";
309  case_stream << "number of steps: " << std::setw(10) << _time_steps.size() << "\n";
310  case_stream << "filename start number: " << std::setw(10) << 0 << "\n";
311  case_stream << "filename increment: " << std::setw(10) << 1 << "\n";
312  case_stream << "time values:\n";
313  for (const auto & time : _time_steps)
314  case_stream << std::setw(12) << std::setprecision(5) << std::scientific << time << "\n";
315  }
316  }
317 }
318 
319 
320 // Write scalar and vector solution
322 {
323  for (const auto & pr : _system_vars_map)
324  {
325  for (const auto & scalar : pr.second.EnsightScalars)
326  this->write_scalar_ascii(pr.first,
327  scalar.scalar_name);
328 
329  for (const auto & vec : pr.second.EnsightVectors)
330  this->write_vector_ascii(pr.first,
331  vec.components,
332  vec.description);
333  }
334 }
335 
336 
337 void EnsightIO::write_scalar_ascii(const std::string & sys,
338  const std::string & var_name)
339 {
340  // Construct scalar variable filename
341  std::ostringstream scl_file;
342  scl_file << _ensight_file_name
343  << "_"
344  << var_name
345  << ".scl"
346  << std::setw(3)
347  << std::setprecision(0)
348  << std::setfill('0')
349  << std::right
350  << _time_steps.size()-1;
351 
352  // Open a stream and start writing scalar variable info.
353  std::ofstream scl_stream(scl_file.str().c_str());
354  scl_stream << "Per node scalar value\n";
355  scl_stream << "part\n";
356  scl_stream << std::setw(10) << 1 << "\n";
357  scl_stream << "coordinates\n";
358 
359  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
360  const unsigned int dim = the_mesh.mesh_dimension();
361  const System & system = _equation_systems.get_system(sys);
362  const DofMap & dof_map = system.get_dof_map();
363  int var = system.variable_number(var_name);
364 
365  std::vector<dof_id_type> dof_indices_scl;
366 
367  // Map from node id -> solution value. We end up just writing this
368  // map out in order, not sure what would happen if there were holes
369  // in the numbering...
370  std::map<int, Real> local_soln;
371 
372  std::vector<Number> elem_soln;
373  std::vector<Number> nodal_soln;
374 
375  // Loop over active local elements, construct the nodal solution, and write it to file.
376  for (const auto & elem : the_mesh.active_local_element_ptr_range())
377  {
378  const FEType & fe_type = system.variable_type(var);
379 
380  dof_map.dof_indices (elem, dof_indices_scl, var);
381 
382  elem_soln.resize(dof_indices_scl.size());
383 
384  for (auto i : index_range(dof_indices_scl))
385  elem_soln[i] = system.current_solution(dof_indices_scl[i]);
386 
387  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln, nodal_soln);
388 
389  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
390 
391 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
392  libmesh_error_msg("Complex-valued Ensight output not yet supported");
393 #endif
394 
395  for (unsigned int n=0; n<elem->n_nodes(); n++)
396  local_soln[elem->node_id(n)] = libmesh_real(nodal_soln[n]);
397  }
398 
399  for (const auto & pr : local_soln)
400  scl_stream << std::setw(12)
401  << std::setprecision(5)
402  << std::scientific
403  << pr.second
404  << "\n";
405 }
406 
407 
408 void EnsightIO::write_vector_ascii(const std::string & sys,
409  const std::vector<std::string> & vec,
410  const std::string & var_name)
411 {
412  // Construct vector variable filename
413  std::ostringstream vec_file;
414  vec_file << _ensight_file_name
415  << "_"
416  << var_name
417  << ".vec"
418  << std::setw(3)
419  << std::setprecision(0)
420  << std::setfill('0')
421  << std::right
422  << _time_steps.size()-1;
423 
424  // Open a stream and start writing vector variable info.
425  std::ofstream vec_stream(vec_file.str().c_str());
426  vec_stream << "Per vector per value\n";
427  vec_stream << "part\n";
428  vec_stream << std::setw(10) << 1 << "\n";
429  vec_stream << "coordinates\n";
430 
431  // Get a constant reference to the mesh object.
432  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
433 
434  // The dimension that we are running
435  const unsigned int dim = the_mesh.mesh_dimension();
436 
437  const System & system = _equation_systems.get_system(sys);
438 
439  const DofMap & dof_map = system.get_dof_map();
440 
441  const unsigned int u_var = system.variable_number(vec[0]);
442  const unsigned int v_var = system.variable_number(vec[1]);
443  const unsigned int w_var = (dim==3) ? system.variable_number(vec[2]) : 0;
444 
445  std::vector<dof_id_type> dof_indices_u;
446  std::vector<dof_id_type> dof_indices_v;
447  std::vector<dof_id_type> dof_indices_w;
448 
449  // Map from node id -> solution value. We end up just writing this
450  // map out in order, not sure what would happen if there were holes
451  // in the numbering...
452  std::map<int,std::vector<Real>> local_soln;
453 
454  // Now we will loop over all the elements in the mesh.
455  for (const auto & elem : the_mesh.active_local_element_ptr_range())
456  {
457  const FEType & fe_type = system.variable_type(u_var);
458 
459  dof_map.dof_indices (elem, dof_indices_u, u_var);
460  dof_map.dof_indices (elem, dof_indices_v, v_var);
461  if (dim==3)
462  dof_map.dof_indices (elem, dof_indices_w, w_var);
463 
464  std::vector<Number> elem_soln_u;
465  std::vector<Number> elem_soln_v;
466  std::vector<Number> elem_soln_w;
467 
468  std::vector<Number> nodal_soln_u;
469  std::vector<Number> nodal_soln_v;
470  std::vector<Number> nodal_soln_w;
471 
472  elem_soln_u.resize(dof_indices_u.size());
473  elem_soln_v.resize(dof_indices_v.size());
474  if (dim == 3)
475  elem_soln_w.resize(dof_indices_w.size());
476 
477  for (auto i : index_range(dof_indices_u))
478  {
479  elem_soln_u[i] = system.current_solution(dof_indices_u[i]);
480  elem_soln_v[i] = system.current_solution(dof_indices_v[i]);
481  if (dim==3)
482  elem_soln_w[i] = system.current_solution(dof_indices_w[i]);
483  }
484 
485  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_u, nodal_soln_u);
486  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_v, nodal_soln_v);
487  if (dim == 3)
488  FEInterface::nodal_soln (dim, fe_type, elem, elem_soln_w, nodal_soln_w);
489 
490  libmesh_assert_equal_to (nodal_soln_u.size(), elem->n_nodes());
491  libmesh_assert_equal_to (nodal_soln_v.size(), elem->n_nodes());
492 
493 #ifdef LIBMESH_ENABLE_COMPLEX
494  libmesh_error_msg("Complex-valued Ensight output not yet supported");
495 #endif
496 
497  for (const auto & n : elem->node_index_range())
498  {
499  std::vector<Real> node_vec(3);
500  node_vec[0] = libmesh_real(nodal_soln_u[n]);
501  node_vec[1] = libmesh_real(nodal_soln_v[n]);
502  node_vec[2] = 0.0;
503  if (dim==3)
504  node_vec[2] = libmesh_real(nodal_soln_w[n]);
505  local_soln[elem->node_id(n)] = node_vec;
506  }
507  }
508 
509  for (unsigned dir=0; dir<3; ++dir)
510  {
511  for (const auto & pr : local_soln)
512  vec_stream << std::setw(12)
513  << std::scientific
514  << std::setprecision(5)
515  << pr.second[dir]
516  << "\n";
517  }
518 }
519 
520 } // namespace libMesh
T libmesh_real(T a)
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
const MeshBase & mesh() const
Definition: mesh_output.h:234
void write_solution_ascii()
Definition: ensight_io.C:321
void write_geometry_ascii()
Definition: ensight_io.C:169
std::string _ensight_file_name
Definition: ensight_io.h:152
Manages multiples systems of equations.
EnsightIO(const std::string &filename, const EquationSystems &eq)
Definition: ensight_io.C:64
static std::map< ElemType, std::string > build_element_map()
Definition: ensight_io.C:44
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1930
const T_sys & get_system(const std::string &name) const
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
PetscBool eq
std::map< std::string, SystemVars > _system_vars_map
Definition: ensight_io.h:156
void write_scalar_ascii(const std::string &sys, const std::string &var)
Definition: ensight_io.C:337
void add_vector(const std::string &system, const std::string &vec_description, const std::string &u, const std::string &v)
Definition: ensight_io.C:81
Base class for Mesh.
Definition: mesh_base.h:77
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:194
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
processor_id_type n_processors() const
static std::map< ElemType, std::string > _element_map
Definition: ensight_io.h:162
void write_vector_ascii(const std::string &sys, const std::vector< std::string > &vec, const std::string &var_name)
Definition: ensight_io.C:408
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1243
void write(Real time=0)
Definition: ensight_io.C:151
const EquationSystems & _equation_systems
Definition: ensight_io.h:159
void add_scalar(const std::string &system, const std::string &scalar_description, const std::string &s)
Definition: ensight_io.C:121
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2183
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Temporarily serializes a DistributedMesh for output.
unsigned int mesh_dimension() const
Definition: mesh_base.C:126
bool has_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)
Definition: fe_interface.C:550
std::vector< Real > _time_steps
Definition: ensight_io.h:153
processor_id_type processor_id() const
const DofMap & get_dof_map() const
Definition: system.h:2049
void write_ascii(Real time=0)
Definition: ensight_io.C:159