libMesh
exodusII_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 // Local includes
20 #include "libmesh/exodusII_io.h"
21 
22 #include "libmesh/boundary_info.h"
23 #include "libmesh/dof_map.h"
24 #include "libmesh/dyna_io.h" // ElementDefinition for BEX
25 #include "libmesh/enum_elem_type.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/enum_to_string.h"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/exodusII_io_helper.h"
30 #include "libmesh/fpe_disabler.h"
31 #include "libmesh/int_range.h"
32 #include "libmesh/libmesh_logging.h"
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/mesh_communication.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/parallel_mesh.h"
37 #include "libmesh/parallel.h"
38 #include "libmesh/system.h"
39 #include "libmesh/utility.h"
40 
41 // TIMPI includes
42 #include "timpi/parallel_sync.h"
43 
44 // C++ includes
45 #include <cmath> // llround
46 #include <cstring>
47 #include <fstream>
48 #include <map>
49 #include <memory>
50 #include <sstream>
51 
52 #ifdef LIBMESH_HAVE_EXODUS_API
53 namespace
54 {
55  using namespace libMesh;
56 
57  const std::vector<Real> & bex_constraint_vec(std::size_t i,
58  const ExodusII_IO_Helper & helper)
59  {
60  std::size_t vec_offset = 0;
61  for (auto block_num : index_range(helper.bex_dense_constraint_vecs))
62  {
63  const auto & vecblock = helper.bex_dense_constraint_vecs[block_num];
64  libmesh_assert_greater_equal(i, vec_offset);
65  if (i - vec_offset < vecblock.size())
66  return vecblock[i - vec_offset];
67 
68  vec_offset += vecblock.size();
69  }
70 
71  libmesh_error_msg("Requested BEX coefficient vector " << i << " not found");
72  }
73 
74 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
75  std::vector<Real>
76  complex_soln_components (const std::vector<Number> & soln,
77  const unsigned int num_vars,
78  const bool write_complex_abs)
79  {
80  unsigned int num_values = soln.size();
81  unsigned int num_elems = num_values / num_vars;
82 
83  // This will contain the real and imaginary parts and the magnitude
84  // of the values in soln
85  int nco = write_complex_abs ? 3 : 2;
86  std::vector<Real> complex_soln(nco * num_values);
87 
88  for (unsigned i=0; i<num_vars; ++i)
89  {
90  for (unsigned int j=0; j<num_elems; ++j)
91  {
92  Number value = soln[i*num_vars + j];
93  complex_soln[nco*i*num_elems + j] = value.real();
94  }
95  for (unsigned int j=0; j<num_elems; ++j)
96  {
97  Number value = soln[i*num_vars + j];
98  complex_soln[nco*i*num_elems + num_elems + j] = value.imag();
99  }
100  if (write_complex_abs)
101  {
102  for (unsigned int j=0; j<num_elems; ++j)
103  {
104  Number value = soln[i*num_vars + j];
105  complex_soln[3*i*num_elems + 2*num_elems + j] = std::abs(value);
106  }
107  }
108  }
109 
110  return complex_soln;
111  }
112 #endif // LIBMESH_USE_COMPLEX_NUMBERS
113 
114 }
115 #endif
116 
117 namespace libMesh
118 {
119 
120 // ------------------------------------------------------------
121 // ExodusII_IO class members
123  bool single_precision) :
126  /* is_parallel_format = */ false,
127  /* serial_only_needed_on_proc_0 = */ true),
129 #ifdef LIBMESH_HAVE_EXODUS_API
130  exio_helper(std::make_unique<ExodusII_IO_Helper>(*this, false, true, single_precision)),
131  _timestep(1),
132  _verbose(false),
133  _append(false),
134 #endif
135  _allow_empty_variables(false),
136  _write_complex_abs(true),
137  _disc_bex(false)
138 {
139  // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
140  libmesh_ignore(single_precision);
141 }
142 
143 
144 
146  bool single_precision) :
147  MeshInput<MeshBase> (),
149  /* is_parallel_format = */ false,
150  /* serial_only_needed_on_proc_0 = */ true),
152 #ifdef LIBMESH_HAVE_EXODUS_API
153  exio_helper(std::make_unique<ExodusII_IO_Helper>(*this, false, true, single_precision)),
154  _timestep(1),
155  _verbose(false),
156  _append(false),
157 #endif
158  _allow_empty_variables(false),
159  _write_complex_abs(true),
160  _disc_bex(false)
161 {
162  // if !LIBMESH_HAVE_EXODUS_API, we didn't use this
163  libmesh_ignore(single_precision);
164 }
165 
166 
167 
169 {
170 #ifdef LIBMESH_HAVE_EXODUS_API
172 #else
173  return 0;
174 #endif
175 }
176 
177 
178 
179 void ExodusII_IO::set_extra_integer_vars(const std::vector<std::string> & extra_integer_vars)
180 {
181  _extra_integer_vars = extra_integer_vars;
182 }
183 
184 void ExodusII_IO::set_output_variables(const std::vector<std::string> & output_variables,
185  bool allow_empty)
186 {
187  _output_variables = output_variables;
188  _allow_empty_variables = allow_empty;
189 }
190 
191 
192 
194  const EquationSystems & es,
195  const std::set<std::string> * system_names)
196 {
197  std::vector<std::string> solution_names;
198  std::vector<Number> v;
199 
200  es.build_variable_names (solution_names, nullptr, system_names);
201  es.build_discontinuous_solution_vector (v, system_names,
202  nullptr, false, /* defaults */
203  this->get_add_sides());
204  this->write_nodal_data_discontinuous(name, v, solution_names);
205 }
206 
207 
208 #ifdef LIBMESH_HAVE_EXODUS_API
209 void ExodusII_IO::write_timestep_discontinuous (const std::string &fname,
210  const EquationSystems &es,
211  const int timestep,
212  const Real time,
213  const std::set<std::string> * system_names)
214 {
215  _timestep = timestep;
216  write_discontinuous_equation_systems (fname,es,system_names);
217 
219  return;
220 
221  exio_helper->write_timestep(timestep, time);
222 }
223 
224 #else
225 void ExodusII_IO::write_timestep_discontinuous (const std::string & /* fname */,
226  const EquationSystems & /* es */,
227  const int /* timestep */,
228  const Real /* time */,
229  const std::set<std::string> * /*system_names*/)
230 { libmesh_error(); }
231 #endif
232 
233 
234 // ------------------------------------------------------------
235 // When the Exodus API is present...
236 #ifdef LIBMESH_HAVE_EXODUS_API
237 
239 {
240  exio_helper->close();
241 }
242 
243 
244 void ExodusII_IO::read (const std::string & fname)
245 {
246  LOG_SCOPE("read()", "ExodusII_IO");
247 
248  // Get a reference to the mesh we are reading
250 
251  // Add extra integers into the mesh
252  std::vector<unsigned int> extra_ids;
253  for (auto & name : _extra_integer_vars)
254  extra_ids.push_back(mesh.add_elem_integer(name));
255 
256  // Clear any existing mesh data
257  mesh.clear();
258 
259  // Keep track of what kinds of elements this file contains
260  elems_of_dimension.clear();
261  elems_of_dimension.resize(4, false);
262 
263  // Open the exodus file in EX_READ mode
264  exio_helper->open(fname.c_str(), /*read_only=*/true);
265 
266  // Get header information from exodus file
267  exio_helper->read_and_store_header_info();
268 
269  // Read the QA records
270  exio_helper->read_qa_records();
271 
272  // Print header information
273  exio_helper->print_header();
274 
275  // Read nodes from the exodus file
276  exio_helper->read_nodes();
277 
278  // Reserve space for the nodes.
279  mesh.reserve_nodes(exio_helper->num_nodes);
280 
281  // Read the node number map from the Exodus file. This is
282  // required if we want to preserve the numbering of nodes as it
283  // exists in the Exodus file. If the Exodus file does not contain
284  // a node_num_map, the identity map is returned by this call.
285  exio_helper->read_node_num_map();
286 
287  // Read the element number map from the Exodus file. This is
288  // required if we want to preserve the numbering of elements as it
289  // exists in the Exodus file. If the Exodus file does not contain
290  // an elem_num_map, the identity map is returned by this call.
291  //
292  // We now do this before creating nodes, so if we have any spline
293  // nodes that need a NodeElem attached we can give them unused elem
294  // ids.
295  exio_helper->read_elem_num_map();
296 
297  // Read any Bezier Extraction coefficient vectors from the file,
298  // such as might occur in an IsoGeometric Analysis (IGA) mesh.
299  exio_helper->read_bex_cv_blocks();
300 
301  // If we have Rational Bezier weights, we'll need to
302  // store them.
303  unsigned char weight_index = 0;
304  const bool weights_exist = !exio_helper->w.empty();
305 
306  // If we have Bezier extraction coefficients, we'll need to put
307  // NODEELEM elements on spline nodes, since our Rational Bezier
308  // elements will be connected to nodes derived from those; nothing
309  // else will be directly connected to the spline nodes.
310  const bool bex_cv_exist = !exio_helper->bex_dense_constraint_vecs.empty();
311 
312  // Even if weights don't exist, we still use RATIONAL_BERNSTEIN for
313  // Bezier-Bernstein BEX elements, we just use 1.0 as the weight on
314  // every node.
315  if (bex_cv_exist)
316  {
317  const Real default_weight = 1.0;
318  weight_index = cast_int<unsigned char>
319  (mesh.add_node_datum<Real>("rational_weight", true,
320  &default_weight));
322  mesh.set_default_mapping_data(weight_index);
323  }
324 
325  std::unordered_map<const Node *, Elem *> spline_nodeelem_ptrs;
326 
327  // Loop over the nodes, create Nodes with local processor_id 0.
328  for (int i=0; i<exio_helper->num_nodes; i++)
329  {
330  // Use the node_num_map to get the correct ID for Exodus
331  int exodus_id = exio_helper->node_num_map[i];
332 
333  // Catch the node that was added to the mesh
334  Node * added_node = mesh.add_point (Point(exio_helper->x[i], exio_helper->y[i], exio_helper->z[i]), exodus_id-1);
335 
336  // If the Mesh assigned an ID different from what is in the
337  // Exodus file, we should probably error.
338  libmesh_error_msg_if(added_node->id() != static_cast<unsigned>(exodus_id-1),
339  "Error! Mesh assigned node ID "
340  << added_node->id()
341  << " which is different from the (zero-based) Exodus ID "
342  << exodus_id-1
343  << "!");
344 
345  // If we have a set of spline weights, these nodes are going to
346  // be used as control points for Bezier elements, and we need
347  // to attach a NodeElem to each to make sure it doesn't get
348  // flagged as an unused node.
349  if (weights_exist)
350  {
351  const auto w = exio_helper->w[i];
352  Point & p = *added_node;
353  p /= w; // Exodus Bezier Extraction stores spline nodes in projective space
354 
355  added_node->set_extra_datum<Real>(weight_index, exio_helper->w[i]);
356  }
357 
358  if (bex_cv_exist)
359  {
360  std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
361 
362  // Give the NodeElem ids at the end, so we can match any
363  // existing ids in the file for other elements
364  elem->set_id() = exio_helper->end_elem_id() + i;
365 
366  elem->set_node(0, added_node);
367  Elem * added_elem = mesh.add_elem(std::move(elem));
368  spline_nodeelem_ptrs[added_node] = added_elem;
369  }
370  }
371 
372  // This assert is no longer valid if the nodes are not numbered
373  // sequentially starting from 1 in the Exodus file.
374  // libmesh_assert_equal_to (static_cast<unsigned int>(exio_helper->num_nodes), mesh.n_nodes());
375 
376  // Get information about all the element and edge blocks
377  exio_helper->read_block_info();
378 
379  // Reserve space for the elements. Account for any NodeElem that
380  // have already been attached to spline control nodes.
381  mesh.reserve_elem(exio_helper->w.size() + exio_helper->num_elem);
382 
383  // Read variables for extra integer IDs
384  std::vector<std::map<dof_id_type, Real>> elem_ids(extra_ids.size());
385  // We use the last time step to load the IDs
386  exio_helper->read_num_time_steps();
387  unsigned int last_step = exio_helper->num_time_steps;
388  for (auto i : index_range(extra_ids))
389  exio_helper->read_elemental_var_values(_extra_integer_vars[i], last_step, elem_ids[i]);
390 
391  // Read in the element connectivity for each block.
392  int nelem_last_block = 0;
393 
394  // If we're building Bezier elements from spline nodes, we need to
395  // calculate those elements' local nodes on the fly, and we'll be
396  // calculating them from constraint matrix columns, and we'll need
397  // to make sure that the same node is found each time it's
398  // calculated from multiple neighboring elements.
399  std::map<std::vector<std::pair<dof_id_type, Real>>, Node *> local_nodes;
400 
401  // We'll set any spline NodeElem subdomain_id() values to exceed the
402  // maximum of subdomain_id() values set via Exodus block ids.
403  int max_subdomain_id = std::numeric_limits<int>::min();
404 
405  // We've already added all the nodes explicitly specified in the
406  // file, but if we have spline nodes we may need to add assembly
407  // element nodes based on them. Add them contiguously so we're
408  // compatible with any subsequent code paths (including our
409  // ExodusII_IO::write()!) that don't support sparse ids.
411 
412  // Loop over all the element blocks
413  for (int i=0; i<exio_helper->num_elem_blk; i++)
414  {
415  // Read the information for block i
416  exio_helper->read_elem_in_block (i);
417  const int subdomain_id = exio_helper->get_block_id(i);
418  max_subdomain_id = std::max(max_subdomain_id, subdomain_id);
419 
420  // populate the map of names
421  std::string subdomain_name = exio_helper->get_block_name(i);
422  if (!subdomain_name.empty())
423  mesh.subdomain_name(static_cast<subdomain_id_type>(subdomain_id)) = subdomain_name;
424 
425  // Set any relevant node/edge maps for this element
426  const std::string type_str (exio_helper->get_elem_type());
427  const auto & conv = exio_helper->get_conversion(type_str);
428 
429  // Loop over all the faces in this block
430  int jmax = nelem_last_block+exio_helper->num_elem_this_blk;
431  for (int j=nelem_last_block; j<jmax; j++)
432  {
433  auto uelem = Elem::build(conv.libmesh_elem_type());
434 
435  const int elem_num = j - nelem_last_block;
436 
437  // Make sure that Exodus's number of nodes per Elem matches
438  // the number of Nodes for this type of Elem. We only check
439  // this for the first Elem in each block, since these values
440  // are the same for every Elem in the block.
441  if (!elem_num)
442  libmesh_error_msg_if(exio_helper->num_nodes_per_elem != static_cast<int>(uelem->n_nodes()),
443  "Error: Exodus file says "
444  << exio_helper->num_nodes_per_elem
445  << " nodes per Elem, but Elem type "
446  << Utility::enum_to_string(uelem->type())
447  << " has " << uelem->n_nodes() << " nodes.");
448 
449  // Assign the current subdomain to this Elem
450  uelem->subdomain_id() = static_cast<subdomain_id_type>(subdomain_id);
451 
452  // Use the elem_num_map to obtain the ID of this element in
453  // the Exodus file. Make sure we aren't reading garbage if
454  // the file is corrupt.
455  libmesh_error_msg_if(std::size_t(j) >= exio_helper->elem_num_map.size(),
456  "Error: Trying to read Exodus file with more elements than elem_num_map entries.\n");
457  int exodus_id = exio_helper->elem_num_map[j];
458 
459  // Assign this element the same ID it had in the Exodus
460  // file, but make it zero-based by subtracting 1. Note:
461  // some day we could use 1-based numbering in libmesh and
462  // thus match the Exodus numbering exactly, but at the
463  // moment libmesh is zero-based.
464  uelem->set_id(exodus_id-1);
465 
466  // Record that we have seen an element of dimension uelem->dim()
467  elems_of_dimension[uelem->dim()] = true;
468 
469  // Catch the Elem pointer that the Mesh throws back
470  Elem * elem = mesh.add_elem(std::move(uelem));
471 
472  // If the Mesh assigned an ID different from what is in the
473  // Exodus file, we should probably error.
474  libmesh_error_msg_if(elem->id() != static_cast<unsigned>(exodus_id-1),
475  "Error! Mesh assigned ID "
476  << elem->id()
477  << " which is different from the (zero-based) Exodus ID "
478  << exodus_id-1
479  << "!");
480 
481  // Assign extra integer IDs
482  for (auto & id : extra_ids)
483  {
484  const Real v = elem_ids[id][elem->id()];
485 
486  if (v == Real(-1))
487  {
489  continue;
490  }
491 
492  // Ignore FE_INVALID here even if we've enabled FPEs; a
493  // thrown exception is preferred over an FPE signal.
494  FPEDisabler disable_fpes;
495  const long long iv = std::llround(v);
496 
497  // Check if the real number is outside of the range we can
498  // convert exactly
499 
500  long long max_representation = 1;
501  max_representation = (max_representation << std::min(std::numeric_limits<Real>::digits,
502  std::numeric_limits<double>::digits));
503  libmesh_error_msg_if(iv > max_representation,
504  "Error! An element integer value higher than "
505  << max_representation
506  << " was found! Exodus uses real numbers for storing element "
507  " integers, which can only represent integers from 0 to "
508  << max_representation
509  << ".");
510 
511  libmesh_error_msg_if(iv < 0,
512  "Error! An element integer value less than -1"
513  << " was found! Exodus uses real numbers for storing element "
514  " integers, which can only represent integers from 0 to "
515  << max_representation
516  << ".");
517 
518 
519  elem->set_extra_integer(id, cast_int<dof_id_type>(iv));
520  }
521 
522  // Set all the nodes for this element
523  //
524  // If we don't have any Bezier extraction operators, this
525  // is easy: we've already built all our nodes and just need
526  // to link to them.
527  if (exio_helper->bex_cv_conn.empty())
528  {
529  for (int k=0; k<exio_helper->num_nodes_per_elem; k++)
530  {
531  // global index
532  int gi = (elem_num)*exio_helper->num_nodes_per_elem + conv.get_node_map(k);
533 
534  // The entries in 'connect' are actually (1-based)
535  // indices into the node_num_map, so to get the right
536  // node ID we:
537  // 1.) Subtract 1 from connect[gi]
538  // 2.) Pass it through node_num_map to get the corresponding Exodus ID
539  // 3.) Subtract 1 from that, since libmesh node numbering is "zero"-based,
540  // even when the Exodus node numbering doesn't start with 1.
541  int libmesh_node_id = exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1;
542 
543  // Set the node pointer in the Elem
544  elem->set_node(k, mesh.node_ptr(libmesh_node_id));
545  }
546  }
547  else // We have Bezier Extraction data
548  {
549  auto & constraint_rows = mesh.get_constraint_rows();
550 
551  const DynaIO::ElementDefinition & dyna_elem_defn =
553  elem->dim(),
554  int(elem->default_order()));
555 
556  std::vector<std::vector<Real>>
557  my_constraint_mat(exio_helper->bex_num_elem_cvs);
558  for (auto spline_node_index :
559  make_range(exio_helper->bex_num_elem_cvs))
560  {
561  my_constraint_mat[spline_node_index].resize(elem->n_nodes());
562 
563  const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];
564  const unsigned long elem_coef_vec_index =
565  my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
566  const auto & my_vec = bex_constraint_vec(elem_coef_vec_index, *exio_helper);
567  for (auto elem_node_index :
568  make_range(elem->n_nodes()))
569  {
570  my_constraint_mat[spline_node_index][elem_node_index] =
571  my_vec[elem_node_index];
572  }
573 
574  }
575 
576  // The tailing entries in each element's connectivity
577  // vector are indices to Exodus constraint coefficient
578  // rows.
579 
580  // Concatenating these rows gives a matrix with
581  // all the constraints for the element nodes: each
582  // column of that matrix is the constraint coefficients
583  // for the node associated with that column (via the
584  // Exodus numbering, not the libMesh numbering).
585  const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];
586 
587  for (auto elem_node_index :
588  make_range(elem->n_nodes()))
589  {
590  // New finite element node data = dot product of
591  // constraint matrix columns with spline node data.
592  // Store each column's non-zero entries, along with
593  // the global spline node indices, as a key to
594  // identify shared finite element nodes.
595  std::vector<std::pair<dof_id_type, Real>> key;
596 
597  for (auto spline_node_index :
598  make_range(exio_helper->bex_num_elem_cvs))
599  {
600  // Pick out a row of the element constraint matrix
601  const unsigned long elem_coef_vec_index =
602  my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based
603 
604  auto & coef_vec =
605  bex_constraint_vec(elem_coef_vec_index, *exio_helper);
606 
607  // Get coef from this node's column intersect that row
608  const Real coef =
609  libmesh_vector_at(coef_vec, elem_node_index);
610 
611  // Get the libMesh node corresponding to that row
612  const int gi = (elem_num)*exio_helper->bex_num_elem_cvs +
613  spline_node_index;
614  const dof_id_type libmesh_node_id =
615  exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1;
616 
617  if (coef != 0) // Ignore irrelevant spline nodes
618  key.emplace_back(libmesh_node_id, coef);
619  }
620 
621  // Have we already created this node? Connect it.
622  if (const auto local_node_it = local_nodes.find(key);
623  local_node_it != local_nodes.end())
624  elem->set_node(dyna_elem_defn.nodes[elem_node_index], local_node_it->second);
625  // Have we not yet created this node? Construct it,
626  // along with its weight and libMesh constraint row,
627  // then connect it.
628  else
629  {
630  Point p(0);
631  Real w = 0;
632  std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;
633 
634  for (auto [libmesh_spline_node_id, coef] : key)
635  {
636  const Node & spline_node = mesh.node_ref(libmesh_spline_node_id);
637 
638  p.add_scaled(spline_node, coef);
639  const Real spline_w = weights_exist ?
640  spline_node.get_extra_datum<Real>(weight_index) : 1;
641  w += coef * spline_w;
642 
643  const Elem * nodeelem =
644  libmesh_map_find(spline_nodeelem_ptrs, &spline_node);
645  constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
646  }
647 
648  Node *n = mesh.add_point(p, n_nodes++);
649  if (weights_exist)
650  n->set_extra_datum<Real>(weight_index, w);
651 
652  // If we're building disconnected Bezier
653  // extraction elements then we don't want to
654  // find the new nodes to reuse later; each
655  // finite element node will connect to only one
656  // element.
657  if (!_disc_bex)
658  local_nodes[key] = n;
659  elem->set_node(dyna_elem_defn.nodes[elem_node_index], n);
660 
661  constraint_rows[n] = constraint_row;
662  }
663  }
664  }
665  }
666 
667  // running sum of # of elements per block,
668  // (should equal total number of elements in the end)
669  nelem_last_block += exio_helper->num_elem_this_blk;
670  }
671 
672  // Now we know enough to fix any spline NodeElem subdomains
673  max_subdomain_id++;
674  for (auto p : spline_nodeelem_ptrs)
675  p.second->subdomain_id() = max_subdomain_id;
676 
677  // Read in edge blocks, storing information in the BoundaryInfo object.
678  // Edge blocks are treated as BCs.
679  exio_helper->read_edge_blocks(mesh);
680 
681  // Set the mesh dimension to the largest encountered for an element
682  for (unsigned char i=0; i!=4; ++i)
683  if (elems_of_dimension[i])
685 
686  // Read in sideset information -- this is useful for applying boundary conditions
687  {
688  // Get basic information about all sidesets
689  exio_helper->read_sideset_info();
690  int offset=0;
691  for (int i=0; i<exio_helper->num_side_sets; i++)
692  {
693  // Compute new offset
694  offset += (i > 0 ? exio_helper->num_sides_per_set[i-1] : 0);
695  exio_helper->read_sideset (i, offset);
696 
697  std::string sideset_name = exio_helper->get_side_set_name(i);
698  if (!sideset_name.empty())
700  (cast_int<boundary_id_type>(exio_helper->get_side_set_id(i)))
701  = sideset_name;
702  }
703 
704  for (auto e : index_range(exio_helper->elem_list))
705  {
706  // The numbers in the Exodus file sidesets should be thought
707  // of as (1-based) indices into the elem_num_map array. So,
708  // to get the right element ID we have to:
709  // 1.) Subtract 1 from elem_list[e] (to get a zero-based index)
710  // 2.) Pass it through elem_num_map (to get the corresponding Exodus ID)
711  // 3.) Subtract 1 from that, since libmesh is "zero"-based,
712  // even when the Exodus numbering doesn't start with 1.
713  dof_id_type libmesh_elem_id =
714  cast_int<dof_id_type>(exio_helper->elem_num_map[exio_helper->elem_list[e] - 1] - 1);
715 
716  // Set any relevant node/edge maps for this element
717  Elem & elem = mesh.elem_ref(libmesh_elem_id);
718 
719  const auto & conv = exio_helper->get_conversion(elem.type());
720 
721  // Map the zero-based Exodus side numbering to the libmesh side numbering
722  unsigned int raw_side_index = exio_helper->side_list[e]-1;
723  std::size_t side_index_offset = conv.get_shellface_index_offset();
724 
725  if (raw_side_index < side_index_offset)
726  {
727  // We assume this is a "shell face"
728  int mapped_shellface = raw_side_index;
729 
730  // Check for errors
731  libmesh_error_msg_if(mapped_shellface < 0 || mapped_shellface >= 2,
732  "Bad 0-based shellface id: "
733  << mapped_shellface
734  << " detected in Exodus file "
735  << exio_helper->current_filename);
736 
737  // Add this (elem,shellface,id) triplet to the BoundaryInfo object.
738  mesh.get_boundary_info().add_shellface (libmesh_elem_id,
739  cast_int<unsigned short>(mapped_shellface),
740  cast_int<boundary_id_type>(exio_helper->id_list[e]));
741  }
742  else
743  {
744  unsigned int side_index = static_cast<unsigned int>(raw_side_index - side_index_offset);
745  int mapped_side = conv.get_side_map(side_index);
746 
747  // Check for errors
748  libmesh_error_msg_if(mapped_side == ExodusII_IO_Helper::Conversion::invalid_id,
749  "Invalid 1-based side id: "
750  << side_index
751  << " detected for "
752  << Utility::enum_to_string(elem.type())
753  << " in Exodus file "
754  << exio_helper->current_filename);
755 
756  libmesh_error_msg_if(mapped_side < 0 ||
757  cast_int<unsigned int>(mapped_side) >= elem.n_sides(),
758  "Bad 0-based side id: "
759  << mapped_side
760  << " detected for "
761  << Utility::enum_to_string(elem.type())
762  << " in Exodus file "
763  << exio_helper->current_filename);
764 
765  // Add this (elem,side,id) triplet to the BoundaryInfo object.
766  mesh.get_boundary_info().add_side (libmesh_elem_id,
767  cast_int<unsigned short>(mapped_side),
768  cast_int<boundary_id_type>(exio_helper->id_list[e]));
769  }
770  } // end for (elem_list)
771  } // end read sideset info
772 
773  // Read in elemset information and apply to Mesh elements if present
774  {
775  exio_helper->read_elemset_info();
776 
777  // Mimic behavior of sideset case where we store all the set
778  // information in a single array with offsets.
779  int offset=0;
780  for (int i=0; i<exio_helper->num_elem_sets; i++)
781  {
782  // Compute new offset
783  offset += (i > 0 ? exio_helper->num_elems_per_set[i-1] : 0);
784  exio_helper->read_elemset (i, offset);
785 
786  // TODO: add support for elemset names
787  // std::string elemset_name = exio_helper->get_elem_set_name(i);
788  // if (!elemset_name.empty())
789  // mesh.get_boundary_info().elemset_name(cast_int<boundary_id_type>(exio_helper->get_elem_set_id(i))) = elemset_name;
790  }
791 
792  // Debugging: print the concatenated list of elemset ids
793  // libMesh::out << "Concatenated list of elemset Elem ids (Exodus numbering):" << std::endl;
794  // for (const auto & id : exio_helper->elemset_list)
795  // libMesh::out << id << " ";
796  // libMesh::out << std::endl;
797 
798  // Next we need to assign the elemset ids to the mesh using the
799  // Elem's "extra_integers" support, if we have any.
800  if (exio_helper->num_elem_all_elemsets)
801  {
802  // Build map from Elem -> {elemsets}. This is needed only
803  // temporarily to determine a unique set of elemset codes.
804  std::map<Elem *, MeshBase::elemset_type> elem_to_elemsets;
805  for (auto e : index_range(exio_helper->elemset_list))
806  {
807  // Follow standard (see sideset case above) approach for
808  // converting the ids stored in the elemset_list to
809  // libmesh Elem ids.
810  //
811  // TODO: this should be moved to a helper function so we
812  // don't duplicate the code.
813  dof_id_type libmesh_elem_id =
814  cast_int<dof_id_type>(exio_helper->elem_num_map[exio_helper->elemset_list[e] - 1] - 1);
815 
816  // Get a pointer to this Elem
817  Elem * elem = mesh.elem_ptr(libmesh_elem_id);
818 
819  // Debugging:
820  // libMesh::out << "Elem " << elem->id() << " is in elemset " << exio_helper->elemset_id_list[e] << std::endl;
821 
822  // Store elemset id in the map
823  elem_to_elemsets[elem].insert(exio_helper->elemset_id_list[e]);
824  }
825 
826  // Create a set of unique elemsets
827  std::set<MeshBase::elemset_type> unique_elemsets;
828  for (const auto & pr : elem_to_elemsets)
829  unique_elemsets.insert(pr.second);
830 
831  // Debugging: print the unique elemsets
832  // libMesh::out << "The set of unique elemsets which exist on the Mesh:" << std::endl;
833  // for (const auto & s : unique_elemsets)
834  // {
835  // for (const auto & elemset_id : s)
836  // libMesh::out << elemset_id << " ";
837  // libMesh::out << std::endl;
838  // }
839 
840  // Enumerate the unique_elemsets and tell the mesh about them
841  dof_id_type code = 0;
842  for (const auto & s : unique_elemsets)
843  mesh.add_elemset_code(code++, s);
844 
845  // Sanity check: make sure that MeshBase::n_elemsets() reports
846  // the expected value after calling MeshBase::add_elemset_code()
847  // one or more times.
848  libmesh_assert_msg(exio_helper->num_elem_sets == cast_int<int>(mesh.n_elemsets()),
849  "Error: mesh.n_elemsets() is " << mesh.n_elemsets()
850  << ", but mesh should have " << exio_helper->num_elem_sets << " elemsets.");
851 
852  // Create storage for the extra integer on all Elems. Elems which
853  // are not in any set will use the default value of DofObject::invalid_id
854  unsigned int elemset_index =
855  mesh.add_elem_integer("elemset_code",
856  /*allocate_data=*/true);
857 
858  // Store the appropriate extra_integer value on all Elems that need it.
859  for (const auto & [elem, s] : elem_to_elemsets)
860  elem->set_extra_integer(elemset_index, mesh.get_elemset_code(s));
861  }
862  } // done reading elemset info
863 
864  // Read nodeset info
865  {
866  // This fills in the following fields of the helper for later use:
867  // nodeset_ids
868  // num_nodes_per_set
869  // num_node_df_per_set
870  // node_sets_node_index
871  // node_sets_dist_index
872  // node_sets_node_list
873  // node_sets_dist_fact
874  exio_helper->read_all_nodesets();
875 
876  for (int nodeset=0; nodeset<exio_helper->num_node_sets; nodeset++)
877  {
878  boundary_id_type nodeset_id =
879  cast_int<boundary_id_type>(exio_helper->nodeset_ids[nodeset]);
880 
881  std::string nodeset_name = exio_helper->get_node_set_name(nodeset);
882  if (!nodeset_name.empty())
883  mesh.get_boundary_info().nodeset_name(nodeset_id) = nodeset_name;
884 
885  // Get starting index of node ids for current nodeset.
886  unsigned int offset = exio_helper->node_sets_node_index[nodeset];
887 
888  for (int i=0; i<exio_helper->num_nodes_per_set[nodeset]; ++i)
889  {
890  int exodus_id = exio_helper->node_sets_node_list[i + offset];
891 
892  // It's possible for nodesets to have invalid ids in them
893  // by accident. Instead of possibly accessing past the
894  // end of node_num_map, let's make sure we have that many
895  // entries.
896  libmesh_error_msg_if(static_cast<std::size_t>(exodus_id - 1) >= exio_helper->node_num_map.size(),
897  "Invalid Exodus node id " << exodus_id
898  << " found in nodeset " << nodeset_id);
899 
900  // As before, the entries in 'node_list' are 1-based
901  // indices into the node_num_map array, so we have to map
902  // them. See comment above.
903  int libmesh_node_id = exio_helper->node_num_map[exodus_id - 1] - 1;
904  mesh.get_boundary_info().add_node(cast_int<dof_id_type>(libmesh_node_id),
905  nodeset_id);
906  }
907  }
908  }
909 
910 #if LIBMESH_DIM < 3
911  libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
912  "Cannot open dimension "
913  << mesh.mesh_dimension()
914  << " mesh file when configured without "
915  << mesh.mesh_dimension()
916  << "D support.");
917 #endif
918 }
919 
920 
921 
923 ExodusII_IO::read_header (const std::string & fname)
924 {
925  // We will need the Communicator of the Mesh we were created with.
927 
928  // Eventual return value
929  ExodusHeaderInfo header_info;
930 
931  // File I/O is done on processor 0, then broadcast to other procs
932  if (mesh.processor_id() == 0)
933  {
934  // Open the exodus file in EX_READ mode
935  exio_helper->open(fname.c_str(), /*read_only=*/true);
936 
937  // Get header information from exodus file without updating the
938  // Helper object's internal data structures.
939  header_info = exio_helper->read_header();
940 
941  // Close the file, we are now done with it. The goal is to keep the
942  // exio_helper object unchanged while calling this function,
943  // although it can't quite be marked "const" because we do have to
944  // actually open/close the file. This way, it should be possible to
945  // use the same ExodusII_IO object to read the headers of multiple
946  // different mesh files.
947  exio_helper->close();
948  }
949 
950  // Broadcast header_info to other procs before returning
951  header_info.broadcast(mesh.comm());
952 
953  // Return the information we read back to the user.
954  return header_info;
955 }
956 
957 
958 
959 void ExodusII_IO::verbose (bool set_verbosity)
960 {
961  _verbose = set_verbosity;
962 
963  // Set the verbose flag in the helper object as well.
964  exio_helper->verbose = _verbose;
965 }
966 
967 
968 
970 {
971  _write_complex_abs = val;
972 }
973 
974 
975 
977 {
978  exio_helper->use_mesh_dimension_instead_of_spatial_dimension(val);
979 }
980 
981 
982 
984 {
985  exio_helper->write_as_dimension(dim);
986 }
987 
988 
989 
991 {
992  libmesh_warning("This method may be deprecated in the future");
993  exio_helper->set_coordinate_offset(p);
994 }
995 
996 
997 
998 void ExodusII_IO::append(bool val)
999 {
1000  _append = val;
1001 }
1002 
1003 
1004 
1006 {
1007  exio_helper->set_add_sides(val);
1008 }
1009 
1010 
1011 
1013 {
1014  return exio_helper->get_add_sides();
1015 }
1016 
1017 
1018 
1019 
1020 const std::vector<Real> & ExodusII_IO::get_time_steps()
1021 {
1022  libmesh_error_msg_if
1023  (!exio_helper->opened_for_reading,
1024  "ERROR, ExodusII file must be opened for reading before calling ExodusII_IO::get_time_steps()!");
1025 
1026  exio_helper->read_time_steps();
1027  return exio_helper->time_steps;
1028 }
1029 
1030 
1031 
1033 {
1034  libmesh_error_msg_if(!exio_helper->opened_for_reading && !exio_helper->opened_for_writing,
1035  "ERROR, ExodusII file must be opened for reading or writing before calling ExodusII_IO::get_num_time_steps()!");
1036 
1037  exio_helper->read_num_time_steps();
1038  return exio_helper->num_time_steps;
1039 }
1040 
1041 
1042 
1044  std::string system_var_name,
1045  std::string exodus_var_name,
1046  unsigned int timestep)
1047 {
1048  LOG_SCOPE("copy_nodal_solution()", "ExodusII_IO");
1049 
1050  const unsigned int var_num = system.variable_number(system_var_name);
1051 
1053 
1054  // With Exodus files we only open them on processor 0, so that's the
1055  // where we have to do the data read too.
1056  if (system.comm().rank() == 0)
1057  {
1058  libmesh_error_msg_if(!exio_helper->opened_for_reading,
1059  "ERROR, ExodusII file must be opened for reading before copying a nodal solution!");
1060 
1061  exio_helper->read_nodal_var_values(exodus_var_name, timestep);
1062  }
1063 
1064  auto & node_var_value_map = exio_helper->nodal_var_values;
1065 
1066  const bool serial_on_zero = mesh.is_serial_on_zero();
1067 
1068  // If our mesh isn't serial, then non-root processors need to
1069  // request the data for their parts of the mesh and insert it
1070  // themselves.
1071  if (!serial_on_zero)
1072  {
1073  std::unordered_map<processor_id_type, std::vector<dof_id_type>> node_ids_to_request;
1074  if (this->processor_id() != 0)
1075  {
1076  std::vector<dof_id_type> node_ids;
1077  for (auto & node : mesh.local_node_ptr_range())
1078  node_ids.push_back(node->id());
1079  if (!node_ids.empty())
1080  node_ids_to_request[0] = std::move(node_ids);
1081  }
1082 
1083  auto value_gather_functor =
1084  [& node_var_value_map]
1086  const std::vector<dof_id_type> & ids,
1087  std::vector<Real> & values)
1088  {
1089  const std::size_t query_size = ids.size();
1090  values.resize(query_size);
1091  for (std::size_t i=0; i != query_size; ++i)
1092  {
1093  if (const auto it = node_var_value_map.find(ids[i]);
1094  it != node_var_value_map.end())
1095  {
1096  values[i] = it->second;
1097  node_var_value_map.erase(it);
1098  }
1099  else
1100  values[i] = std::numeric_limits<Real>::quiet_NaN();
1101  }
1102  };
1103 
1104  auto value_action_functor =
1105  [& node_var_value_map]
1107  const std::vector<dof_id_type> & ids,
1108  const std::vector<Real> & values)
1109  {
1110  const std::size_t query_size = ids.size();
1111  for (std::size_t i=0; i != query_size; ++i)
1112  if (!libmesh_isnan(values[i]))
1113  node_var_value_map[ids[i]] = values[i];
1114  };
1115 
1116  Real * value_ex = nullptr;
1117  Parallel::pull_parallel_vector_data
1118  (system.comm(), node_ids_to_request, value_gather_functor,
1119  value_action_functor, value_ex);
1120  }
1121 
1122  // Everybody inserts the data they've received. If we're
1123  // serial_on_zero then proc 0 inserts everybody's data and other
1124  // procs have empty map ranges.
1125  for (auto p : exio_helper->nodal_var_values)
1126  {
1127  dof_id_type i = p.first;
1128  const Node * node = MeshInput<MeshBase>::mesh().query_node_ptr(i);
1129 
1130  if (node &&
1131  (serial_on_zero || node->processor_id() == system.processor_id()) &&
1132  node->n_comp(system.number(), var_num) > 0)
1133  {
1134  dof_id_type dof_index = node->dof_number(system.number(), var_num, 0);
1135 
1136  // If the dof_index is local to this processor, set the value
1137  system.solution->set (dof_index, p.second);
1138  }
1139  }
1140 
1141  system.solution->close();
1142  system.update();
1143 }
1144 
1145 
1146 
1148  std::string system_var_name,
1149  std::string exodus_var_name,
1150  unsigned int timestep)
1151 {
1152  LOG_SCOPE("copy_elemental_solution()", "ExodusII_IO");
1153 
1154  const unsigned int var_num = system.variable_number(system_var_name);
1155  // Assert that variable is an elemental one.
1156  //
1157  // NOTE: Currently, this reader is capable of reading only individual components of MONOMIAL_VEC
1158  // types, and each must be written out to its own CONSTANT MONOMIAL variable
1159  libmesh_error_msg_if((system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL))
1160  && (system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL_VEC)),
1161  "Error! Trying to copy elemental solution into a variable that is not of CONSTANT MONOMIAL nor CONSTANT MONOMIAL_VEC type.");
1162 
1164  const DofMap & dof_map = system.get_dof_map();
1165 
1166  // Map from element ID to elemental variable value. We need to use
1167  // a map here rather than a vector (e.g. elem_var_values) since the
1168  // libmesh element numbering can contain "holes". This is the case
1169  // if we are reading elemental var values from an adaptively refined
1170  // mesh that has not been sequentially renumbered.
1171  std::map<dof_id_type, Real> elem_var_value_map;
1172 
1173  // With Exodus files we only open them on processor 0, so that's the
1174  // where we have to do the data read too.
1175  if (system.comm().rank() == 0)
1176  {
1177  libmesh_error_msg_if(!exio_helper->opened_for_reading,
1178  "ERROR, ExodusII file must be opened for reading before copying an elemental solution!");
1179 
1180  exio_helper->read_elemental_var_values(exodus_var_name, timestep, elem_var_value_map);
1181  }
1182 
1183  const bool serial_on_zero = mesh.is_serial_on_zero();
1184 
1185  // If our mesh isn't serial, then non-root processors need to
1186  // request the data for their parts of the mesh and insert it
1187  // themselves.
1188  if (!serial_on_zero)
1189  {
1190  std::unordered_map<processor_id_type, std::vector<dof_id_type>> elem_ids_to_request;
1191  if (this->processor_id() != 0)
1192  {
1193  std::vector<dof_id_type> elem_ids;
1194  for (auto & elem : mesh.active_local_element_ptr_range())
1195  elem_ids.push_back(elem->id());
1196 
1197  if (!elem_ids.empty())
1198  elem_ids_to_request[0] = std::move(elem_ids);
1199  }
1200 
1201  auto value_gather_functor =
1202  [& elem_var_value_map]
1204  const std::vector<dof_id_type> & ids,
1205  std::vector<Real> & values)
1206  {
1207  const std::size_t query_size = ids.size();
1208  values.resize(query_size);
1209  for (std::size_t i=0; i != query_size; ++i)
1210  {
1211  if (const auto it = elem_var_value_map.find(ids[i]);
1212  it != elem_var_value_map.end())
1213  {
1214  values[i] = it->second;
1215  elem_var_value_map.erase(it);
1216  }
1217  else
1218  values[i] = std::numeric_limits<Real>::quiet_NaN();
1219  }
1220  };
1221 
1222  auto value_action_functor =
1223  [& elem_var_value_map]
1225  const std::vector<dof_id_type> & ids,
1226  const std::vector<Real> & values)
1227  {
1228  const std::size_t query_size = ids.size();
1229  for (std::size_t i=0; i != query_size; ++i)
1230  if (!libmesh_isnan(values[i]))
1231  elem_var_value_map[ids[i]] = values[i];
1232  };
1233 
1234  Real * value_ex = nullptr;
1235  Parallel::pull_parallel_vector_data
1236  (system.comm(), elem_ids_to_request, value_gather_functor,
1237  value_action_functor, value_ex);
1238  }
1239 
1240  std::map<dof_id_type, Real>::iterator
1241  it = elem_var_value_map.begin(),
1242  end = elem_var_value_map.end();
1243 
1244  // Everybody inserts the data they've received. If we're
1245  // serial_on_zero then proc 0 inserts everybody's data and other
1246  // procs have empty map ranges.
1247  for (; it!=end; ++it)
1248  {
1249  const Elem * elem = mesh.query_elem_ptr(it->first);
1250 
1251  if (elem && elem->n_comp(system.number(), var_num) > 0)
1252  {
1253  dof_id_type dof_index = elem->dof_number(system.number(), var_num, 0);
1254  if (serial_on_zero || dof_map.local_index(dof_index ))
1255  system.solution->set (dof_index, it->second);
1256  }
1257  }
1258 
1259  system.solution->close();
1260  system.update();
1261 }
1262 
1264  std::vector<std::string> system_var_names,
1265  std::vector<std::string> exodus_var_names,
1266  unsigned int timestep)
1267 {
1268  LOG_SCOPE("copy_scalar_solution()", "ExodusII_IO");
1269 
1270  libmesh_error_msg_if(!exio_helper->opened_for_reading,
1271  "ERROR, ExodusII file must be opened for reading before copying a scalar solution!");
1272 
1273  libmesh_error_msg_if(system_var_names.size() != exodus_var_names.size(),
1274  "ERROR, the number of system_var_names must match exodus_var_names.");
1275 
1276  std::vector<Real> values_from_exodus;
1277  read_global_variable(exodus_var_names, timestep, values_from_exodus);
1278 
1279 #ifdef LIBMESH_HAVE_MPI
1280  if (this->n_processors() > 1)
1281  {
1282  const Parallel::MessageTag tag = this->comm().get_unique_tag(1);
1283  if (this->processor_id() == this->n_processors()-1)
1284  this->comm().receive(0, values_from_exodus, tag);
1285  if (this->processor_id() == 0)
1286  this->comm().send(this->n_processors()-1, values_from_exodus, tag);
1287  }
1288 #endif
1289 
1290  if (system.processor_id() == (system.n_processors()-1))
1291  {
1292  const DofMap & dof_map = system.get_dof_map();
1293 
1294  for (auto i : index_range(system_var_names))
1295  {
1296  const unsigned int var_num = system.variable_scalar_number(system_var_names[i], 0);
1297 
1298  std::vector<dof_id_type> SCALAR_dofs;
1299  dof_map.SCALAR_dof_indices(SCALAR_dofs, var_num);
1300 
1301  system.solution->set (SCALAR_dofs[0], values_from_exodus[i]);
1302  }
1303  }
1304 
1305  system.solution->close();
1306  system.update();
1307 }
1308 
1309 void ExodusII_IO::read_elemental_variable(std::string elemental_var_name,
1310  unsigned int timestep,
1311  std::map<unsigned int, Real> & unique_id_to_value_map)
1312 {
1313  LOG_SCOPE("read_elemental_variable()", "ExodusII_IO");
1314 
1315  // Note that this function MUST be called before renumbering
1316  std::map<dof_id_type, Real> elem_var_value_map;
1317 
1318  exio_helper->read_elemental_var_values(elemental_var_name, timestep, elem_var_value_map);
1319  for (auto & pr : elem_var_value_map)
1320  {
1321  const Elem * elem = MeshInput<MeshBase>::mesh().query_elem_ptr(pr.first);
1322  unique_id_to_value_map.emplace(elem->top_parent()->unique_id(), pr.second);
1323  }
1324 }
1325 
1326 void ExodusII_IO::read_global_variable(std::vector<std::string> global_var_names,
1327  unsigned int timestep,
1328  std::vector<Real> & global_values)
1329 {
1330  LOG_SCOPE("read_global_variable()", "ExodusII_IO");
1331 
1332  std::size_t size = global_var_names.size();
1333  libmesh_error_msg_if(size == 0, "ERROR, empty list of global variables to read from the Exodus file.");
1334 
1335  // read the values for all global variables
1336  std::vector<Real> values_from_exodus;
1337  exio_helper->read_var_names(ExodusII_IO_Helper::GLOBAL);
1338  exio_helper->read_global_values(values_from_exodus, timestep);
1339  std::vector<std::string> global_var_names_exodus = exio_helper->global_var_names;
1340 
1341  if (values_from_exodus.size() == 0)
1342  return; // This will happen in parallel on procs that are not 0
1343 
1344  global_values.clear();
1345  for (std::size_t i = 0; i != size; ++i)
1346  {
1347  // for each global variable in global_var_names, look the corresponding one in global_var_names_from_exodus
1348  // and fill global_values accordingly
1349  auto it = find(global_var_names_exodus.begin(), global_var_names_exodus.end(), global_var_names[i]);
1350  if (it != global_var_names_exodus.end())
1351  global_values.push_back(values_from_exodus[it - global_var_names_exodus.begin()]);
1352  else
1353  libmesh_error_msg("ERROR, Global variable " << global_var_names[i] << \
1354  " not found in Exodus file.");
1355  }
1356 
1357 }
1358 
1360 {
1361  LOG_SCOPE("write_element_data()", "ExodusII_IO");
1362 
1363  // Be sure the file has been opened for writing!
1364  libmesh_error_msg_if(MeshOutput<MeshBase>::mesh().processor_id() == 0 && !exio_helper->opened_for_writing,
1365  "ERROR, ExodusII file must be initialized before outputting element variables.");
1366 
1367  // This function currently only works on serialized meshes. We rely
1368  // on having a reference to a non-const MeshBase object from our
1369  // MeshInput parent class to construct a MeshSerializer object,
1370  // similar to what is done in ExodusII_IO::write(). Note that
1371  // calling ExodusII_IO::write_timestep() followed by
1372  // ExodusII_IO::write_element_data() when the underlying Mesh is a
1373  // DistributedMesh will result in an unnecessary additional
1374  // serialization/re-parallelization step.
1375  // The "true" specifies that we only need the mesh serialized to processor 0
1377 
1378  // To be (possibly) filled with a filtered list of variable names to output.
1379  std::vector<std::string> names;
1380 
1381  // If _output_variables is populated, only output the monomials which are
1382  // also in the _output_variables vector.
1383  if (_output_variables.size() > 0)
1384  {
1385  // Create a list of CONSTANT MONOMIAL variable names
1386  std::vector<std::string> monomials;
1387  FEType type(CONSTANT, MONOMIAL);
1388  es.build_variable_names(monomials, &type);
1389 
1390  // Now concatenate a list of CONSTANT MONOMIAL_VEC variable names
1391  type = FEType(CONSTANT, MONOMIAL_VEC);
1392  es.build_variable_names(monomials, &type);
1393 
1394  // Filter that list against the _output_variables list. Note: if names is still empty after
1395  // all this filtering, all the monomial variables will be gathered
1396  for (const auto & var : monomials)
1397  if (std::find(_output_variables.begin(), _output_variables.end(), var) != _output_variables.end())
1398  names.push_back(var);
1399  }
1400 
1401  // If we pass in a list of names to "build_elemental_solution_vector()"
1402  // it'll filter the variables coming back.
1403  std::vector<Number> soln;
1404  es.build_elemental_solution_vector(soln, names);
1405 
1406  // Also, store the list of subdomains on which each variable is active
1407  std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
1408  es.get_vars_active_subdomains(names, vars_active_subdomains);
1409 
1410  if (soln.empty()) // If there is nothing to write just return
1411  return;
1412 
1413  // The data must ultimately be written block by block. This means that this data
1414  // must be sorted appropriately.
1415  if (MeshOutput<MeshBase>::mesh().processor_id())
1416  return;
1417 
1419 
1420 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1421 
1422  std::vector<std::string> complex_names =
1423  exio_helper->get_complex_names(names, _write_complex_abs);
1424 
1425  std::vector<std::set<subdomain_id_type>>
1426  complex_vars_active_subdomains =
1427  exio_helper->get_complex_vars_active_subdomains(vars_active_subdomains,
1429  exio_helper->initialize_element_variables(complex_names, complex_vars_active_subdomains);
1430 
1431  const std::vector<Real> complex_soln =
1432  complex_soln_components(soln, names.size(), _write_complex_abs);
1433 
1434  exio_helper->write_element_values(mesh, complex_soln, _timestep, complex_vars_active_subdomains);
1435 
1436 #else
1437  exio_helper->initialize_element_variables(names, vars_active_subdomains);
1438  exio_helper->write_element_values(mesh, soln, _timestep, vars_active_subdomains);
1439 #endif
1440 }
1441 
1442 
1443 
1444 void
1446 (const EquationSystems & es,
1447  const std::set<std::string> * system_names,
1448  const std::string & var_suffix)
1449 {
1450  LOG_SCOPE("write_element_data_from_discontinuous_nodal_data()", "ExodusII_IO");
1451 
1452  // Be sure that some other function has already opened the file and prepared it
1453  // for writing. This is the same behavior as the write_element_data() function
1454  // which we are trying to mimic.
1455  libmesh_error_msg_if(MeshOutput<MeshBase>::mesh().processor_id() == 0 && !exio_helper->opened_for_writing,
1456  "ERROR, ExodusII file must be initialized before outputting element variables.");
1457 
1458  // This function currently only works on serialized meshes. The
1459  // "true" flag specifies that we only need the mesh serialized to
1460  // processor 0
1463  true);
1464 
1465  // Note: in general we want to respect the contents of
1466  // _output_variables, only building a solution vector with values
1467  // from the requested variables. First build a list of all variable
1468  // names, then throw out ones that aren't in _output_variables, if
1469  // any.
1470  std::vector<std::string> var_names;
1471  es.build_variable_names (var_names, /*fetype=*/nullptr, system_names);
1472 
1473  // Get a subset of all variable names that are CONSTANT,
1474  // MONOMIALs. We treat those slightly differently since they can
1475  // truly only have a single value per Elem.
1476  //
1477  // Should the same apply here for CONSTANT MONOMIAL_VECs? [CW]
1478  // That is, get rid of 'const' on 'fe_type' and rerun:
1479  // fe_type = FEType(CONSTANT, MONOMIAL_VEC);
1480  // es.build_variable_names(monomial_var_names, &fe_type);
1481  // Then, es.find_variable_numbers() can be used without a type
1482  // (since we know for sure they're monomials) like:
1483  // var_nums = es.find_variable_numbers(monomial_var_names)
1484  // for which the DOF indices for 'var_nums' have to be resolved
1485  // manually like in build_elemental_solution_vector()
1486  std::vector<std::string> monomial_var_names;
1487  const FEType fe_type(CONSTANT, MONOMIAL);
1488  es.build_variable_names(monomial_var_names, &fe_type);
1489 
1490  // Remove all names from var_names that are not in _output_variables.
1491  // Note: This approach avoids errors when the user provides invalid
1492  // variable names in _output_variables, as the code will not try to
1493  // write a variable that doesn't exist.
1494  if (!_output_variables.empty())
1495  {
1496  var_names.erase
1497  (std::remove_if
1498  (var_names.begin(),
1499  var_names.end(),
1500  [this](const std::string & name)
1501  {return !std::count(_output_variables.begin(),
1502  _output_variables.end(),
1503  name);}),
1504  var_names.end());
1505 
1506  // Also filter the monomial variable names.
1507  monomial_var_names.erase
1508  (std::remove_if
1509  (monomial_var_names.begin(),
1510  monomial_var_names.end(),
1511  [this](const std::string & name)
1512  {return !std::count(_output_variables.begin(),
1513  _output_variables.end(),
1514  name);}),
1515  monomial_var_names.end());
1516  }
1517 
1518  // Build a solution vector, limiting the results to the variables in
1519  // var_names and the Systems in system_names, and only computing values
1520  // at the vertices.
1521  std::vector<Number> v;
1523  (v, system_names, &var_names, /*vertices_only=*/true);
1524 
1525  // Get active subdomains for each variable in var_names.
1526  std::vector<std::set<subdomain_id_type>> vars_active_subdomains;
1527  es.get_vars_active_subdomains(var_names, vars_active_subdomains);
1528 
1529  // Determine names of variables to write based on the number of
1530  // nodes/vertices the elements in different subdomains have.
1532  std::map<subdomain_id_type, unsigned int> subdomain_id_to_vertices_per_elem;
1533  for (const auto & elem : mesh.active_element_ptr_range())
1534  {
1535  // Try to insert key/value pair into the map. If this returns
1536  // false, check the returned iterator's value to make sure it
1537  // matches. It shouldn't actually be possible for this to fail
1538  // (since if the Mesh was like this it would have already
1539  // failed) but it doesn't hurt to be on the safe side.
1540  auto pr2 = subdomain_id_to_vertices_per_elem.emplace
1541  (elem->subdomain_id(), elem->n_vertices());
1542  libmesh_error_msg_if(!pr2.second && pr2.first->second != elem->n_vertices(),
1543  "Elem with different number of vertices found.");
1544  }
1545 
1546  // Determine "derived" variable names. These names are created by
1547  // starting with the base variable name and appending the user's
1548  // variable_suffix (default: "_elem_node_") followed by a node id.
1549  //
1550  // Not every derived variable will be active on every subdomain,
1551  // even if the original variable _is_ active. Subdomains can have
1552  // different geometric element types (with differing numbers of
1553  // nodes), so some of the derived variable names will be inactive on
1554  // those subdomains.
1555  //
1556  // Since we would otherwise generate the same name once per
1557  // subdomain, we keep the list of names unique as we are creating
1558  // it. We can't use a std::set for this because we don't want the
1559  // variables names to be in a different order from the order
1560  // they were written in the call to: build_discontinuous_solution_vector()
1561  //
1562  // The list of derived variable names includes one for each vertex,
1563  // for higher-order elements we currently only write out vertex
1564  // values, but this could be changed in the future without too much
1565  // trouble.
1566  std::vector<std::string> derived_var_names;
1567 
1568  // Keep track of mapping from derived_name to (orig_name, node_id)
1569  // pair. We will use this later to determine whether a given
1570  // variable is active on a given subdomain.
1571  std::map<std::string, std::pair<std::string, unsigned int>>
1572  derived_name_to_orig_name_and_node_id;
1573 
1574  for (const auto & pr : subdomain_id_to_vertices_per_elem)
1575  {
1576  const subdomain_id_type sbd_id = pr.first;
1577  const unsigned int vertices_per_elem =
1578  subdomain_id_to_vertices_per_elem[sbd_id];
1579 
1580  std::ostringstream oss;
1581  for (unsigned int n=0; n<vertices_per_elem; ++n)
1582  for (const auto & orig_var_name : var_names)
1583  {
1584  oss.str("");
1585  oss.clear();
1586  oss << orig_var_name << var_suffix << n;
1587  std::string derived_name = oss.str();
1588 
1589  // Only add this var name if it's not already in the list.
1590  if (!std::count(derived_var_names.begin(), derived_var_names.end(), derived_name))
1591  {
1592  derived_var_names.push_back(derived_name);
1593  // Add entry for derived_name -> (orig_name, node_id) mapping.
1594  derived_name_to_orig_name_and_node_id[derived_name] =
1595  std::make_pair(orig_var_name, n);
1596  }
1597  }
1598  }
1599 
1600  // For each derived variable name, determine whether it is active
1601  // based on how many nodes/vertices the elements in a given subdomain have,
1602  // and whether they were active on the subdomain to begin with.
1603  std::vector<std::set<subdomain_id_type>>
1604  derived_vars_active_subdomains(derived_var_names.size());
1605 
1606  // A new data structure for keeping track of a list of variable names
1607  // that are in the discontinuous solution vector on each subdomain. Used
1608  // for indexing. Note: if a variable was inactive at the System level,
1609  // an entry for it will still be in the discontinuous solution vector,
1610  // but it will just have a value of zero. On the other hand, when we
1611  // create the derived variable names some of them are "inactive" on
1612  // different subdomains in the sense that they don't exist at all, i.e.
1613  // there is no zero padding for them. We need to be able to distinguish
1614  // between these two types in order to do the indexing into this vector
1615  // correctly.
1616  std::map<subdomain_id_type, std::vector<std::string>>
1617  subdomain_to_var_names;
1618 
1619  for (auto derived_var_id : index_range(derived_var_names))
1620  {
1621  const auto & derived_name = derived_var_names[derived_var_id];
1622  const auto & [orig_name, node_id] =
1623  libmesh_map_find (derived_name_to_orig_name_and_node_id,
1624  derived_name);
1625 
1626  // For each subdomain, determine whether the current variable
1627  // should be active on that subdomain.
1628  for (const auto & pr : subdomain_id_to_vertices_per_elem)
1629  {
1630  // Convenience variables for the current subdomain and the
1631  // number of nodes elements in this subdomain have.
1632  subdomain_id_type sbd_id = pr.first;
1633  unsigned int vertices_per_elem_this_sbd =
1634  subdomain_id_to_vertices_per_elem[sbd_id];
1635 
1636  // Check whether variable orig_name was active on this
1637  // subdomain to begin with by looking in the
1638  // vars_active_subdomains container. We assume that the
1639  // location of orig_name in the var_names vector matches its
1640  // index in the vars_active_subdomains container.
1641  auto var_loc = std::find(var_names.begin(), var_names.end(), orig_name);
1642  libmesh_error_msg_if(var_loc == var_names.end(),
1643  "Variable " << orig_name << " somehow not found in var_names array.");
1644  auto var_id = std::distance(var_names.begin(), var_loc);
1645 
1646  // The derived_var will only be active if this subdomain has
1647  // enough vertices for that to be the case.
1648  if (node_id < vertices_per_elem_this_sbd)
1649  {
1650  // Regardless of whether the original variable was not active on this subdomain,
1651  // the discontinuous solution vector will have zero padding for it, and
1652  // we will need to account for it. Therefore it should still be added to
1653  // the subdomain_to_var_names data structure!
1654  subdomain_to_var_names[sbd_id].push_back(derived_name);
1655 
1656  // If the original variable was not active on the
1657  // current subdomain, it should not be added to the
1658  // derived_vars_active_subdomains data structure, since
1659  // it will not be written to the Exodus file.
1660 
1661  // Determine if the original variable was active on the
1662  // current subdomain.
1663  bool orig_var_active =
1664  (vars_active_subdomains[var_id].empty() ||
1665  vars_active_subdomains[var_id].count(sbd_id));
1666 
1667  // And only if it was, add it to the
1668  // derived_vars_active_subdomains data structure.
1669  if (orig_var_active)
1670  derived_vars_active_subdomains[derived_var_id].insert(sbd_id);
1671  }
1672  } // end loop over subdomain_id_to_vertices_per_elem
1673  } // end loop over derived_var_names
1674 
1675  // At this point we've built the "true" list of derived names, but
1676  // if there are any CONSTANT MONOMIALS in this list, we now want to
1677  // remove all but one copy of them from the derived_var_names list,
1678  // and rename them in (but not remove them from) the
1679  // subdomain_to_var_names list, and then update the
1680  // derived_vars_active_subdomains containers before finally calling
1681  // the Exodus helper functions.
1682  for (auto & derived_var_name : derived_var_names)
1683  {
1684  // Get the original name associated with this derived name.
1685  const auto & name_and_id =
1686  libmesh_map_find (derived_name_to_orig_name_and_node_id,
1687  derived_var_name);
1688 
1689  // Convenience variables for the map entry's contents.
1690  const std::string & orig_name = name_and_id.first;
1691 
1692  // Was the original name a constant monomial?
1693  if (std::count(monomial_var_names.begin(),
1694  monomial_var_names.end(),
1695  orig_name))
1696  {
1697  // Rename this variable in the subdomain_to_var_names vectors.
1698  for (auto & pr : subdomain_to_var_names)
1699  {
1700  // Reference to ordered list of variable names on this subdomain.
1701  auto & name_vec = pr.second;
1702 
1703  auto name_vec_it =
1704  std::find(name_vec.begin(),
1705  name_vec.end(),
1706  derived_var_name);
1707 
1708  if (name_vec_it != name_vec.end())
1709  {
1710  // Actually rename it back to the orig_name, dropping
1711  // the "_elem_corner_" stuff.
1712  *name_vec_it = orig_name;
1713  }
1714  }
1715 
1716  // Finally, rename the variable in the derived_var_names vector itself.
1717  derived_var_name = orig_name;
1718  } // if (monomial)
1719  } // end loop over derived names
1720 
1721  // Now remove duplicate entries from derived_var_names after the first.
1722  // Also update the derived_vars_active_subdomains container in a consistent way.
1723  {
1724  std::vector<std::string> derived_var_names_edited;
1725  std::vector<std::set<subdomain_id_type>> derived_vars_active_subdomains_edited;
1726  std::vector<unsigned int> found_first(monomial_var_names.size());
1727 
1728  for (auto i : index_range(derived_var_names))
1729  {
1730  const auto & derived_var_name = derived_var_names[i];
1731  const auto & active_set = derived_vars_active_subdomains[i];
1732 
1733  // Determine whether we will keep this derived variable name in
1734  // the final container.
1735  bool keep = true;
1736  for (auto j : index_range(monomial_var_names))
1737  if (derived_var_name == monomial_var_names[j])
1738  {
1739  if (!found_first[j])
1740  found_first[j] = 1;
1741 
1742  else
1743  keep = false;
1744  }
1745 
1746  // We also don't keep variables that are not active on any subdomains.
1747  // Contrary to other uses of the var_active_subdomains container where
1748  // the empty set means "all" subdomains, here it really means "none".
1749  if (active_set.empty())
1750  keep = false;
1751 
1752  if (keep)
1753  {
1754  derived_var_names_edited.push_back(derived_var_name);
1755  derived_vars_active_subdomains_edited.push_back(active_set);
1756  }
1757  }
1758 
1759  // We built the filtered ranges, now swap them with the originals.
1760  derived_var_names.swap(derived_var_names_edited);
1761  derived_vars_active_subdomains.swap(derived_vars_active_subdomains_edited);
1762  }
1763 
1764 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1765  // Build complex variable names "r_foo", "i_foo", "a_foo" and the lists of
1766  // subdomains on which they are active.
1767  auto complex_var_names =
1768  exio_helper->get_complex_names(derived_var_names,
1769  _write_complex_abs);
1770  auto complex_vars_active_subdomains =
1771  exio_helper->get_complex_vars_active_subdomains(derived_vars_active_subdomains,
1772  _write_complex_abs);
1773  auto complex_subdomain_to_var_names =
1774  exio_helper->get_complex_subdomain_to_var_names(subdomain_to_var_names,
1775  _write_complex_abs);
1776 
1777  // Make expanded version of vector "v" in which each entry in the
1778  // original expands to an ("r_", "i_", "a_") triple.
1779  // "nco" is the number of complex outputs, which depends on whether
1780  // or not we are writing out the complex magnitudes.
1781  std::vector<Real> complex_v;
1782  int nco = _write_complex_abs ? 3 : 2;
1783  complex_v.reserve(nco * v.size());
1784  for (const auto & val : v)
1785  {
1786  complex_v.push_back(val.real());
1787  complex_v.push_back(val.imag());
1788  if (_write_complex_abs)
1789  complex_v.push_back(std::abs(val));
1790  }
1791 
1792  // Finally, initialize storage for the variables and write them to file.
1793  exio_helper->initialize_element_variables
1794  (complex_var_names, complex_vars_active_subdomains);
1795  exio_helper->write_element_values_element_major
1796  (mesh, complex_v, _timestep,
1797  complex_vars_active_subdomains,
1798  complex_var_names,
1799  complex_subdomain_to_var_names);
1800 #else
1801 
1802  // Call function which writes the derived variable names to the
1803  // Exodus file.
1804  exio_helper->initialize_element_variables(derived_var_names, derived_vars_active_subdomains);
1805 
1806  // ES::build_discontinuous_solution_vector() creates a vector with
1807  // an element-major ordering, so call Helper::write_element_values()
1808  // passing false for the last argument.
1809  exio_helper->write_element_values_element_major
1810  (mesh, v, _timestep,
1811  derived_vars_active_subdomains,
1812  derived_var_names,
1813  subdomain_to_var_names);
1814 #endif
1815 }
1816 
1817 
1818 
1819 void ExodusII_IO::write_nodal_data (const std::string & fname,
1820  const std::vector<Number> & soln,
1821  const std::vector<std::string> & names)
1822 {
1823  LOG_SCOPE("write_nodal_data()", "ExodusII_IO");
1824 
1826 
1827  int num_vars = cast_int<int>(names.size());
1828  dof_id_type num_nodes = mesh.n_nodes();
1829 
1830  // The names of the variables to be output
1831  std::vector<std::string> output_names;
1832 
1833  if (_allow_empty_variables || !_output_variables.empty())
1834  output_names = _output_variables;
1835  else
1836  output_names = names;
1837 
1838 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1839  std::vector<std::string> complex_names =
1840  exio_helper->get_complex_names(output_names,
1842 
1843  // Call helper function for opening/initializing data, giving it the
1844  // complex variable names
1845  this->write_nodal_data_common(fname, complex_names, /*continuous=*/true);
1846 #else
1847  // Call helper function for opening/initializing data
1848  this->write_nodal_data_common(fname, output_names, /*continuous=*/true);
1849 #endif
1850 
1851  if (mesh.processor_id())
1852  return;
1853 
1854  // This will count the number of variables actually output
1855  for (int c=0; c<num_vars; c++)
1856  {
1857  std::stringstream name_to_find;
1858 
1859  std::vector<std::string>::iterator pos =
1860  std::find(output_names.begin(), output_names.end(), names[c]);
1861  if (pos == output_names.end())
1862  continue;
1863 
1864  unsigned int variable_name_position =
1865  cast_int<unsigned int>(pos - output_names.begin());
1866 
1867  // Set up temporary vectors to be passed to Exodus to write the
1868  // nodal values for a single variable at a time.
1869 #ifdef LIBMESH_USE_REAL_NUMBERS
1870  std::vector<Number> cur_soln;
1871 
1872  // num_nodes is either exactly how much space we will need for
1873  // each vector, or a safe upper bound for the amount of memory
1874  // we will require when there are gaps in the numbering.
1875  cur_soln.reserve(num_nodes);
1876 #else
1877  std::vector<Real> real_parts;
1878  std::vector<Real> imag_parts;
1879  std::vector<Real> magnitudes;
1880  real_parts.reserve(num_nodes);
1881  imag_parts.reserve(num_nodes);
1882  if (_write_complex_abs)
1883  magnitudes.reserve(num_nodes);
1884 #endif
1885 
1886  // There could be gaps in soln based on node numbering, but in
1887  // serial the empty numbers are left empty.
1888  // There could also be offsets in soln based on "fake" nodes
1889  // inserted on each processor (because NumericVector indices
1890  // have to be contiguous); the helper keeps track of those.
1891  // We now copy the proper solution values contiguously into
1892  // "cur_soln", removing the gaps.
1893  for (const auto & node : mesh.node_ptr_range())
1894  {
1895  const dof_id_type idx =
1896  (exio_helper->node_id_to_vec_id(node->id()))
1897  * num_vars + c;
1898 #ifdef LIBMESH_USE_REAL_NUMBERS
1899  cur_soln.push_back(soln[idx]);
1900 #else
1901  real_parts.push_back(soln[idx].real());
1902  imag_parts.push_back(soln[idx].imag());
1903  if (_write_complex_abs)
1904  magnitudes.push_back(std::abs(soln[idx]));
1905 #endif
1906  }
1907 
1908  // If we're adding extra sides, we need to add their data too.
1909  //
1910  // Because soln was created from a parallel NumericVector, its
1911  // numbering was contiguous on each processor; we need to use
1912  // the same offsets here, and we need to loop through elements
1913  // from earlier ranks first.
1914  if (exio_helper->get_add_sides())
1915  {
1916  std::vector<std::vector<const Elem *>>
1917  elems_by_pid(mesh.n_processors());
1918 
1919  for (const auto & elem : mesh.active_element_ptr_range())
1920  elems_by_pid[elem->processor_id()].push_back(elem);
1921 
1922  for (auto p : index_range(elems_by_pid))
1923  {
1924  dof_id_type global_idx =
1925  exio_helper->added_node_offset_on(p) * num_vars + c;
1926  for (const Elem * elem : elems_by_pid[p])
1927  {
1928  for (auto s : elem->side_index_range())
1929  {
1931  continue;
1932 
1933  const std::vector<unsigned int> side_nodes =
1934  elem->nodes_on_side(s);
1935 
1936  for (auto n : index_range(side_nodes))
1937  {
1938  libmesh_ignore(n);
1939  libmesh_assert_less(global_idx, soln.size());
1940 #ifdef LIBMESH_USE_REAL_NUMBERS
1941  cur_soln.push_back(soln[global_idx]);
1942 #else
1943  real_parts.push_back(soln[global_idx].real());
1944  imag_parts.push_back(soln[global_idx].imag());
1945  if (_write_complex_abs)
1946  magnitudes.push_back(std::abs(soln[global_idx]));
1947 #endif
1948  global_idx += num_vars;
1949  }
1950  }
1951  }
1952  }
1953  }
1954 
1955  // Finally, actually call the Exodus API to write to file.
1956 #ifdef LIBMESH_USE_REAL_NUMBERS
1957  exio_helper->write_nodal_values(variable_name_position+1, cur_soln, _timestep);
1958 #else
1959  int nco = _write_complex_abs ? 3 : 2;
1960  exio_helper->write_nodal_values(nco*variable_name_position+1, real_parts, _timestep);
1961  exio_helper->write_nodal_values(nco*variable_name_position+2, imag_parts, _timestep);
1962  if (_write_complex_abs)
1963  exio_helper->write_nodal_values(3*variable_name_position+3, magnitudes, _timestep);
1964 #endif
1965 
1966  }
1967 }
1968 
1969 
1970 
1971 
1972 void ExodusII_IO::write_information_records (const std::vector<std::string> & records)
1973 {
1975  return;
1976 
1977  libmesh_error_msg_if(!exio_helper->opened_for_writing,
1978  "ERROR, ExodusII file must be initialized before outputting information records.");
1979 
1980  exio_helper->write_information_records(records);
1981 }
1982 
1983 
1984 
1985 void ExodusII_IO::write_global_data (const std::vector<Number> & soln,
1986  const std::vector<std::string> & names)
1987 {
1988  LOG_SCOPE("write_global_data()", "ExodusII_IO");
1989 
1991  return;
1992 
1993  libmesh_error_msg_if(!exio_helper->opened_for_writing,
1994  "ERROR, ExodusII file must be initialized before outputting global variables.");
1995 
1996 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1997 
1998  std::vector<std::string> complex_names =
1999  exio_helper->get_complex_names(names,
2001 
2002  exio_helper->initialize_global_variables(complex_names);
2003 
2004  const std::vector<Real> complex_soln =
2005  complex_soln_components(soln, names.size(), _write_complex_abs);
2006 
2007  exio_helper->write_global_values(complex_soln, _timestep);
2008 
2009 #else
2010  exio_helper->initialize_global_variables(names);
2011  exio_helper->write_global_values(soln, _timestep);
2012 #endif
2013 }
2014 
2015 
2016 
2017 void ExodusII_IO::write_timestep (const std::string & fname,
2018  const EquationSystems & es,
2019  const int timestep,
2020  const Real time,
2021  const std::set<std::string> * system_names)
2022 {
2023  _timestep = timestep;
2024  MeshOutput<MeshBase>::write_equation_systems(fname,es,system_names);
2025 
2027  return;
2028 
2029  exio_helper->write_timestep(timestep, time);
2030 }
2031 
2032 
2033 void ExodusII_IO::write_equation_systems (const std::string & fname,
2034  const EquationSystems & es,
2035  const std::set<std::string> * system_names)
2036 {
2037  write_timestep(fname, es, 1, 0, system_names);
2038 }
2039 
2040 
2042 {
2043  libmesh_error_msg_if(!exio_helper->opened_for_writing,
2044  "ERROR, ExodusII file must be opened for writing "
2045  "before calling ExodusII_IO::write_elemsets()!");
2046 
2048  exio_helper->write_elemsets(mesh);
2049 }
2050 
2051 void
2053 write_sideset_data(int timestep,
2054  const std::vector<std::string> & var_names,
2055  const std::vector<std::set<boundary_id_type>> & side_ids,
2056  const std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
2057 {
2058  libmesh_error_msg_if(!exio_helper->opened_for_writing,
2059  "ERROR, ExodusII file must be opened for writing "
2060  "before calling ExodusII_IO::write_sideset_data()!");
2061 
2063  exio_helper->write_sideset_data(mesh, timestep, var_names, side_ids, bc_vals);
2064 }
2065 
2066 
2067 
2068 void
2070 read_sideset_data(int timestep,
2071  std::vector<std::string> & var_names,
2072  std::vector<std::set<boundary_id_type>> & side_ids,
2073  std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
2074 {
2075  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2076  "ERROR, ExodusII file must be opened for reading "
2077  "before calling ExodusII_IO::read_sideset_data()!");
2078 
2080  exio_helper->read_sideset_data(mesh, timestep, var_names, side_ids, bc_vals);
2081 }
2082 
2083 
2084 
2085 void
2087 get_sideset_data_indices (std::map<BoundaryInfo::BCTuple, unsigned int> & bc_array_indices)
2088 
2089 {
2090  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2091  "ERROR, ExodusII file must be opened for reading "
2092  "before calling ExodusII_IO::get_sideset_data_indices()!");
2093 
2095  exio_helper->get_sideset_data_indices(mesh, bc_array_indices);
2096 }
2097 
2098 void
2100 get_nodeset_data_indices (std::map<BoundaryInfo::NodeBCTuple, unsigned int> & bc_array_indices)
2101 {
2102  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2103  "ERROR, ExodusII file must be opened for reading "
2104  "before calling ExodusII_IO::get_nodeset_data_indices()!");
2105 
2106  exio_helper->get_nodeset_data_indices(bc_array_indices);
2107 }
2108 
2109 void
2111 write_nodeset_data (int timestep,
2112  const std::vector<std::string> & var_names,
2113  const std::vector<std::set<boundary_id_type>> & node_boundary_ids,
2114  const std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
2115 {
2116  libmesh_error_msg_if(!exio_helper->opened_for_writing,
2117  "ERROR, ExodusII file must be opened for writing "
2118  "before calling ExodusII_IO::write_nodeset_data()!");
2119 
2120  exio_helper->write_nodeset_data(timestep, var_names, node_boundary_ids, bc_vals);
2121 }
2122 
2123 
2124 
2125 void
2127 read_nodeset_data (int timestep,
2128  std::vector<std::string> & var_names,
2129  std::vector<std::set<boundary_id_type>> & node_boundary_ids,
2130  std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
2131 {
2132  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2133  "ERROR, ExodusII file must be opened for reading "
2134  "before calling ExodusII_IO::read_nodeset_data()!");
2135 
2136  exio_helper->read_nodeset_data(timestep, var_names, node_boundary_ids, bc_vals);
2137 }
2138 
2139 void
2141 write_elemset_data (int timestep,
2142  const std::vector<std::string> & var_names,
2143  const std::vector<std::set<elemset_id_type>> & elemset_ids_in,
2144  const std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
2145 {
2146  libmesh_error_msg_if(!exio_helper->opened_for_writing,
2147  "ERROR, ExodusII file must be opened for writing "
2148  "before calling ExodusII_IO::write_elemset_data()!");
2149 
2150  exio_helper->write_elemset_data(timestep, var_names, elemset_ids_in, elemset_vals);
2151 }
2152 
2153 
2154 
2155 void
2157 read_elemset_data (int timestep,
2158  std::vector<std::string> & var_names,
2159  std::vector<std::set<elemset_id_type>> & elemset_ids_in,
2160  std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
2161 {
2162  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2163  "ERROR, ExodusII file must be opened for reading "
2164  "before calling ExodusII_IO::read_elemset_data()!");
2165 
2166  exio_helper->read_elemset_data(timestep, var_names, elemset_ids_in, elemset_vals);
2167 }
2168 
2169 void
2170 ExodusII_IO::get_elemset_data_indices (std::map<std::pair<dof_id_type, elemset_id_type>, unsigned int> & elemset_array_indices)
2171 {
2172  libmesh_error_msg_if(!exio_helper->opened_for_reading,
2173  "ERROR, ExodusII file must be opened for reading "
2174  "before calling ExodusII_IO::get_elemset_data_indices()!");
2175 
2176  exio_helper->get_elemset_data_indices(elemset_array_indices);
2177 }
2178 
2179 
2180 void ExodusII_IO::write (const std::string & fname)
2181 {
2182  LOG_SCOPE("write()", "ExodusII_IO");
2183 
2185 
2186  // We may need to gather a DistributedMesh to output it, making that
2187  // const qualifier in our constructor a dirty lie
2188  // The "true" specifies that we only need the mesh serialized to processor 0
2189  MeshSerializer serialize
2190  (const_cast<MeshBase &>(mesh),
2192 
2193  libmesh_assert( !exio_helper->opened_for_writing );
2194 
2195  // If the user has set the append flag here, it doesn't really make
2196  // sense: the intent of this function is to write a Mesh with no
2197  // data, while "appending" is really intended to add data to an
2198  // existing file. If we're verbose, print a message to this effect.
2199  if (_append && _verbose)
2200  libmesh_warning("Warning: Appending in ExodusII_IO::write() does not make sense.\n"
2201  "Creating a new file instead!");
2202 
2203  exio_helper->create(fname);
2204  exio_helper->initialize(fname,mesh);
2205  exio_helper->write_nodal_coordinates(mesh);
2206  exio_helper->write_elements(mesh);
2207  exio_helper->write_sidesets(mesh);
2208  exio_helper->write_nodesets(mesh);
2209  exio_helper->write_elemsets(mesh);
2210 
2211  if ((mesh.get_boundary_info().n_edge_conds() > 0) && _verbose)
2212  libmesh_warning("Warning: Mesh contains edge boundary IDs, but these "
2213  "are not supported by the ExodusII format.");
2214 }
2215 
2216 
2217 
2218 void ExodusII_IO::write_nodal_data_discontinuous (const std::string & fname,
2219  const std::vector<Number> & soln,
2220  const std::vector<std::string> & names)
2221 {
2222  LOG_SCOPE("write_nodal_data_discontinuous()", "ExodusII_IO");
2223 
2225 
2226 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2227 
2228  std::vector<std::string> complex_names =
2229  exio_helper->get_complex_names(names,
2231 
2232  // Call helper function for opening/initializing data, giving it the
2233  // complex variable names
2234  this->write_nodal_data_common(fname, complex_names, /*continuous=*/false);
2235 #else
2236  // Call helper function for opening/initializing data
2237  this->write_nodal_data_common(fname, names, /*continuous=*/false);
2238 #endif
2239 
2240  if (mesh.processor_id())
2241  return;
2242 
2243  int num_vars = cast_int<int>(names.size());
2244  libmesh_assert_equal_to(soln.size() % num_vars, 0);
2245  int num_nodes = soln.size() / num_vars;
2246  libmesh_assert_equal_to(exio_helper->num_nodes, num_nodes);
2247 
2248 #ifndef NDEBUG
2249  if (!this->get_add_sides())
2250  {
2251  int num_real_nodes = 0;
2252  for (const auto & elem : mesh.active_element_ptr_range())
2253  num_real_nodes += elem->n_nodes();
2254  libmesh_assert_equal_to(num_real_nodes, num_nodes);
2255  }
2256 #endif
2257 
2258  for (int c=0; c<num_vars; c++)
2259  {
2260 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2261  std::vector<Real> real_parts(num_nodes);
2262  std::vector<Real> imag_parts(num_nodes);
2263  std::vector<Real> magnitudes;
2264  if (_write_complex_abs)
2265  magnitudes.resize(num_nodes);
2266 
2267  // The number of complex outputs depends on whether or not we are
2268  // writing out the absolute values.
2269  int nco = _write_complex_abs ? 3 : 2;
2270 
2271  for (int i=0; i<num_nodes; ++i)
2272  {
2273  real_parts[i] = soln[i*num_vars + c].real();
2274  imag_parts[i] = soln[i*num_vars + c].imag();
2275  if (_write_complex_abs)
2276  magnitudes[i] = std::abs(soln[i*num_vars + c]);
2277  }
2278  exio_helper->write_nodal_values(nco*c+1, real_parts, _timestep);
2279  exio_helper->write_nodal_values(nco*c+2, imag_parts, _timestep);
2280  if (_write_complex_abs)
2281  exio_helper->write_nodal_values(3*c+3, magnitudes, _timestep);
2282 #else
2283  // Copy out this variable's solution
2284  std::vector<Number> cur_soln(num_nodes);
2285 
2286  for (int i=0; i<num_nodes; i++)
2287  cur_soln[i] = soln[i*num_vars + c];
2288 
2289  exio_helper->write_nodal_values(c+1,cur_soln,_timestep);
2290 #endif
2291  }
2292 }
2293 
2294 
2295 
2297  const std::vector<std::string> & names,
2298  bool continuous)
2299 {
2301 
2302  // This function can be called multiple times, we only want to open
2303  // the ExodusII file the first time it's called.
2304  if (!exio_helper->opened_for_writing)
2305  {
2306  // If we're appending, open() the file with read_only=false,
2307  // otherwise create() it and write the contents of the mesh to
2308  // it.
2309  if (_append)
2310  {
2311  // We do our writing only from proc 0, to avoid race
2312  // conditions with Exodus 8
2314  {
2315  exio_helper->open(fname.c_str(), /*read_only=*/false);
2316  // If we're appending, it's not valid to call exio_helper->initialize()
2317  // or exio_helper->initialize_nodal_variables(), but we do need to set up
2318  // certain aspects of the Helper object itself, such as the number of nodes
2319  // and elements. We do that by reading the header...
2320  exio_helper->read_and_store_header_info();
2321 
2322  // ...and reading the block info
2323  exio_helper->read_block_info();
2324  }
2325  // Keep other processors aware of what we've done on root
2326  else
2327  {
2328  exio_helper->opened_for_writing = true;
2329  exio_helper->current_filename = fname;
2330  }
2331  }
2332  else
2333  {
2334  exio_helper->create(fname);
2335 
2336  // But some of our write calls are parallel-only, due to
2337  // calls to parallel-only getter functions.
2338  exio_helper->initialize(fname, mesh, !continuous);
2339 
2340  exio_helper->write_nodal_coordinates(mesh, !continuous);
2341  exio_helper->write_elements(mesh, !continuous);
2342 
2343  exio_helper->write_sidesets(mesh);
2344  exio_helper->write_nodesets(mesh);
2345  exio_helper->write_elemsets(mesh);
2346 
2347  exio_helper->initialize_nodal_variables(names);
2348  }
2349  }
2350  else
2351  {
2352  // We are already open for writing, so check that the filename
2353  // passed to this function matches the filename currently in use
2354  // by the helper.
2355  libmesh_error_msg_if(fname != exio_helper->current_filename,
2356  "Error! This ExodusII_IO object is already associated with file: "
2357  << exio_helper->current_filename
2358  << ", cannot use it with requested file: "
2359  << fname);
2360  }
2361 }
2362 
2363 const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
2364 {
2365  exio_helper->read_var_names(ExodusII_IO_Helper::NODAL);
2366  return exio_helper->nodal_var_names;
2367 }
2368 
2369 const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
2370 {
2371  exio_helper->read_var_names(ExodusII_IO_Helper::ELEMENTAL);
2372  return exio_helper->elem_var_names;
2373 }
2374 
2375 const std::vector<std::string> & ExodusII_IO::get_global_var_names()
2376 {
2377  exio_helper->read_var_names(ExodusII_IO_Helper::GLOBAL);
2378  return exio_helper->global_var_names;
2379 }
2380 
2381 const std::vector<int> & ExodusII_IO::get_elem_num_map() const
2382 {
2383  // We could make this function non-const and have it call
2384  // exio_helper->read_elem_num_map() before returning a reference...
2385  // but the intention is that this will be called some time after a
2386  // mesh is read in, in which case it would be doing extra work to
2387  // read in the elem_num_map twice.
2388  return exio_helper->elem_num_map;
2389 }
2390 
2391 const std::vector<int> & ExodusII_IO::get_node_num_map() const
2392 {
2393  return exio_helper->node_num_map;
2394 }
2395 
2397 {
2398  // Provide a warning when accessing the helper object
2399  // since it is a non-public API and is likely to see
2400  // future API changes
2401  libmesh_experimental();
2402 
2403  return *exio_helper;
2404 }
2405 
2406 
2407 void ExodusII_IO::set_hdf5_writing(bool write_hdf5)
2408 {
2409  exio_helper->set_hdf5_writing(write_hdf5);
2410 }
2411 
2412 
2414 {
2415  _disc_bex = disc_bex;
2416 }
2417 
2418 
2419 
2420 // LIBMESH_HAVE_EXODUS_API is not defined, declare error() versions of functions...
2421 #else
2422 
2423 
2424 
2425 ExodusII_IO::~ExodusII_IO () = default;
2426 
2427 
2428 
2429 void ExodusII_IO::read (const std::string &)
2430 {
2431  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2432 }
2433 
2434 
2435 
2436 ExodusHeaderInfo ExodusII_IO::read_header (const std::string &)
2437 {
2438  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2439 }
2440 
2441 
2442 
2443 void ExodusII_IO::verbose (bool)
2444 {
2445  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2446 }
2447 
2448 
2449 
2451 {
2452  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2453 }
2454 
2455 
2456 
2457 void ExodusII_IO::write_as_dimension(unsigned)
2458 {
2459  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2460 }
2461 
2462 
2463 
2465 {
2466  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2467 }
2468 
2469 
2470 
2471 void ExodusII_IO::append(bool)
2472 {
2473  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2474 }
2475 
2476 
2477 
2479 {
2480  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2481 }
2482 
2483 
2484 
2486 {
2487  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2488  return false;
2489 }
2490 
2491 
2492 
2493 const std::vector<Real> & ExodusII_IO::get_time_steps()
2494 {
2495  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2496 }
2497 
2498 
2499 
2501 {
2502  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2503 }
2504 
2505 
2507  std::string,
2508  std::string,
2509  unsigned int)
2510 {
2511  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2512 }
2513 
2514 
2515 
2517  std::string,
2518  std::string,
2519  unsigned int)
2520 {
2521  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2522 }
2523 
2524 
2525 
2527  std::vector<std::string>,
2528  std::vector<std::string>,
2529  unsigned int)
2530 {
2531  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2532 }
2533 
2534 
2535 
2537 {
2538  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2539 }
2540 
2541 
2542 
2543 void
2545 (const EquationSystems &,
2546  const std::set<std::string> *,
2547  const std::string & )
2548 {
2549  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2550 }
2551 
2552 
2553 
2554 void ExodusII_IO::write_nodal_data (const std::string &,
2555  const std::vector<Number> &,
2556  const std::vector<std::string> &)
2557 {
2558  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2559 }
2560 
2561 
2562 
2563 void ExodusII_IO::write_information_records (const std::vector<std::string> &)
2564 {
2565  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2566 }
2567 
2568 
2569 
2570 void ExodusII_IO::write_global_data (const std::vector<Number> &,
2571  const std::vector<std::string> &)
2572 {
2573  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2574 }
2575 
2576 
2577 
2578 void ExodusII_IO::write_timestep (const std::string &,
2579  const EquationSystems &,
2580  const int,
2581  const Real,
2582  const std::set<std::string> *)
2583 {
2584  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2585 }
2586 
2587 
2588 void ExodusII_IO::write_equation_systems (const std::string &,
2589  const EquationSystems &,
2590  const std::set<std::string> *)
2591 {
2592  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2593 }
2594 
2595 
2597 {
2598  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2599 }
2600 
2601 void
2603 write_sideset_data (int,
2604  const std::vector<std::string> &,
2605  const std::vector<std::set<boundary_id_type>> &,
2606  const std::vector<std::map<BoundaryInfo::BCTuple, Real>> &)
2607 {
2608  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2609 }
2610 
2611 
2612 
2613 void
2615 read_sideset_data (int,
2616  std::vector<std::string> &,
2617  std::vector<std::set<boundary_id_type>> &,
2618  std::vector<std::map<BoundaryInfo::BCTuple, Real>> &)
2619 {
2620  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2621 }
2622 
2623 void
2625 get_sideset_data_indices (std::map<BoundaryInfo::BCTuple, unsigned int> &)
2626 
2627 {
2628  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2629 }
2630 
2631 void
2633 write_nodeset_data (int,
2634  const std::vector<std::string> &,
2635  const std::vector<std::set<boundary_id_type>> &,
2636  const std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> &)
2637 {
2638  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2639 }
2640 
2641 void
2643 read_nodeset_data (int,
2644  std::vector<std::string> &,
2645  std::vector<std::set<boundary_id_type>> &,
2646  std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> &)
2647 {
2648  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2649 }
2650 
2651 void
2653 get_nodeset_data_indices (std::map<BoundaryInfo::NodeBCTuple, unsigned int> &)
2654 {
2655  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2656 }
2657 
2658 void
2660 write_elemset_data (int,
2661  const std::vector<std::string> &,
2662  const std::vector<std::set<elemset_id_type>> &,
2663  const std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> &)
2664 {
2665  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2666 }
2667 
2668 void
2670 read_elemset_data (int,
2671  std::vector<std::string> &,
2672  std::vector<std::set<elemset_id_type>> &,
2673  std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> &)
2674 {
2675  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2676 }
2677 
2678 void
2680 get_elemset_data_indices (std::map<std::pair<dof_id_type, elemset_id_type>, unsigned int> &)
2681 {
2682  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2683 }
2684 
2685 void ExodusII_IO::write (const std::string &)
2686 {
2687  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2688 }
2689 
2690 
2691 
2692 void ExodusII_IO::write_nodal_data_discontinuous (const std::string &,
2693  const std::vector<Number> &,
2694  const std::vector<std::string> &)
2695 {
2696  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2697 }
2698 
2699 
2700 
2701 void ExodusII_IO::write_nodal_data_common(std::string,
2702  const std::vector<std::string> &,
2703  bool)
2704 {
2705  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2706 }
2707 
2708 
2709 const std::vector<std::string> & ExodusII_IO::get_elem_var_names()
2710 {
2711  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2712 }
2713 
2714 const std::vector<std::string> & ExodusII_IO::get_nodal_var_names()
2715 {
2716  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2717 }
2718 
2719 const std::vector<std::string> & ExodusII_IO::get_global_var_names()
2720 {
2721  libmesh_error_msg("ERROR, ExodusII API is not defined.");
2722 }
2723 
2724 void ExodusII_IO::set_hdf5_writing(bool) {}
2725 
2726 #endif // LIBMESH_HAVE_EXODUS_API
2727 } // namespace libMesh
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 MT & mesh() const
Definition: mesh_output.h:259
The FPEDisabler class puts Floating-Point Exception (FPE) trapping on hold during its lifetime...
Definition: fpe_disabler.h:49
const std::vector< std::string > & get_global_var_names()
Return list of the global variable names.
Definition: exodusII_io.C:2375
void use_mesh_dimension_instead_of_spatial_dimension(bool val)
In the general case, meshes containing 2D elements can be manifolds living in 3D space, thus by default we write all meshes with the Exodus dimension set to LIBMESH_DIM = mesh.spatial_dimension().
Definition: exodusII_io.C:976
This is the EquationSystems class.
virtual void reserve_nodes(const dof_id_type nn)=0
Reserves space for a known number of nodes.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
unsigned int variable_scalar_number(std::string_view var, unsigned int component) const
Definition: system.h:2489
void build_variable_names(std::vector< std::string > &var_names, const FEType *type=nullptr, const std::set< std::string > *system_names=nullptr) const
Fill the input vector var_names with the names of the variables for each system.
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
Definition: mesh_base.h:1703
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
A Node is like a Point, but with more information.
Definition: node.h:52
std::string & nodeset_name(boundary_id_type id)
boost::multiprecision::float128 real(const boost::multiprecision::float128 in)
void write_as_dimension(unsigned dim)
Directly control the num_dim which is written to the Exodus file.
Definition: exodusII_io.C:983
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
ExodusII_IO_Helper & get_exio_helper()
Return a reference to the ExodusII_IO_Helper object.
Definition: exodusII_io.C:2396
void write_elemset_data(int timestep, const std::vector< std::string > &var_names, const std::vector< std::set< elemset_id_type >> &elemset_ids_in, const std::vector< std::map< std::pair< dof_id_type, elemset_id_type >, Real >> &elemset_vals)
The Exodus format can also store values on elemsets.
Definition: exodusII_io.C:2141
This class is used as both an external data structure for passing around Exodus file header informati...
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
std::size_t n_edge_conds() const
std::vector< std::string > _extra_integer_vars
An optional list of variables in the EXODUS file that are to be used to set extra integers when loadi...
Definition: exodusII_io.h:644
unsigned int dim
static const int invalid_id
An invalid_id that can be returned to signal failure in case something goes wrong.
void write_sideset_data(int timestep, const std::vector< std::string > &var_names, const std::vector< std::set< boundary_id_type >> &side_ids, const std::vector< std::map< BoundaryInfo::BCTuple, Real >> &bc_vals)
The Exodus format can also store values on sidesets.
Definition: exodusII_io.C:2053
void build_discontinuous_solution_vector(std::vector< Number > &soln, const std::set< std::string > *system_names=nullptr, const std::vector< std::string > *var_names=nullptr, bool vertices_only=false, bool add_sides=false) const
Fill the input vector soln with solution values.
const Elem * top_parent() const
Definition: elem.h:3056
std::vector< bool > elems_of_dimension
A vector of bools describing what dimension elements have been encountered when reading a mesh...
Definition: mesh_input.h:107
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
void read_elemset_data(int timestep, std::vector< std::string > &var_names, std::vector< std::set< elemset_id_type >> &elemset_ids_in, std::vector< std::map< std::pair< dof_id_type, elemset_id_type >, Real >> &elemset_vals)
Read all the elemset data at a particular timestep.
Definition: exodusII_io.C:2157
unsigned int add_elem_integer(std::string name, bool allocate_data=true, dof_id_type default_value=DofObject::invalid_id)
Register an integer datum (of type dof_id_type) to be added to each element in the mesh...
Definition: mesh_base.C:560
void add_elemset_code(dof_id_type code, MeshBase::elemset_type id_set)
Tabulate a user-defined "code" for elements which belong to the element sets specified in id_set...
Definition: mesh_base.C:398
Defines mapping from libMesh element types to LS-DYNA element types or vice-versa.
Definition: dyna_io.h:125
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
processor_id_type rank() const
unique_id_type unique_id() const
Definition: dof_object.h:844
const std::vector< int > & get_node_num_map() const
Identical to the behavior of get_elem_num_map(), but for the node_num_map instead.
Definition: exodusII_io.C:2391
const Parallel::Communicator & comm() const
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
void get_sideset_data_indices(std::map< BoundaryInfo::BCTuple, unsigned int > &bc_array_indices)
Similar to read_sideset_data(), but instead of creating one std::map per sideset per variable...
Definition: exodusII_io.C:2087
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type get_elemset_code(const MeshBase::elemset_type &id_set) const
Definition: mesh_base.C:439
ExodusII_IO(MeshBase &mesh, bool single_precision=false)
Constructor.
Definition: exodusII_io.C:122
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
bool libmesh_isnan(T x)
void append(bool val)
If true, this flag will cause the ExodusII_IO object to attempt to open an existing file for writing...
Definition: exodusII_io.C:998
Real distance(const Point &p)
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
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
void get_nodeset_data_indices(std::map< BoundaryInfo::NodeBCTuple, unsigned int > &bc_array_indices)
Similar to read_nodeset_data(), but instead of creating one std::map per nodeset per variable...
Definition: exodusII_io.C:2100
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn...
Definition: dof_map.C:2547
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
virtual bool is_serial_on_zero() const
Definition: mesh_base.h:218
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void get_vars_active_subdomains(const std::vector< std::string > &names, std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) const
Retrieve vars_active_subdomains, which indicates the active subdomains for each variable in names...
int _timestep
Stores the current value of the timestep when calling ExodusII_IO::write_timestep().
Definition: exodusII_io.h:625
static const ElementDefinition & find_elem_definition(dyna_int_type dyna_elem, int dim, int p)
Finds the ElementDefinition corresponding to a particular element type.
Definition: dyna_io.C:755
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)
Writes a discontinuous solution at a specific timestep.
Definition: exodusII_io.C:209
processor_id_type n_processors() const
void write_nodeset_data(int timestep, const std::vector< std::string > &var_names, const std::vector< std::set< boundary_id_type >> &node_boundary_ids, const std::vector< std::map< BoundaryInfo::NodeBCTuple, Real >> &bc_vals)
The Exodus format can also store values on nodesets.
Definition: exodusII_io.C:2111
void libmesh_ignore(const Args &...)
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
const dof_id_type n_nodes
Definition: tecplot_io.C:67
unsigned int number() const
Definition: system.h:2350
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
void read_nodeset_data(int timestep, std::vector< std::string > &var_names, std::vector< std::set< boundary_id_type >> &node_boundary_ids, std::vector< std::map< BoundaryInfo::NodeBCTuple, Real >> &bc_vals)
Read all the nodeset data at a particular timestep.
Definition: exodusII_io.C:2127
virtual bool get_add_sides() override
Definition: exodusII_io.C:1012
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)
This method implements writing a mesh with discontinuous data to a specified file where the data is t...
Definition: mesh_output.C:89
dof_id_type id() const
Definition: dof_object.h:828
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:57
void set_extra_integer_vars(const std::vector< std::string > &extra_integer_vars)
Set the elemental variables in the Exodus file to be read into extra element integers.
Definition: exodusII_io.C:179
void set_output_variables(const std::vector< std::string > &output_variables, bool allow_empty=true)
Sets the list of variable names to be included in the output.
Definition: exodusII_io.C:184
virtual unsigned int n_nodes() const =0
bool _verbose
should we be verbose?
Definition: exodusII_io.h:630
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
Writes out the solution for no specific time or timestep.
Definition: exodusII_io.C:2033
bool _append
Default false.
Definition: exodusII_io.h:636
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
void copy_nodal_solution(System &system, std::string system_var_name, std::string exodus_var_name, unsigned int timestep=1)
If we read in a nodal solution while reading in a mesh, we can attempt to copy that nodal solution in...
Definition: exodusII_io.C:1043
void set_default_mapping_type(const ElemMappingType type)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:821
void copy_elemental_solution(System &system, std::string system_var_name, std::string exodus_var_name, unsigned int timestep=1)
If we read in a elemental solution while reading in a mesh, we can attempt to copy that elemental sol...
Definition: exodusII_io.C:1147
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
const std::vector< Real > & get_time_steps()
Definition: exodusII_io.C:1020
unsigned int n_elemsets() const
Returns the number of unique elemset ids which have been added via add_elemset_code(), which is the size of the _all_elemset_ids set.
Definition: mesh_base.C:424
std::vector< std::string > _output_variables
The names of the variables to be output.
Definition: exodusII_io.h:650
This is the ExodusII_IO_Helper class.
void set_extra_datum(const unsigned int index, const T value)
Sets the value on this object of the extra datum associated with index, which should have been obtain...
Definition: dof_object.h:1126
void broadcast(const Parallel::Communicator &comm)
Broadcasts data from processor 0 to other procs using the provided Communicator.
void write_information_records(const std::vector< std::string > &)
Write out information records.
Definition: exodusII_io.C:1972
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:1692
void read_elemental_variable(std::string elemental_var_name, unsigned int timestep, std::map< unsigned int, Real > &unique_id_to_value_map)
Given an elemental variable and a time step, returns a mapping from the elements (top parent) unique ...
Definition: exodusII_io.C:1309
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
void verbose(bool set_verbosity)
Set the flag indicating if we should be verbose.
Definition: exodusII_io.C:959
void write_elemsets()
Write elemsets stored on the Mesh to file.
Definition: exodusII_io.C:2041
void write_nodal_data_common(std::string fname, const std::vector< std::string > &names, bool continuous=true)
This function factors out a bunch of code which is common to the write_nodal_data() and write_nodal_d...
Definition: exodusII_io.C:2296
An object whose state is distributed along a set of processors.
virtual void read(const std::string &name) override
This method implements reading a mesh from a specified file.
Definition: exodusII_io.C:244
bool _disc_bex
Set to true (false is the default) to generate independent nodes for every Bezier Extraction element...
Definition: exodusII_io.h:677
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
void read_global_variable(std::vector< std::string > global_var_names, unsigned int timestep, std::vector< Real > &global_values)
Given a vector of global variables and a time step, returns the values of the global variable at the ...
Definition: exodusII_io.C:1326
unsigned int add_node_datum(const std::string &name, bool allocate_data=true, const T *default_value=nullptr)
Register a datum (of type T) to be added to each node in the mesh.
Definition: mesh_base.h:2341
std::string & sideset_name(boundary_id_type id)
void write_global_data(const std::vector< Number > &, const std::vector< std::string > &)
Write out global variables.
Definition: exodusII_io.C:1985
void write_nodal_data_discontinuous(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
Write out a discontinuous nodal solution.
Definition: exodusII_io.C:2218
void copy_scalar_solution(System &system, std::vector< std::string > system_var_names, std::vector< std::string > exodus_var_names, unsigned int timestep=1)
Copy global variables into scalar variables of a System object.
Definition: exodusII_io.C:1263
std::string enum_to_string(const T e)
virtual const Elem * elem_ptr(const dof_id_type i) const =0
bool _write_complex_abs
By default, when complex numbers are enabled, for each variable we write out three values: the real p...
Definition: exodusII_io.h:671
std::unique_ptr< ExodusII_IO_Helper > exio_helper
Only attempt to instantiate an ExodusII helper class if the Exodus API is defined.
Definition: exodusII_io.h:619
const std::vector< int > & get_elem_num_map() const
Returns a const reference to the elem_num_map, which is a vector that is created when a Mesh is read ...
Definition: exodusII_io.C:2381
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
void write_discontinuous_exodusII(const std::string &name, const EquationSystems &es, const std::set< std::string > *system_names=nullptr)
Writes a exodusII file with discontinuous data.
Definition: exodusII_io.C:193
static int get_exodus_version()
Definition: exodusII_io.C:168
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2508
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
void write_complex_magnitude(bool val)
Set the flag indicating whether the complex modulus should be written when complex numbers are enable...
Definition: exodusII_io.C:969
virtual unsigned short dim() const =0
Temporarily serialize a DistributedMesh for non-distributed-mesh capable code paths.
virtual void write(const std::string &fname) override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:2180
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Writes out the solution at a specific timestep.
Definition: exodusII_io.C:2017
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
virtual ~ExodusII_IO()
Definition: exodusII_io.C:238
const std::vector< std::string > & get_elem_var_names()
Return list of the elemental variable names.
Definition: exodusII_io.C:2369
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
Add shell face shellface of element number elem with boundary id id to the boundary information data ...
static const bool value
Definition: xdr_io.C:54
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void write_element_data(const EquationSystems &es)
Write out element solution.
Definition: exodusII_io.C:1359
ExodusHeaderInfo read_header(const std::string &name)
Read only the header information, instead of the entire mesh.
Definition: exodusII_io.C:923
bool _allow_empty_variables
Flag which controls the behavior of _output_variables: .) If true, _output_variables is allowed to re...
Definition: exodusII_io.h:660
void set_default_mapping_data(const unsigned char data)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:839
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
void write_element_data_from_discontinuous_nodal_data(const EquationSystems &es, const std::set< std::string > *system_names=nullptr, const std::string &var_suffix="_elem_node_")
Similar to the function above, but instead of only handling (CONSTANT, MONOMIAL) data, writes out a general discontinuous solution field, e.g.
Definition: exodusII_io.C:1446
void set_discontinuous_bex(bool disc_bex)
Set to true (false is the default) to generate independent nodes for every Bezier Extraction element ...
Definition: exodusII_io.C:2413
boost::multiprecision::float128 imag(const boost::multiprecision::float128)
T get_extra_datum(const unsigned int index) const
Gets the value on this object of the extra datum associated with index, which should have been obtain...
Definition: dof_object.h:1146
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
Write out a nodal solution.
Definition: exodusII_io.C:1819
std::vector< unsigned int > nodes
Definition: dyna_io.h:142
processor_id_type processor_id() const
virtual Order default_order() const =0
void write_added_sides(bool val)
By default, we only write out the elements physically stored in the mesh.
Definition: exodusII_io.C:1005
const DofMap & get_dof_map() const
Definition: system.h:2374
void set_hdf5_writing(bool write_hdf5)
Set to true (the default) to write files in an HDF5-based file format (when HDF5 is available)...
Definition: exodusII_io.C:2407
processor_id_type processor_id() const
Definition: dof_object.h:905
const std::vector< std::string > & get_nodal_var_names()
Return list of the nodal variable names.
Definition: exodusII_io.C:2363
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual void reserve_elem(const dof_id_type ne)=0
Reserves space for a known number of elements.
void read_sideset_data(int timestep, std::vector< std::string > &var_names, std::vector< std::set< boundary_id_type >> &side_ids, std::vector< std::map< BoundaryInfo::BCTuple, Real >> &bc_vals)
Similar to write_sideset_data(), this function is used to read the data at a particular timestep...
Definition: exodusII_io.C:2070
void get_elemset_data_indices(std::map< std::pair< dof_id_type, elemset_id_type >, unsigned int > &elemset_array_indices)
Similar to read_elemset_data(), but instead of creating one std::map per nodeset per variable...
Definition: exodusII_io.C:2170
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
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
void set_extra_integer(const unsigned int index, const dof_id_type value)
Sets the value on this object of the extra integer associated with index, which should have been obta...
Definition: dof_object.h:1086
virtual dof_id_type n_nodes() const =0
void set_coordinate_offset(Point p)
Allows you to set a vector that is added to the coordinates of all of the nodes.
Definition: exodusII_io.C:990
void build_elemental_solution_vector(std::vector< Number > &soln, std::vector< std::string > &names) const
Retrieve the solution data for CONSTANT MONOMIALs and/or components of CONSTANT MONOMIAL_VECs.
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.
uint8_t dof_id_type
Definition: id_types.h:67
bool local_index(dof_id_type dof_index) const
Definition: dof_map.h:839
static bool redundant_added_side(const Elem &elem, unsigned int side)
std::vector< std::vector< std::vector< Real > > > bex_dense_constraint_vecs