exodusII_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 // C++ includes
20 #include <fstream>
21 #include <cstring>
22 #include <sstream>
23 #include <map>
24 
25 // Local includes
26 #include "libmesh/exodusII_io.h"
27 #include "libmesh/boundary_info.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/enum_elem_type.h"
30 #include "libmesh/elem.h"
33 #include "libmesh/system.h"
34 #include "libmesh/numeric_vector.h"
36 #include "libmesh/string_to_enum.h"
38 #include "libmesh/parallel_mesh.h"
39 
40 namespace libMesh
41 {
42 
43 // ------------------------------------------------------------
44 // ExodusII_IO class members
46 #ifdef LIBMESH_HAVE_EXODUS_API
47  bool single_precision
48 #else
49  bool
50 #endif
51  ) :
54  /* is_parallel_format = */ false,
55  /* serial_only_needed_on_proc_0 = */ true),
57 #ifdef LIBMESH_HAVE_EXODUS_API
58  exio_helper(new ExodusII_IO_Helper(*this, false, true, single_precision)),
59  _timestep(1),
60  _verbose(false),
61  _append(false),
62 #endif
63  _allow_empty_variables(false)
64 {
65 }
66 
67 
68 void ExodusII_IO::set_output_variables(const std::vector<std::string> & output_variables,
69  bool allow_empty)
70 {
71  _output_variables = output_variables;
72  _allow_empty_variables = allow_empty;
73 }
74 
75 
76 
77 #ifdef LIBMESH_ENABLE_DEPRECATED
79  std::string var_name,
80  unsigned int timestep)
81 {
82  libmesh_deprecated();
83  copy_nodal_solution(system, var_name, var_name, timestep);
84 }
85 #endif
86 
87 
88 
90  const EquationSystems & es,
91  const std::set<std::string> * system_names)
92 {
93  std::vector<std::string> solution_names;
94  std::vector<Number> v;
95 
96  es.build_variable_names (solution_names, nullptr, system_names);
97  es.build_discontinuous_solution_vector (v, system_names);
98 
99  this->write_nodal_data_discontinuous(name, v, solution_names);
100 }
101 
102 
103 #ifdef LIBMESH_HAVE_EXODUS_API
104 void ExodusII_IO::write_timestep_discontinuous (const std::string &fname,
105  const EquationSystems &es,
106  const int timestep,
107  const Real time,
108  const std::set<std::string> * system_names)
109 {
110  _timestep = timestep;
111  write_discontinuous_equation_systems (fname,es,system_names);
112 
114  return;
115 
116  exio_helper->write_timestep(timestep, time);
117 }
118 
119 #else
120 void ExodusII_IO::write_timestep_discontinuous (const std::string & /* fname */,
121  const EquationSystems & /* es */,
122  const int /* timestep */,
123  const Real /* time */,
124  const std::set<std::string> * /*system_names*/)
125 { libmesh_error(); }
126 #endif
127 
128 
129 // ------------------------------------------------------------
130 // When the Exodus API is present...
131 #ifdef LIBMESH_HAVE_EXODUS_API
132 
134 {
135  exio_helper->close();
136 }
137 
138 
139 
140 void ExodusII_IO::read (const std::string & fname)
141 {
142  // Get a reference to the mesh we are reading
144 
145  // Clear any existing mesh data
146  mesh.clear();
147 
148  // Keep track of what kinds of elements this file contains
149  elems_of_dimension.clear();
150  elems_of_dimension.resize(4, false);
151 
152  // Instantiate the ElementMaps interface
154 
155  // Open the exodus file in EX_READ mode
156  exio_helper->open(fname.c_str(), /*read_only=*/true);
157 
158  // Get header information from exodus file
159  exio_helper->read_header();
160 
161  // Read the QA records
162  exio_helper->read_qa_records();
163 
164  // Print header information
165  exio_helper->print_header();
166 
167  // Read nodes from the exodus file
168  exio_helper->read_nodes();
169 
170  // Reserve space for the nodes.
171  mesh.reserve_nodes(exio_helper->num_nodes);
172 
173  // Read the node number map from the Exodus file. This is
174  // required if we want to preserve the numbering of nodes as it
175  // exists in the Exodus file. If the Exodus file does not contain
176  // a node_num_map, the identity map is returned by this call.
177  exio_helper->read_node_num_map();
178 
179  // Loop over the nodes, create Nodes with local processor_id 0.
180  for (int i=0; i<exio_helper->num_nodes; i++)
181  {
182  // Use the node_num_map to get the correct ID for Exodus
183  int exodus_id = exio_helper->node_num_map[i];
184 
185  // Catch the node that was added to the mesh
186  Node * added_node = mesh.add_point (Point(exio_helper->x[i], exio_helper->y[i], exio_helper->z[i]), exodus_id-1);
187 
188  // If the Mesh assigned an ID different from what is in the
189  // Exodus file, we should probably error.
190  if (added_node->id() != static_cast<unsigned>(exodus_id-1))
191  libmesh_error_msg("Error! Mesh assigned node ID " \
192  << added_node->id() \
193  << " which is different from the (zero-based) Exodus ID " \
194  << exodus_id-1 \
195  << "!");
196  }
197 
198  // This assert is no longer valid if the nodes are not numbered
199  // sequentially starting from 1 in the Exodus file.
200  // libmesh_assert_equal_to (static_cast<unsigned int>(exio_helper->num_nodes), mesh.n_nodes());
201 
202  // Get information about all the blocks
203  exio_helper->read_block_info();
204 
205  // Reserve space for the elements
206  mesh.reserve_elem(exio_helper->num_elem);
207 
208  // Read the element number map from the Exodus file. This is
209  // required if we want to preserve the numbering of elements as it
210  // exists in the Exodus file. If the Exodus file does not contain
211  // an elem_num_map, the identity map is returned by this call.
212  exio_helper->read_elem_num_map();
213 
214  // Read in the element connectivity for each block.
215  int nelem_last_block = 0;
216 
217  // Loop over all the blocks
218  for (int i=0; i<exio_helper->num_elem_blk; i++)
219  {
220  // Read the information for block i
221  exio_helper->read_elem_in_block (i);
222  int subdomain_id = exio_helper->get_block_id(i);
223 
224  // populate the map of names
225  std::string subdomain_name = exio_helper->get_block_name(i);
226  if (!subdomain_name.empty())
227  mesh.subdomain_name(static_cast<subdomain_id_type>(subdomain_id)) = subdomain_name;
228 
229  // Set any relevant node/edge maps for this element
230  const std::string type_str (exio_helper->get_elem_type());
231  const ExodusII_IO_Helper::Conversion conv = em.assign_conversion(type_str);
232 
233  // Loop over all the faces in this block
234  int jmax = nelem_last_block+exio_helper->num_elem_this_blk;
235  for (int j=nelem_last_block; j<jmax; j++)
236  {
237  Elem * elem = Elem::build (conv.get_canonical_type()).release();
238  libmesh_assert (elem);
239  elem->subdomain_id() = static_cast<subdomain_id_type>(subdomain_id) ;
240 
241  // Use the elem_num_map to obtain the ID of this element in the Exodus file
242  int exodus_id = exio_helper->elem_num_map[j];
243 
244  // Assign this element the same ID it had in the Exodus
245  // file, but make it zero-based by subtracting 1. Note:
246  // some day we could use 1-based numbering in libmesh and
247  // thus match the Exodus numbering exactly, but at the
248  // moment libmesh is zero-based.
249  elem->set_id(exodus_id-1);
250 
251  // Record that we have seen an element of dimension elem->dim()
252  elems_of_dimension[elem->dim()] = true;
253 
254  // Catch the Elem pointer that the Mesh throws back
255  elem = mesh.add_elem (elem);
256 
257  // If the Mesh assigned an ID different from what is in the
258  // Exodus file, we should probably error.
259  if (elem->id() != static_cast<unsigned>(exodus_id-1))
260  libmesh_error_msg("Error! Mesh assigned ID " \
261  << elem->id() \
262  << " which is different from the (zero-based) Exodus ID " \
263  << exodus_id-1 \
264  << "!");
265 
266  // Set all the nodes for this element
267  for (int k=0; k<exio_helper->num_nodes_per_elem; k++)
268  {
269  // global index
270  int gi = (j-nelem_last_block)*exio_helper->num_nodes_per_elem + conv.get_node_map(k);
271 
272  // The entries in 'connect' are actually (1-based)
273  // indices into the node_num_map, so to get the right
274  // node ID we:
275  // 1.) Subtract 1 from connect[gi]
276  // 2.) Pass it through node_num_map to get the corresponding Exodus ID
277  // 3.) Subtract 1 from that, since libmesh node numbering is "zero"-based,
278  // even when the Exodus node numbering doesn't start with 1.
279  int libmesh_node_id = exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1;
280 
281  // Set the node pointer in the Elem
282  elem->set_node(k) = mesh.node_ptr(libmesh_node_id);
283  }
284  }
285 
286  // running sum of # of elements per block,
287  // (should equal total number of elements in the end)
288  nelem_last_block += exio_helper->num_elem_this_blk;
289  }
290 
291  // This assert isn't valid if the Exodus file's numbering doesn't
292  // start with 1! For example, if Exodus's elem_num_map is 21, 22,
293  // 23, 24, 25, 26, 27, 28, 29, 30, ... 84, then by the time you are
294  // done with the loop above, mesh.n_elem() will report 84 and
295  // nelem_last_block will be 64.
296  // libmesh_assert_equal_to (static_cast<unsigned>(nelem_last_block), mesh.n_elem());
297 
298  // Set the mesh dimension to the largest encountered for an element
299  for (unsigned char i=0; i!=4; ++i)
300  if (elems_of_dimension[i])
302 
303  // Read in sideset information -- this is useful for applying boundary conditions
304  {
305  // Get basic information about all sidesets
306  exio_helper->read_sideset_info();
307  int offset=0;
308  for (int i=0; i<exio_helper->num_side_sets; i++)
309  {
310  // Compute new offset
311  offset += (i > 0 ? exio_helper->num_sides_per_set[i-1] : 0);
312  exio_helper->read_sideset (i, offset);
313 
314  std::string sideset_name = exio_helper->get_side_set_name(i);
315  if (!sideset_name.empty())
317  (cast_int<boundary_id_type>(exio_helper->get_side_set_id(i)))
318  = sideset_name;
319  }
320 
321  for (auto e : index_range(exio_helper->elem_list))
322  {
323  // The numbers in the Exodus file sidesets should be thought
324  // of as (1-based) indices into the elem_num_map array. So,
325  // to get the right element ID we have to:
326  // 1.) Subtract 1 from elem_list[e] (to get a zero-based index)
327  // 2.) Pass it through elem_num_map (to get the corresponding Exodus ID)
328  // 3.) Subtract 1 from that, since libmesh is "zero"-based,
329  // even when the Exodus numbering doesn't start with 1.
330  dof_id_type libmesh_elem_id =
331  cast_int<dof_id_type>(exio_helper->elem_num_map[exio_helper->elem_list[e] - 1] - 1);
332 
333  // Set any relevant node/edge maps for this element
334  Elem & elem = mesh.elem_ref(libmesh_elem_id);
335 
336  const ExodusII_IO_Helper::Conversion conv = em.assign_conversion(elem.type());
337 
338  // Map the zero-based Exodus side numbering to the libmesh side numbering
339  unsigned int raw_side_index = exio_helper->side_list[e]-1;
340  std::size_t side_index_offset = conv.get_shellface_index_offset();
341 
342  if (raw_side_index < side_index_offset)
343  {
344  // We assume this is a "shell face"
345  int mapped_shellface = raw_side_index;
346 
347  // Check for errors
348  if (mapped_shellface == ExodusII_IO_Helper::Conversion::invalid_id)
349  libmesh_error_msg("Invalid 1-based side id: " \
350  << mapped_shellface \
351  << " detected for " \
352  << Utility::enum_to_string(elem.type()));
353 
354  // Add this (elem,shellface,id) triplet to the BoundaryInfo object.
355  mesh.get_boundary_info().add_shellface (libmesh_elem_id,
356  cast_int<unsigned short>(mapped_shellface),
357  cast_int<boundary_id_type>(exio_helper->id_list[e]));
358  }
359  else
360  {
361  unsigned int side_index = static_cast<unsigned int>(raw_side_index - side_index_offset);
362  int mapped_side = conv.get_side_map(side_index);
363 
364  // Check for errors
366  libmesh_error_msg("Invalid 1-based side id: " \
367  << side_index \
368  << " detected for " \
369  << Utility::enum_to_string(elem.type()));
370 
371  // Add this (elem,side,id) triplet to the BoundaryInfo object.
372  mesh.get_boundary_info().add_side (libmesh_elem_id,
373  cast_int<unsigned short>(mapped_side),
374  cast_int<boundary_id_type>(exio_helper->id_list[e]));
375  }
376  }
377  }
378 
379  // Read nodeset info
380  {
381  exio_helper->read_nodeset_info();
382 
383  for (int nodeset=0; nodeset<exio_helper->num_node_sets; nodeset++)
384  {
385  boundary_id_type nodeset_id =
386  cast_int<boundary_id_type>(exio_helper->nodeset_ids[nodeset]);
387 
388  std::string nodeset_name = exio_helper->get_node_set_name(nodeset);
389  if (!nodeset_name.empty())
390  mesh.get_boundary_info().nodeset_name(nodeset_id) = nodeset_name;
391 
392  exio_helper->read_nodeset(nodeset);
393 
394  for (const auto & exodus_id : exio_helper->node_list)
395  {
396  // As before, the entries in 'node_list' are 1-based
397  // indices into the node_num_map array, so we have to map
398  // them. See comment above.
399  int libmesh_node_id = exio_helper->node_num_map[exodus_id - 1] - 1;
400  mesh.get_boundary_info().add_node(cast_int<dof_id_type>(libmesh_node_id),
401  nodeset_id);
402  }
403  }
404  }
405 
406 #if LIBMESH_DIM < 3
407  if (mesh.mesh_dimension() > LIBMESH_DIM)
408  libmesh_error_msg("Cannot open dimension " \
409  << mesh.mesh_dimension() \
410  << " mesh file when configured without " \
411  << mesh.mesh_dimension() \
412  << "D support.");
413 #endif
414 }
415 
416 
417 
418 void ExodusII_IO::verbose (bool set_verbosity)
419 {
420  _verbose = set_verbosity;
421 
422  // Set the verbose flag in the helper object as well.
423  exio_helper->verbose = _verbose;
424 }
425 
426 
427 
429 {
430  exio_helper->use_mesh_dimension_instead_of_spatial_dimension(val);
431 }
432 
433 
434 
436 {
437  exio_helper->write_as_dimension(dim);
438 }
439 
440 
441 
443 {
444  libmesh_warning("This method may be deprecated in the future");
445  exio_helper->set_coordinate_offset(p);
446 }
447 
448 
449 
450 void ExodusII_IO::append(bool val)
451 {
452  _append = val;
453 }
454 
455 
456 
457 const std::vector<Real> & ExodusII_IO::get_time_steps()
458 {
459  if (!exio_helper->opened_for_reading)
460  libmesh_error_msg("ERROR, ExodusII file must be opened for reading before calling ExodusII_IO::get_time_steps()!");
461 
462  exio_helper->read_time_steps();
463  return exio_helper->time_steps;
464 }
465 
466 
467 
469 {
470  if (!exio_helper->opened_for_reading && !exio_helper->opened_for_writing)
471  libmesh_error_msg("ERROR, ExodusII file must be opened for reading or writing before calling ExodusII_IO::get_num_time_steps()!");
472 
473  exio_helper->read_num_time_steps();
474  return exio_helper->num_time_steps;
475 }
476 
477 
478 
480  std::string system_var_name,
481  std::string exodus_var_name,
482  unsigned int timestep)
483 {
484  if (!exio_helper->opened_for_reading)
485  libmesh_error_msg("ERROR, ExodusII file must be opened for reading before copying a nodal solution!");
486 
487  exio_helper->read_nodal_var_values(exodus_var_name, timestep);
488 
489  const unsigned int var_num = system.variable_number(system_var_name);
490 
491  for (dof_id_type i=0,
492  n_nodal = cast_int<dof_id_type>(exio_helper->nodal_var_values.size());
493  i != n_nodal; ++i)
494  {
495  const Node * node = MeshInput<MeshBase>::mesh().query_node_ptr(i);
496 
497  if (node && node->n_comp(system.number(), var_num) > 0)
498  {
499  dof_id_type dof_index = node->dof_number(system.number(), var_num, 0);
500 
501  // If the dof_index is local to this processor, set the value
502  if ((dof_index >= system.solution->first_local_index()) && (dof_index < system.solution->last_local_index()))
503  system.solution->set (dof_index, exio_helper->nodal_var_values[i]);
504  }
505  }
506 
507  system.solution->close();
508  system.update();
509 }
510 
511 
512 
514  std::string system_var_name,
515  std::string exodus_var_name,
516  unsigned int timestep)
517 {
518  if (system.comm().rank() == 0)
519  {
520  if (!exio_helper->opened_for_reading)
521  libmesh_error_msg("ERROR, ExodusII file must be opened for reading before copying an elemental solution!");
522 
523  // Map from element ID to elemental variable value. We need to use
524  // a map here rather than a vector (e.g. elem_var_values) since the
525  // libmesh element numbering can contain "holes". This is the case
526  // if we are reading elemental var values from an adaptively refined
527  // mesh that has not been sequentially renumbered.
528  std::map<dof_id_type, Real> elem_var_value_map;
529  exio_helper->read_elemental_var_values(exodus_var_name, timestep, elem_var_value_map);
530 
531  const unsigned int var_num = system.variable_number(system_var_name);
532  if (system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL))
533  libmesh_error_msg("Error! Trying to copy elemental solution into a variable that is not of CONSTANT MONOMIAL type.");
534 
535  std::map<dof_id_type, Real>::iterator
536  it = elem_var_value_map.begin(),
537  end = elem_var_value_map.end();
538 
539  for (; it!=end; ++it)
540  {
541  const Elem * elem = MeshInput<MeshBase>::mesh().query_elem_ptr(it->first);
542 
543  if (elem && elem->n_comp(system.number(), var_num) > 0)
544  {
545  dof_id_type dof_index = elem->dof_number(system.number(), var_num, 0);
546  system.solution->set (dof_index, it->second);
547  }
548  }
549  }
550 
551  system.solution->close();
552  system.update();
553 }
554 
555 void ExodusII_IO::read_elemental_variable(std::string elemental_var_name,
556  unsigned int timestep,
557  std::map<unsigned int, Real> & unique_id_to_value_map)
558 {
559  // Note that this function MUST be called before renumbering
560  std::map<dof_id_type, Real> elem_var_value_map;
561 
562  exio_helper->read_elemental_var_values(elemental_var_name, timestep, elem_var_value_map);
563  for (auto & pr : elem_var_value_map)
564  {
565  const Elem * elem = MeshInput<MeshBase>::mesh().query_elem_ptr(pr.first);
566  unique_id_to_value_map.insert(std::make_pair(elem->top_parent()->unique_id(), pr.second));
567  }
568 }
569 
570 void ExodusII_IO::read_global_variable(std::vector<std::string> global_var_names,
571  unsigned int timestep,
572  std::vector<Real> & global_values)
573 {
574  std::size_t size = global_var_names.size();
575  if (size == 0)
576  libmesh_error_msg("ERROR, empty list of global variables to read from the Exodus file.");
577 
578  // read the values for all global variables
579  std::vector<Real> values_from_exodus;
580  exio_helper->read_var_names(ExodusII_IO_Helper::GLOBAL);
581  exio_helper->read_global_values(values_from_exodus, timestep);
582  std::vector<std::string> global_var_names_exodus = exio_helper->global_var_names;
583 
584  global_values.clear();
585  for (std::size_t i = 0; i != size; ++i)
586  {
587  // for each global variable in global_var_names, look the corresponding one in global_var_names_from_exodus
588  // and fill global_values accordingly
589  auto it = find(global_var_names_exodus.begin(), global_var_names_exodus.end(), global_var_names[i]);
590  if (it != global_var_names_exodus.end())
591  global_values.push_back(values_from_exodus[it - global_var_names_exodus.begin()]);
592  else
593  libmesh_error_msg("ERROR, Global variable " << global_var_names[i] << \
594  " not found in Exodus file.");
595  }
596 
597 }
598 
600 {
601  // Be sure the file has been opened for writing!
602  if (MeshOutput<MeshBase>::mesh().processor_id() == 0 && !exio_helper->opened_for_writing)
603  libmesh_error_msg("ERROR, ExodusII file must be initialized before outputting element variables.");
604 
605  // This function currently only works on serialized meshes. We rely
606  // on having a reference to a non-const MeshBase object from our
607  // MeshInput parent class to construct a MeshSerializer object,
608  // similar to what is done in ExodusII_IO::write(). Note that
609  // calling ExodusII_IO::write_timestep() followed by
610  // ExodusII_IO::write_element_data() when the underlying Mesh is a
611  // DistributedMesh will result in an unnecessary additional
612  // serialization/re-parallelization step.
613  // The "true" specifies that we only need the mesh serialized to processor 0
615 
616  // To be (possibly) filled with a filtered list of variable names to output.
617  std::vector<std::string> names;
618 
619  // If _output_variables is populated, only output the monomials which are
620  // also in the _output_variables vector.
621  if (_output_variables.size() > 0)
622  {
623  std::vector<std::string> monomials;
624  const FEType type(CONSTANT, MONOMIAL);
625 
626  // Create a list of monomial variable names
627  es.build_variable_names(monomials, &type);
628 
629  // Filter that list against the _output_variables list. Note: if names is still empty after
630  // all this filtering, all the monomial variables will be gathered
631  for (const auto & var : monomials)
632  if (std::find(_output_variables.begin(), _output_variables.end(), var) != _output_variables.end())
633  names.push_back(var);
634  }
635 
636  // If we pass in a list of names to "get_solution" it'll filter the variables coming back
637  std::vector<Number> soln;
638  es.build_elemental_solution_vector(soln, names);
639 
640  // Also, store the list of subdomains on which each variable is active
641  std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
642  es.get_vars_active_subdomains(names, vars_active_subdomains);
643 
644  if (soln.empty()) // If there is nothing to write just return
645  return;
646 
647  // The data must ultimately be written block by block. This means that this data
648  // must be sorted appropriately.
649  if (MeshOutput<MeshBase>::mesh().processor_id())
650  return;
651 
653 
654 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
655 
656  std::vector<std::string> complex_names = exio_helper->get_complex_names(names);
657 
658  std::vector<std::set<subdomain_id_type>> complex_vars_active_subdomains =
659  exio_helper->get_complex_vars_active_subdomains(vars_active_subdomains);
660  exio_helper->initialize_element_variables(complex_names, complex_vars_active_subdomains);
661 
662  unsigned int num_values = soln.size();
663  unsigned int num_vars = names.size();
664  unsigned int num_elems = num_values / num_vars;
665 
666  // This will contain the real and imaginary parts and the magnitude
667  // of the values in soln
668  std::vector<Real> complex_soln(3*num_values);
669 
670  for (unsigned i=0; i<num_vars; ++i)
671  {
672 
673  for (unsigned int j=0; j<num_elems; ++j)
674  {
675  Number value = soln[i*num_vars + j];
676  complex_soln[3*i*num_elems + j] = value.real();
677  }
678  for (unsigned int j=0; j<num_elems; ++j)
679  {
680  Number value = soln[i*num_vars + j];
681  complex_soln[3*i*num_elems + num_elems +j] = value.imag();
682  }
683  for (unsigned int j=0; j<num_elems; ++j)
684  {
685  Number value = soln[i*num_vars + j];
686  complex_soln[3*i*num_elems + 2*num_elems + j] = std::abs(value);
687  }
688  }
689 
690  exio_helper->write_element_values(mesh, complex_soln, _timestep, complex_vars_active_subdomains);
691 
692 #else
693  exio_helper->initialize_element_variables(names, vars_active_subdomains);
694  exio_helper->write_element_values(mesh, soln, _timestep, vars_active_subdomains);
695 #endif
696 }
697 
698 
699 
700 void ExodusII_IO::write_nodal_data (const std::string & fname,
701  const std::vector<Number> & soln,
702  const std::vector<std::string> & names)
703 {
704  LOG_SCOPE("write_nodal_data()", "ExodusII_IO");
705 
707 
708  int num_vars = cast_int<int>(names.size());
709  dof_id_type num_nodes = mesh.n_nodes();
710 
711  // The names of the variables to be output
712  std::vector<std::string> output_names;
713 
715  output_names = _output_variables;
716  else
717  output_names = names;
718 
719 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
720 
721  std::vector<std::string> complex_names = exio_helper->get_complex_names(names);
722 
723  // Call helper function for opening/initializing data, giving it the
724  // complex variable names
725  this->write_nodal_data_common(fname, complex_names, /*continuous=*/true);
726 #else
727  // Call helper function for opening/initializing data
728  this->write_nodal_data_common(fname, output_names, /*continuous=*/true);
729 #endif
730 
731  if (mesh.processor_id())
732  return;
733 
734  // This will count the number of variables actually output
735  for (int c=0; c<num_vars; c++)
736  {
737  std::stringstream name_to_find;
738 
739  std::vector<std::string>::iterator pos =
740  std::find(output_names.begin(), output_names.end(), names[c]);
741  if (pos == output_names.end())
742  continue;
743 
744  unsigned int variable_name_position =
745  cast_int<unsigned int>(pos - output_names.begin());
746 
747  // Set up temporary vectors to be passed to Exodus to write the
748  // nodal values for a single variable at a time.
749 #ifdef LIBMESH_USE_REAL_NUMBERS
750  std::vector<Number> cur_soln;
751 
752  // num_nodes is either exactly how much space we will need for
753  // each vector, or a safe upper bound for the amount of memory
754  // we will require when there are gaps in the numbering.
755  cur_soln.reserve(num_nodes);
756 #else
757  std::vector<Real> real_parts;
758  std::vector<Real> imag_parts;
759  std::vector<Real> magnitudes;
760  real_parts.reserve(num_nodes);
761  imag_parts.reserve(num_nodes);
762  magnitudes.reserve(num_nodes);
763 #endif
764 
765  // There could be gaps in "soln", but it will always be in the
766  // order of [num_vars * node_id + var_id]. We now copy the
767  // proper solution values contiguously into "cur_soln",
768  // removing the gaps.
769  for (const auto & node : mesh.node_ptr_range())
770  {
771  dof_id_type idx = node->id()*num_vars + c;
772 #ifdef LIBMESH_USE_REAL_NUMBERS
773  cur_soln.push_back(soln[idx]);
774 #else
775  real_parts.push_back(soln[idx].real());
776  imag_parts.push_back(soln[idx].imag());
777  magnitudes.push_back(std::abs(soln[idx]));
778 #endif
779  }
780 
781  // Finally, actually call the Exodus API to write to file.
782 #ifdef LIBMESH_USE_REAL_NUMBERS
783  exio_helper->write_nodal_values(variable_name_position+1, cur_soln, _timestep);
784 #else
785  exio_helper->write_nodal_values(3*variable_name_position+1, real_parts, _timestep);
786  exio_helper->write_nodal_values(3*variable_name_position+2, imag_parts, _timestep);
787  exio_helper->write_nodal_values(3*variable_name_position+3, magnitudes, _timestep);
788 #endif
789 
790  }
791 }
792 
793 
794 
795 
796 void ExodusII_IO::write_information_records (const std::vector<std::string> & records)
797 {
799  return;
800 
801  if (!exio_helper->opened_for_writing)
802  libmesh_error_msg("ERROR, ExodusII file must be initialized before outputting information records.");
803 
804  exio_helper->write_information_records(records);
805 }
806 
807 
808 
809 void ExodusII_IO::write_global_data (const std::vector<Number> & soln,
810  const std::vector<std::string> & names)
811 {
813  return;
814 
815  if (!exio_helper->opened_for_writing)
816  libmesh_error_msg("ERROR, ExodusII file must be initialized before outputting global variables.");
817 
818 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
819 
820  std::vector<std::string> complex_names = exio_helper->get_complex_names(names);
821 
822  exio_helper->initialize_global_variables(complex_names);
823 
824  unsigned int num_values = soln.size();
825  unsigned int num_vars = names.size();
826  unsigned int num_elems = num_values / num_vars;
827 
828  // This will contain the real and imaginary parts and the magnitude
829  // of the values in soln
830  std::vector<Real> complex_soln(3*num_values);
831 
832  for (unsigned i=0; i<num_vars; ++i)
833  {
834 
835  for (unsigned int j=0; j<num_elems; ++j)
836  {
837  Number value = soln[i*num_vars + j];
838  complex_soln[3*i*num_elems + j] = value.real();
839  }
840  for (unsigned int j=0; j<num_elems; ++j)
841  {
842  Number value = soln[i*num_vars + j];
843  complex_soln[3*i*num_elems + num_elems +j] = value.imag();
844  }
845  for (unsigned int j=0; j<num_elems; ++j)
846  {
847  Number value = soln[i*num_vars + j];
848  complex_soln[3*i*num_elems + 2*num_elems + j] = std::abs(value);
849  }
850  }
851 
852  exio_helper->write_global_values(complex_soln, _timestep);
853 
854 #else
855  exio_helper->initialize_global_variables(names);
856  exio_helper->write_global_values(soln, _timestep);
857 #endif
858 }
859 
860 
861 
862 void ExodusII_IO::write_timestep (const std::string & fname,
863  const EquationSystems & es,
864  const int timestep,
865  const Real time,
866  const std::set<std::string> * system_names)
867 {
868  _timestep = timestep;
869  write_equation_systems(fname,es,system_names);
870 
872  return;
873 
874  exio_helper->write_timestep(timestep, time);
875 }
876 
877 
878 
879 void ExodusII_IO::write (const std::string & fname)
880 {
882 
883  // We may need to gather a DistributedMesh to output it, making that
884  // const qualifier in our constructor a dirty lie
885  // The "true" specifies that we only need the mesh serialized to processor 0
887 
888  libmesh_assert( !exio_helper->opened_for_writing );
889 
890  // If the user has set the append flag here, it doesn't really make
891  // sense: the intent of this function is to write a Mesh with no
892  // data, while "appending" is really intended to add data to an
893  // existing file. If we're verbose, print a message to this effect.
894  if (_append && _verbose)
895  libmesh_warning("Warning: Appending in ExodusII_IO::write() does not make sense.\n"
896  "Creating a new file instead!");
897 
898  exio_helper->create(fname);
899  exio_helper->initialize(fname,mesh);
900  exio_helper->write_nodal_coordinates(mesh);
901  exio_helper->write_elements(mesh);
902  exio_helper->write_sidesets(mesh);
903  exio_helper->write_nodesets(mesh);
904 
905  if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
906  libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
907  "are not supported by the ExodusII format.");
908 }
909 
910 
911 
912 void ExodusII_IO::write_nodal_data_discontinuous (const std::string & fname,
913  const std::vector<Number> & soln,
914  const std::vector<std::string> & names)
915 {
916  LOG_SCOPE("write_nodal_data_discontinuous()", "ExodusII_IO");
917 
919 
920  int num_vars = cast_int<int>(names.size());
921  int num_nodes = 0;
922  for (const auto & elem : mesh.active_element_ptr_range())
923  num_nodes += elem->n_nodes();
924 
925 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
926 
927  std::vector<std::string> complex_names = exio_helper->get_complex_names(names);
928 
929  // Call helper function for opening/initializing data, giving it the
930  // complex variable names
931  this->write_nodal_data_common(fname, complex_names, /*continuous=*/false);
932 #else
933  // Call helper function for opening/initializing data
934  this->write_nodal_data_common(fname, names, /*continuous=*/false);
935 #endif
936 
937  if (mesh.processor_id())
938  return;
939 
940  for (int c=0; c<num_vars; c++)
941  {
942 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
943  std::vector<Real> real_parts(num_nodes);
944  std::vector<Real> imag_parts(num_nodes);
945  std::vector<Real> magnitudes(num_nodes);
946 
947  for (int i=0; i<num_nodes; ++i)
948  {
949  real_parts[i] = soln[i*num_vars + c].real();
950  imag_parts[i] = soln[i*num_vars + c].imag();
951  magnitudes[i] = std::abs(soln[i*num_vars + c]);
952  }
953  exio_helper->write_nodal_values(3*c+1,real_parts,_timestep);
954  exio_helper->write_nodal_values(3*c+2,imag_parts,_timestep);
955  exio_helper->write_nodal_values(3*c+3,magnitudes,_timestep);
956 #else
957  // Copy out this variable's solution
958  std::vector<Number> cur_soln(num_nodes);
959 
960  for (int i=0; i<num_nodes; i++)
961  cur_soln[i] = soln[i*num_vars + c];
962 
963  exio_helper->write_nodal_values(c+1,cur_soln,_timestep);
964 #endif
965  }
966 }
967 
968 
969 
971  const std::vector<std::string> & names,
972  bool continuous)
973 {
975 
976  // This function can be called multiple times, we only want to open
977  // the ExodusII file the first time it's called.
978  if (!exio_helper->opened_for_writing)
979  {
980  // If we're appending, open() the file with read_only=false,
981  // otherwise create() it and write the contents of the mesh to
982  // it.
983  if (_append)
984  {
985  exio_helper->open(fname.c_str(), /*read_only=*/false);
986  // If we're appending, it's not valid to call exio_helper->initialize()
987  // or exio_helper->initialize_nodal_variables(), but we do need to set up
988  // certain aspects of the Helper object itself, such as the number of nodes
989  // and elements. We do that by reading the header...
990  exio_helper->read_header();
991 
992  // ...and reading the block info
993  exio_helper->read_block_info();
994  }
995  else
996  {
997  exio_helper->create(fname);
998 
999  exio_helper->initialize(fname, mesh, !continuous);
1000  exio_helper->write_nodal_coordinates(mesh, !continuous);
1001  exio_helper->write_elements(mesh, !continuous);
1002 
1003  exio_helper->write_sidesets(mesh);
1004  exio_helper->write_nodesets(mesh);
1005 
1006  if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
1007  libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
1008  "are not supported by the ExodusII format.");
1009 
1010  exio_helper->initialize_nodal_variables(names);
1011  }
1012  }
1013  else
1014  {
1015  // We are already open for writing, so check that the filename
1016  // passed to this function matches the filename currently in use
1017  // by the helper.
1018  if (fname != exio_helper->current_filename)
1019  libmesh_error_msg("Error! This ExodusII_IO object is already associated with file: " \
1020  << exio_helper->current_filename \
1021  << ", cannot use it with requested file: " \
1022  << fname);
1023  }
1024 }
1025 
1026 const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
1027 {
1028  exio_helper->read_var_names(ExodusII_IO_Helper::NODAL);
1029  return exio_helper->nodal_var_names;
1030 }
1031 
1032 const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
1033 {
1034  exio_helper->read_var_names(ExodusII_IO_Helper::ELEMENTAL);
1035  return exio_helper->elem_var_names;
1036 }
1037 
1039 {
1040  // Provide a warning when accessing the helper object
1041  // since it is a non-public API and is likely to see
1042  // future API changes
1043  libmesh_experimental();
1044 
1045  return *exio_helper;
1046 }
1047 
1048 
1049 // LIBMESH_HAVE_EXODUS_API is not defined, declare error() versions of functions...
1050 #else
1051 
1052 
1053 
1055 {
1056 }
1057 
1058 
1059 
1060 void ExodusII_IO::read (const std::string &)
1061 {
1062  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1063 }
1064 
1065 
1066 
1067 void ExodusII_IO::verbose (bool)
1068 {
1069  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1070 }
1071 
1072 
1073 
1075 {
1076  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1077 }
1078 
1079 
1080 
1081 void ExodusII_IO::write_as_dimension(unsigned)
1082 {
1083  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1084 }
1085 
1086 
1087 
1089 {
1090  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1091 }
1092 
1093 
1094 
1095 void ExodusII_IO::append(bool)
1096 {
1097  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1098 }
1099 
1100 
1101 
1102 const std::vector<Real> & ExodusII_IO::get_time_steps()
1103 {
1104  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1105 }
1106 
1107 
1108 
1110 {
1111  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1112 }
1113 
1114 
1115 void ExodusII_IO::copy_nodal_solution(System &,
1116  std::string,
1117  std::string,
1118  unsigned int)
1119 {
1120  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1121 }
1122 
1123 
1124 
1126  std::string,
1127  std::string,
1128  unsigned int)
1129 {
1130  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1131 }
1132 
1133 
1134 
1135 void ExodusII_IO::write_element_data (const EquationSystems &)
1136 {
1137  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1138 }
1139 
1140 
1141 
1142 void ExodusII_IO::write_nodal_data (const std::string &,
1143  const std::vector<Number> &,
1144  const std::vector<std::string> &)
1145 {
1146  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1147 }
1148 
1149 
1150 
1151 void ExodusII_IO::write_information_records (const std::vector<std::string> &)
1152 {
1153  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1154 }
1155 
1156 
1157 
1158 void ExodusII_IO::write_global_data (const std::vector<Number> &,
1159  const std::vector<std::string> &)
1160 {
1161  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1162 }
1163 
1164 
1165 
1166 void ExodusII_IO::write_timestep (const std::string &,
1167  const EquationSystems &,
1168  const int,
1169  const Real,
1170  const std::set<std::string> *)
1171 {
1172  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1173 }
1174 
1175 
1176 
1177 void ExodusII_IO::write (const std::string &)
1178 {
1179  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1180 }
1181 
1182 
1183 
1184 void ExodusII_IO::write_nodal_data_discontinuous (const std::string &,
1185  const std::vector<Number> &,
1186  const std::vector<std::string> &)
1187 {
1188  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1189 }
1190 
1191 
1192 
1193 void ExodusII_IO::write_nodal_data_common(std::string,
1194  const std::vector<std::string> &,
1195  bool)
1196 {
1197  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1198 }
1199 
1200 
1201 const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
1202 {
1203  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1204 }
1205 
1206 const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
1207 {
1208  libmesh_error_msg("ERROR, ExodusII API is not defined.");
1209 }
1210 
1211 #endif // LIBMESH_HAVE_EXODUS_API
1212 } // namespace libMesh
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 MT & mesh() const
Definition: mesh_output.h:234
void use_mesh_dimension_instead_of_spatial_dimension(bool val)
Definition: exodusII_io.C:428
double abs(double a)
Manages multiples systems of equations.
virtual void reserve_nodes(const dof_id_type nn)=0
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:833
void build_variable_names(std::vector< std::string > &var_names, const FEType *type=nullptr, const std::set< std::string > *system_names=nullptr) const
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2024
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
std::string & nodeset_name(boundary_id_type id)
void write_as_dimension(unsigned dim)
Definition: exodusII_io.C:435
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:803
ExodusII_IO_Helper & get_exio_helper()
Definition: exodusII_io.C:1038
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
Definition: mesh_output.C:31
std::size_t n_edge_conds() const
const Elem * top_parent() const
Definition: elem.h:2503
std::vector< bool > elems_of_dimension
Definition: mesh_input.h:97
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
ExodusII_IO_Helper::Conversion assign_conversion(std::string type_str)
unique_id_type unique_id() const
Definition: dof_object.h:672
const Parallel::Communicator & comm() const
void copy_nodal_solution(System &system, std::string var_name, unsigned int timestep=1)
Definition: exodusII_io.C:78
IterBase * end
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
ExodusII_IO(MeshBase &mesh, bool single_precision=false)
Definition: exodusII_io.C:45
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:131
void append(bool val)
Definition: exodusII_io.C:450
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Base class for Mesh.
Definition: mesh_base.h:77
dof_id_type & set_id()
Definition: dof_object.h:664
void get_vars_active_subdomains(const std::vector< std::string > &names, std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) const
void write_timestep_discontinuous(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Definition: exodusII_io.C:104
void add_node(const Node *node, const boundary_id_type id)
unsigned int number() const
Definition: system.h:2025
int8_t boundary_id_type
Definition: id_types.h:51
virtual void write_discontinuous_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
Definition: mesh_output.C:92
dof_id_type id() const
Definition: dof_object.h:655
processor_id_type rank() const
Definition: communicator.h:173
void set_output_variables(const std::vector< std::string > &output_variables, bool allow_empty=true)
Definition: exodusII_io.C:68
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
virtual Elem * add_elem(Elem *e)=0
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:245
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1243
virtual SimpleRange< node_iterator > node_ptr_range()=0
void copy_elemental_solution(System &system, std::string system_var_name, std::string exodus_var_name, unsigned int timestep=1)
Definition: exodusII_io.C:513
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
const std::vector< Real > & get_time_steps()
Definition: exodusII_io.C:457
std::vector< std::string > _output_variables
Definition: exodusII_io.h:347
void write_information_records(const std::vector< std::string > &)
Definition: exodusII_io.C:796
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:538
void read_elemental_variable(std::string elemental_var_name, unsigned int timestep, std::map< unsigned int, Real > &unique_id_to_value_map)
Definition: exodusII_io.C:555
void set_mesh_dimension(unsigned char d)
Definition: mesh_base.h:213
void verbose(bool set_verbosity)
Definition: exodusII_io.C:418
void write_nodal_data_common(std::string fname, const std::vector< std::string > &names, bool continuous=true)
Definition: exodusII_io.C:970
An object whose state is distributed along a set of processors.
virtual void read(const std::string &name) override
Definition: exodusII_io.C:140
virtual void clear()
Definition: mesh_base.C:260
void read_global_variable(std::vector< std::string > global_var_names, unsigned int timestep, std::vector< Real > &global_values)
Definition: exodusII_io.C:570
std::string & sideset_name(boundary_id_type id)
void write_global_data(const std::vector< Number > &, const std::vector< std::string > &)
Definition: exodusII_io.C:809
void write_nodal_data_discontinuous(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
Definition: exodusII_io.C:912
std::string enum_to_string(const T e)
std::unique_ptr< ExodusII_IO_Helper > exio_helper
Definition: exodusII_io.h:323
virtual void update()
Definition: system.C:408
void write_discontinuous_exodusII(const std::string &name, const EquationSystems &es, const std::set< std::string > *system_names=nullptr)
Definition: exodusII_io.C:89
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2183
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
virtual unsigned short dim() const =0
Temporarily serializes a DistributedMesh for output.
virtual void write(const std::string &fname) override
Definition: exodusII_io.C:879
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Definition: exodusII_io.C:862
virtual ~ExodusII_IO()
Definition: exodusII_io.C:133
const std::vector< std::string > & get_elem_var_names()
Definition: exodusII_io.C:1032
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
static const bool value
Definition: xdr_io.C:109
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:504
unsigned int mesh_dimension() const
Definition: mesh_base.C:126
void write_element_data(const EquationSystems &es)
Definition: exodusII_io.C:599
virtual const Node * node_ptr(const dof_id_type i) const =0
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
Definition: exodusII_io.C:700
processor_id_type processor_id() const
void build_discontinuous_solution_vector(std::vector< Number > &soln, const std::set< std::string > *system_names=nullptr) const
const std::vector< std::string > & get_nodal_var_names()
Definition: exodusII_io.C:1026
A geometric point in (x,y,z) space.
Definition: point.h:38
virtual void reserve_elem(const dof_id_type ne)=0
virtual dof_id_type n_nodes() const =0
void set_coordinate_offset(Point p)
Definition: exodusII_io.C:442
void build_elemental_solution_vector(std::vector< Number > &soln, std::vector< std::string > &names) const
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
uint8_t dof_id_type
Definition: id_types.h:64