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