libMesh
exodusII_io_helper.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 #include "libmesh/exodusII_io_helper.h"
20 
21 
22 #ifdef LIBMESH_HAVE_EXODUS_API
23 
24 // libMesh includes
25 #include "libmesh/boundary_info.h"
26 #include "libmesh/enum_elem_type.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/fpe_disabler.h"
30 #include "libmesh/remote_elem.h"
31 #include "libmesh/system.h"
32 #include "libmesh/numeric_vector.h"
33 #include "libmesh/enum_to_string.h"
34 #include "libmesh/enum_elem_type.h"
35 #include "libmesh/int_range.h"
36 #include "libmesh/utility.h"
37 #include "libmesh/libmesh_logging.h"
38 
39 #ifdef DEBUG
40 #include "libmesh/mesh_tools.h" // for elem_types warning
41 #endif
42 
43 #include <libmesh/ignore_warnings.h>
44 namespace exII {
45 extern "C" {
46 #include "exodusII.h" // defines MAX_LINE_LENGTH, MAX_STR_LENGTH used later
47 }
48 }
49 #include <libmesh/restore_warnings.h>
50 
51 // C++ includes
52 #include <algorithm>
53 #include <cfenv> // workaround for HDF5 bug
54 #include <cstdlib> // std::strtol
55 #include <sstream>
56 #include <unordered_map>
57 
58 // Anonymous namespace for file local data and helper functions
59 namespace
60 {
61 
62 using namespace libMesh;
63 
64 // File scope constant node/edge/face mapping arrays.
65 // 2D inverse face map definitions.
66 // These take a libMesh ID and turn it into an Exodus ID
67 const std::vector<int> trishell3_inverse_edge_map = {3, 4, 5};
68 const std::vector<int> quadshell4_inverse_edge_map = {3, 4, 5, 6};
69 
70 // 3D node map definitions
71 // The hex27, prism20-21, and tet14 appear to be the only elements
72 // with a non-identity mapping between Exodus' node numbering and
73 // libmesh's. Exodus doesn't even number prisms hierarchically!
74 const std::vector<int> hex27_node_map = {
75  // Vertex and mid-edge nodes
76  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
77  // Mid-face nodes and center node
78  21, 25, 24, 26, 23, 22, 20};
79 //20 21 22 23 24 25 26 // LibMesh indices
80 
81 const std::vector<int> hex27_inverse_node_map = {
82  // Vertex and mid-edge nodes
83  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
84  // Mid-face nodes and center node
85  26, 20, 25, 24, 22, 21, 23};
86 //20 21 22 23 24 25 26
87 
88 const std::vector<int> prism20_node_map = {
89  // Vertices
90  0, 1, 2, 3, 4, 5,
91  // Matching mid-edge nodes
92  6, 7, 8, 9, 10, 11, 12, 13, 14,
93  // Non-matching nodes
94  19, 17, 18, 15, 16};
95 //15 16 17 18 19 // LibMesh indices
96 
97 const std::vector<int> prism20_inverse_node_map = {
98  // Vertices
99  0, 1, 2, 3, 4, 5,
100  // Matching mid-edge nodes
101  6, 7, 8, 9, 10, 11, 12, 13, 14,
102  // Non-matching nodes
103  18, 19, 16, 17, 15};
104 //15 16 17 18 19
105 
106 const std::vector<int> prism21_node_map = {
107  // Vertices
108  0, 1, 2, 3, 4, 5,
109  // Matching mid-edge nodes
110  6, 7, 8, 9, 10, 11, 12, 13, 14,
111  // Non-matching nodes
112  20, 18, 19, 16, 17, 15};
113 //15 16 17 18 19 20 // LibMesh indices
114 
115 const std::vector<int> prism21_inverse_node_map = {
116  // Vertices
117  0, 1, 2, 3, 4, 5,
118  // Matching mid-edge nodes
119  6, 7, 8, 9, 10, 11, 12, 13, 14,
120  // Non-matching nodes
121  20, 18, 19, 16, 17, 15};
122 //15 16 17 18 19 20
123 
124 const std::vector<int> tet14_node_map = {
125  // Vertex and mid-edge nodes
126  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
127  // Mid-face nodes
128  10, 13, 11, 12};
129 //10 11 12 13 // LibMesh indices
130 
131 const std::vector<int> tet14_inverse_node_map = {
132  // Vertex and mid-edge nodes
133  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
134  // Mid-face nodes
135  10, 12, 13, 11};
136 //10 11 12 13
137 
138 
139 // 3D face map definitions
140 const std::vector<int> tet_face_map = {1, 2, 3, 0};
141 const std::vector<int> hex_face_map = {1, 2, 3, 4, 0, 5};
142 const std::vector<int> prism_face_map = {1, 2, 3, 0, 4};
143 
144 // These take a libMesh ID and turn it into an Exodus ID
145 const std::vector<int> tet_inverse_face_map = {4, 1, 2, 3};
146 const std::vector<int> hex_inverse_face_map = {5, 1, 2, 3, 4, 6};
147 const std::vector<int> prism_inverse_face_map = {4, 1, 2, 3, 5};
148 
149 // 3D element edge maps. Map 0-based Exodus id -> libMesh id.
150 // Commented out until we have code that needs it, to keep compiler
151 // warnings happy.
152 // const std::vector<int> hex_edge_map =
153  // {0,1,2,3,8,9,10,11,4,5,7,6};
154 
155 // 3D inverse element edge maps. Map libmesh edge ids to 1-based Exodus edge ids.
156 // Commented out until we have code that needs it, to keep compiler
157 // warnings happy.
158 // const std::vector<int> hex_inverse_edge_map =
159  // {1,2,3,4,9,10,12,11,5,6,7,8};
160 
164  int inquire(libMesh::ExodusII_IO_Helper & e2h, exII::ex_inquiry req_info_in, std::string error_msg="")
165  {
166  int ret_int = 0;
167  char ret_char = 0;
168  float ret_float = 0.;
169 
170  e2h.ex_err = exII::ex_inquire(e2h.ex_id,
171  req_info_in,
172  &ret_int,
173  &ret_float,
174  &ret_char);
175 
176  EX_CHECK_ERR(e2h.ex_err, error_msg);
177 
178  return ret_int;
179  }
180 
181  // Bezier Extraction test: if we see BEx data we had better be in a
182  // Bezier element block
183  inline bool is_bezier_elem(const char * elem_type_str)
184  {
185  // Reading Bezier Extraction from Exodus files requires ExodusII v8
186 #if EX_API_VERS_NODOT < 800
187  libMesh::libmesh_ignore(elem_type_str);
188  return false;
189 #else
190  if (strlen(elem_type_str) <= 4)
191  return false;
192  return (std::string(elem_type_str, elem_type_str+4) == "BEX_");
193 #endif
194  }
195 
196 
197  std::map<subdomain_id_type, std::vector<unsigned int>>
198  build_subdomain_map(const MeshBase & mesh, bool add_sides, subdomain_id_type & subdomain_id_end)
199  {
200  std::map<subdomain_id_type, std::vector<unsigned int>> subdomain_map;
201 
202  // If we've been asked to add side elements, those will go in
203  // their own blocks.
204  if (add_sides)
205  {
206  std::set<subdomain_id_type> sbd_ids;
207  mesh.subdomain_ids(sbd_ids);
208  if (!sbd_ids.empty())
209  subdomain_id_end = *sbd_ids.rbegin()+1;
210  }
211 
212  // Loop through element and map between block and element vector.
213  for (const auto & elem : mesh.active_element_ptr_range())
214  {
215  // We skip writing infinite elements to the Exodus file, so
216  // don't put them in the subdomain_map. That way the number of
217  // blocks should be correct.
218  if (elem->infinite())
219  continue;
220 
221  subdomain_map[ elem->subdomain_id() ].push_back(elem->id());
222 
223  // If we've been asked to add side elements, those will go in their own
224  // blocks. We don't have any ids to list for elements that don't
225  // explicitly exist in the mesh, but we do an entry to keep
226  // track of the number of elements we'll add in each new block.
227  if (add_sides)
228  for (auto s : elem->side_index_range())
229  {
231  continue;
232 
233  auto & marker =
234  subdomain_map[subdomain_id_end + elem->side_type(s)];
235  if (marker.empty())
236  marker.push_back(1);
237  else
238  ++marker[0];
239  }
240  }
241 
242  if (!add_sides && !subdomain_map.empty())
243  subdomain_id_end = subdomain_map.rbegin()->first + 1;
244 
245  return subdomain_map;
246  }
247 } // end anonymous namespace
248 
249 
250 
251 namespace libMesh
252 {
253 
254 // ExodusII_IO_Helper::Conversion static data
255 const int ExodusII_IO_Helper::Conversion::invalid_id = std::numeric_limits<int>::max();
256 
258  bool v,
259  bool run_only_on_proc0,
260  bool single_precision) :
261  ParallelObject(parent),
262  ex_id(0),
263  ex_err(0),
264  header_info(), // zero-initialize
265  title(header_info.title),
266  num_dim(header_info.num_dim),
267  num_nodes(header_info.num_nodes),
268  num_elem(header_info.num_elem),
269  num_elem_blk(header_info.num_elem_blk),
270  num_edge(header_info.num_edge),
271  num_edge_blk(header_info.num_edge_blk),
272  num_node_sets(header_info.num_node_sets),
273  num_side_sets(header_info.num_side_sets),
274  num_elem_sets(header_info.num_elem_sets),
275  num_global_vars(0),
276  num_sideset_vars(0),
277  num_nodeset_vars(0),
278  num_elemset_vars(0),
279  num_elem_this_blk(0),
280  num_nodes_per_elem(0),
281  num_attr(0),
282  num_elem_all_sidesets(0),
283  num_elem_all_elemsets(0),
284  bex_num_elem_cvs(0),
285  num_time_steps(0),
286  num_nodal_vars(0),
287  num_elem_vars(0),
288  verbose(v),
289  opened_for_writing(false),
290  opened_for_reading(false),
291  _run_only_on_proc0(run_only_on_proc0),
292  _opened_by_create(false),
293  _elem_vars_initialized(false),
294  _global_vars_initialized(false),
295  _nodal_vars_initialized(false),
296  _use_mesh_dimension_instead_of_spatial_dimension(false),
297  _write_hdf5(true),
298  _end_elem_id(0),
299  _write_as_dimension(0),
300  _single_precision(single_precision)
301 {
302  title.resize(MAX_LINE_LENGTH+1);
303  elem_type.resize(MAX_STR_LENGTH);
306 }
307 
308 
309 
311 
312 
313 
315 {
316  return EX_API_VERS_NODOT;
317 }
318 
319 
320 
321 // Initialization function for conversion_map object
323 {
324  auto convert_type = [this](ElemType type,
325  std::string_view exodus_type,
326  const std::vector<int> * node_map = nullptr,
327  const std::vector<int> * inverse_node_map = nullptr,
328  const std::vector<int> * side_map = nullptr,
329  const std::vector<int> * inverse_side_map = nullptr,
330  const std::vector<int> * shellface_map = nullptr,
331  const std::vector<int> * inverse_shellface_map = nullptr,
332  size_t shellface_index_offset = 0)
333  {
334  std::unique_ptr<Elem> elem = Elem::build(type);
335  auto & conv = conversion_map[elem->dim()][type];
336  conv.libmesh_type = type;
337  conv.exodus_type = exodus_type;
338  conv.node_map = node_map;
339  conv.inverse_node_map = inverse_node_map;
340  conv.side_map = side_map;
341  conv.inverse_side_map = inverse_side_map;
342  conv.shellface_map = shellface_map;
343  conv.inverse_shellface_map = inverse_shellface_map;
344  conv.shellface_index_offset = shellface_index_offset;
345  conv.n_nodes = elem->n_nodes();
346  for (int d = elem->dim()+1; d <= 3; ++d)
347  conversion_map[d][type] = conv;
348  };
349 
350  convert_type(NODEELEM, "SPHERE");
351  convert_type(EDGE2, "EDGE2");
352  convert_type(EDGE3, "EDGE3");
353  convert_type(EDGE4, "EDGE4");
354  convert_type(QUAD4, "QUAD4");
355  convert_type(QUAD8, "QUAD8");
356  convert_type(QUAD9, "QUAD9");
357  convert_type(QUADSHELL4, "SHELL4", nullptr, nullptr, nullptr,
358  /* inverse_side_map = */ &quadshell4_inverse_edge_map,
359  nullptr, nullptr, /* shellface_index_offset = */ 2);
360  convert_type(QUADSHELL8, "SHELL8", nullptr, nullptr, nullptr,
361  /* inverse_side_map = */ &quadshell4_inverse_edge_map,
362  nullptr, nullptr, /* shellface_index_offset = */ 2);
363  convert_type(QUADSHELL9, "SHELL9", nullptr, nullptr, nullptr,
364  /* inverse_side_map = */ &quadshell4_inverse_edge_map,
365  nullptr, nullptr, /* shellface_index_offset = */ 2);
366 
367  convert_type(TRI3, "TRI3");
368  convert_type(TRI6, "TRI6");
369  convert_type(TRI7, "TRI7");
370  // Exodus does weird things to triangle side mapping in 3D. See
371  // https://sandialabs.github.io/seacas-docs/html/element_types.html#tri
372  conversion_map[3][TRI3].inverse_side_map = &trishell3_inverse_edge_map;
373  conversion_map[3][TRI3].shellface_index_offset = 2;
374  conversion_map[3][TRI6].inverse_side_map = &trishell3_inverse_edge_map;
375  conversion_map[3][TRI6].shellface_index_offset = 2;
376  conversion_map[3][TRI7].inverse_side_map = &trishell3_inverse_edge_map;
377  conversion_map[3][TRI7].shellface_index_offset = 2;
378 
379  convert_type(TRISHELL3, "TRISHELL3", nullptr, nullptr, nullptr,
380  /* inverse_side_map = */ &trishell3_inverse_edge_map,
381  nullptr, nullptr, /* shellface_index_offset = */ 2);
382  convert_type(TRI3SUBDIVISION, "TRI3");
383  convert_type(HEX8, "HEX8", nullptr, nullptr,
384  &hex_face_map, &hex_inverse_face_map);
385  convert_type(HEX20, "HEX20", nullptr, nullptr,
386  &hex_face_map, &hex_inverse_face_map);
387  convert_type(HEX27, "HEX27", &hex27_node_map,
388  &hex27_inverse_node_map,
389  &hex_face_map, &hex_inverse_face_map);
390  convert_type(TET4, "TETRA4", nullptr, nullptr,
391  &tet_face_map, &tet_inverse_face_map);
392  convert_type(TET10, "TETRA10", nullptr, nullptr,
393  &tet_face_map, &tet_inverse_face_map);
394  convert_type(TET14, "TETRA14", &tet14_node_map,
395  &tet14_inverse_node_map,
396  &tet_face_map, &tet_inverse_face_map);
397  convert_type(PRISM6, "WEDGE", nullptr, nullptr,
398  &prism_face_map, &prism_inverse_face_map);
399  convert_type(PRISM15, "WEDGE15", nullptr, nullptr,
400  &prism_face_map, &prism_inverse_face_map);
401  convert_type(PRISM18, "WEDGE18", nullptr, nullptr,
402  &prism_face_map, &prism_inverse_face_map);
403  convert_type(PRISM20, "WEDGE20", &prism20_node_map,
404  &prism20_inverse_node_map,
405  &prism_face_map, &prism_inverse_face_map);
406  convert_type(PRISM21, "WEDGE21", &prism21_node_map,
407  &prism21_inverse_node_map,
408  &prism_face_map, &prism_inverse_face_map);
409  convert_type(PYRAMID5, "PYRAMID5");
410  convert_type(PYRAMID13, "PYRAMID13");
411  convert_type(PYRAMID14, "PYRAMID14");
412  convert_type(PYRAMID18, "PYRAMID18");
413 }
414 
415 
416 
417 // This function initializes the element_equivalence_map the first time it
418 // is called, and returns early all other times.
420 {
421  // We use an ExodusII SPHERE element to represent a NodeElem
422  element_equivalence_map["SPHERE"] = NODEELEM;
423 
424  // EDGE2 equivalences
425  element_equivalence_map["EDGE"] = EDGE2;
426  element_equivalence_map["EDGE2"] = EDGE2;
427  element_equivalence_map["TRUSS"] = EDGE2;
428  element_equivalence_map["BEAM"] = EDGE2;
430  element_equivalence_map["TRUSS2"] = EDGE2;
431  element_equivalence_map["BEAM2"] = EDGE2;
432  element_equivalence_map["BAR2"] = EDGE2;
433 
434  // EDGE3 equivalences
435  element_equivalence_map["EDGE3"] = EDGE3;
436  element_equivalence_map["TRUSS3"] = EDGE3;
437  element_equivalence_map["BEAM3"] = EDGE3;
438  element_equivalence_map["BAR3"] = EDGE3;
439 
440  // EDGE4 equivalences
441  element_equivalence_map["EDGE4"] = EDGE4;
442  element_equivalence_map["TRUSS4"] = EDGE4;
443  element_equivalence_map["BEAM4"] = EDGE4;
444  element_equivalence_map["BAR4"] = EDGE4;
445 
446  // This whole design is going to need to be refactored whenever we
447  // support higher-order IGA, with one element type having variable
448  // polynomiaal degree...
449  element_equivalence_map["BEX_CURVE"] = EDGE3;
450 
451  // QUAD4 equivalences
452  element_equivalence_map["QUAD"] = QUAD4;
453  element_equivalence_map["QUAD4"] = QUAD4;
454 
455  // QUADSHELL4 equivalences
458 
459  // QUAD8 equivalences
460  element_equivalence_map["QUAD8"] = QUAD8;
461 
462  // QUADSHELL8 equivalences
464 
465  // QUAD9 equivalences
466  element_equivalence_map["QUAD9"] = QUAD9;
467  // This only supports p==2 IGA:
468  element_equivalence_map["BEX_QUAD"] = QUAD9;
469 
470  // QUADSHELL9 equivalences
472 
473  // TRI3 equivalences
474  element_equivalence_map["TRI"] = TRI3;
475  element_equivalence_map["TRI3"] = TRI3;
476  element_equivalence_map["TRIANGLE"] = TRI3;
477 
478  // TRISHELL3 equivalences
479  element_equivalence_map["TRISHELL"] = TRISHELL3;
480  element_equivalence_map["TRISHELL3"] = TRISHELL3;
481 
482  // TRI6 equivalences
483  element_equivalence_map["TRI6"] = TRI6;
484  // element_equivalence_map["TRISHELL6"] = TRI6;
485  // This only supports p==2 IGA:
486  element_equivalence_map["BEX_TRIANGLE"] = TRI6;
487 
488  // TRI7 equivalences
489  element_equivalence_map["TRI7"] = TRI7;
490 
491  // HEX8 equivalences
492  element_equivalence_map["HEX"] = HEX8;
493  element_equivalence_map["HEX8"] = HEX8;
494 
495  // HEX20 equivalences
496  element_equivalence_map["HEX20"] = HEX20;
497 
498  // HEX27 equivalences
499  element_equivalence_map["HEX27"] = HEX27;
500  // This only supports p==2 IGA:
501  element_equivalence_map["BEX_HEX"] = HEX27;
502 
503  // TET4 equivalences
504  element_equivalence_map["TETRA"] = TET4;
505  element_equivalence_map["TETRA4"] = TET4;
506 
507  // TET10 equivalences
508  element_equivalence_map["TETRA10"] = TET10;
509  // This only supports p==2 IGA:
510  element_equivalence_map["BEX_TETRA"] = TET10;
511 
512  // TET14 (in Exodus 8) equivalence
513  element_equivalence_map["TETRA14"] = TET14;
514 
515  // PRISM6 equivalences
516  element_equivalence_map["WEDGE"] = PRISM6;
517  element_equivalence_map["WEDGE6"] = PRISM6;
518 
519  // PRISM15 equivalences
520  element_equivalence_map["WEDGE15"] = PRISM15;
521 
522  // PRISM18 equivalences
523  element_equivalence_map["WEDGE18"] = PRISM18;
524  // This only supports p==2 IGA:
525  element_equivalence_map["BEX_WEDGE"] = PRISM18;
526 
527  // PRISM20 equivalences
528  element_equivalence_map["WEDGE20"] = PRISM20;
529 
530  // PRISM21 equivalences
531  element_equivalence_map["WEDGE21"] = PRISM21;
532 
533  // PYRAMID equivalences
534  element_equivalence_map["PYRAMID"] = PYRAMID5;
535  element_equivalence_map["PYRAMID5"] = PYRAMID5;
536  element_equivalence_map["PYRAMID13"] = PYRAMID13;
537  element_equivalence_map["PYRAMID14"] = PYRAMID14;
538  element_equivalence_map["PYRAMID18"] = PYRAMID18;
539 }
540 
543 {
544  auto & maps_for_dim = libmesh_map_find(conversion_map, this->num_dim);
545  return libmesh_map_find(maps_for_dim, type);
546 }
547 
549 ExodusII_IO_Helper::get_conversion(std::string type_str) const
550 {
551  // Do only upper-case comparisons
552  std::transform(type_str.begin(), type_str.end(), type_str.begin(), ::toupper);
553  return get_conversion (libmesh_map_find(element_equivalence_map, type_str));
554 }
555 
557 {
558  return elem_type.data();
559 }
560 
561 
562 
563 void ExodusII_IO_Helper::message(std::string_view msg)
564 {
565  if (verbose) libMesh::out << msg << std::endl;
566 }
567 
568 
569 
570 void ExodusII_IO_Helper::message(std::string_view msg, int i)
571 {
572  if (verbose) libMesh::out << msg << i << "." << std::endl;
573 }
574 
575 
577 MappedOutputVector(const std::vector<Real> & our_data_in,
578  bool single_precision_in)
579  : our_data(our_data_in),
580  single_precision(single_precision_in)
581 {
582  if (single_precision)
583  {
584  if (sizeof(Real) != sizeof(float))
585  {
586  float_vec.resize(our_data.size());
587  // boost float128 demands explicit downconversions
588  for (std::size_t i : index_range(our_data))
589  float_vec[i] = float(our_data[i]);
590  }
591  }
592 
593  else if (sizeof(Real) != sizeof(double))
594  {
595  double_vec.resize(our_data.size());
596  // boost float128 demands explicit downconversions
597  for (std::size_t i : index_range(our_data))
598  double_vec[i] = double(our_data[i]);
599  }
600 }
601 
602 void *
604 {
605  if (single_precision)
606  {
607  if (sizeof(Real) != sizeof(float))
608  return static_cast<void*>(float_vec.data());
609  }
610 
611  else if (sizeof(Real) != sizeof(double))
612  return static_cast<void*>(double_vec.data());
613 
614  // Otherwise return a (suitably casted) pointer to the original underlying data.
615  return const_cast<void *>(static_cast<const void *>(our_data.data()));
616 }
617 
619 MappedInputVector(std::vector<Real> & our_data_in,
620  bool single_precision_in)
621  : our_data(our_data_in),
622  single_precision(single_precision_in)
623 {
624  // Allocate temporary space to store enough floats/doubles, if required.
625  if (single_precision)
626  {
627  if (sizeof(Real) != sizeof(float))
628  float_vec.resize(our_data.size());
629  }
630  else if (sizeof(Real) != sizeof(double))
631  double_vec.resize(our_data.size());
632 }
633 
636 {
637  if (single_precision)
638  {
639  if (sizeof(Real) != sizeof(float))
640  our_data.assign(float_vec.begin(), float_vec.end());
641  }
642  else if (sizeof(Real) != sizeof(double))
643  our_data.assign(double_vec.begin(), double_vec.end());
644 }
645 
646 void *
648 {
649  if (single_precision)
650  {
651  if (sizeof(Real) != sizeof(float))
652  return static_cast<void*>(float_vec.data());
653  }
654 
655  else if (sizeof(Real) != sizeof(double))
656  return static_cast<void*>(double_vec.data());
657 
658  // Otherwise return a (suitably casted) pointer to the original underlying data.
659  return static_cast<void *>(our_data.data());
660 }
661 
662 void ExodusII_IO_Helper::open(const char * filename, bool read_only)
663 {
664  // Version of Exodus you are using
665  float ex_version = 0.;
666 
667  int comp_ws = 0;
668 
669  if (_single_precision)
670  comp_ws = cast_int<int>(sizeof(float));
671 
672  // Fall back on double precision when necessary since ExodusII
673  // doesn't seem to support long double
674  else
675  comp_ws = cast_int<int>(std::min(sizeof(Real), sizeof(double)));
676 
677  // Word size in bytes of the floating point data as they are stored
678  // in the ExodusII file. "If this argument is 0, the word size of the
679  // floating point data already stored in the file is returned"
680  int io_ws = 0;
681 
682  {
683  FPEDisabler disable_fpes;
684  ex_id = exII::ex_open(filename,
685  read_only ? EX_READ : EX_WRITE,
686  &comp_ws,
687  &io_ws,
688  &ex_version);
689  }
690 
691  std::string err_msg = std::string("Error opening ExodusII mesh file: ") + std::string(filename);
692  EX_CHECK_ERR(ex_id, err_msg);
693  if (verbose) libMesh::out << "File opened successfully." << std::endl;
694 
695  if (read_only)
696  opened_for_reading = true;
697  else
698  opened_for_writing = true;
699 
700  current_filename = std::string(filename);
701 }
702 
703 
704 
707 {
708  // Read init params using newer API that reads into a struct. For
709  // backwards compatibility, assign local member values from struct
710  // afterwards. Note: using the new API allows us to automatically
711  // read edge and face block/set information if it's present in the
712  // file.
713  exII::ex_init_params params = {};
714  int err_flag = exII::ex_get_init_ext(ex_id, &params);
715  EX_CHECK_ERR(err_flag, "Error retrieving header info.");
716 
717  // Extract required data into our struct
719  h.title.assign(params.title, params.title + MAX_LINE_LENGTH);
720  h.num_dim = params.num_dim;
721  h.num_nodes = params.num_nodes;
722  h.num_elem = params.num_elem;
723  h.num_elem_blk = params.num_elem_blk;
724  h.num_node_sets = params.num_node_sets;
725  h.num_side_sets = params.num_side_sets;
726  h.num_elem_sets = params.num_elem_sets;
727  h.num_edge_blk = params.num_edge_blk;
728  h.num_edge = params.num_edge;
729 
730  // And return it
731  return h;
732 }
733 
734 
735 
737 {
738  // Read header params from file, storing them in this class's
739  // ExodusHeaderInfo struct. This automatically updates the local
740  // num_dim, num_elem, etc. references.
741  this->header_info = this->read_header();
742 
743  // Read the number of timesteps which are present in the file
744  this->read_num_time_steps();
745 
746  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_NODAL, &num_nodal_vars);
747  EX_CHECK_ERR(ex_err, "Error reading number of nodal variables.");
748 
749  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_BLOCK, &num_elem_vars);
750  EX_CHECK_ERR(ex_err, "Error reading number of elemental variables.");
751 
752  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_GLOBAL, &num_global_vars);
753  EX_CHECK_ERR(ex_err, "Error reading number of global variables.");
754 
755  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_SIDE_SET, &num_sideset_vars);
756  EX_CHECK_ERR(ex_err, "Error reading number of sideset variables.");
757 
758  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_NODE_SET, &num_nodeset_vars);
759  EX_CHECK_ERR(ex_err, "Error reading number of nodeset variables.");
760 
761  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_SET, &num_elemset_vars);
762  EX_CHECK_ERR(ex_err, "Error reading number of elemset variables.");
763 
764  message("Exodus header info retrieved successfully.");
765 }
766 
767 
768 
769 
771 {
772  // The QA records are four MAX_STR_LENGTH-byte character strings.
773  int num_qa_rec =
774  inquire(*this, exII::EX_INQ_QA, "Error retrieving number of QA records");
775 
776  if (verbose)
777  libMesh::out << "Found "
778  << num_qa_rec
779  << " QA record(s) in the Exodus file."
780  << std::endl;
781 
782  if (num_qa_rec > 0)
783  {
784  // Actual (num_qa_rec x 4) storage for strings. The object we
785  // pass to the Exodus API will just contain pointers into the
786  // qa_storage object, which will have all automatic memory
787  // management.
788  std::vector<std::vector<std::vector<char>>> qa_storage(num_qa_rec);
789  for (auto i : make_range(num_qa_rec))
790  {
791  qa_storage[i].resize(4);
792  for (auto j : make_range(4))
793  qa_storage[i][j].resize(MAX_STR_LENGTH+1);
794  }
795 
796  // inner_array_t is a fixed-size array of 4 strings
797  typedef char * inner_array_t[4];
798 
799  // There is at least one compiler (Clang 12.0.1) that complains about
800  // "a non-scalar type used in a pseudo-destructor expression" when
801  // we try to instantiate a std::vector of inner_array_t objects as in:
802  // std::vector<inner_array_t> qa_record(num_qa_rec);
803  // So, we instead attempt to achieve the same effect with a std::unique_ptr.
804  auto qa_record = std::make_unique<inner_array_t[]>(num_qa_rec);
805 
806  // Create data structure to be passed to Exodus API by setting
807  // pointers to the actual strings which are in qa_storage.
808  for (auto i : make_range(num_qa_rec))
809  for (auto j : make_range(4))
810  qa_record[i][j] = qa_storage[i][j].data();
811 
812  ex_err = exII::ex_get_qa (ex_id, qa_record.get());
813  EX_CHECK_ERR(ex_err, "Error reading the QA records.");
814 
815  // Print the QA records
816  if (verbose)
817  {
818  for (auto i : make_range(num_qa_rec))
819  {
820  libMesh::out << "QA Record: " << i << std::endl;
821  for (auto j : make_range(4))
822  libMesh::out << qa_record[i][j] << std::endl;
823  }
824  }
825  }
826 }
827 
828 
829 
830 
832 {
833  if (verbose)
834  libMesh::out << "Title: \t" << title.data() << std::endl
835  << "Mesh Dimension: \t" << num_dim << std::endl
836  << "Number of Nodes: \t" << num_nodes << std::endl
837  << "Number of elements: \t" << num_elem << std::endl
838  << "Number of elt blocks: \t" << num_elem_blk << std::endl
839  << "Number of node sets: \t" << num_node_sets << std::endl
840  << "Number of side sets: \t" << num_side_sets << std::endl
841  << "Number of elem sets: \t" << num_elem_sets << std::endl;
842 }
843 
844 
845 
847 {
848  LOG_SCOPE("read_nodes()", "ExodusII_IO_Helper");
849 
850  x.resize(num_nodes);
851  y.resize(num_nodes);
852  z.resize(num_nodes);
853 
854  if (num_nodes)
855  {
856  ex_err = exII::ex_get_coord
857  (ex_id,
861 
862  EX_CHECK_ERR(ex_err, "Error retrieving nodal data.");
863  message("Nodal data retrieved successfully.");
864  }
865 
866  // If a nodal attribute bex_weight exists, we get spline weights
867  // from it
868  int n_nodal_attr = 0;
869  ex_err = exII::ex_get_attr_param(ex_id, exII::EX_NODAL, 0, & n_nodal_attr);
870  EX_CHECK_ERR(ex_err, "Error getting number of nodal attributes.");
871 
872  if (n_nodal_attr > 0)
873  {
874  std::vector<std::vector<char>> attr_name_data
875  (n_nodal_attr, std::vector<char>(MAX_STR_LENGTH + 1));
876  std::vector<char *> attr_names(n_nodal_attr);
877  for (auto i : index_range(attr_names))
878  attr_names[i] = attr_name_data[i].data();
879 
880  ex_err = exII::ex_get_attr_names(ex_id, exII::EX_NODAL, 0, attr_names.data());
881  EX_CHECK_ERR(ex_err, "Error getting nodal attribute names.");
882 
883  for (auto i : index_range(attr_names))
884  if (std::string("bex_weight") == attr_names[i])
885  {
886  w.resize(num_nodes);
887  ex_err =
888  exII::ex_get_one_attr (ex_id, exII::EX_NODAL, 0, i+1,
890  EX_CHECK_ERR(ex_err, "Error getting Bezier Extraction nodal weights");
891  }
892  }
893 }
894 
895 
896 
898 {
899  node_num_map.resize(num_nodes);
900 
901  // Note: we cannot use the exII::ex_get_num_map() here because it
902  // (apparently) does not behave like ex_get_node_num_map() when
903  // there is no node number map in the file: it throws an error
904  // instead of returning a default identity array (1,2,3,...).
905  ex_err = exII::ex_get_node_num_map
906  (ex_id, node_num_map.empty() ? nullptr : node_num_map.data());
907 
908  EX_CHECK_ERR(ex_err, "Error retrieving nodal number map.");
909  message("Nodal numbering map retrieved successfully.");
910 
911  if (verbose)
912  {
913  libMesh::out << "[" << this->processor_id() << "] node_num_map[i] = ";
914  for (unsigned int i=0; i<static_cast<unsigned int>(std::min(10, num_nodes-1)); ++i)
915  libMesh::out << node_num_map[i] << ", ";
916  libMesh::out << "... " << node_num_map.back() << std::endl;
917  }
918 }
919 
920 
922 {
923  // If a bex blob exists, we look for Bezier Extraction coefficient
924  // data there.
925 
926  // These APIs require newer Exodus than 5.22
927 #if EX_API_VERS_NODOT >= 800
928  int n_blobs = exII::ex_inquire_int(ex_id, exII::EX_INQ_BLOB);
929 
930  if (n_blobs > 0)
931  {
932  std::vector<exII::ex_blob> blobs(n_blobs);
933  std::vector<std::vector<char>> blob_names(n_blobs);
934  for (auto i : make_range(n_blobs))
935  {
936  blob_names[i].resize(MAX_STR_LENGTH+1);
937  blobs[i].name = blob_names[i].data();
938  }
939 
940  ex_err = exII::ex_get_blobs(ex_id, blobs.data());
941  EX_CHECK_ERR(ex_err, "Error getting blobs.");
942 
943  bool found_blob = false;
944  const exII::ex_blob * my_blob = &blobs[0];
945  for (const auto & blob : blobs)
946  {
947  if (std::string("bex_cv_blob") == blob.name)
948  {
949  found_blob = true;
950  my_blob = &blob;
951  }
952  }
953 
954  if (!found_blob)
955  libmesh_error_msg("Found no bex_cv_blob for bezier elements");
956 
957  const int n_blob_attr =
958  exII::ex_get_attribute_count(ex_id, exII::EX_BLOB,
959  my_blob->id);
960 
961  std::vector<exII::ex_attribute> attributes(n_blob_attr);
962  ex_err = exII::ex_get_attribute_param(ex_id, exII::EX_BLOB,
963  my_blob->id,
964  attributes.data());
965  EX_CHECK_ERR(ex_err, "Error getting bex blob attribute parameters.");
966 
967  int bex_num_dense_cv_blocks = 0;
968  std::vector<int> bex_dense_cv_info;
969  for (auto & attr : attributes)
970  {
971  if (std::string("bex_dense_cv_info") == attr.name)
972  {
973  const std::size_t value_count = attr.value_count;
974  if (value_count % 2)
975  libmesh_error_msg("Found odd number of bex_dense_cv_info");
976 
977  bex_dense_cv_info.resize(value_count);
978  attr.values = bex_dense_cv_info.data();
979  exII::ex_get_attribute(ex_id, &attr);
980 
981  bex_num_dense_cv_blocks = value_count / 2;
982 
983  libmesh_error_msg_if(bex_num_dense_cv_blocks > 1,
984  "Found more than 1 dense bex CV block; unsure how to handle that");
985  }
986  }
987 
988  if (bex_dense_cv_info.empty())
989  libmesh_error_msg("No bex_dense_cv_info found");
990 
991  int n_blob_vars;
992  exII::ex_get_variable_param(ex_id, exII::EX_BLOB, &n_blob_vars);
993  std::vector<char> var_name (MAX_STR_LENGTH + 1);
994  for (auto v_id : make_range(1,n_blob_vars+1))
995  {
996  ex_err = exII::ex_get_variable_name(ex_id, exII::EX_BLOB, v_id, var_name.data());
997  EX_CHECK_ERR(ex_err, "Error reading bex blob var name.");
998 
999  if (std::string("bex_dense_cv_blocks") == var_name.data())
1000  {
1001  std::vector<double> bex_dense_cv_blocks(my_blob->num_entry);
1002 
1003  ex_err = exII::ex_get_var(ex_id, 1, exII::EX_BLOB, v_id,
1004  my_blob->id, my_blob->num_entry,
1005  bex_dense_cv_blocks.data());
1006  EX_CHECK_ERR(ex_err, "Error reading bex_dense_cv_blocks.");
1007 
1008  bex_dense_constraint_vecs.clear();
1009  bex_dense_constraint_vecs.resize(bex_num_dense_cv_blocks);
1010 
1011  std::size_t offset = 0;
1012  for (auto i : IntRange<std::size_t>(0, bex_num_dense_cv_blocks))
1013  {
1014  bex_dense_constraint_vecs[i].resize(bex_dense_cv_info[2*i]);
1015  const int vecsize = bex_dense_cv_info[2*i+1];
1016  for (auto & vec : bex_dense_constraint_vecs[i])
1017  {
1018  vec.resize(vecsize);
1019  std::copy(std::next(bex_dense_cv_blocks.begin(), offset),
1020  std::next(bex_dense_cv_blocks.begin(), offset + vecsize),
1021  vec.begin());
1022  offset += vecsize;
1023  }
1024  }
1025  libmesh_assert(offset == bex_dense_cv_blocks.size());
1026  }
1027  }
1028  }
1029 #endif // EX_API_VERS_NODOT >= 800
1030 }
1031 
1032 
1033 void ExodusII_IO_Helper::print_nodes(std::ostream & out_stream)
1034 {
1035  for (int i=0; i<num_nodes; i++)
1036  out_stream << "(" << x[i] << ", " << y[i] << ", " << z[i] << ")" << std::endl;
1037 }
1038 
1039 
1040 
1042 {
1043  if (num_elem_blk)
1044  {
1045  // Read all element block IDs.
1046  block_ids.resize(num_elem_blk);
1047  ex_err = exII::ex_get_ids(ex_id,
1048  exII::EX_ELEM_BLOCK,
1049  block_ids.data());
1050 
1051  EX_CHECK_ERR(ex_err, "Error getting block IDs.");
1052  message("All block IDs retrieved successfully.");
1053 
1054  char name_buffer[MAX_STR_LENGTH+1];
1055  for (int i=0; i<num_elem_blk; ++i)
1056  {
1057  ex_err = exII::ex_get_name(ex_id, exII::EX_ELEM_BLOCK,
1058  block_ids[i], name_buffer);
1059  EX_CHECK_ERR(ex_err, "Error getting block name.");
1060  id_to_block_names[block_ids[i]] = name_buffer;
1061  }
1062  message("All block names retrieved successfully.");
1063  }
1064 
1065  if (num_edge_blk)
1066  {
1067  // Read all edge block IDs.
1068  edge_block_ids.resize(num_edge_blk);
1069  ex_err = exII::ex_get_ids(ex_id,
1070  exII::EX_EDGE_BLOCK,
1071  edge_block_ids.data());
1072 
1073  EX_CHECK_ERR(ex_err, "Error getting edge block IDs.");
1074  message("All edge block IDs retrieved successfully.");
1075 
1076  // Read in edge block names
1077  char name_buffer[MAX_STR_LENGTH+1];
1078  for (int i=0; i<num_edge_blk; ++i)
1079  {
1080  ex_err = exII::ex_get_name(ex_id, exII::EX_EDGE_BLOCK,
1081  edge_block_ids[i], name_buffer);
1082  EX_CHECK_ERR(ex_err, "Error getting block name.");
1083  id_to_edge_block_names[edge_block_ids[i]] = name_buffer;
1084  }
1085  message("All edge block names retrieved successfully.");
1086  }
1087 }
1088 
1089 
1090 
1092 {
1093  libmesh_assert_less (index, block_ids.size());
1094 
1095  return block_ids[index];
1096 }
1097 
1098 
1099 
1101 {
1102  libmesh_assert_less (index, block_ids.size());
1103 
1104  return id_to_block_names[block_ids[index]];
1105 }
1106 
1107 
1108 
1110 {
1111  libmesh_assert_less (index, ss_ids.size());
1112 
1113  return ss_ids[index];
1114 }
1115 
1116 
1117 
1119 {
1120  libmesh_assert_less (index, ss_ids.size());
1121 
1122  return id_to_ss_names[ss_ids[index]];
1123 }
1124 
1125 
1126 
1128 {
1129  libmesh_assert_less (index, nodeset_ids.size());
1130 
1131  return nodeset_ids[index];
1132 }
1133 
1134 
1135 
1137 {
1138  libmesh_assert_less (index, nodeset_ids.size());
1139 
1140  return id_to_ns_names[nodeset_ids[index]];
1141 }
1142 
1143 
1144 
1145 
1147 {
1148  LOG_SCOPE("read_elem_in_block()", "ExodusII_IO_Helper");
1149 
1150  libmesh_assert_less (block, block_ids.size());
1151 
1152  // Unlike the other "extended" APIs, this one does not use a parameter struct.
1153  int num_edges_per_elem = 0;
1154  int num_faces_per_elem = 0;
1155  int num_node_data_per_elem = 0;
1156  ex_err = exII::ex_get_block(ex_id,
1157  exII::EX_ELEM_BLOCK,
1158  block_ids[block],
1159  elem_type.data(),
1161  &num_node_data_per_elem,
1162  &num_edges_per_elem, // 0 or -1 if no "extended" block info
1163  &num_faces_per_elem, // 0 or -1 if no "extended" block info
1164  &num_attr);
1165 
1166  EX_CHECK_ERR(ex_err, "Error getting block info.");
1167  message("Info retrieved successfully for block: ", block);
1168 
1169  // Warn that we don't currently support reading blocks with extended info.
1170  // Note: the docs say -1 will be returned for this but I found that it was
1171  // actually 0, so not sure which it will be in general.
1172  if (!(num_edges_per_elem == 0) && !(num_edges_per_elem == -1))
1173  libmesh_warning("Exodus files with extended edge connectivity not currently supported.");
1174  if (!(num_faces_per_elem == 0) && !(num_faces_per_elem == -1))
1175  libmesh_warning("Exodus files with extended face connectivity not currently supported.");
1176 
1177  // If we have a Bezier element here, then we've packed constraint
1178  // vector connectivity at the end of the nodal connectivity, and
1179  // num_nodes_per_elem reflected both.
1180  const bool is_bezier = is_bezier_elem(elem_type.data());
1181  if (is_bezier)
1182  {
1183  const auto & conv = get_conversion(std::string(elem_type.data()));
1184  num_nodes_per_elem = conv.n_nodes;
1185  }
1186  else
1187  num_nodes_per_elem = num_node_data_per_elem;
1188 
1189  if (verbose)
1190  libMesh::out << "Read a block of " << num_elem_this_blk
1191  << " " << elem_type.data() << "(s)"
1192  << " having " << num_nodes_per_elem
1193  << " nodes per element." << std::endl;
1194 
1195  // Read in the connectivity of the elements of this block,
1196  // watching out for the case where we actually have no
1197  // elements in this block (possible with parallel files)
1198  connect.resize(num_node_data_per_elem*num_elem_this_blk);
1199 
1200  if (!connect.empty())
1201  {
1202  ex_err = exII::ex_get_conn(ex_id,
1203  exII::EX_ELEM_BLOCK,
1204  block_ids[block],
1205  connect.data(), // node_conn
1206  nullptr, // elem_edge_conn (unused)
1207  nullptr); // elem_face_conn (unused)
1208 
1209  EX_CHECK_ERR(ex_err, "Error reading block connectivity.");
1210  message("Connectivity retrieved successfully for block: ", block);
1211  }
1212 
1213  // If we had any attributes for this block, check to see if some of
1214  // them were Bezier-extension attributes.
1215 
1216  // num_attr above is zero, not actually the number of block attributes?
1217  // ex_get_attr_param *also* gives me zero? Really, Exodus?
1218 #if EX_API_VERS_NODOT >= 800
1219  int real_n_attr = exII::ex_get_attribute_count(ex_id, exII::EX_ELEM_BLOCK, block_ids[block]);
1220  EX_CHECK_ERR(real_n_attr, "Error getting number of element block attributes.");
1221 
1222  if (real_n_attr > 0)
1223  {
1224  std::vector<exII::ex_attribute> attributes(real_n_attr);
1225 
1226  ex_err = exII::ex_get_attribute_param(ex_id, exII::EX_ELEM_BLOCK, block_ids[block], attributes.data());
1227  EX_CHECK_ERR(ex_err, "Error getting element block attribute parameters.");
1228 
1229  ex_err = exII::ex_get_attributes(ex_id, real_n_attr, attributes.data());
1230  EX_CHECK_ERR(ex_err, "Error getting element block attribute values.");
1231 
1232  for (auto attr : attributes)
1233  {
1234  if (std::string("bex_elem_degrees") == attr.name)
1235  {
1236  if (attr.type != exII::EX_INTEGER)
1237  libmesh_error_msg("Found non-integer bex_elem_degrees");
1238 
1239  if (attr.value_count > 3)
1240  libmesh_error_msg("Looking for at most 3 bex_elem_degrees; found " << attr.value_count);
1241 
1242  libmesh_assert(is_bezier);
1243 
1244  std::vector<int> bex_elem_degrees(3); // max dim
1245 
1246  const int * as_int = static_cast<int *>(attr.values);
1247  std::copy(as_int, as_int+attr.value_count, bex_elem_degrees.begin());
1248 
1249 
1250  // Right now Bezier extraction elements aren't possible
1251  // for p>2 and aren't useful for p<2, and we don't
1252  // support anisotropic p...
1253 #ifndef NDEBUG
1254  const auto & conv = get_conversion(std::string(elem_type.data()));
1255 
1256  for (auto d : IntRange<int>(0, conv.dim))
1257  libmesh_assert_equal_to(bex_elem_degrees[d], 2);
1258 #endif
1259  }
1260  // ex_get_attributes did a values=calloc(); free() is our job.
1261  if (attr.values)
1262  free(attr.values);
1263  }
1264  }
1265 
1266  if (is_bezier)
1267  {
1268  // We'd better have the number of cvs we expect
1269  if( num_node_data_per_elem > num_nodes_per_elem )
1270  bex_num_elem_cvs = num_node_data_per_elem / 2;
1271  else
1273  libmesh_assert_greater_equal(bex_num_elem_cvs, 0);
1274 
1275  // The old connect vector is currently a mix of the expected
1276  // connectivity and any Bezier extraction connectivity;
1277  // disentangle that, if necessary.
1279  if (num_node_data_per_elem > num_nodes_per_elem)
1280  {
1281  std::vector<int> old_connect(bex_num_elem_cvs * num_elem_this_blk);
1282  old_connect.swap(connect);
1283  auto src = old_connect.data();
1284  auto dst = connect.data();
1285  for (auto e : IntRange<std::size_t>(0, num_elem_this_blk))
1286  {
1287  std::copy(src, src + bex_num_elem_cvs, dst);
1288  src += bex_num_elem_cvs;
1289  dst += bex_num_elem_cvs;
1290 
1291  bex_cv_conn[e].resize(bex_num_elem_cvs);
1292  std::copy(src, src + bex_num_elem_cvs,
1293  bex_cv_conn[e].begin());
1294  src += bex_num_elem_cvs;
1295  }
1296  }
1297  }
1298 
1299 #endif // EX_API_VERS_NODOT >= 800
1300 }
1301 
1302 
1303 
1305 {
1306  LOG_SCOPE("read_edge_blocks()", "ExodusII_IO_Helper");
1307 
1308  // Check for quick return if there are no edge blocks.
1309  if (num_edge_blk == 0)
1310  return;
1311 
1312  // Build data structure that we can quickly search for edges
1313  // and then add required BoundaryInfo information. This is a
1314  // map from edge->key() to a list of (elem_id, edge_id) pairs
1315  // for the Edge in question. Since edge->key() is edge orientation
1316  // invariant, this map does not distinguish different orientations
1317  // of the same Edge. Since edge->key() is also not guaranteed to be
1318  // unique (though it is very unlikely for two distinct edges to have
1319  // the same key()), when we later look up an (elem_id, edge_id) pair
1320  // in the edge_map, we need to verify that the edge indeed matches
1321  // the searched edge by doing some further checks.
1322  typedef std::pair<dof_id_type, unsigned int> ElemEdgePair;
1323  std::unordered_map<dof_id_type, std::vector<ElemEdgePair>> edge_map;
1324  std::unique_ptr<Elem> edge_ptr;
1325  for (const auto & elem : mesh.element_ptr_range())
1326  for (auto e : elem->edge_index_range())
1327  {
1328  elem->build_edge_ptr(edge_ptr, e);
1329  dof_id_type edge_key = edge_ptr->key();
1330 
1331  // Creates vector if not already there
1332  auto & vec = edge_map[edge_key];
1333  vec.emplace_back(elem->id(), e);
1334 
1335  // If edge_ptr is a higher-order Elem (EDGE3 or higher) then also add
1336  // a map entry for the lower-order (EDGE2) element which has matching
1337  // vertices. This allows us to match lower-order edge blocks to edges
1338  // of higher-order 3D elems (e.g. HEX20, TET10) and simplifies the
1339  // definition of edge blocks.
1340  if (edge_ptr->default_order() != FIRST)
1341  {
1342  // Construct a temporary low-order edge so that we can compute its key()
1343  auto low_order_edge =
1345 
1346  // Assign node pointers to low-order edge
1347  for (unsigned int v=0; v<edge_ptr->n_vertices(); ++v)
1348  low_order_edge->set_node(v, edge_ptr->node_ptr(v));
1349 
1350  // Compute the key for the temporary low-order edge we just built
1351  dof_id_type low_order_edge_key = low_order_edge->key();
1352 
1353  // Add this key to the map associated with the same (elem,
1354  // edge) pair as the higher-order edge
1355  auto & low_order_vec = edge_map[low_order_edge_key];
1356  low_order_vec.emplace_back(elem->id(), e);
1357  }
1358  }
1359 
1360  // Get reference to the mesh's BoundaryInfo object, as we will be
1361  // adding edges to this below.
1363 
1364  for (const auto & edge_block_id : edge_block_ids)
1365  {
1366  // exII::ex_get_block() output parameters. Unlike the other
1367  // "extended" APIs, exII::ex_get_block() does not use a
1368  // parameter struct.
1369  int num_edge_this_blk = 0;
1370  int num_nodes_per_edge = 0;
1371  int num_edges_per_edge = 0;
1372  int num_faces_per_edge = 0;
1373  int num_attr_per_edge = 0;
1374  ex_err = exII::ex_get_block(ex_id,
1375  exII::EX_EDGE_BLOCK,
1376  edge_block_id,
1377  elem_type.data(),
1378  &num_edge_this_blk,
1379  &num_nodes_per_edge,
1380  &num_edges_per_edge, // 0 or -1 for edge blocks
1381  &num_faces_per_edge, // 0 or -1 for edge blocks
1382  &num_attr_per_edge);
1383 
1384  EX_CHECK_ERR(ex_err, "Error getting edge block info.");
1385  message("Info retrieved successfully for block: ", edge_block_id);
1386 
1387  // Read in the connectivity of the edges of this block,
1388  // watching out for the case where we actually have no
1389  // elements in this block (possible with parallel files)
1390  connect.resize(num_nodes_per_edge * num_edge_this_blk);
1391 
1392  if (!connect.empty())
1393  {
1394  ex_err = exII::ex_get_conn(ex_id,
1395  exII::EX_EDGE_BLOCK,
1396  edge_block_id,
1397  connect.data(), // node_conn
1398  nullptr, // elem_edge_conn (unused)
1399  nullptr); // elem_face_conn (unused)
1400 
1401  EX_CHECK_ERR(ex_err, "Error reading block connectivity.");
1402  message("Connectivity retrieved successfully for block: ", edge_block_id);
1403 
1404  // All edge types have an identity mapping from the corresponding
1405  // Exodus type, so we don't need to bother with mapping ids, but
1406  // we do need to know what kind of elements to build.
1407  const auto & conv = get_conversion(std::string(elem_type.data()));
1408 
1409  // Loop over indices in connectivity array, build edge elements,
1410  // look them up in the edge_map.
1411  for (unsigned int i=0, sz=connect.size(); i<sz; i+=num_nodes_per_edge)
1412  {
1413  auto edge = Elem::build(conv.libmesh_elem_type());
1414  for (int n=0; n<num_nodes_per_edge; ++n)
1415  {
1416  int exodus_node_id = connect[i+n];
1417  int exodus_node_id_zero_based = exodus_node_id - 1;
1418  int libmesh_node_id = node_num_map[exodus_node_id_zero_based] - 1;
1419 
1420  edge->set_node(n, mesh.node_ptr(libmesh_node_id));
1421  }
1422 
1423  // Compute key for the edge Elem we just built.
1424  dof_id_type edge_key = edge->key();
1425 
1426  // If this key is not found in the edge_map, which is
1427  // supposed to include every edge in the Mesh, then we
1428  // will throw an error now.
1429  auto & elem_edge_pair_vec =
1430  libmesh_map_find(edge_map, edge_key);
1431 
1432  for (const auto & elem_edge_pair : elem_edge_pair_vec)
1433  {
1434  // We only want to match edges which have the same
1435  // nodes (possibly with different orientation) to the one in the
1436  // Exodus file, otherwise we ignore this elem_edge_pair.
1437  //
1438  // Note: this also handles the situation where two
1439  // edges have the same key (hash collision) as then
1440  // this check avoids a false positive.
1441 
1442  // Build edge indicated by elem_edge_pair
1443  mesh.elem_ptr(elem_edge_pair.first)->
1444  build_edge_ptr(edge_ptr, elem_edge_pair.second);
1445 
1446  // Determine whether this candidate edge is a "real" match,
1447  // i.e. has the same nodes with a possibly different
1448  // orientation. Note that here we only check that
1449  // the vertices match regardless of how many nodes
1450  // the edge has, which allows us to match a
1451  // lower-order edge to a higher-order Elem.
1452  bool is_match =
1453  ((edge_ptr->node_id(0) == edge->node_id(0)) && (edge_ptr->node_id(1) == edge->node_id(1))) ||
1454  ((edge_ptr->node_id(0) == edge->node_id(1)) && (edge_ptr->node_id(1) == edge->node_id(0)));
1455 
1456  if (is_match)
1457  {
1458  // Add this (elem, edge, id) combo to the BoundaryInfo object.
1459  bi.add_edge(elem_edge_pair.first,
1460  elem_edge_pair.second,
1461  edge_block_id);
1462  }
1463  } // end loop over elem_edge_pairs
1464  } // end loop over connectivity array
1465 
1466  // Set edgeset name in the BoundaryInfo object.
1467  bi.edgeset_name(edge_block_id) = id_to_edge_block_names[edge_block_id];
1468  } // end if !connect.empty()
1469  } // end for edge_block_id : edge_block_ids
1470 }
1471 
1472 
1473 
1475 {
1476  elem_num_map.resize(num_elem);
1477 
1478  // Note: we cannot use the exII::ex_get_num_map() here because it
1479  // (apparently) does not behave like ex_get_elem_num_map() when
1480  // there is no elem number map in the file: it throws an error
1481  // instead of returning a default identity array (1,2,3,...).
1482  ex_err = exII::ex_get_elem_num_map
1483  (ex_id, elem_num_map.empty() ? nullptr : elem_num_map.data());
1484 
1485  EX_CHECK_ERR(ex_err, "Error retrieving element number map.");
1486  message("Element numbering map retrieved successfully.");
1487 
1488  if (num_elem)
1489  {
1490  auto it = std::max_element(elem_num_map.begin(), elem_num_map.end());
1491  _end_elem_id = *it;
1492  }
1493  else
1494  _end_elem_id = 0;
1495 
1496  if (verbose)
1497  {
1498  libMesh::out << "[" << this->processor_id() << "] elem_num_map[i] = ";
1499  for (unsigned int i=0; i<static_cast<unsigned int>(std::min(10, num_elem-1)); ++i)
1500  libMesh::out << elem_num_map[i] << ", ";
1501  libMesh::out << "... " << elem_num_map.back() << std::endl;
1502  }
1503 }
1504 
1505 
1506 
1508 {
1509  ss_ids.resize(num_side_sets);
1510  if (num_side_sets > 0)
1511  {
1512  ex_err = exII::ex_get_ids(ex_id,
1513  exII::EX_SIDE_SET,
1514  ss_ids.data());
1515  EX_CHECK_ERR(ex_err, "Error retrieving sideset information.");
1516  message("All sideset information retrieved successfully.");
1517 
1518  // Resize appropriate data structures -- only do this once outside the loop
1520  num_df_per_set.resize(num_side_sets);
1521 
1522  // Inquire about the length of the concatenated side sets element list
1523  num_elem_all_sidesets = inquire(*this, exII::EX_INQ_SS_ELEM_LEN, "Error retrieving length of the concatenated side sets element list!");
1524 
1527  id_list.resize (num_elem_all_sidesets);
1528  }
1529 
1530  char name_buffer[MAX_STR_LENGTH+1];
1531  for (int i=0; i<num_side_sets; ++i)
1532  {
1533  ex_err = exII::ex_get_name(ex_id, exII::EX_SIDE_SET,
1534  ss_ids[i], name_buffer);
1535  EX_CHECK_ERR(ex_err, "Error getting side set name.");
1536  id_to_ss_names[ss_ids[i]] = name_buffer;
1537  }
1538  message("All side set names retrieved successfully.");
1539 }
1540 
1541 
1543 {
1544  nodeset_ids.resize(num_node_sets);
1545  if (num_node_sets > 0)
1546  {
1547  ex_err = exII::ex_get_ids(ex_id,
1548  exII::EX_NODE_SET,
1549  nodeset_ids.data());
1550  EX_CHECK_ERR(ex_err, "Error retrieving nodeset information.");
1551  message("All nodeset information retrieved successfully.");
1552 
1553  // Resize appropriate data structures -- only do this once outside the loop
1556  }
1557 
1558  char name_buffer[MAX_STR_LENGTH+1];
1559  for (int i=0; i<num_node_sets; ++i)
1560  {
1561  ex_err = exII::ex_get_name(ex_id, exII::EX_NODE_SET,
1562  nodeset_ids[i], name_buffer);
1563  EX_CHECK_ERR(ex_err, "Error getting node set name.");
1564  id_to_ns_names[nodeset_ids[i]] = name_buffer;
1565  }
1566  message("All node set names retrieved successfully.");
1567 }
1568 
1569 
1570 
1572 {
1573  elemset_ids.resize(num_elem_sets);
1574  if (num_elem_sets > 0)
1575  {
1576  ex_err = exII::ex_get_ids(ex_id,
1577  exII::EX_ELEM_SET,
1578  elemset_ids.data());
1579  EX_CHECK_ERR(ex_err, "Error retrieving elemset information.");
1580  message("All elemset information retrieved successfully.");
1581 
1582  // Resize appropriate data structures -- only do this once outside the loop
1585 
1586  // Inquire about the length of the concatenated elemset list
1588  inquire(*this, exII::EX_INQ_ELS_LEN,
1589  "Error retrieving length of the concatenated elem sets element list!");
1590 
1593 
1594  // Debugging
1595  // libMesh::out << "num_elem_all_elemsets = " << num_elem_all_elemsets << std::endl;
1596  }
1597 
1598  char name_buffer[MAX_STR_LENGTH+1];
1599  for (int i=0; i<num_elem_sets; ++i)
1600  {
1601  ex_err = exII::ex_get_name(ex_id, exII::EX_ELEM_SET,
1602  elemset_ids[i], name_buffer);
1603  EX_CHECK_ERR(ex_err, "Error getting node set name.");
1604  id_to_elemset_names[elemset_ids[i]] = name_buffer;
1605  }
1606  message("All elem set names retrieved successfully.");
1607 }
1608 
1609 
1610 
1611 void ExodusII_IO_Helper::read_sideset(int id, int offset)
1612 {
1613  LOG_SCOPE("read_sideset()", "ExodusII_IO_Helper");
1614 
1615  libmesh_assert_less (id, ss_ids.size());
1616  libmesh_assert_less (id, num_sides_per_set.size());
1617  libmesh_assert_less (id, num_df_per_set.size());
1618  libmesh_assert_less_equal (offset, elem_list.size());
1619  libmesh_assert_less_equal (offset, side_list.size());
1620 
1621  ex_err = exII::ex_get_set_param(ex_id,
1622  exII::EX_SIDE_SET,
1623  ss_ids[id],
1624  &num_sides_per_set[id],
1625  &num_df_per_set[id]);
1626  EX_CHECK_ERR(ex_err, "Error retrieving sideset parameters.");
1627  message("Parameters retrieved successfully for sideset: ", id);
1628 
1629 
1630  // It's OK for offset==elem_list.size() as long as num_sides_per_set[id]==0
1631  // because in that case we don't actually read anything...
1632 #ifdef DEBUG
1633  if (static_cast<unsigned int>(offset) == elem_list.size() ||
1634  static_cast<unsigned int>(offset) == side_list.size() )
1635  libmesh_assert_equal_to (num_sides_per_set[id], 0);
1636 #endif
1637 
1638 
1639  // Don't call ex_get_set unless there are actually sides there to get.
1640  // Exodus prints an annoying warning in DEBUG mode otherwise...
1641  if (num_sides_per_set[id] > 0)
1642  {
1643  ex_err = exII::ex_get_set(ex_id,
1644  exII::EX_SIDE_SET,
1645  ss_ids[id],
1646  &elem_list[offset],
1647  &side_list[offset]);
1648  EX_CHECK_ERR(ex_err, "Error retrieving sideset data.");
1649  message("Data retrieved successfully for sideset: ", id);
1650 
1651  for (int i=0; i<num_sides_per_set[id]; i++)
1652  id_list[i+offset] = ss_ids[id];
1653  }
1654 }
1655 
1656 
1657 
1658 void ExodusII_IO_Helper::read_elemset(int id, int offset)
1659 {
1660  LOG_SCOPE("read_elemset()", "ExodusII_IO_Helper");
1661 
1662  libmesh_assert_less (id, elemset_ids.size());
1663  libmesh_assert_less (id, num_elems_per_set.size());
1664  libmesh_assert_less (id, num_elem_df_per_set.size());
1665  libmesh_assert_less_equal (offset, elemset_list.size());
1666 
1667  ex_err = exII::ex_get_set_param(ex_id,
1668  exII::EX_ELEM_SET,
1669  elemset_ids[id],
1670  &num_elems_per_set[id],
1671  &num_elem_df_per_set[id]);
1672  EX_CHECK_ERR(ex_err, "Error retrieving elemset parameters.");
1673  message("Parameters retrieved successfully for elemset: ", id);
1674 
1675 
1676  // It's OK for offset==elemset_list.size() as long as num_elems_per_set[id]==0
1677  // because in that case we don't actually read anything...
1678  #ifdef DEBUG
1679  if (static_cast<unsigned int>(offset) == elemset_list.size())
1680  libmesh_assert_equal_to (num_elems_per_set[id], 0);
1681  #endif
1682 
1683  // Don't call ex_get_set() unless there are actually elems there to get.
1684  // Exodus prints an annoying warning in DEBUG mode otherwise...
1685  if (num_elems_per_set[id] > 0)
1686  {
1687  ex_err = exII::ex_get_set(ex_id,
1688  exII::EX_ELEM_SET,
1689  elemset_ids[id],
1690  &elemset_list[offset],
1691  /*set_extra_list=*/nullptr);
1692  EX_CHECK_ERR(ex_err, "Error retrieving elemset data.");
1693  message("Data retrieved successfully for elemset: ", id);
1694 
1695  // Create vector containing elemset ids for each element in the set
1696  for (int i=0; i<num_elems_per_set[id]; i++)
1697  elemset_id_list[i+offset] = elemset_ids[id];
1698  }
1699 }
1700 
1701 
1702 
1704 {
1705  LOG_SCOPE("read_all_nodesets()", "ExodusII_IO_Helper");
1706 
1707  // Figure out how many nodesets there are in the file so we can
1708  // properly resize storage as necessary.
1709  num_node_sets =
1710  inquire
1711  (*this, exII::EX_INQ_NODE_SETS,
1712  "Error retrieving number of node sets");
1713 
1714  // Figure out how many nodes there are in all the nodesets.
1715  int total_nodes_in_all_sets =
1716  inquire
1717  (*this, exII::EX_INQ_NS_NODE_LEN,
1718  "Error retrieving number of nodes in all node sets.");
1719 
1720  // Figure out how many distribution factors there are in all the nodesets.
1721  int total_df_in_all_sets =
1722  inquire
1723  (*this, exII::EX_INQ_NS_DF_LEN,
1724  "Error retrieving number of distribution factors in all node sets.");
1725 
1726  // If there are no nodesets, there's nothing to read in.
1727  if (num_node_sets == 0)
1728  return;
1729 
1730  // Allocate space to read all the nodeset data.
1731  // Use existing class members where possible to avoid shadowing
1732  nodeset_ids.clear(); nodeset_ids.resize(num_node_sets);
1737  node_sets_node_list.clear(); node_sets_node_list.resize(total_nodes_in_all_sets);
1738  node_sets_dist_fact.clear(); node_sets_dist_fact.resize(total_df_in_all_sets);
1739 
1740  // Handle single-precision files
1741  MappedInputVector mapped_node_sets_dist_fact(node_sets_dist_fact, _single_precision);
1742 
1743  // Build exII::ex_set_spec struct
1744  exII::ex_set_specs set_specs = {};
1745  set_specs.sets_ids = nodeset_ids.data();
1746  set_specs.num_entries_per_set = num_nodes_per_set.data();
1747  set_specs.num_dist_per_set = num_node_df_per_set.data();
1748  set_specs.sets_entry_index = node_sets_node_index.data();
1749  set_specs.sets_dist_index = node_sets_dist_index.data();
1750  set_specs.sets_entry_list = node_sets_node_list.data();
1751  set_specs.sets_extra_list = nullptr;
1752  set_specs.sets_dist_fact = total_df_in_all_sets ? mapped_node_sets_dist_fact.data() : nullptr;
1753 
1754  ex_err = exII::ex_get_concat_sets(ex_id, exII::EX_NODE_SET, &set_specs);
1755  EX_CHECK_ERR(ex_err, "Error reading concatenated nodesets");
1756 
1757  // Read the nodeset names from file!
1758  char name_buffer[MAX_STR_LENGTH+1];
1759  for (int i=0; i<num_node_sets; ++i)
1760  {
1761  ex_err = exII::ex_get_name
1762  (ex_id,
1763  exII::EX_NODE_SET,
1764  nodeset_ids[i],
1765  name_buffer);
1766  EX_CHECK_ERR(ex_err, "Error getting node set name.");
1767  id_to_ns_names[nodeset_ids[i]] = name_buffer;
1768  }
1769 }
1770 
1771 
1772 
1774 {
1775  // Call ex_close on every processor that did ex_open or ex_create;
1776  // newer Exodus versions error if we try to reopen a file that
1777  // hasn't been officially closed. Don't close the file if we didn't
1778  // open it; this also raises an Exodus error.
1779 
1780  // We currently do read-only ex_open on every proc (to do read
1781  // operations on every proc), but we do ex_open and ex_create for
1782  // writes on every proc only with Nemesis files.
1784  (this->processor_id() == 0) ||
1785  (!_run_only_on_proc0))
1786  {
1788  {
1789  ex_err = exII::ex_close(ex_id);
1790  // close() is called from the destructor, so it may be called e.g.
1791  // during stack unwinding while processing an exception. In that case
1792  // we don't want to throw another exception or immediately terminate
1793  // the code, since that would prevent any possible recovery from the
1794  // exception in question. So we just log the error closing the file
1795  // and continue.
1796  if (ex_err < 0)
1797  message("Error closing Exodus file.");
1798  else
1799  message("Exodus file closed successfully.");
1800  }
1801  }
1802 
1803  // Now that the file is closed, it's no longer opened for
1804  // reading or writing.
1805  opened_for_writing = false;
1806  opened_for_reading = false;
1807  _opened_by_create = false;
1808 }
1809 
1810 
1811 
1813 {
1814  // Make sure we have an up-to-date count of the number of time steps in the file.
1815  this->read_num_time_steps();
1816 
1817  if (num_time_steps > 0)
1818  {
1819  time_steps.resize(num_time_steps);
1820  ex_err = exII::ex_get_all_times
1821  (ex_id,
1823  EX_CHECK_ERR(ex_err, "Error reading timesteps!");
1824  }
1825 }
1826 
1827 
1828 
1830 {
1831  num_time_steps =
1832  inquire(*this, exII::EX_INQ_TIME, "Error retrieving number of time steps");
1833 }
1834 
1835 
1836 
1837 void ExodusII_IO_Helper::read_nodal_var_values(std::string nodal_var_name, int time_step)
1838 {
1839  LOG_SCOPE("read_nodal_var_values()", "ExodusII_IO_Helper");
1840 
1841  // Read the nodal variable names from file, so we can see if we have the one we're looking for
1842  this->read_var_names(NODAL);
1843 
1844  // See if we can find the variable we are looking for
1845  unsigned int var_index = 0;
1846  bool found = false;
1847 
1848  // Do a linear search for nodal_var_name in nodal_var_names
1849  for (; var_index<nodal_var_names.size(); ++var_index)
1850  {
1851  found = (nodal_var_names[var_index] == nodal_var_name);
1852  if (found)
1853  break;
1854  }
1855 
1856  if (!found)
1857  {
1858  libMesh::err << "Available variables: " << std::endl;
1859  for (const auto & var_name : nodal_var_names)
1860  libMesh::err << var_name << std::endl;
1861 
1862  libmesh_error_msg("Unable to locate variable named: " << nodal_var_name);
1863  }
1864 
1865  // Clear out any previously read nodal variable values
1866  nodal_var_values.clear();
1867 
1868  std::vector<Real> unmapped_nodal_var_values(num_nodes);
1869 
1870  // Call the Exodus API to read the nodal variable values
1871  ex_err = exII::ex_get_var
1872  (ex_id,
1873  time_step,
1874  exII::EX_NODAL,
1875  var_index+1,
1876  1, // exII::ex_entity_id, not sure exactly what this is but in the ex_get_nodal_var.c shim, they pass 1
1877  num_nodes,
1878  MappedInputVector(unmapped_nodal_var_values, _single_precision).data());
1879  EX_CHECK_ERR(ex_err, "Error reading nodal variable values!");
1880 
1881  for (unsigned i=0; i<static_cast<unsigned>(num_nodes); i++)
1882  {
1883  libmesh_assert_less(i, this->node_num_map.size());
1884 
1885  // Use the node_num_map to obtain the ID of this node in the Exodus file,
1886  // and remember to subtract 1 since libmesh is zero-based and Exodus is 1-based.
1887  const unsigned mapped_node_id = this->node_num_map[i] - 1;
1888 
1889  libmesh_assert_less(i, unmapped_nodal_var_values.size());
1890 
1891  // Store the nodal value in the map.
1892  nodal_var_values[mapped_node_id] = unmapped_nodal_var_values[i];
1893  }
1894 }
1895 
1896 
1897 
1899 {
1900  switch (type)
1901  {
1902  case NODAL:
1904  break;
1905  case ELEMENTAL:
1907  break;
1908  case GLOBAL:
1910  break;
1911  case SIDESET:
1913  break;
1914  case NODESET:
1916  break;
1917  case ELEMSET:
1919  break;
1920  default:
1921  libmesh_error_msg("Unrecognized ExodusVarType " << type);
1922  }
1923 }
1924 
1925 
1926 
1927 void ExodusII_IO_Helper::read_var_names_impl(const char * var_type,
1928  int & count,
1929  std::vector<std::string> & result)
1930 {
1931  // First read and store the number of names we have
1932  ex_err = exII::ex_get_var_param(ex_id, var_type, &count);
1933  EX_CHECK_ERR(ex_err, "Error reading number of variables.");
1934 
1935  // Do nothing if no variables are detected
1936  if (count == 0)
1937  return;
1938 
1939  // Second read the actual names and convert them into a format we can use
1940  NamesData names_table(count, MAX_STR_LENGTH);
1941 
1942  ex_err = exII::ex_get_var_names(ex_id,
1943  var_type,
1944  count,
1945  names_table.get_char_star_star()
1946  );
1947  EX_CHECK_ERR(ex_err, "Error reading variable names!");
1948 
1949  if (verbose)
1950  {
1951  libMesh::out << "Read the variable(s) from the file:" << std::endl;
1952  for (int i=0; i<count; i++)
1953  libMesh::out << names_table.get_char_star(i) << std::endl;
1954  }
1955 
1956  // Allocate enough space for our variable name strings.
1957  result.resize(count);
1958 
1959  // Copy the char buffers into strings.
1960  for (int i=0; i<count; i++)
1961  result[i] = names_table.get_char_star(i); // calls string::op=(const char *)
1962 }
1963 
1964 
1965 
1966 
1967 void
1969  const std::vector<std::string> & names)
1970 {
1971  switch (type)
1972  {
1973  case NODAL:
1974  this->write_var_names_impl("n", num_nodal_vars, names);
1975  break;
1976  case ELEMENTAL:
1977  this->write_var_names_impl("e", num_elem_vars, names);
1978  break;
1979  case GLOBAL:
1980  this->write_var_names_impl("g", num_global_vars, names);
1981  break;
1982  case SIDESET:
1983  {
1984  // Note: calling this function *sets* num_sideset_vars to the
1985  // number of entries in the 'names' vector, num_sideset_vars
1986  // does not already need to be set before calling this.
1987  this->write_var_names_impl("s", num_sideset_vars, names);
1988  break;
1989  }
1990  case NODESET:
1991  {
1992  this->write_var_names_impl("m", num_nodeset_vars, names);
1993  break;
1994  }
1995  case ELEMSET:
1996  {
1997  this->write_var_names_impl("t", num_elemset_vars, names);
1998  break;
1999  }
2000  default:
2001  libmesh_error_msg("Unrecognized ExodusVarType " << type);
2002  }
2003 }
2004 
2005 
2006 
2007 void
2009  int & count,
2010  const std::vector<std::string> & names)
2011 {
2012  // Update the count variable so that it's available to other parts of the class.
2013  count = cast_int<int>(names.size());
2014 
2015  // Write that number of variables to the file.
2016  ex_err = exII::ex_put_var_param(ex_id, var_type, count);
2017  EX_CHECK_ERR(ex_err, "Error setting number of vars.");
2018 
2019  // Nemesis doesn't like trying to write nodal variable names in
2020  // files with no nodes.
2021  if (!this->num_nodes)
2022  return;
2023 
2024  if (count > 0)
2025  {
2026  NamesData names_table(count, MAX_STR_LENGTH);
2027 
2028  // Store the input names in the format required by Exodus.
2029  for (int i=0; i != count; ++i)
2030  {
2031  if(names[i].length() > MAX_STR_LENGTH)
2032  libmesh_warning(
2033  "*** Warning, Exodus variable name \""
2034  << names[i] << "\" too long (max " << MAX_STR_LENGTH
2035  << " characters). Name will be truncated. ");
2036  names_table.push_back_entry(names[i]);
2037  }
2038 
2039  if (verbose)
2040  {
2041  libMesh::out << "Writing variable name(s) to file: " << std::endl;
2042  for (int i=0; i != count; ++i)
2043  libMesh::out << names_table.get_char_star(i) << std::endl;
2044  }
2045 
2046  ex_err = exII::ex_put_var_names(ex_id,
2047  var_type,
2048  count,
2049  names_table.get_char_star_star()
2050  );
2051 
2052  EX_CHECK_ERR(ex_err, "Error writing variable names.");
2053  }
2054 }
2055 
2056 
2057 
2058 void ExodusII_IO_Helper::read_elemental_var_values(std::string elemental_var_name,
2059  int time_step,
2060  std::map<dof_id_type, Real> & elem_var_value_map)
2061 {
2062  LOG_SCOPE("read_elemental_var_values()", "ExodusII_IO_Helper");
2063 
2064  this->read_var_names(ELEMENTAL);
2065 
2066  // See if we can find the variable we are looking for
2067  unsigned int var_index = 0;
2068  bool found = false;
2069 
2070  // Do a linear search for elem_var_name in elemental_var_names
2071  for (; var_index != elem_var_names.size(); ++var_index)
2072  if (elem_var_names[var_index] == elemental_var_name)
2073  {
2074  found = true;
2075  break;
2076  }
2077 
2078  if (!found)
2079  {
2080  libMesh::err << "Available variables: " << std::endl;
2081  for (const auto & var_name : elem_var_names)
2082  libMesh::err << var_name << std::endl;
2083 
2084  libmesh_error_msg("Unable to locate variable named: " << elemental_var_name);
2085  }
2086 
2087  // Sequential index which we can use to look up the element ID in the elem_num_map.
2088  unsigned ex_el_num = 0;
2089 
2090  // Element variable truth table
2091  std::vector<int> var_table(block_ids.size() * elem_var_names.size());
2092  exII::ex_get_truth_table(ex_id, exII::EX_ELEM_BLOCK, block_ids.size(), elem_var_names.size(), var_table.data());
2093 
2094  for (unsigned i=0; i<static_cast<unsigned>(num_elem_blk); i++)
2095  {
2096  ex_err = exII::ex_get_block(ex_id,
2097  exII::EX_ELEM_BLOCK,
2098  block_ids[i],
2099  /*elem_type=*/nullptr,
2101  /*num_nodes_per_entry=*/nullptr,
2102  /*num_edges_per_entry=*/nullptr,
2103  /*num_faces_per_entry=*/nullptr,
2104  /*num_attr=*/nullptr);
2105  EX_CHECK_ERR(ex_err, "Error getting number of elements in block.");
2106 
2107  // If the current variable isn't active on this subdomain, advance
2108  // the index by the number of elements on this block and go to the
2109  // next loop iteration.
2110  if (!var_table[elem_var_names.size()*i + var_index])
2111  {
2112  ex_el_num += num_elem_this_blk;
2113  continue;
2114  }
2115 
2116  std::vector<Real> block_elem_var_values(num_elem_this_blk);
2117 
2118  ex_err = exII::ex_get_var
2119  (ex_id,
2120  time_step,
2121  exII::EX_ELEM_BLOCK,
2122  var_index+1,
2123  block_ids[i],
2125  MappedInputVector(block_elem_var_values, _single_precision).data());
2126  EX_CHECK_ERR(ex_err, "Error getting elemental values.");
2127 
2128  for (unsigned j=0; j<static_cast<unsigned>(num_elem_this_blk); j++)
2129  {
2130  // Use the elem_num_map to obtain the ID of this element in the Exodus file,
2131  // and remember to subtract 1 since libmesh is zero-based and Exodus is 1-based.
2132  unsigned mapped_elem_id = this->elem_num_map[ex_el_num] - 1;
2133 
2134  // Store the elemental value in the map.
2135  elem_var_value_map[mapped_elem_id] = block_elem_var_values[j];
2136 
2137  // Go to the next sequential element ID.
2138  ex_el_num++;
2139  }
2140  }
2141 }
2142 
2143 
2144 // For Writing Solutions
2145 
2146 void ExodusII_IO_Helper::create(std::string filename)
2147 {
2148  // If we're processor 0, always create the file.
2149  // If we running on all procs, e.g. as one of several Nemesis files, also
2150  // call create there.
2151  if ((this->processor_id() == 0) || (!_run_only_on_proc0))
2152  {
2153  int
2154  comp_ws = 0,
2155  io_ws = 0;
2156 
2157  if (_single_precision)
2158  {
2159  comp_ws = cast_int<int>(sizeof(float));
2160  io_ws = cast_int<int>(sizeof(float));
2161  }
2162  // Fall back on double precision when necessary since ExodusII
2163  // doesn't seem to support long double
2164  else
2165  {
2166  comp_ws = cast_int<int>
2167  (std::min(sizeof(Real), sizeof(double)));
2168  io_ws = cast_int<int>
2169  (std::min(sizeof(Real), sizeof(double)));
2170  }
2171 
2172  // By default we just open the Exodus file in "EX_CLOBBER" mode,
2173  // which, according to "ncdump -k", writes the file in "64-bit
2174  // offset" mode, which is a NETCDF3 file format.
2175  int mode = EX_CLOBBER;
2176 
2177  // If HDF5 is available, by default we will write Exodus files
2178  // in a more modern NETCDF4-compatible format. For this file
2179  // type, "ncdump -k" will report "netCDF-4".
2180 #ifdef LIBMESH_HAVE_HDF5
2181  if (this->_write_hdf5)
2182  {
2183  mode |= EX_NETCDF4;
2184  mode |= EX_NOCLASSIC;
2185  }
2186 #endif
2187 
2188  {
2189  FPEDisabler disable_fpes;
2190  ex_id = exII::ex_create(filename.c_str(), mode, &comp_ws, &io_ws);
2191  }
2192 
2193  EX_CHECK_ERR(ex_id, "Error creating ExodusII/Nemesis mesh file.");
2194 
2195  if (verbose)
2196  libMesh::out << "File created successfully." << std::endl;
2197  }
2198 
2199  opened_for_writing = true;
2200  _opened_by_create = true;
2201  current_filename = filename;
2202 }
2203 
2204 
2205 
2206 void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh, bool use_discontinuous)
2207 {
2208  // The majority of this function only executes on processor 0, so any functions
2209  // which are collective, like n_active_elem() or n_edge_conds() must be called
2210  // before the processors' execution paths diverge.
2211  libmesh_parallel_only(mesh.comm());
2212 
2213  unsigned int n_active_elem = mesh.n_active_elem();
2214  const BoundaryInfo & bi = mesh.get_boundary_info();
2215  num_edge = bi.n_edge_conds();
2216 
2217  // We need to know about all processors' subdomains
2218  subdomain_id_type subdomain_id_end = 0;
2219  auto subdomain_map = build_subdomain_map(mesh, _add_sides, subdomain_id_end);
2220 
2221  num_elem = n_active_elem;
2222  num_nodes = 0;
2223 
2224  // If we're adding face elements they'll need copies of their nodes.
2225  // We also have to count of how many nodes (and gaps between nodes!)
2226  // are on each processor, to calculate offsets for any nodal data
2227  // writing later.
2228  _added_side_node_offsets.clear();
2229  if (_add_sides)
2230  {
2231  dof_id_type num_side_elem = 0;
2232  dof_id_type num_local_side_nodes = 0;
2233 
2234  for (const auto & elem : mesh.active_local_element_ptr_range())
2235  {
2236  for (auto s : elem->side_index_range())
2237  {
2239  continue;
2240 
2241  num_side_elem++;
2242  num_local_side_nodes += elem->nodes_on_side(s).size();
2243  }
2244  }
2245 
2246  mesh.comm().sum(num_side_elem);
2247  num_elem += num_side_elem;
2248 
2249  mesh.comm().allgather(num_local_side_nodes, _added_side_node_offsets);
2250  const processor_id_type n_proc = mesh.n_processors();
2251  libmesh_assert_equal_to(n_proc, _added_side_node_offsets.size());
2252 
2253  for (auto p : make_range(n_proc-1))
2255 
2256  num_nodes = _added_side_node_offsets[n_proc-1];
2257 
2258  dof_id_type n_local_nodes = cast_int<dof_id_type>
2259  (std::distance(mesh.local_nodes_begin(),
2260  mesh.local_nodes_end()));
2261  dof_id_type n_total_nodes = n_local_nodes;
2262  mesh.comm().sum(n_total_nodes);
2263 
2264  const dof_id_type max_nn = mesh.max_node_id();
2265  const dof_id_type n_gaps = max_nn - n_total_nodes;
2266  const dof_id_type gaps_per_processor = n_gaps / n_proc;
2267  const dof_id_type remainder_gaps = n_gaps % n_proc;
2268 
2269  n_local_nodes = n_local_nodes + // Actual nodes
2270  gaps_per_processor + // Our even share of gaps
2271  (mesh.processor_id() < remainder_gaps); // Leftovers
2272 
2273  mesh.comm().allgather(n_local_nodes, _true_node_offsets);
2274  for (auto p : make_range(n_proc-1))
2276  libmesh_assert_equal_to(_true_node_offsets[n_proc-1], mesh.max_node_id());
2277  }
2278 
2279  // If _write_as_dimension is nonzero, use it to set num_dim in the Exodus file.
2280  if (_write_as_dimension)
2284  else
2286 
2287  if ((_run_only_on_proc0) && (this->processor_id() != 0))
2288  return;
2289 
2290  if (!use_discontinuous)
2291  {
2292  // Don't rely on mesh.n_nodes() here. If ReplicatedMesh nodes
2293  // have been deleted without renumbering after, it will be
2294  // incorrect.
2295  num_nodes += cast_int<int>(std::distance(mesh.nodes_begin(),
2296  mesh.nodes_end()));
2297  }
2298  else
2299  {
2300  for (const auto & elem : mesh.active_element_ptr_range())
2301  num_nodes += elem->n_nodes();
2302  }
2303 
2304  std::set<boundary_id_type> unique_side_boundaries;
2305  std::vector<boundary_id_type> unique_node_boundaries;
2306 
2307  // Build set of unique sideset (+shellface) ids
2308  {
2309  // Start with "side" boundaries (i.e. of 3D elements)
2310  std::vector<boundary_id_type> side_boundaries;
2311  bi.build_side_boundary_ids(side_boundaries);
2312  unique_side_boundaries.insert(side_boundaries.begin(), side_boundaries.end());
2313 
2314  // Add shell face boundaries to the list of side boundaries, since ExodusII
2315  // treats these the same way.
2316  std::vector<boundary_id_type> shellface_boundaries;
2317  bi.build_shellface_boundary_ids(shellface_boundaries);
2318  unique_side_boundaries.insert(shellface_boundaries.begin(), shellface_boundaries.end());
2319 
2320  // Add any empty-but-named side boundary ids
2321  for (const auto & pr : bi.get_sideset_name_map())
2322  unique_side_boundaries.insert(pr.first);
2323  }
2324 
2325  // Build set of unique nodeset ids
2326  bi.build_node_boundary_ids(unique_node_boundaries);
2327  for (const auto & pair : bi.get_nodeset_name_map())
2328  {
2329  const boundary_id_type id = pair.first;
2330 
2331  if (std::find(unique_node_boundaries.begin(),
2332  unique_node_boundaries.end(), id)
2333  == unique_node_boundaries.end())
2334  unique_node_boundaries.push_back(id);
2335  }
2336 
2337  num_side_sets = cast_int<int>(unique_side_boundaries.size());
2338  num_node_sets = cast_int<int>(unique_node_boundaries.size());
2339 
2340  num_elem_blk = cast_int<int>(subdomain_map.size());
2341 
2342  if (str_title.size() > MAX_LINE_LENGTH)
2343  {
2344  libMesh::err << "Warning, Exodus files cannot have titles longer than "
2345  << MAX_LINE_LENGTH
2346  << " characters. Your title will be truncated."
2347  << std::endl;
2348  str_title.resize(MAX_LINE_LENGTH);
2349  }
2350 
2351  // Edge BCs are handled a bit differently than sidesets and nodesets.
2352  // They are written as separate "edge blocks", and then edge variables
2353  // can be defined on those blocks. That is, they are not written as
2354  // edge sets, since edge sets must refer to edges stored elsewhere.
2355  // We write a separate edge block for each unique boundary id that
2356  // we have.
2357  num_edge_blk = bi.get_edge_boundary_ids().size();
2358 
2359  // Check whether the Mesh Elems have an extra_integer called "elemset_code".
2360  // If so, this means that the mesh defines elemsets via the
2361  // extra_integers capability of Elems.
2362  if (mesh.has_elem_integer("elemset_code"))
2363  {
2364  // unsigned int elemset_index =
2365  // mesh.get_elem_integer_index("elemset_code");
2366 
2367  // Debugging
2368  // libMesh::out << "Mesh defines an elemset_code at index " << elemset_index << std::endl;
2369 
2370  // Store the number of elemsets in the exo file header.
2372  }
2373 
2374  // Build an ex_init_params() structure that is to be passed to the
2375  // newer ex_put_init_ext() API. The new API will eventually allow us
2376  // to store edge and face data in the Exodus file.
2377  //
2378  // Notes:
2379  // * We use C++11 zero initialization syntax to make sure that all
2380  // members of the struct (including ones we aren't using) are
2381  // given sensible values.
2382  // * For the "title" field, we manually do a null-terminated string
2383  // copy since std::string does not null-terminate but it does
2384  // return the number of characters successfully copied.
2385  exII::ex_init_params params = {};
2386  params.title[str_title.copy(params.title, MAX_LINE_LENGTH)] = '\0';
2387  params.num_dim = num_dim;
2388  params.num_nodes = num_nodes;
2389  params.num_elem = num_elem;
2390  params.num_elem_blk = num_elem_blk;
2391  params.num_node_sets = num_node_sets;
2392  params.num_side_sets = num_side_sets;
2393  params.num_elem_sets = num_elem_sets;
2394  params.num_edge_blk = num_edge_blk;
2395  params.num_edge = num_edge;
2396 
2397  ex_err = exII::ex_put_init_ext(ex_id, &params);
2398  EX_CHECK_ERR(ex_err, "Error initializing new Exodus file.");
2399 }
2400 
2401 
2402 
2403 void ExodusII_IO_Helper::write_nodal_coordinates(const MeshBase & mesh, bool use_discontinuous)
2404 {
2405  if ((_run_only_on_proc0) && (this->processor_id() != 0))
2406  return;
2407 
2408  // Clear existing data from any previous calls.
2409  x.clear();
2410  y.clear();
2411  z.clear();
2412  node_num_map.clear();
2413 
2414  // Reserve space in the nodal coordinate vectors. num_nodes is
2415  // exact, this just allows us to do away with one potentially
2416  // error-inducing loop index.
2417  x.reserve(num_nodes);
2418  y.reserve(num_nodes);
2419  z.reserve(num_nodes);
2420 
2421  auto push_node = [this](const Point & p) {
2422  x.push_back(p(0) + _coordinate_offset(0));
2423 
2424 #if LIBMESH_DIM > 1
2425  y.push_back(p(1) + _coordinate_offset(1));
2426 #else
2427  y.push_back(0.);
2428 #endif
2429 #if LIBMESH_DIM > 2
2430  z.push_back(p(2) + _coordinate_offset(2));
2431 #else
2432  z.push_back(0.);
2433 #endif
2434  };
2435 
2436  // And in the node_num_map - since the nodes aren't organized in
2437  // blocks, libmesh will always write out the identity map
2438  // here... unless there has been some refinement and coarsening, or
2439  // node deletion without a corresponding call to contract(). You
2440  // need to write this any time there could be 'holes' in the node
2441  // numbering, so we write it every time.
2442 
2443  // Let's skip the node_num_map in the discontinuous and add_sides
2444  // cases, since we're effectively duplicating nodes for the sake of
2445  // discontinuous visualization, so it isn't clear how to deal with
2446  // node_num_map here. This means that writing meshes in such a way
2447  // won't work with element numberings that have id "holes".
2448 
2449  if (!use_discontinuous && !_add_sides)
2450  node_num_map.reserve(num_nodes);
2451 
2452  // Clear out any previously-mapped node IDs.
2454 
2455  if (!use_discontinuous)
2456  {
2457  for (const auto & node_ptr : mesh.node_ptr_range())
2458  {
2459  const Node & node = *node_ptr;
2460 
2461  push_node(node);
2462 
2463  // Fill in node_num_map entry with the proper (1-based) node
2464  // id, unless we're not going to be able to keep the map up
2465  // later.
2466  if (!_add_sides)
2467  node_num_map.push_back(node.id() + 1);
2468 
2469  // Also map the zero-based libmesh node id to the 1-based
2470  // Exodus ID it will be assigned (this is equivalent to the
2471  // current size of the x vector).
2472  libmesh_node_num_to_exodus[ cast_int<int>(node.id()) ] = cast_int<int>(x.size());
2473  }
2474  }
2475  else
2476  {
2477  for (const auto & elem : mesh.active_element_ptr_range())
2478  for (const Node & node : elem->node_ref_range())
2479  {
2480  push_node(node);
2481 
2482  // Let's skip the node_num_map in the discontinuous
2483  // case, since we're effectively duplicating nodes for
2484  // the sake of discontinuous visualization, so it isn't
2485  // clear how to deal with node_num_map here. This means
2486  // that writing discontinuous meshes won't work with
2487  // element numberings that have "holes".
2488  }
2489  }
2490 
2491  if (_add_sides)
2492  {
2493  // To match the numbering of parallel-generated nodal solutions
2494  // on fake side nodes, we need to loop through elements from
2495  // earlier ranks first.
2496  std::vector<std::vector<const Elem *>>
2497  elems_by_pid(mesh.n_processors());
2498 
2499  for (const auto & elem : mesh.active_element_ptr_range())
2500  elems_by_pid[elem->processor_id()].push_back(elem);
2501 
2502  for (auto p : index_range(elems_by_pid))
2503  for (const Elem * elem : elems_by_pid[p])
2504  for (auto s : elem->side_index_range())
2505  {
2507  continue;
2508 
2509  const std::vector<unsigned int> side_nodes =
2510  elem->nodes_on_side(s);
2511 
2512  for (auto n : side_nodes)
2513  push_node(elem->point(n));
2514  }
2515 
2516  // Node num maps just don't make sense if we're adding a bunch
2517  // of visualization nodes that are independent copies of the
2518  // same libMesh node.
2519  node_num_map.clear();
2520  }
2521 
2522  ex_err = exII::ex_put_coord
2523  (ex_id,
2524  x.empty() ? nullptr : MappedOutputVector(x, _single_precision).data(),
2525  y.empty() ? nullptr : MappedOutputVector(y, _single_precision).data(),
2526  z.empty() ? nullptr : MappedOutputVector(z, _single_precision).data());
2527 
2528  EX_CHECK_ERR(ex_err, "Error writing coordinates to Exodus file.");
2529 
2530  if (!use_discontinuous && !_add_sides)
2531  {
2532  // Also write the (1-based) node_num_map to the file.
2533  ex_err = exII::ex_put_node_num_map(ex_id, node_num_map.data());
2534  EX_CHECK_ERR(ex_err, "Error writing node_num_map");
2535  }
2536 }
2537 
2538 
2539 
2540 void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_discontinuous)
2541 {
2542  LOG_SCOPE("write_elements()", "ExodusII_IO_Helper");
2543 
2544  // Map from block ID to a vector of element IDs in that block. Element
2545  // IDs are now of type dof_id_type, subdomain IDs are of type subdomain_id_type.
2546  subdomain_id_type subdomain_id_end = 0;
2547  auto subdomain_map = build_subdomain_map(mesh, _add_sides, subdomain_id_end);
2548 
2549  if ((_run_only_on_proc0) && (this->processor_id() != 0))
2550  return;
2551 
2552  // element map vector
2553  num_elem_blk = cast_int<int>(subdomain_map.size());
2554  block_ids.resize(num_elem_blk);
2555 
2556  std::vector<int> elem_blk_id;
2557  std::vector<int> num_elem_this_blk_vec;
2558  std::vector<int> num_nodes_per_elem_vec;
2559  std::vector<int> num_edges_per_elem_vec;
2560  std::vector<int> num_faces_per_elem_vec;
2561  std::vector<int> num_attr_vec;
2562  NamesData elem_type_table(num_elem_blk, MAX_STR_LENGTH);
2563 
2564  // Note: It appears that there is a bug in exodusII::ex_put_name where
2565  // the index returned from the ex_id_lkup is erroneously used. For now
2566  // the work around is to use the alternative function ex_put_names, but
2567  // this function requires a char ** data structure.
2568  NamesData names_table(num_elem_blk, MAX_STR_LENGTH);
2569 
2570  num_elem = 0;
2571 
2572  // counter indexes into the block_ids vector
2573  unsigned int counter = 0;
2574  for (auto & [subdomain_id, element_id_vec] : subdomain_map)
2575  {
2576  block_ids[counter] = subdomain_id;
2577 
2578  const ElemType elem_t = (subdomain_id >= subdomain_id_end) ?
2579  ElemType(subdomain_id - subdomain_id_end) :
2580  mesh.elem_ref(element_id_vec[0]).type();
2581 
2582  if (subdomain_id >= subdomain_id_end)
2583  {
2585  libmesh_assert(element_id_vec.size() == 1);
2586  num_elem_this_blk_vec.push_back
2587  (cast_int<int>(element_id_vec[0]));
2588  names_table.push_back_entry
2589  (Utility::enum_to_string<ElemType>(elem_t));
2590  }
2591  else
2592  {
2593  libmesh_assert(!element_id_vec.empty());
2594  num_elem_this_blk_vec.push_back
2595  (cast_int<int>(element_id_vec.size()));
2596  names_table.push_back_entry
2597  (mesh.subdomain_name(subdomain_id));
2598  }
2599 
2600  num_elem += num_elem_this_blk_vec.back();
2601 
2602  // Use the first element in this block to get representative information.
2603  // Note that Exodus assumes all elements in a block are of the same type!
2604  // We are using that same assumption here!
2605  const auto & conv = get_conversion(elem_t);
2607  if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint)
2608  libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
2609 
2610  elem_blk_id.push_back(subdomain_id);
2611  elem_type_table.push_back_entry(conv.exodus_elem_type().c_str());
2612  num_nodes_per_elem_vec.push_back(num_nodes_per_elem);
2613  num_attr_vec.push_back(0); // we don't currently use elem block attributes.
2614  num_edges_per_elem_vec.push_back(0); // We don't currently store any edge blocks
2615  num_faces_per_elem_vec.push_back(0); // We don't currently store any face blocks
2616  ++counter;
2617  }
2618 
2619  elem_num_map.resize(num_elem);
2620  std::vector<int>::iterator curr_elem_map_end = elem_num_map.begin();
2621 
2622  // In the case of discontinuous plotting we initialize a map from
2623  // (element, node) pairs to the corresponding discontinuous node index.
2624  // This ordering must match the ordering used in write_nodal_coordinates.
2625  //
2626  // Note: This map takes the place of the libmesh_node_num_to_exodus map in
2627  // the discontinuous case.
2628  std::map<std::pair<dof_id_type, unsigned int>, dof_id_type> discontinuous_node_indices;
2629  dof_id_type node_counter = 1; // Exodus numbering is 1-based
2630  if (use_discontinuous)
2631  {
2632  for (const auto & elem : mesh.active_element_ptr_range())
2633  for (auto n : elem->node_index_range())
2634  discontinuous_node_indices[std::make_pair(elem->id(),n)] =
2635  node_counter++;
2636  }
2637  else
2638  node_counter = mesh.max_node_id() + 1; // Exodus numbering is 1-based
2639 
2640  if (_add_sides)
2641  {
2642  for (const Elem * elem : mesh.active_element_ptr_range())
2643  {
2644  // We'll use "past-the-end" indices to indicate side node
2645  // copies
2646  unsigned int local_node_index = elem->n_nodes();
2647 
2648  for (auto s : elem->side_index_range())
2649  {
2651  continue;
2652 
2653  const std::vector<unsigned int> side_nodes =
2654  elem->nodes_on_side(s);
2655 
2656  for (auto n : index_range(side_nodes))
2657  {
2658  libmesh_ignore(n);
2659  discontinuous_node_indices
2660  [std::make_pair(elem->id(),local_node_index++)] =
2661  node_counter++;
2662  }
2663  }
2664  }
2665  }
2666 
2667  // Reference to the BoundaryInfo object for convenience.
2668  const BoundaryInfo & bi = mesh.get_boundary_info();
2669 
2670  // Build list of (elem, edge, id) triples
2671  std::vector<BoundaryInfo::BCTuple> edge_tuples = bi.build_edge_list();
2672 
2673  // Build the connectivity array for each edge block. The connectivity array
2674  // is a vector<int> with "num_edges * num_nodes_per_edge" entries. We write
2675  // the Exodus node numbers to the connectivity arrays so that they can
2676  // be used directly in the calls to exII::ex_put_conn() below. We also keep
2677  // track of the ElemType and the number of nodes for each boundary_id. All
2678  // edges with a given boundary_id must be of the same type.
2679  std::map<boundary_id_type, std::vector<int>> edge_id_to_conn;
2680  std::map<boundary_id_type, std::pair<ElemType, unsigned int>> edge_id_to_elem_type;
2681 
2682  std::unique_ptr<const Elem> edge;
2683  for (const auto & t : edge_tuples)
2684  {
2685  dof_id_type elem_id = std::get<0>(t);
2686  unsigned int edge_id = std::get<1>(t);
2687  boundary_id_type b_id = std::get<2>(t);
2688 
2689  // Build the edge in question
2690  mesh.elem_ptr(elem_id)->build_edge_ptr(edge, edge_id);
2691 
2692  // Error checking: make sure that all edges in this block are
2693  // the same geometric type.
2694  if (const auto check_it = edge_id_to_elem_type.find(b_id);
2695  check_it == edge_id_to_elem_type.end())
2696  {
2697  // Keep track of the ElemType and number of nodes in this boundary id.
2698  edge_id_to_elem_type[b_id] = std::make_pair(edge->type(), edge->n_nodes());
2699  }
2700  else
2701  {
2702  // Make sure the existing data is consistent
2703  const auto & val_pair = check_it->second;
2704  libmesh_error_msg_if(val_pair.first != edge->type() || val_pair.second != edge->n_nodes(),
2705  "All edges in a block must have same geometric type.");
2706  }
2707 
2708  // Get reference to the connectivity array for this block
2709  auto & conn = edge_id_to_conn[b_id];
2710 
2711  // For each node on the edge, look up the exodus node id and
2712  // store it in the conn array. Note: all edge types have
2713  // identity node mappings so we don't bother with Conversion
2714  // objects here.
2715  for (auto n : edge->node_index_range())
2716  {
2717  // We look up Exodus node numbers differently if we are
2718  // writing a discontinuous Exodus file.
2719  int exodus_node_id = -1;
2720 
2721  if (!use_discontinuous)
2722  {
2723  dof_id_type libmesh_node_id = edge->node_ptr(n)->id();
2724  exodus_node_id = libmesh_map_find
2725  (libmesh_node_num_to_exodus, cast_int<int>(libmesh_node_id));
2726  }
2727  else
2728  {
2729  // Get the node on the element containing this edge
2730  // which corresponds to edge node n. Then use that id to look up
2731  // the exodus_node_id in the discontinuous_node_indices map.
2732  unsigned int pn = mesh.elem_ptr(elem_id)->local_edge_node(edge_id, n);
2733  exodus_node_id = libmesh_map_find
2734  (discontinuous_node_indices, std::make_pair(elem_id, pn));
2735  }
2736 
2737  conn.push_back(exodus_node_id);
2738  }
2739  }
2740 
2741  // Make sure we have the same number of edge ids that we thought we would.
2742  libmesh_assert(static_cast<int>(edge_id_to_conn.size()) == num_edge_blk);
2743 
2744  // Build data structures describing edge blocks. This information must be
2745  // be passed to exII::ex_put_concat_all_blocks() at the same time as the
2746  // information about elem blocks.
2747  std::vector<int> edge_blk_id;
2748  NamesData edge_type_table(num_edge_blk, MAX_STR_LENGTH);
2749  std::vector<int> num_edge_this_blk_vec;
2750  std::vector<int> num_nodes_per_edge_vec;
2751  std::vector<int> num_attr_edge_vec;
2752 
2753  // We also build a data structure of edge block names which can
2754  // later be passed to exII::ex_put_names().
2755  NamesData edge_block_names_table(num_edge_blk, MAX_STR_LENGTH);
2756 
2757  // Note: We are going to use the edge **boundary** ids as **block** ids.
2758  for (const auto & pr : edge_id_to_conn)
2759  {
2760  // Store the edge block id in the array to be passed to Exodus.
2761  boundary_id_type id = pr.first;
2762  edge_blk_id.push_back(id);
2763 
2764  // Set Exodus element type and number of nodes for this edge block.
2765  const auto & elem_type_node_count = edge_id_to_elem_type[id];
2766  const auto & conv = get_conversion(elem_type_node_count.first);
2767  edge_type_table.push_back_entry(conv.exodus_type.c_str());
2768  num_nodes_per_edge_vec.push_back(elem_type_node_count.second);
2769 
2770  // The number of edges is the number of entries in the connectivity
2771  // array divided by the number of nodes per edge.
2772  num_edge_this_blk_vec.push_back(pr.second.size() / elem_type_node_count.second);
2773 
2774  // We don't store any attributes currently
2775  num_attr_edge_vec.push_back(0);
2776 
2777  // Store the name of this edge block
2778  edge_block_names_table.push_back_entry(bi.get_edgeset_name(id));
2779  }
2780 
2781  // Zero-initialize and then fill in an exII::ex_block_params struct
2782  // with the data we have collected. This new API replaces the old
2783  // exII::ex_put_concat_elem_block() API, and will eventually allow
2784  // us to also allocate space for edge/face blocks if desired.
2785  //
2786  // TODO: It seems like we should be able to take advantage of the
2787  // optimization where you set define_maps==1, but when I tried this
2788  // I got the error: "failed to find node map size". I think the
2789  // problem is that we need to first specify a nonzero number of
2790  // node/elem maps during the call to ex_put_init_ext() in order for
2791  // this to work correctly.
2792  exII::ex_block_params params = {};
2793 
2794  // Set pointers for information about elem blocks.
2795  params.elem_blk_id = elem_blk_id.data();
2796  params.elem_type = elem_type_table.get_char_star_star();
2797  params.num_elem_this_blk = num_elem_this_blk_vec.data();
2798  params.num_nodes_per_elem = num_nodes_per_elem_vec.data();
2799  params.num_edges_per_elem = num_edges_per_elem_vec.data();
2800  params.num_faces_per_elem = num_faces_per_elem_vec.data();
2801  params.num_attr_elem = num_attr_vec.data();
2802  params.define_maps = 0;
2803 
2804  // Set pointers to edge block information only if we actually have some.
2805  if (num_edge_blk)
2806  {
2807  params.edge_blk_id = edge_blk_id.data();
2808  params.edge_type = edge_type_table.get_char_star_star();
2809  params.num_edge_this_blk = num_edge_this_blk_vec.data();
2810  params.num_nodes_per_edge = num_nodes_per_edge_vec.data();
2811  params.num_attr_edge = num_attr_edge_vec.data();
2812  }
2813 
2814  ex_err = exII::ex_put_concat_all_blocks(ex_id, &params);
2815  EX_CHECK_ERR(ex_err, "Error writing element blocks.");
2816 
2817  // This counter is used to fill up the libmesh_elem_num_to_exodus map in the loop below.
2818  unsigned libmesh_elem_num_to_exodus_counter = 0;
2819 
2820  // We need these later if we're adding fake sides, but we don't need
2821  // to recalculate it.
2822  auto num_elem_this_blk_it = num_elem_this_blk_vec.begin();
2823  auto next_fake_id = mesh.max_elem_id() + 1; // 1-based numbering in Exodus
2824 
2825  for (auto & [subdomain_id, element_id_vec] : subdomain_map)
2826  {
2827  // Use the first element in the block to get representative
2828  // information for a "real" block. Note that Exodus assumes all
2829  // elements in a block are of the same type! We are using that
2830  // same assumption here!
2831  const ElemType elem_t = (subdomain_id >= subdomain_id_end) ?
2832  ElemType(subdomain_id - subdomain_id_end) :
2833  mesh.elem_ref(element_id_vec[0]).type();
2834 
2835  const auto & conv = get_conversion(elem_t);
2837  if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint)
2838  libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
2839 
2840  // If this is a *real* block, we just loop over vectors of
2841  // element ids to add.
2842  if (subdomain_id < subdomain_id_end)
2843  {
2844  connect.resize(element_id_vec.size()*num_nodes_per_elem);
2845 
2846  for (auto i : index_range(element_id_vec))
2847  {
2848  unsigned int elem_id = element_id_vec[i];
2849  libmesh_elem_num_to_exodus[elem_id] = ++libmesh_elem_num_to_exodus_counter; // 1-based indexing for Exodus
2850 
2851  const Elem & elem = mesh.elem_ref(elem_id);
2852 
2853  // We *might* be able to get away with writing mixed element
2854  // types which happen to have the same number of nodes, but
2855  // do we actually *want* to get away with that?
2856  // .) No visualization software would be able to handle it.
2857  // .) There'd be no way for us to read it back in reliably.
2858  // .) Even elements with the same number of nodes may have different connectivities (?)
2859 
2860  // This needs to be more than an assert so we don't fail
2861  // with a mysterious segfault while trying to write mixed
2862  // element meshes in optimized mode.
2863  libmesh_error_msg_if(elem.type() != conv.libmesh_elem_type(),
2864  "Error: Exodus requires all elements with a given subdomain ID to be the same type.\n"
2865  << "Can't write both "
2866  << Utility::enum_to_string(elem.type())
2867  << " and "
2868  << Utility::enum_to_string(conv.libmesh_elem_type())
2869  << " in the same block!");
2870 
2871  for (unsigned int j=0; j<static_cast<unsigned int>(num_nodes_per_elem); ++j)
2872  {
2873  unsigned int connect_index = cast_int<unsigned int>((i*num_nodes_per_elem)+j);
2874  unsigned elem_node_index = conv.get_inverse_node_map(j); // inverse node map is for writing.
2875  if (!use_discontinuous)
2876  {
2877  // The global id for the current node in libmesh.
2878  dof_id_type libmesh_node_id = elem.node_id(elem_node_index);
2879 
2880  // Write the Exodus global node id associated with
2881  // this libmesh node number to the connectivity
2882  // array, or throw an error if it's not found.
2883  connect[connect_index] =
2884  libmesh_map_find(libmesh_node_num_to_exodus,
2885  cast_int<int>(libmesh_node_id));
2886  }
2887  else
2888  {
2889  // Look up the (elem_id, elem_node_index) pair in the map.
2890  connect[connect_index] =
2891  libmesh_map_find(discontinuous_node_indices,
2892  std::make_pair(elem_id, elem_node_index));
2893  }
2894  }
2895  }
2896 
2897  // This transform command stores its result in a range that
2898  // begins at the third argument, so this command is adding
2899  // values to the elem_num_map vector starting from
2900  // curr_elem_map_end. Here we add 1 to each id to make a
2901  // 1-based exodus file.
2902  curr_elem_map_end = std::transform
2903  (element_id_vec.begin(),
2904  element_id_vec.end(),
2905  curr_elem_map_end,
2906  [](dof_id_type id){return id+1;});
2907  }
2908  // If this is a "fake" block of added sides, we build those as
2909  // we go.
2910  else
2911  {
2913 
2914  libmesh_assert(num_elem_this_blk_it != num_elem_this_blk_vec.end());
2915  num_elem_this_blk = *num_elem_this_blk_it;
2916 
2918 
2919  std::size_t connect_index = 0;
2920  for (const auto & elem : mesh.active_element_ptr_range())
2921  {
2922  unsigned int local_node_index = elem->n_nodes();
2923 
2924  for (auto s : elem->side_index_range())
2925  {
2927  continue;
2928 
2929  if (elem->side_type(s) != elem_t)
2930  continue;
2931 
2932  const std::vector<unsigned int> side_nodes =
2933  elem->nodes_on_side(s);
2934 
2935  for (auto n : index_range(side_nodes))
2936  {
2937  libmesh_ignore(n);
2938  const int exodus_node_id = libmesh_map_find
2939  (discontinuous_node_indices,
2940  std::make_pair(elem->id(), local_node_index++));
2941  libmesh_assert_less(connect_index, connect.size());
2942  connect[connect_index++] = exodus_node_id;
2943  }
2944  }
2945  }
2946 
2947  auto old_curr_map_end = curr_elem_map_end;
2948  curr_elem_map_end += num_elem_this_blk;
2949 
2950  std::generate
2951  (old_curr_map_end, curr_elem_map_end,
2952  [&next_fake_id](){return next_fake_id++;});
2953  }
2954 
2955  ++num_elem_this_blk_it;
2956 
2957  ex_err = exII::ex_put_conn
2958  (ex_id,
2959  exII::EX_ELEM_BLOCK,
2960  subdomain_id,
2961  connect.data(), // node_conn
2962  nullptr, // elem_edge_conn (unused)
2963  nullptr); // elem_face_conn (unused)
2964  EX_CHECK_ERR(ex_err, "Error writing element connectivities");
2965  }
2966 
2967  // write out the element number map that we created
2968  ex_err = exII::ex_put_elem_num_map(ex_id, elem_num_map.data());
2969  EX_CHECK_ERR(ex_err, "Error writing element map");
2970 
2971  // Write out the block names
2972  if (num_elem_blk > 0)
2973  {
2974  ex_err = exII::ex_put_names(ex_id, exII::EX_ELEM_BLOCK, names_table.get_char_star_star());
2975  EX_CHECK_ERR(ex_err, "Error writing element block names");
2976  }
2977 
2978  // Write out edge blocks if we have any
2979  for (const auto & pr : edge_id_to_conn)
2980  {
2981  ex_err = exII::ex_put_conn
2982  (ex_id,
2983  exII::EX_EDGE_BLOCK,
2984  pr.first,
2985  pr.second.data(), // node_conn
2986  nullptr, // elem_edge_conn (unused)
2987  nullptr); // elem_face_conn (unused)
2988  EX_CHECK_ERR(ex_err, "Error writing element connectivities");
2989  }
2990 
2991  // Write out the edge block names, if any.
2992  if (num_edge_blk > 0)
2993  {
2994  ex_err = exII::ex_put_names
2995  (ex_id,
2996  exII::EX_EDGE_BLOCK,
2997  edge_block_names_table.get_char_star_star());
2998  EX_CHECK_ERR(ex_err, "Error writing edge block names");
2999  }
3000 }
3001 
3002 
3003 
3004 
3006 {
3007  LOG_SCOPE("write_sidesets()", "ExodusII_IO_Helper");
3008 
3009  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3010  return;
3011 
3012  // Maps from sideset id to lists of corresponding element ids and side ids
3013  std::map<int, std::vector<int>> elem_lists;
3014  std::map<int, std::vector<int>> side_lists;
3015  std::set<boundary_id_type> side_boundary_ids;
3016 
3017  {
3018  // Accumulate the vectors to pass into ex_put_side_set
3019  // build_side_lists() returns a vector of (elem, side, bc) tuples.
3020  for (const auto & t : mesh.get_boundary_info().build_side_list())
3021  {
3022  std::vector<const Elem *> family;
3023 #ifdef LIBMESH_ENABLE_AMR
3024 
3028  mesh.elem_ref(std::get<0>(t)).active_family_tree_by_side(family, std::get<1>(t), false);
3029 #else
3030  family.push_back(mesh.elem_ptr(std::get<0>(t)));
3031 #endif
3032 
3033  for (const auto & f : family)
3034  {
3035  const auto & conv = get_conversion(mesh.elem_ptr(f->id())->type());
3036 
3037  // Use the libmesh to exodus data structure map to get the proper sideset IDs
3038  // The data structure contains the "collapsed" contiguous ids
3039  elem_lists[std::get<2>(t)].push_back(libmesh_elem_num_to_exodus[f->id()]);
3040  side_lists[std::get<2>(t)].push_back(conv.get_inverse_side_map(std::get<1>(t)));
3041  }
3042  }
3043 
3044  std::vector<boundary_id_type> tmp;
3046  side_boundary_ids.insert(tmp.begin(), tmp.end());
3047  }
3048 
3049  {
3050  // add data for shell faces, if needed
3051 
3052  // Accumulate the vectors to pass into ex_put_side_set
3053  for (const auto & t : mesh.get_boundary_info().build_shellface_list())
3054  {
3055  std::vector<const Elem *> family;
3056 #ifdef LIBMESH_ENABLE_AMR
3057 
3061  mesh.elem_ref(std::get<0>(t)).active_family_tree_by_side(family, std::get<1>(t), false);
3062 #else
3063  family.push_back(mesh.elem_ptr(std::get<0>(t)));
3064 #endif
3065 
3066  for (const auto & f : family)
3067  {
3068  const auto & conv = get_conversion(mesh.elem_ptr(f->id())->type());
3069 
3070  // Use the libmesh to exodus data structure map to get the proper sideset IDs
3071  // The data structure contains the "collapsed" contiguous ids
3072  elem_lists[std::get<2>(t)].push_back(libmesh_elem_num_to_exodus[f->id()]);
3073  side_lists[std::get<2>(t)].push_back(conv.get_inverse_shellface_map(std::get<1>(t)));
3074  }
3075  }
3076 
3077  std::vector<boundary_id_type> tmp;
3079  side_boundary_ids.insert(tmp.begin(), tmp.end());
3080  }
3081 
3082  // Add any empty-but-named side boundary ids
3083  for (const auto & pr : mesh.get_boundary_info().get_sideset_name_map())
3084  side_boundary_ids.insert(pr.first);
3085 
3086  // Write out the sideset names, but only if there is something to write
3087  if (side_boundary_ids.size() > 0)
3088  {
3089  NamesData names_table(side_boundary_ids.size(), MAX_STR_LENGTH);
3090 
3091  std::vector<exII::ex_set> sets(side_boundary_ids.size());
3092 
3093  // Loop over "side_boundary_ids" and "sets" simultaneously
3094  for (auto [i, it] = std::tuple{0u, side_boundary_ids.begin()}; i<sets.size(); ++i, ++it)
3095  {
3096  boundary_id_type ss_id = *it;
3098 
3099  sets[i].id = ss_id;
3100  sets[i].type = exII::EX_SIDE_SET;
3101  sets[i].num_distribution_factor = 0;
3102  sets[i].distribution_factor_list = nullptr;
3103 
3104  if (const auto elem_it = elem_lists.find(ss_id);
3105  elem_it == elem_lists.end())
3106  {
3107  sets[i].num_entry = 0;
3108  sets[i].entry_list = nullptr;
3109  sets[i].extra_list = nullptr;
3110  }
3111  else
3112  {
3113  sets[i].num_entry = elem_it->second.size();
3114  sets[i].entry_list = elem_it->second.data();
3115  sets[i].extra_list = libmesh_map_find(side_lists, ss_id).data();
3116  }
3117  }
3118 
3119  ex_err = exII::ex_put_sets(ex_id, side_boundary_ids.size(), sets.data());
3120  EX_CHECK_ERR(ex_err, "Error writing sidesets");
3121 
3122  ex_err = exII::ex_put_names(ex_id, exII::EX_SIDE_SET, names_table.get_char_star_star());
3123  EX_CHECK_ERR(ex_err, "Error writing sideset names");
3124  }
3125 }
3126 
3127 
3128 
3130 {
3131  LOG_SCOPE("write_nodesets()", "ExodusII_IO_Helper");
3132 
3133  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3134  return;
3135 
3136  // build_node_list() builds a sorted list of (node-id, bc-id) tuples
3137  // that is sorted by node-id, but we actually want it to be sorted
3138  // by bc-id, i.e. the second argument of the tuple.
3139  auto bc_tuples =
3141 
3142  // We use std::stable_sort() here so that the entries within a
3143  // single nodeset remain sorted in node-id order, but now the
3144  // smallest boundary id's nodes appear first in the list, followed
3145  // by the second smallest, etc. That is, we are purposely doing two
3146  // different sorts here, with the first one being within the
3147  // build_node_list() call itself.
3148  std::stable_sort(bc_tuples.begin(), bc_tuples.end(),
3149  [](const BoundaryInfo::NodeBCTuple & t1,
3150  const BoundaryInfo::NodeBCTuple & t2)
3151  { return std::get<1>(t1) < std::get<1>(t2); });
3152 
3153  std::vector<boundary_id_type> node_boundary_ids;
3154  mesh.get_boundary_info().build_node_boundary_ids(node_boundary_ids);
3155 
3156  // Add any empty-but-named node boundary ids
3157  for (const auto & pair : mesh.get_boundary_info().get_nodeset_name_map())
3158  {
3159  const boundary_id_type id = pair.first;
3160 
3161  if (std::find(node_boundary_ids.begin(),
3162  node_boundary_ids.end(), id)
3163  == node_boundary_ids.end())
3164  node_boundary_ids.push_back(id);
3165  }
3166 
3167  // Write out the nodeset names, but only if there is something to write
3168  if (node_boundary_ids.size() > 0)
3169  {
3170  NamesData names_table(node_boundary_ids.size(), MAX_STR_LENGTH);
3171 
3172  // Vectors to be filled and passed to exII::ex_put_concat_sets()
3173  // Use existing class members and avoid variable shadowing.
3174  nodeset_ids.clear();
3175  num_nodes_per_set.clear();
3176  num_node_df_per_set.clear();
3177  node_sets_node_index.clear();
3178  node_sets_node_list.clear();
3179 
3180  // Pre-allocate space
3181  nodeset_ids.reserve(node_boundary_ids.size());
3182  num_nodes_per_set.reserve(node_boundary_ids.size());
3183  num_node_df_per_set.resize(node_boundary_ids.size()); // all zeros
3184  node_sets_node_index.reserve(node_boundary_ids.size());
3185  node_sets_node_list.reserve(bc_tuples.size());
3186 
3187  // Assign entries to node_sets_node_list, keeping track of counts as we go.
3188  std::map<boundary_id_type, unsigned int> nodeset_counts;
3189  for (auto id : node_boundary_ids)
3190  nodeset_counts[id] = 0;
3191 
3192  for (const auto & t : bc_tuples)
3193  {
3194  const dof_id_type & node_id = std::get<0>(t) + 1; // Note: we use 1-based node ids in Exodus!
3195  const boundary_id_type & nodeset_id = std::get<1>(t);
3196  node_sets_node_list.push_back(node_id);
3197  nodeset_counts[nodeset_id] += 1;
3198  }
3199 
3200  // Fill in other indexing vectors needed by Exodus
3201  unsigned int running_sum = 0;
3202  for (const auto & pr : nodeset_counts)
3203  {
3204  nodeset_ids.push_back(pr.first);
3205  num_nodes_per_set.push_back(pr.second);
3206  node_sets_node_index.push_back(running_sum);
3207  names_table.push_back_entry(mesh.get_boundary_info().get_nodeset_name(pr.first));
3208  running_sum += pr.second;
3209  }
3210 
3211  // Fill in an exII::ex_set_specs object which can then be passed to
3212  // the ex_put_concat_sets() function.
3213  exII::ex_set_specs set_data = {};
3214  set_data.sets_ids = nodeset_ids.data();
3215  set_data.num_entries_per_set = num_nodes_per_set.data();
3216  set_data.num_dist_per_set = num_node_df_per_set.data(); // zeros
3217  set_data.sets_entry_index = node_sets_node_index.data();
3218  set_data.sets_dist_index = node_sets_node_index.data(); // dummy value
3219  set_data.sets_entry_list = node_sets_node_list.data();
3220 
3221  // Write all nodesets together.
3222  ex_err = exII::ex_put_concat_sets(ex_id, exII::EX_NODE_SET, &set_data);
3223  EX_CHECK_ERR(ex_err, "Error writing concatenated nodesets");
3224 
3225  // Write out the nodeset names
3226  ex_err = exII::ex_put_names(ex_id, exII::EX_NODE_SET, names_table.get_char_star_star());
3227  EX_CHECK_ERR(ex_err, "Error writing nodeset names");
3228  }
3229 }
3230 
3231 
3232 
3233 void ExodusII_IO_Helper::initialize_element_variables(std::vector<std::string> names,
3234  const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains)
3235 {
3236  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3237  return;
3238 
3239  // Quick return if there are no element variables to write
3240  if (names.size() == 0)
3241  return;
3242 
3243  // Be sure that variables in the file match what we are asking for
3244  if (num_elem_vars > 0)
3245  {
3246  this->check_existing_vars(ELEMENTAL, names, this->elem_var_names);
3247  return;
3248  }
3249 
3250  // Quick return if we have already called this function
3252  return;
3253 
3254  // Set the flag so we can skip this stuff on subsequent calls to
3255  // initialize_element_variables()
3256  _elem_vars_initialized = true;
3257 
3258  this->write_var_names(ELEMENTAL, names);
3259 
3260  // Use the truth table to indicate which subdomain/variable pairs are
3261  // active according to vars_active_subdomains.
3262  std::vector<int> truth_tab(num_elem_blk*num_elem_vars, 0);
3263  for (auto var_num : index_range(vars_active_subdomains))
3264  {
3265  // If the list of active subdomains is empty, it is interpreted as being
3266  // active on *all* subdomains.
3267  std::set<subdomain_id_type> current_set;
3268  if (vars_active_subdomains[var_num].empty())
3269  for (auto block_id : block_ids)
3270  current_set.insert(cast_int<subdomain_id_type>(block_id));
3271  else
3272  current_set = vars_active_subdomains[var_num];
3273 
3274  // Find index into the truth table for each id in current_set.
3275  for (auto block_id : current_set)
3276  {
3277  auto it = std::find(block_ids.begin(), block_ids.end(), block_id);
3278  libmesh_error_msg_if(it == block_ids.end(),
3279  "ExodusII_IO_Helper: block id " << block_id << " not found in block_ids.");
3280 
3281  std::size_t block_index =
3282  std::distance(block_ids.begin(), it);
3283 
3284  std::size_t truth_tab_index = block_index*num_elem_vars + var_num;
3285  truth_tab[truth_tab_index] = 1;
3286  }
3287  }
3288 
3289  ex_err = exII::ex_put_truth_table
3290  (ex_id,
3291  exII::EX_ELEM_BLOCK,
3292  num_elem_blk,
3293  num_elem_vars,
3294  truth_tab.data());
3295  EX_CHECK_ERR(ex_err, "Error writing element truth table.");
3296 }
3297 
3298 
3299 
3300 void ExodusII_IO_Helper::initialize_nodal_variables(std::vector<std::string> names)
3301 {
3302  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3303  return;
3304 
3305  // Quick return if there are no nodal variables to write
3306  if (names.size() == 0)
3307  return;
3308 
3309  // Quick return if we have already called this function
3311  return;
3312 
3313  // Be sure that variables in the file match what we are asking for
3314  if (num_nodal_vars > 0)
3315  {
3316  this->check_existing_vars(NODAL, names, this->nodal_var_names);
3317  return;
3318  }
3319 
3320  // Set the flag so we can skip the rest of this function on subsequent calls.
3321  _nodal_vars_initialized = true;
3322 
3323  this->write_var_names(NODAL, names);
3324 }
3325 
3326 
3327 
3328 void ExodusII_IO_Helper::initialize_global_variables(std::vector<std::string> names)
3329 {
3330  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3331  return;
3332 
3333  // Quick return if there are no global variables to write
3334  if (names.size() == 0)
3335  return;
3336 
3338  return;
3339 
3340  // Be sure that variables in the file match what we are asking for
3341  if (num_global_vars > 0)
3342  {
3343  this->check_existing_vars(GLOBAL, names, this->global_var_names);
3344  return;
3345  }
3346 
3347  _global_vars_initialized = true;
3348 
3349  this->write_var_names(GLOBAL, names);
3350 }
3351 
3352 
3353 
3355  std::vector<std::string> & names,
3356  std::vector<std::string> & names_from_file)
3357 {
3358  // There may already be global variables in the file (for example,
3359  // if we're appending) and in that case, we
3360  // 1.) Cannot initialize them again.
3361  // 2.) Should check to be sure that the global variable names are the same.
3362 
3363  // Fills up names_from_file for us
3364  this->read_var_names(type);
3365 
3366  // Both the number of variables and their names (up to the first
3367  // MAX_STR_LENGTH characters) must match for the names we are
3368  // planning to write and the names already in the file.
3369  bool match =
3370  std::equal(names.begin(), names.end(),
3371  names_from_file.begin(),
3372  [](const std::string & a,
3373  const std::string & b) -> bool
3374  {
3375  return a.compare(/*pos=*/0, /*len=*/MAX_STR_LENGTH, b) == 0;
3376  });
3377 
3378  if (!match)
3379  {
3380  libMesh::err << "Error! The Exodus file already contains the variables:" << std::endl;
3381  for (const auto & name : names_from_file)
3382  libMesh::err << name << std::endl;
3383 
3384  libMesh::err << "And you asked to write:" << std::endl;
3385  for (const auto & name : names)
3386  libMesh::err << name << std::endl;
3387 
3388  libmesh_error_msg("Cannot overwrite existing variables in Exodus II file.");
3389  }
3390 }
3391 
3392 
3393 
3395 {
3396  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3397  return;
3398 
3399  if (_single_precision)
3400  {
3401  float cast_time = float(time);
3402  ex_err = exII::ex_put_time(ex_id, timestep, &cast_time);
3403  }
3404  else
3405  {
3406  double cast_time = double(time);
3407  ex_err = exII::ex_put_time(ex_id, timestep, &cast_time);
3408  }
3409  EX_CHECK_ERR(ex_err, "Error writing timestep.");
3410 
3411  this->update();
3412 }
3413 
3414 
3415 
3416 void
3418 {
3419  LOG_SCOPE("write_elemsets()", "ExodusII_IO_Helper");
3420 
3421  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3422  return;
3423 
3424  // TODO: Add support for named elemsets
3425  // NamesData names_table(elemsets.size(), MAX_STR_LENGTH);
3426 
3427  // We only need to write elemsets if the Mesh has an extra elem
3428  // integer called "elemset_code" defined on it.
3429  if (mesh.has_elem_integer("elemset_code"))
3430  {
3431  std::map<elemset_id_type, std::vector<int>> exodus_elemsets;
3432 
3433  unsigned int elemset_index =
3434  mesh.get_elem_integer_index("elemset_code");
3435 
3436  // Catch ids returned from MeshBase::get_elemsets() calls
3437  MeshBase::elemset_type set_ids;
3438  for (const auto & elem : mesh.element_ptr_range())
3439  {
3440  dof_id_type elemset_code =
3441  elem->get_extra_integer(elemset_index);
3442 
3443  // Look up which element set ids (if any) this elemset_code corresponds to.
3444  mesh.get_elemsets(elemset_code, set_ids);
3445 
3446  // Debugging
3447  // libMesh::out << "elemset_code = " << elemset_code << std::endl;
3448  // for (const auto & set_id : set_ids)
3449  // libMesh::out << set_id << " ";
3450  // libMesh::out << std::endl;
3451 
3452  // Store this Elem id in every set to which it belongs.
3453  for (const auto & set_id : set_ids)
3454  exodus_elemsets[set_id].push_back(libmesh_elem_num_to_exodus[elem->id()]);
3455  }
3456 
3457  // Debugging: print contents of exodus_elemsets map
3458  // for (const auto & [set_id, elem_ids] : exodus_elemsets)
3459  // {
3460  // libMesh::out << "elemset " << set_id << ": ";
3461  // for (const auto & elem_id : elem_ids)
3462  // libMesh::out << elem_id << " ";
3463  // libMesh::out << std::endl;
3464  // }
3465 
3466  // Only continue if we actually had some elements in sets
3467  if (!exodus_elemsets.empty())
3468  {
3469  // Reserve space, loop over newly-created map, construct
3470  // exII::ex_set objects to be passed to exII::ex_put_sets(). Note:
3471  // we do non-const iteration since Exodus requires non-const pointers
3472  // to be passed to its APIs.
3473  std::vector<exII::ex_set> sets;
3474  sets.reserve(exodus_elemsets.size());
3475 
3476  for (auto & [elem_set_id, ids_vec] : exodus_elemsets)
3477  {
3478  // TODO: Add support for named elemsets
3479  // names_table.push_back_entry(mesh.get_elemset_name(elem_set_id));
3480 
3481  exII::ex_set & current_set = sets.emplace_back();
3482  current_set.id = elem_set_id;
3483  current_set.type = exII::EX_ELEM_SET;
3484  current_set.num_entry = ids_vec.size();
3485  current_set.num_distribution_factor = 0;
3486  current_set.entry_list = ids_vec.data();
3487  current_set.extra_list = nullptr; // extra_list is used for sidesets, not needed for elemsets
3488  current_set.distribution_factor_list = nullptr; // not used for elemsets
3489  }
3490 
3491  // Sanity check: make sure the number of elemsets we already wrote to the header
3492  // matches the number of elemsets we just constructed by looping over the Mesh.
3493  libmesh_assert_msg(num_elem_sets == cast_int<int>(exodus_elemsets.size()),
3494  "Mesh has " << exodus_elemsets.size()
3495  << " elemsets, but header was written with num_elem_sets == " << num_elem_sets);
3496  libmesh_assert_msg(num_elem_sets == cast_int<int>(mesh.n_elemsets()),
3497  "mesh.n_elemsets() == " << mesh.n_elemsets()
3498  << ", but header was written with num_elem_sets == " << num_elem_sets);
3499 
3500  ex_err = exII::ex_put_sets(ex_id, exodus_elemsets.size(), sets.data());
3501  EX_CHECK_ERR(ex_err, "Error writing elemsets");
3502 
3503  // TODO: Add support for named elemsets
3504  // ex_err = exII::ex_put_names(ex_id, exII::EX_ELEM_SET, names_table.get_char_star_star());
3505  // EX_CHECK_ERR(ex_err, "Error writing elemset names");
3506  } // end if (!exodus_elemsets.empty())
3507  } // end if (mesh.has_elem_integer("elemset_code"))
3508 }
3509 
3510 
3511 
3512 void
3515  int timestep,
3516  const std::vector<std::string> & var_names,
3517  const std::vector<std::set<boundary_id_type>> & side_ids,
3518  const std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
3519 {
3520  LOG_SCOPE("write_sideset_data()", "ExodusII_IO_Helper");
3521 
3522  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3523  return;
3524 
3525  // Write the sideset variable names to file. This function should
3526  // only be called once for SIDESET variables, repeated calls to
3527  // write_var_names overwrites/changes the order of names that were
3528  // there previously, and will mess up any data that has already been
3529  // written.
3530  this->write_var_names(SIDESET, var_names);
3531 
3532  // I hope that we are allowed to call read_sideset_info() even
3533  // though we are in the middle of writing? It seems to work provided
3534  // that you have already written the mesh itself... read_sideset_info()
3535  // fills in the following data members:
3536  // .) num_side_sets
3537  // .) ss_ids
3538  this->read_sideset_info();
3539 
3540  // Write "truth" table for sideset variables. The function
3541  // exII::ex_put_variable_param() must be called before
3542  // exII::ex_put_truth_table(). For us, this happens during the call
3543  // to ExodusII_IO_Helper::write_var_names(). sset_var_tab is a logically
3544  // (num_side_sets x num_sset_var) integer array of 0s and 1s
3545  // indicating which sidesets a given sideset variable is defined on.
3546  std::vector<int> sset_var_tab(num_side_sets * var_names.size());
3547 
3548  // We now call read_sideset() once per sideset and write any sideset
3549  // variable values which are defined there.
3550  int offset=0;
3551  for (int ss=0; ss<num_side_sets; ++ss)
3552  {
3553  // We don't know num_sides_per_set for each set until we call
3554  // read_sideset(). The values for each sideset are stored (using
3555  // the offsets) into the 'elem_list' and 'side_list' arrays of
3556  // this class.
3557  offset += (ss > 0 ? num_sides_per_set[ss-1] : 0);
3558  this->read_sideset(ss, offset);
3559 
3560  // For each variable in var_names, write the values for the
3561  // current sideset, if any.
3562  for (auto var : index_range(var_names))
3563  {
3564  // If this var has no values on this sideset, go to the next one.
3565  if (!side_ids[var].count(ss_ids[ss]))
3566  continue;
3567 
3568  // Otherwise, fill in this entry of the sideset truth table.
3569  sset_var_tab[ss*var_names.size() + var] = 1;
3570 
3571  // Data vector that will eventually be passed to exII::ex_put_var().
3572  std::vector<Real> sset_var_vals(num_sides_per_set[ss]);
3573 
3574  // Get reference to the BCTuple -> Real map for this variable.
3575  const auto & data_map = bc_vals[var];
3576 
3577  // Loop over elem_list, side_list entries in current sideset.
3578  for (int i=0; i<num_sides_per_set[ss]; ++i)
3579  {
3580  // Get elem_id and side_id from the respective lists that
3581  // are filled in by calling read_sideset().
3582  //
3583  // Note: these are Exodus-specific ids, so we have to convert them
3584  // to libmesh ids, as that is what will be in the bc_tuples.
3585  //
3586  // TODO: we should probably consult the exodus_elem_num_to_libmesh
3587  // mapping in order to figure out which libmesh element id 'elem_id'
3588  // actually corresponds to here, instead of just assuming it will be
3589  // off by one. Unfortunately that data structure does not seem to
3590  // be used at the moment. If we assume that write_sideset_data() is
3591  // always called following write(), then this should be a fairly safe
3592  // assumption...
3593  dof_id_type elem_id = elem_list[i + offset] - 1;
3594  unsigned int side_id = side_list[i + offset] - 1;
3595 
3596  // Sanity check: make sure that the "off by one"
3597  // assumption we used above to set 'elem_id' is valid.
3598  libmesh_error_msg_if
3599  (libmesh_map_find(libmesh_elem_num_to_exodus, cast_int<int>(elem_id)) !=
3600  cast_int<dof_id_type>(elem_list[i + offset]),
3601  "Error mapping Exodus elem id to libmesh elem id.");
3602 
3603  // Map from Exodus side ids to libmesh side ids.
3604  const auto & conv = get_conversion(mesh.elem_ptr(elem_id)->type());
3605 
3606  // Map from Exodus side ids to libmesh side ids.
3607  unsigned int converted_side_id = conv.get_side_map(side_id);
3608 
3609  // Construct a key so we can quickly see whether there is any
3610  // data for this variable in the map.
3611  BoundaryInfo::BCTuple key = std::make_tuple
3612  (elem_id,
3613  converted_side_id,
3614  ss_ids[ss]);
3615 
3616  // Find the data for this (elem,side,id) tuple. Throw an
3617  // error if not found. Then store value in vector which
3618  // will be passed to Exodus.
3619  sset_var_vals[i] = libmesh_map_find(data_map, key);
3620  } // end for (i)
3621 
3622  // As far as I can tell, there is no "concat" version of writing
3623  // sideset data, you have to call ex_put_sset_var() once per (variable,
3624  // sideset) pair.
3625  if (sset_var_vals.size() > 0)
3626  {
3627  ex_err = exII::ex_put_var
3628  (ex_id,
3629  timestep,
3630  exII::EX_SIDE_SET,
3631  var + 1, // 1-based variable index of current variable
3632  ss_ids[ss],
3633  num_sides_per_set[ss],
3634  MappedOutputVector(sset_var_vals, _single_precision).data());
3635  EX_CHECK_ERR(ex_err, "Error writing sideset vars.");
3636  }
3637  } // end for (var)
3638  } // end for (ss)
3639 
3640  // Finally, write the sideset truth table.
3641  ex_err =
3642  exII::ex_put_truth_table(ex_id,
3643  exII::EX_SIDE_SET,
3644  num_side_sets,
3645  cast_int<int>(var_names.size()),
3646  sset_var_tab.data());
3647  EX_CHECK_ERR(ex_err, "Error writing sideset var truth table.");
3648 }
3649 
3650 
3651 
3652 void
3655  int timestep,
3656  std::vector<std::string> & var_names,
3657  std::vector<std::set<boundary_id_type>> & side_ids,
3658  std::vector<std::map<BoundaryInfo::BCTuple, Real>> & bc_vals)
3659 {
3660  LOG_SCOPE("read_sideset_data()", "ExodusII_IO_Helper");
3661 
3662  // This reads the sideset variable names into the local
3663  // sideset_var_names data structure.
3664  this->read_var_names(SIDESET);
3665 
3666  if (num_sideset_vars)
3667  {
3668  // Read the sideset data truth table
3669  std::vector<int> sset_var_tab(num_side_sets * num_sideset_vars);
3670  ex_err = exII::ex_get_truth_table
3671  (ex_id,
3672  exII::EX_SIDE_SET,
3673  num_side_sets,
3675  sset_var_tab.data());
3676  EX_CHECK_ERR(ex_err, "Error reading sideset variable truth table.");
3677 
3678  // Set up/allocate space in incoming data structures.
3679  var_names = sideset_var_names;
3680  side_ids.resize(num_sideset_vars);
3681  bc_vals.resize(num_sideset_vars);
3682 
3683  // Read the sideset data.
3684  //
3685  // Note: we assume that read_sideset() has already been called
3686  // for each sideset, so the required values in elem_list and
3687  // side_list are already present.
3688  //
3689  // TODO: As a future optimization, we could read only the values
3690  // requested by the user by looking at the input parameter
3691  // var_names and checking whether it already has entries in
3692  // it. We could do the same thing with the input side_ids
3693  // container and only read values for requested sidesets.
3694  int offset=0;
3695  for (int ss=0; ss<num_side_sets; ++ss)
3696  {
3697  offset += (ss > 0 ? num_sides_per_set[ss-1] : 0);
3698  for (int var=0; var<num_sideset_vars; ++var)
3699  {
3700  int is_present = sset_var_tab[num_sideset_vars*ss + var];
3701 
3702  if (is_present)
3703  {
3704  // Record the fact that this variable is defined on this sideset.
3705  side_ids[var].insert(ss_ids[ss]);
3706 
3707  // Note: the assumption here is that a previous call
3708  // to this->read_sideset_info() has already set the
3709  // values of num_sides_per_set, so we just use those values here.
3710  std::vector<Real> sset_var_vals(num_sides_per_set[ss]);
3711  ex_err = exII::ex_get_var
3712  (ex_id,
3713  timestep,
3714  exII::EX_SIDE_SET,
3715  var + 1, // 1-based sideset variable index!
3716  ss_ids[ss],
3717  num_sides_per_set[ss],
3718  MappedInputVector(sset_var_vals, _single_precision).data());
3719  EX_CHECK_ERR(ex_err, "Error reading sideset variable.");
3720 
3721  for (int i=0; i<num_sides_per_set[ss]; ++i)
3722  {
3723  dof_id_type exodus_elem_id = elem_list[i + offset];
3724  unsigned int exodus_side_id = side_list[i + offset];
3725 
3726  // FIXME: We should use exodus_elem_num_to_libmesh for this,
3727  // but it apparently is never set up, so just
3728  // subtract 1 from the Exodus elem id.
3729  dof_id_type converted_elem_id = exodus_elem_id - 1;
3730 
3731  // Map Exodus side id to libmesh side id.
3732  // Map from Exodus side ids to libmesh side ids.
3733  const auto & conv = get_conversion(mesh.elem_ptr(converted_elem_id)->type());
3734 
3735  // Map from Exodus side id to libmesh side id.
3736  // Note: the mapping is defined on 0-based indices, so subtract
3737  // 1 before doing the mapping.
3738  unsigned int converted_side_id = conv.get_side_map(exodus_side_id - 1);
3739 
3740  // Make a BCTuple key from the converted information.
3741  BoundaryInfo::BCTuple key = std::make_tuple
3742  (converted_elem_id,
3743  converted_side_id,
3744  ss_ids[ss]);
3745 
3746  // Store (elem, side, b_id) tuples in bc_vals[var]
3747  bc_vals[var].emplace(key, sset_var_vals[i]);
3748  } // end for (i)
3749  } // end if (present)
3750  } // end for (var)
3751  } // end for (ss)
3752  } // end if (num_sideset_vars)
3753 }
3754 
3755 
3756 void
3759  std::map<BoundaryInfo::BCTuple, unsigned int> & bc_array_indices)
3760 {
3761  // Clear any existing data, we are going to build this data structure from scratch
3762  bc_array_indices.clear();
3763 
3764  // Store the sideset data array indices.
3765  //
3766  // Note: we assume that read_sideset() has already been called
3767  // for each sideset, so the required values in elem_list and
3768  // side_list are already present.
3769  int offset=0;
3770  for (int ss=0; ss<num_side_sets; ++ss)
3771  {
3772  offset += (ss > 0 ? num_sides_per_set[ss-1] : 0);
3773  for (int i=0; i<num_sides_per_set[ss]; ++i)
3774  {
3775  dof_id_type exodus_elem_id = elem_list[i + offset];
3776  unsigned int exodus_side_id = side_list[i + offset];
3777 
3778  // FIXME: We should use exodus_elem_num_to_libmesh for this,
3779  // but it apparently is never set up, so just
3780  // subtract 1 from the Exodus elem id.
3781  dof_id_type converted_elem_id = exodus_elem_id - 1;
3782 
3783  // Conversion operator for this Elem type
3784  const auto & conv = get_conversion(mesh.elem_ptr(converted_elem_id)->type());
3785 
3786  // Map from Exodus side id to libmesh side id.
3787  // Note: the mapping is defined on 0-based indices, so subtract
3788  // 1 before doing the mapping.
3789  unsigned int converted_side_id = conv.get_side_map(exodus_side_id - 1);
3790 
3791  // Make a BCTuple key from the converted information.
3792  BoundaryInfo::BCTuple key = std::make_tuple
3793  (converted_elem_id,
3794  converted_side_id,
3795  ss_ids[ss]);
3796 
3797  // Store (elem, side, b_id) tuple with corresponding array index
3798  bc_array_indices.emplace(key, cast_int<unsigned int>(i));
3799  } // end for (i)
3800  } // end for (ss)
3801 }
3802 
3803 
3804 
3806 write_nodeset_data (int timestep,
3807  const std::vector<std::string> & var_names,
3808  const std::vector<std::set<boundary_id_type>> & node_boundary_ids,
3809  const std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
3810 {
3811  LOG_SCOPE("write_nodeset_data()", "ExodusII_IO_Helper");
3812 
3813  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3814  return;
3815 
3816  // Write the nodeset variable names to file. This function should
3817  // only be called once for NODESET variables, repeated calls to
3818  // write_var_names() overwrites/changes the order of names that were
3819  // there previously, and will mess up any data that has already been
3820  // written.
3821  this->write_var_names(NODESET, var_names);
3822 
3823  // For all nodesets, reads and fills in the arrays:
3824  // nodeset_ids
3825  // num_nodes_per_set
3826  // node_sets_node_index - starting index for each nodeset in the node_sets_node_list vector
3827  // node_sets_node_list
3828  // Note: we need these arrays so that we know what data to write
3829  this->read_all_nodesets();
3830 
3831  // The "truth" table for nodeset variables. nset_var_tab is a
3832  // logically (num_node_sets x num_nset_var) integer array of 0s and
3833  // 1s indicating which nodesets a given nodeset variable is defined
3834  // on.
3835  std::vector<int> nset_var_tab(num_node_sets * var_names.size());
3836 
3837  for (int ns=0; ns<num_node_sets; ++ns)
3838  {
3839  // The offset into the node_sets_node_list for the current nodeset
3840  int offset = node_sets_node_index[ns];
3841 
3842  // For each variable in var_names, write the values for the
3843  // current nodeset, if any.
3844  for (auto var : index_range(var_names))
3845  {
3846  // If this var has no values on this nodeset, go to the next one.
3847  if (!node_boundary_ids[var].count(nodeset_ids[ns]))
3848  continue;
3849 
3850  // Otherwise, fill in this entry of the nodeset truth table.
3851  nset_var_tab[ns*var_names.size() + var] = 1;
3852 
3853  // Data vector that will eventually be passed to exII::ex_put_var().
3854  std::vector<Real> nset_var_vals(num_nodes_per_set[ns]);
3855 
3856  // Get reference to the NodeBCTuple -> Real map for this variable.
3857  const auto & data_map = bc_vals[var];
3858 
3859  // Loop over entries in current nodeset.
3860  for (int i=0; i<num_nodes_per_set[ns]; ++i)
3861  {
3862  // Here we convert Exodus node ids to libMesh node ids by
3863  // subtracting 1. We should probably use the
3864  // exodus_node_num_to_libmesh data structure for this, but
3865  // I don't think it is set up at the time when
3866  // write_nodeset_data() would normally be called.
3867  dof_id_type libmesh_node_id = node_sets_node_list[i + offset] - 1;
3868 
3869  // Construct a key to look up values in data_map.
3871  std::make_tuple(libmesh_node_id, nodeset_ids[ns]);
3872 
3873  // We require that the user provided either no values for
3874  // this (var, nodeset) combination (in which case we don't
3875  // reach this point) or a value for _every_ node in this
3876  // nodeset for this var, so we use the libmesh_map_find()
3877  // macro to check for this.
3878  nset_var_vals[i] = libmesh_map_find(data_map, key);
3879  } // end for (node in nodeset[ns])
3880 
3881  // Write nodeset values to Exodus file
3882  if (nset_var_vals.size() > 0)
3883  {
3884  ex_err = exII::ex_put_var
3885  (ex_id,
3886  timestep,
3887  exII::EX_NODE_SET,
3888  var + 1, // 1-based variable index of current variable
3889  nodeset_ids[ns],
3890  num_nodes_per_set[ns],
3891  MappedOutputVector(nset_var_vals, _single_precision).data());
3892  EX_CHECK_ERR(ex_err, "Error writing nodeset vars.");
3893  }
3894  } // end for (var in var_names)
3895  } // end for (ns)
3896 
3897  // Finally, write the nodeset truth table.
3898  ex_err =
3899  exII::ex_put_truth_table(ex_id,
3900  exII::EX_NODE_SET,
3901  num_node_sets,
3902  cast_int<int>(var_names.size()),
3903  nset_var_tab.data());
3904  EX_CHECK_ERR(ex_err, "Error writing nodeset var truth table.");
3905 }
3906 
3907 
3908 
3909 void
3911 write_elemset_data (int timestep,
3912  const std::vector<std::string> & var_names,
3913  const std::vector<std::set<elemset_id_type>> & elemset_ids_in,
3914  const std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
3915 {
3916  LOG_SCOPE("write_elemset_data()", "ExodusII_IO_Helper");
3917 
3918  if ((_run_only_on_proc0) && (this->processor_id() != 0))
3919  return;
3920 
3921  // Write the elemset variable names to file. This function should
3922  // only be called once for ELEMSET variables, repeated calls to
3923  // write_var_names() overwrites/changes the order of names that were
3924  // there previously, and will mess up any data that has already been
3925  // written.
3926  this->write_var_names(ELEMSET, var_names);
3927 
3928  // We now call the API to read the elemset info even though we are
3929  // in the middle of writing. This is a bit counter-intuitive, but it
3930  // seems to work provided that you have already written the mesh
3931  // itself... read_elemset_info() fills in the following data
3932  // members:
3933  // .) id_to_elemset_names
3934  // .) num_elems_per_set
3935  // .) num_elem_df_per_set
3936  // .) elemset_list
3937  // .) elemset_id_list
3938  // .) id_to_elemset_names
3939  this->read_elemset_info();
3940 
3941  // The "truth" table for elemset variables. elemset_var_tab is a
3942  // logically (num_elem_sets x num_elemset_vars) integer array of 0s and
3943  // 1s indicating which elemsets a given elemset variable is defined
3944  // on.
3945  std::vector<int> elemset_var_tab(num_elem_sets * var_names.size());
3946 
3947  int offset=0;
3948  for (int es=0; es<num_elem_sets; ++es)
3949  {
3950  // Debugging
3951  // libMesh::out << "Writing elemset variable values for elemset "
3952  // << es << ", elemset_id = " << elemset_ids[es]
3953  // << std::endl;
3954 
3955  // We know num_elems_per_set because we called read_elemset_info() above.
3956  offset += (es > 0 ? num_elems_per_set[es-1] : 0);
3957  this->read_elemset(es, offset);
3958 
3959  // For each variable in var_names, write the values for the
3960  // current elemset, if any.
3961  for (auto var : index_range(var_names))
3962  {
3963  // Debugging
3964  // libMesh::out << "Writing elemset variable values for var " << var << std::endl;
3965 
3966  // If this var has no values on this elemset, go to the next one.
3967  if (!elemset_ids_in[var].count(elemset_ids[es]))
3968  continue;
3969 
3970  // Otherwise, fill in this entry of the nodeset truth table.
3971  elemset_var_tab[es*var_names.size() + var] = 1;
3972 
3973  // Data vector that will eventually be passed to exII::ex_put_var().
3974  std::vector<Real> elemset_var_vals(num_elems_per_set[es]);
3975 
3976  // Get reference to the (elem_id, elemset_id) -> Real map for this variable.
3977  const auto & data_map = elemset_vals[var];
3978 
3979  // Loop over entries in current elemset.
3980  for (int i=0; i<num_elems_per_set[es]; ++i)
3981  {
3982  // Here we convert Exodus elem ids to libMesh node ids
3983  // simply by subtracting 1. We should probably use the
3984  // exodus_elem_num_to_libmesh data structure for this,
3985  // but I don't think it is set up at the time when this
3986  // function is normally called.
3987  dof_id_type libmesh_elem_id = elemset_list[i + offset] - 1;
3988 
3989  // Construct a key to look up values in data_map.
3990  std::pair<dof_id_type, elemset_id_type> key =
3991  std::make_pair(libmesh_elem_id, elemset_ids[es]);
3992 
3993  // Debugging:
3994  // libMesh::out << "Searching for key = (" << key.first << ", " << key.second << ")" << std::endl;
3995 
3996  // We require that the user provided either no values for
3997  // this (var, elemset) combination (in which case we don't
3998  // reach this point) or a value for _every_ elem in this
3999  // elemset for this var, so we use the libmesh_map_find()
4000  // macro to check for this.
4001  elemset_var_vals[i] = libmesh_map_find(data_map, key);
4002  } // end for (node in nodeset[ns])
4003 
4004  // Write elemset values to Exodus file
4005  if (elemset_var_vals.size() > 0)
4006  {
4007  ex_err = exII::ex_put_var
4008  (ex_id,
4009  timestep,
4010  exII::EX_ELEM_SET,
4011  var + 1, // 1-based variable index of current variable
4012  elemset_ids[es],
4013  num_elems_per_set[es],
4014  MappedOutputVector(elemset_var_vals, _single_precision).data());
4015  EX_CHECK_ERR(ex_err, "Error writing elemset vars.");
4016  }
4017  } // end for (var in var_names)
4018  } // end for (ns)
4019 
4020  // Finally, write the elemset truth table to file.
4021  ex_err =
4022  exII::ex_put_truth_table(ex_id,
4023  exII::EX_ELEM_SET, // exII::ex_entity_type
4024  num_elem_sets,
4025  cast_int<int>(var_names.size()),
4026  elemset_var_tab.data());
4027  EX_CHECK_ERR(ex_err, "Error writing elemset var truth table.");
4028 }
4029 
4030 
4031 
4032 void
4034 read_elemset_data (int timestep,
4035  std::vector<std::string> & var_names,
4036  std::vector<std::set<elemset_id_type>> & elemset_ids_in,
4037  std::vector<std::map<std::pair<dof_id_type, elemset_id_type>, Real>> & elemset_vals)
4038 {
4039  LOG_SCOPE("read_elemset_data()", "ExodusII_IO_Helper");
4040 
4041  // This reads the elemset variable names into the local
4042  // elemset_var_names data structure.
4043  this->read_var_names(ELEMSET);
4044 
4045  // Debugging
4046  // libMesh::out << "elmeset variable names:" << std::endl;
4047  // for (const auto & name : elemset_var_names)
4048  // libMesh::out << name << " ";
4049  // libMesh::out << std::endl;
4050 
4051  if (num_elemset_vars)
4052  {
4053  // Debugging
4054  // std::cout << "Reading " << num_elem_sets
4055  // << " elemsets and " << num_elemset_vars
4056  // << " elemset variables." << std::endl;
4057 
4058  // Read the elemset data truth table.
4059  std::vector<int> elemset_var_tab(num_elem_sets * num_elemset_vars);
4060  exII::ex_get_truth_table(ex_id,
4061  exII::EX_ELEM_SET, // exII::ex_entity_type
4062  num_elem_sets,
4064  elemset_var_tab.data());
4065  EX_CHECK_ERR(ex_err, "Error reading elemset variable truth table.");
4066 
4067  // Debugging
4068  // libMesh::out << "Elemset variable truth table:" << std::endl;
4069  // for (const auto & val : elemset_var_tab)
4070  // libMesh::out << val << " ";
4071  // libMesh::out << std::endl;
4072 
4073  // Debugging
4074  // for (auto i : make_range(num_elem_sets))
4075  // {
4076  // for (auto j : make_range(num_elemset_vars))
4077  // libMesh::out << elemset_var_tab[num_elemset_vars*i + j] << " ";
4078  // libMesh::out << std::endl;
4079  // }
4080 
4081  // Set up/allocate space in incoming data structures. All vectors are
4082  // num_elemset_vars in length.
4083  var_names = elemset_var_names;
4084  elemset_ids_in.resize(num_elemset_vars);
4085  elemset_vals.resize(num_elemset_vars);
4086 
4087  // Read the elemset data
4088  int offset=0;
4089  for (int es=0; es<num_elem_sets; ++es)
4090  {
4091  offset += (es > 0 ? num_elems_per_set[es-1] : 0);
4092  for (int var=0; var<num_elemset_vars; ++var)
4093  {
4094  int is_present = elemset_var_tab[num_elemset_vars*es + var];
4095 
4096  if (is_present)
4097  {
4098  // Debugging
4099  // libMesh::out << "Variable " << var << " is present on elemset " << es << std::endl;
4100 
4101  // Record the fact that this variable is defined on this elemset.
4102  elemset_ids_in[var].insert(elemset_ids[es]);
4103 
4104  // Note: the assumption here is that a previous call
4105  // to this->read_elemset_info() has already set the
4106  // values of num_elems_per_set, so we just use those values here.
4107  std::vector<Real> elemset_var_vals(num_elems_per_set[es]);
4108  ex_err = exII::ex_get_var
4109  (ex_id,
4110  timestep,
4111  exII::EX_ELEM_SET, // exII::ex_entity_type
4112  var + 1, // 1-based sideset variable index!
4113  elemset_ids[es],
4114  num_elems_per_set[es],
4115  MappedInputVector(elemset_var_vals, _single_precision).data());
4116  EX_CHECK_ERR(ex_err, "Error reading elemset variable.");
4117 
4118  for (int i=0; i<num_elems_per_set[es]; ++i)
4119  {
4120  dof_id_type exodus_elem_id = elemset_list[i + offset];
4121 
4122  // FIXME: We should use exodus_elem_num_to_libmesh for this,
4123  // but it apparently is never set up, so just
4124  // subtract 1 from the Exodus elem id.
4125  dof_id_type converted_elem_id = exodus_elem_id - 1;
4126 
4127  // Make key based on the elem and set ids
4128  auto key = std::make_pair(converted_elem_id,
4129  static_cast<elemset_id_type>(elemset_ids[es]));
4130 
4131  // Store value in the map
4132  elemset_vals[var].emplace(key, elemset_var_vals[i]);
4133  } // end for (i)
4134  } // end if (present)
4135  } // end for (var)
4136  } // end for (es)
4137  } // end if (num_elemset_vars)
4138 }
4139 
4140 
4141 
4143 get_elemset_data_indices (std::map<std::pair<dof_id_type, elemset_id_type>, unsigned int> & elemset_array_indices)
4144 {
4145  // Clear existing data, we are going to build these data structures from scratch
4146  elemset_array_indices.clear();
4147 
4148  // Read the elemset data.
4149  //
4150  // Note: we assume that the functions
4151  // 1.) this->read_elemset_info() and
4152  // 2.) this->read_elemset()
4153  // have already been called, so that we already know e.g. how
4154  // many elems are in each set, their ids, etc.
4155  int offset=0;
4156  for (int es=0; es<num_elem_sets; ++es)
4157  {
4158  offset += (es > 0 ? num_elems_per_set[es-1] : 0);
4159 
4160  // Note: we don't actually call exII::ex_get_var() here because
4161  // we don't need the values. We only need the indices into that vector
4162  // for each (elem_id, elemset_id) tuple.
4163  for (int i=0; i<num_elems_per_set[es]; ++i)
4164  {
4165  dof_id_type exodus_elem_id = elemset_list[i + offset];
4166 
4167  // FIXME: We should use exodus_elem_num_to_libmesh for this,
4168  // but it apparently is never set up, so just
4169  // subtract 1 from the Exodus elem id.
4170  dof_id_type converted_elem_id = exodus_elem_id - 1;
4171 
4172  // Make key based on the elem and set ids
4173  // Make a NodeBCTuple key from the converted information.
4174  auto key = std::make_pair(converted_elem_id,
4175  static_cast<elemset_id_type>(elemset_ids[es]));
4176 
4177  // Store the array index of this (node, b_id) tuple
4178  elemset_array_indices.emplace(key, cast_int<unsigned int>(i));
4179  } // end for (i)
4180  } // end for (es)
4181 }
4182 
4183 
4184 
4186 read_nodeset_data (int timestep,
4187  std::vector<std::string> & var_names,
4188  std::vector<std::set<boundary_id_type>> & node_boundary_ids,
4189  std::vector<std::map<BoundaryInfo::NodeBCTuple, Real>> & bc_vals)
4190 {
4191  LOG_SCOPE("read_nodeset_data()", "ExodusII_IO_Helper");
4192 
4193  // This reads the sideset variable names into the local
4194  // sideset_var_names data structure.
4195  this->read_var_names(NODESET);
4196 
4197  if (num_nodeset_vars)
4198  {
4199  // Read the nodeset data truth table
4200  std::vector<int> nset_var_tab(num_node_sets * num_nodeset_vars);
4201  ex_err = exII::ex_get_truth_table
4202  (ex_id,
4203  exII::EX_NODE_SET,
4204  num_node_sets,
4206  nset_var_tab.data());
4207  EX_CHECK_ERR(ex_err, "Error reading nodeset variable truth table.");
4208 
4209  // Set up/allocate space in incoming data structures.
4210  var_names = nodeset_var_names;
4211  node_boundary_ids.resize(num_nodeset_vars);
4212  bc_vals.resize(num_nodeset_vars);
4213 
4214  // Read the nodeset data.
4215  //
4216  // Note: we assume that the functions
4217  // 1.) this->read_nodeset_info() and
4218  // 2.) this->read_all_nodesets()
4219  // have already been called, so that we already know e.g. how
4220  // many nodes are in each set, their ids, etc.
4221  //
4222  // TODO: As a future optimization, we could read only the values
4223  // requested by the user by looking at the input parameter
4224  // var_names and checking whether it already has entries in
4225  // it.
4226  int offset=0;
4227  for (int ns=0; ns<num_node_sets; ++ns)
4228  {
4229  offset += (ns > 0 ? num_nodes_per_set[ns-1] : 0);
4230  for (int var=0; var<num_nodeset_vars; ++var)
4231  {
4232  int is_present = nset_var_tab[num_nodeset_vars*ns + var];
4233 
4234  if (is_present)
4235  {
4236  // Record the fact that this variable is defined on this nodeset.
4237  node_boundary_ids[var].insert(nodeset_ids[ns]);
4238 
4239  // Note: the assumption here is that a previous call
4240  // to this->read_nodeset_info() has already set the
4241  // values of num_nodes_per_set, so we just use those values here.
4242  std::vector<Real> nset_var_vals(num_nodes_per_set[ns]);
4243  ex_err = exII::ex_get_var
4244  (ex_id,
4245  timestep,
4246  exII::EX_NODE_SET,
4247  var + 1, // 1-based nodeset variable index!
4248  nodeset_ids[ns],
4249  num_nodes_per_set[ns],
4250  MappedInputVector(nset_var_vals, _single_precision).data());
4251  EX_CHECK_ERR(ex_err, "Error reading nodeset variable.");
4252 
4253  for (int i=0; i<num_nodes_per_set[ns]; ++i)
4254  {
4255  // The read_all_nodesets() function now reads all the node ids into the
4256  // node_sets_node_list vector, which is of length "total_nodes_in_all_sets"
4257  // The old read_nodset() function is no longer called as far as I can tell,
4258  // and should probably be removed? The "offset" that we are using only
4259  // depends on the current nodeset index and the num_nodes_per_set vector,
4260  // which gets filled in by the call to read_all_nodesets().
4261  dof_id_type exodus_node_id = node_sets_node_list[i + offset];
4262 
4263  // FIXME: We should use exodus_node_num_to_libmesh for this,
4264  // but it apparently is never set up, so just
4265  // subtract 1 from the Exodus node id.
4266  dof_id_type converted_node_id = exodus_node_id - 1;
4267 
4268  // Make a NodeBCTuple key from the converted information.
4269  BoundaryInfo::NodeBCTuple key = std::make_tuple
4270  (converted_node_id, nodeset_ids[ns]);
4271 
4272  // Store (node, b_id) tuples in bc_vals[var]
4273  bc_vals[var].emplace(key, nset_var_vals[i]);
4274  } // end for (i)
4275  } // end if (present)
4276  } // end for (var)
4277  } // end for (ns)
4278  } // end if (num_nodeset_vars)
4279 }
4280 
4281 
4282 
4283 void
4285 get_nodeset_data_indices (std::map<BoundaryInfo::NodeBCTuple, unsigned int> & bc_array_indices)
4286 {
4287  // Clear existing data, we are going to build these data structures from scratch
4288  bc_array_indices.clear();
4289 
4290  // Read the nodeset data.
4291  //
4292  // Note: we assume that the functions
4293  // 1.) this->read_nodeset_info() and
4294  // 2.) this->read_all_nodesets()
4295  // have already been called, so that we already know e.g. how
4296  // many nodes are in each set, their ids, etc.
4297  int offset=0;
4298  for (int ns=0; ns<num_node_sets; ++ns)
4299  {
4300  offset += (ns > 0 ? num_nodes_per_set[ns-1] : 0);
4301  // Note: we don't actually call exII::ex_get_var() here because
4302  // we don't need the values. We only need the indices into that vector
4303  // for each (node_id, boundary_id) tuple.
4304  for (int i=0; i<num_nodes_per_set[ns]; ++i)
4305  {
4306  // The read_all_nodesets() function now reads all the node ids into the
4307  // node_sets_node_list vector, which is of length "total_nodes_in_all_sets"
4308  // The old read_nodset() function is no longer called as far as I can tell,
4309  // and should probably be removed? The "offset" that we are using only
4310  // depends on the current nodeset index and the num_nodes_per_set vector,
4311  // which gets filled in by the call to read_all_nodesets().
4312  dof_id_type exodus_node_id = node_sets_node_list[i + offset];
4313 
4314  // FIXME: We should use exodus_node_num_to_libmesh for this,
4315  // but it apparently is never set up, so just
4316  // subtract 1 from the Exodus node id.
4317  dof_id_type converted_node_id = exodus_node_id - 1;
4318 
4319  // Make a NodeBCTuple key from the converted information.
4320  BoundaryInfo::NodeBCTuple key = std::make_tuple
4321  (converted_node_id, nodeset_ids[ns]);
4322 
4323  // Store the array index of this (node, b_id) tuple
4324  bc_array_indices.emplace(key, cast_int<unsigned int>(i));
4325  } // end for (i)
4326  } // end for (ns)
4327 }
4328 
4330 (const MeshBase & mesh,
4331  const std::vector<Real> & values,
4332  int timestep,
4333  const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains)
4334 {
4335  LOG_SCOPE("write_element_values()", "ExodusII_IO_Helper");
4336 
4337  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4338  return;
4339 
4340  // Ask the file how many element vars it has, store it in the num_elem_vars variable.
4341  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_BLOCK, &num_elem_vars);
4342  EX_CHECK_ERR(ex_err, "Error reading number of elemental variables.");
4343 
4344  // We will eventually loop over the element blocks (subdomains) and
4345  // write the data one block at a time. Build a data structure that
4346  // maps each subdomain to a list of element ids it contains.
4347  std::map<subdomain_id_type, std::vector<unsigned int>> subdomain_map;
4348  for (const auto & elem : mesh.active_element_ptr_range())
4349  subdomain_map[elem->subdomain_id()].push_back(elem->id());
4350 
4351  // Use mesh.n_elem() to access into the values vector rather than
4352  // the number of elements the Exodus writer thinks the mesh has,
4353  // which may not include inactive elements.
4355 
4356  // Sanity check: we must have an entry in vars_active_subdomains for
4357  // each variable that we are potentially writing out.
4358  libmesh_assert_equal_to
4359  (vars_active_subdomains.size(),
4360  static_cast<unsigned>(num_elem_vars));
4361 
4362  // For each variable, create a 'data' array which holds all the elemental variable
4363  // values *for a given block* on this processor, then write that data vector to file
4364  // before moving onto the next block.
4365  for (unsigned int var_id=0; var_id<static_cast<unsigned>(num_elem_vars); ++var_id)
4366  {
4367  // The size of the subdomain map is the number of blocks.
4368  auto it = subdomain_map.begin();
4369 
4370  // Reference to the set of active subdomains for the current variable.
4371  const auto & active_subdomains
4372  = vars_active_subdomains[var_id];
4373 
4374  for (unsigned int j=0; it!=subdomain_map.end(); ++it, ++j)
4375  {
4376  // Skip any variable/subdomain pairs that are inactive.
4377  // Note that if active_subdomains is empty, it is interpreted
4378  // as being active on *all* subdomains.
4379  if (!(active_subdomains.empty() || active_subdomains.count(it->first)))
4380  continue;
4381 
4382  // Get reference to list of elem ids which are in the
4383  // current subdomain and count, allocate storage to hold
4384  // data that will be written to file.
4385  const auto & elem_nums = it->second;
4386  const unsigned int num_elems_this_block =
4387  cast_int<unsigned int>(elem_nums.size());
4388  std::vector<Real> data(num_elems_this_block);
4389 
4390  // variable-major ordering is:
4391  // (u1, u2, u3, ..., uN), (v1, v2, v3, ..., vN), ...
4392  // where N is the number of elements.
4393  for (unsigned int k=0; k<num_elems_this_block; ++k)
4394  data[k] = values[var_id*n_elem + elem_nums[k]];
4395 
4396  ex_err = exII::ex_put_var
4397  (ex_id,
4398  timestep,
4399  exII::EX_ELEM_BLOCK,
4400  var_id+1,
4401  this->get_block_id(j),
4402  num_elems_this_block,
4403  MappedOutputVector(data, _single_precision).data());
4404 
4405  EX_CHECK_ERR(ex_err, "Error writing element values.");
4406  }
4407  }
4408 
4409  this->update();
4410 }
4411 
4412 
4413 
4415 (const MeshBase & mesh,
4416  const std::vector<Real> & values,
4417  int timestep,
4418  const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains,
4419  const std::vector<std::string> & derived_var_names,
4420  const std::map<subdomain_id_type, std::vector<std::string>> & subdomain_to_var_names)
4421 {
4422  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4423  return;
4424 
4425  // Ask the file how many element vars it has, store it in the num_elem_vars variable.
4426  ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_BLOCK, &num_elem_vars);
4427  EX_CHECK_ERR(ex_err, "Error reading number of elemental variables.");
4428 
4429  // We will eventually loop over the element blocks (subdomains) and
4430  // write the data one block (subdomain) at a time. Build a data
4431  // structure that keeps track of how many elements are in each
4432  // subdomain. This will allow us to reserve space in the data vector
4433  // we are going to write.
4434  std::map<subdomain_id_type, unsigned int> subdomain_to_n_elem;
4435  for (const auto & elem : mesh.active_element_ptr_range())
4436  subdomain_to_n_elem[elem->subdomain_id()] += 1;
4437 
4438  // Sanity check: we must have an entry in vars_active_subdomains for
4439  // each variable that we are potentially writing out.
4440  libmesh_assert_equal_to
4441  (vars_active_subdomains.size(),
4442  static_cast<unsigned>(num_elem_vars));
4443 
4444  // The size of the subdomain map is the number of blocks.
4445  auto subdomain_to_n_elem_iter = subdomain_to_n_elem.begin();
4446 
4447  // Store range of active Elem pointers. We are going to loop over
4448  // the elements n_vars * n_subdomains times, so let's make sure
4449  // the predicated iterators aren't slowing us down too much.
4450  ConstElemRange elem_range
4451  (mesh.active_elements_begin(),
4452  mesh.active_elements_end());
4453 
4454  for (unsigned int sbd_idx=0;
4455  subdomain_to_n_elem_iter != subdomain_to_n_elem.end();
4456  ++subdomain_to_n_elem_iter, ++sbd_idx)
4457  for (unsigned int var_id=0; var_id<static_cast<unsigned>(num_elem_vars); ++var_id)
4458  {
4459  // Reference to the set of active subdomains for the current variable.
4460  const auto & active_subdomains
4461  = vars_active_subdomains[var_id];
4462 
4463  // If the vars_active_subdomains container passed to this function
4464  // has an empty entry, it means the variable really is not active on
4465  // _any_ subdomains, not that it is active on _all_ subdomains. This
4466  // is just due to the way that we build the vars_active_subdomains
4467  // container.
4468  if (!active_subdomains.count(subdomain_to_n_elem_iter->first))
4469  continue;
4470 
4471  // Vector to hold values that will be written to Exodus file.
4472  std::vector<Real> data;
4473  data.reserve(subdomain_to_n_elem_iter->second);
4474 
4475  unsigned int values_offset = 0;
4476  for (auto & elem : elem_range)
4477  {
4478  // We'll use the Elem's subdomain id in several places below.
4479  subdomain_id_type sbd_id = elem->subdomain_id();
4480 
4481  // Get reference to the list of variable names defining
4482  // the indexing for the current Elem's subdomain.
4483  auto subdomain_to_var_names_iter =
4484  subdomain_to_var_names.find(sbd_id);
4485 
4486  // It's possible, but unusual, for there to be an Elem
4487  // from a subdomain that has no active variables from the
4488  // set of variables we are currently writing. If that
4489  // happens, we can just go to the next Elem because we
4490  // don't need to advance the offset into the values
4491  // vector, etc.
4492  if (subdomain_to_var_names_iter == subdomain_to_var_names.end())
4493  continue;
4494 
4495  const auto & var_names_this_sbd
4496  = subdomain_to_var_names_iter->second;
4497 
4498  // Only extract values if Elem is in the current subdomain.
4499  if (sbd_id == subdomain_to_n_elem_iter->first)
4500  {
4501  // Location of current var_id in the list of all variables on this
4502  // subdomain. FIXME: linear search but it's over a typically relatively
4503  // short vector of active variable names on this subdomain. We could do
4504  // a nested std::map<string,index> instead of a std::vector where the
4505  // location of the string is implicitly the index..
4506  auto pos =
4507  std::find(var_names_this_sbd.begin(),
4508  var_names_this_sbd.end(),
4509  derived_var_names[var_id]);
4510 
4511  libmesh_error_msg_if(pos == var_names_this_sbd.end(),
4512  "Derived name " << derived_var_names[var_id] << " not found!");
4513 
4514  // Find the current variable's location in the list of all variable
4515  // names on the current Elem's subdomain.
4516  auto true_index =
4517  std::distance(var_names_this_sbd.begin(), pos);
4518 
4519  data.push_back(values[values_offset + true_index]);
4520  }
4521 
4522  // The "true" offset is how much we have to advance the index for each Elem
4523  // in this subdomain.
4524  auto true_offset = var_names_this_sbd.size();
4525 
4526  // Increment to the next Elem's values
4527  values_offset += true_offset;
4528  } // for elem
4529 
4530  // Now write 'data' to Exodus file, in single precision if requested.
4531  if (!data.empty())
4532  {
4533  ex_err = exII::ex_put_var
4534  (ex_id,
4535  timestep,
4536  exII::EX_ELEM_BLOCK,
4537  var_id+1,
4538  this->get_block_id(sbd_idx),
4539  data.size(),
4541 
4542  EX_CHECK_ERR(ex_err, "Error writing element values.");
4543  }
4544  } // for each var_id
4545 
4546  this->update();
4547 }
4548 
4549 
4550 
4551 void
4553  const std::vector<Real> & values,
4554  int timestep)
4555 {
4556  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4557  return;
4558 
4559  if (!values.empty())
4560  {
4561  libmesh_assert_equal_to(values.size(), std::size_t(num_nodes));
4562 
4563  ex_err = exII::ex_put_var
4564  (ex_id,
4565  timestep,
4566  exII::EX_NODAL,
4567  var_id,
4568  1, // exII::ex_entity_id, not sure exactly what this is but in the ex_put_nodal_var.c shim, they pass 1
4569  num_nodes,
4570  MappedOutputVector(values, _single_precision).data());
4571 
4572  EX_CHECK_ERR(ex_err, "Error writing nodal values.");
4573 
4574  this->update();
4575  }
4576 }
4577 
4578 
4579 
4580 void ExodusII_IO_Helper::write_information_records(const std::vector<std::string> & records)
4581 {
4582  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4583  return;
4584 
4585  // There may already be information records in the file (for
4586  // example, if we're appending) and in that case, according to the
4587  // Exodus documentation, writing more information records is not
4588  // supported.
4589  int num_info = inquire(*this, exII::EX_INQ_INFO, "Error retrieving the number of information records from file!");
4590  if (num_info > 0)
4591  {
4592  libMesh::err << "Warning! The Exodus file already contains information records.\n"
4593  << "Exodus does not support writing additional records in this situation."
4594  << std::endl;
4595  return;
4596  }
4597 
4598  int num_records = cast_int<int>(records.size());
4599 
4600  if (num_records > 0)
4601  {
4602  NamesData info(num_records, MAX_LINE_LENGTH);
4603 
4604  // If an entry is longer than MAX_LINE_LENGTH characters it's not an error, we just
4605  // write the first MAX_LINE_LENGTH characters to the file.
4606  for (const auto & record : records)
4607  info.push_back_entry(record);
4608 
4609  ex_err = exII::ex_put_info(ex_id, num_records, info.get_char_star_star());
4610  EX_CHECK_ERR(ex_err, "Error writing global values.");
4611 
4612  this->update();
4613  }
4614 }
4615 
4616 
4617 
4618 void ExodusII_IO_Helper::write_global_values(const std::vector<Real> & values, int timestep)
4619 {
4620  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4621  return;
4622 
4623  if (!values.empty())
4624  {
4625  ex_err = exII::ex_put_var
4626  (ex_id,
4627  timestep,
4628  exII::EX_GLOBAL,
4629  1, // var index
4630  0, // obj_id (not used)
4632  MappedOutputVector(values, _single_precision).data());
4633 
4634  EX_CHECK_ERR(ex_err, "Error writing global values.");
4635 
4636  this->update();
4637  }
4638 }
4639 
4640 
4641 
4643 {
4644  ex_err = exII::ex_update(ex_id);
4645  EX_CHECK_ERR(ex_err, "Error flushing buffers to file.");
4646 }
4647 
4648 
4649 
4650 void ExodusII_IO_Helper::read_global_values(std::vector<Real> & values, int timestep)
4651 {
4652  if ((_run_only_on_proc0) && (this->processor_id() != 0))
4653  return;
4654 
4655  values.clear();
4656  values.resize(num_global_vars);
4657  ex_err = exII::ex_get_var
4658  (ex_id,
4659  timestep,
4660  exII::EX_GLOBAL,
4661  1, // var_index
4662  1, // obj_id
4664  MappedInputVector(values, _single_precision).data());
4665 
4666  EX_CHECK_ERR(ex_err, "Error reading global values.");
4667 }
4668 
4669 
4670 
4672 {
4674 }
4675 
4676 
4678 {
4679  _write_hdf5 = write_hdf5;
4680 }
4681 
4682 
4683 
4684 
4686 {
4688 }
4689 
4690 
4691 
4693 {
4694  _coordinate_offset = p;
4695 }
4696 
4697 
4698 std::vector<std::string>
4699 ExodusII_IO_Helper::get_complex_names(const std::vector<std::string> & names,
4700  bool write_complex_abs) const
4701 {
4702  std::vector<std::string> complex_names;
4703 
4704  // This will loop over all names and create new "complex" names
4705  // (i.e. names that start with r_, i_ or a_)
4706  for (const auto & name : names)
4707  {
4708  complex_names.push_back("r_" + name);
4709  complex_names.push_back("i_" + name);
4710  if (write_complex_abs)
4711  complex_names.push_back("a_" + name);
4712  }
4713 
4714  return complex_names;
4715 }
4716 
4717 
4718 
4719 std::vector<std::set<subdomain_id_type>>
4722 (const std::vector<std::set<subdomain_id_type>> & vars_active_subdomains,
4723  bool write_complex_abs) const
4724 {
4725  std::vector<std::set<subdomain_id_type>> complex_vars_active_subdomains;
4726 
4727  for (auto & s : vars_active_subdomains)
4728  {
4729  // Push back the same data enough times for the real, imag, (and
4730  // possibly modulus) for the complex-valued solution.
4731  complex_vars_active_subdomains.push_back(s);
4732  complex_vars_active_subdomains.push_back(s);
4733  if (write_complex_abs)
4734  complex_vars_active_subdomains.push_back(s);
4735  }
4736 
4737  return complex_vars_active_subdomains;
4738 }
4739 
4740 
4741 
4742 std::map<subdomain_id_type, std::vector<std::string>>
4745 (const std::map<subdomain_id_type, std::vector<std::string>> & subdomain_to_var_names,
4746  bool write_complex_abs) const
4747 {
4748  // Eventual return value
4749  std::map<subdomain_id_type, std::vector<std::string>> ret;
4750 
4751  unsigned int num_complex_outputs = write_complex_abs ? 3 : 2;
4752 
4753  for (const auto & pr : subdomain_to_var_names)
4754  {
4755  // Initialize entry for current subdomain
4756  auto & vec = ret[pr.first];
4757 
4758  // Get list of non-complex variable names active on this subdomain.
4759  const auto & varnames = pr.second;
4760 
4761  // Allocate space for the complex-valued entries
4762  vec.reserve(num_complex_outputs * varnames.size());
4763 
4764  // For each varname in the input map, write three variable names
4765  // to the output formed by prepending "r_", "i_", and "a_",
4766  // respectively.
4767  for (const auto & varname : varnames)
4768  {
4769  vec.push_back("r_" + varname);
4770  vec.push_back("i_" + varname);
4771  if (write_complex_abs)
4772  vec.push_back("a_" + varname);
4773  }
4774  }
4775  return ret;
4776 }
4777 
4778 
4779 
4781 {
4782  if (!node_map)
4783  return i;
4784 
4785  libmesh_assert_less (i, node_map->size());
4786  return (*node_map)[i];
4787 }
4788 
4789 
4790 
4792 {
4793  if (!inverse_node_map)
4794  return i;
4795 
4796  libmesh_assert_less (i, inverse_node_map->size());
4797  return (*inverse_node_map)[i];
4798 }
4799 
4800 
4801 
4803 {
4804  if (!side_map)
4805  return i;
4806 
4807  // If we asked for a side that doesn't exist, return an invalid_id
4808  // and allow higher-level code to handle it.
4809  if (static_cast<size_t>(i) >= side_map->size())
4810  return invalid_id;
4811 
4812  return (*side_map)[i];
4813 }
4814 
4815 
4816 
4818 {
4819  // For identity side mappings, we our convention is to return a 1-based index.
4820  if (!inverse_side_map)
4821  return i + 1;
4822 
4823  libmesh_assert_less (i, inverse_side_map->size());
4824  return (*inverse_side_map)[i];
4825 }
4826 
4827 
4828 
4834 {
4835  if (!shellface_map)
4836  return i;
4837 
4838  libmesh_assert_less (i, shellface_map->size());
4839  return (*shellface_map)[i];
4840 }
4841 
4842 
4843 
4845 {
4846  if (!inverse_shellface_map)
4847  return i + 1;
4848 
4849  libmesh_assert_less (i, inverse_shellface_map->size());
4850  return (*inverse_shellface_map)[i];
4851 }
4852 
4853 
4854 
4856 {
4857  return libmesh_type;
4858 }
4859 
4860 
4861 
4863 {
4864  return exodus_type;
4865 }
4866 
4867 
4868 
4873 {
4874  return shellface_index_offset;
4875 }
4876 
4877 ExodusII_IO_Helper::NamesData::NamesData(size_t n_strings, size_t string_length) :
4878  data_table(n_strings),
4879  data_table_pointers(n_strings),
4880  counter(0),
4881  table_size(n_strings)
4882 {
4883  for (size_t i=0; i<n_strings; ++i)
4884  {
4885  data_table[i].resize(string_length + 1);
4886 
4887  // Properly terminate these C-style strings, just to be safe.
4888  data_table[i][0] = '\0';
4889 
4890  // Set pointer into the data_table
4891  data_table_pointers[i] = data_table[i].data();
4892  }
4893 }
4894 
4895 
4896 
4898 {
4899  libmesh_assert_less (counter, table_size);
4900 
4901  // 1.) Copy the C++ string into the vector<char>...
4902  size_t num_copied = name.copy(data_table[counter].data(), data_table[counter].size()-1);
4903 
4904  // 2.) ...And null-terminate it.
4905  data_table[counter][num_copied] = '\0';
4906 
4907  // Go to next row
4908  ++counter;
4909 }
4910 
4911 
4912 
4914 {
4915  return data_table_pointers.data();
4916 }
4917 
4918 
4919 
4921 {
4922  libmesh_error_msg_if(static_cast<unsigned>(i) >= table_size,
4923  "Requested char * " << i << " but only have " << table_size << "!");
4924 
4925  return data_table[i].data();
4926 }
4927 
4928 
4929 
4930 } // namespace libMesh
4931 
4932 
4933 
4934 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
virtual void write_sidesets(const MeshBase &mesh)
Writes the sidesets contained in "mesh".
std::vector< int > num_elems_per_set
void write_sideset_data(const MeshBase &mesh, 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)
Write sideset data for the requested timestep.
OStreamProxy err
void use_mesh_dimension_instead_of_spatial_dimension(bool val)
Sets the underlying value of the boolean flag _use_mesh_dimension_instead_of_spatial_dimension.
void write_var_names(ExodusVarType type, const std::vector< std::string > &names)
Wraps calls to exII::ex_put_var_names() and exII::ex_put_var_param().
The FPEDisabler class puts Floating-Point Exception (FPE) trapping on hold during its lifetime...
Definition: fpe_disabler.h:49
std::tuple< dof_id_type, unsigned short int, boundary_id_type > BCTuple
As above, but the library creates and fills in a vector of (elem-id, side-id, bc-id) triplets and ret...
This class facilitates inline conversion of an input data vector to a different precision level...
ElemType
Defines an enum for geometric element types.
std::vector< std::string > sideset_var_names
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
void read_node_num_map()
Reads the optional node_num_map from the ExodusII mesh file.
std::vector< int > num_sides_per_set
const char * get_elem_type() const
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 elemset per variable...
void read_elemset(int id, int offset)
Reads information about elemset id and inserts it into the global elemset array at the position offse...
std::vector< std::string > elem_var_names
void active_family_tree_by_side(std::vector< const Elem *> &family, unsigned int side, bool reset=true) const
Same as the active_family_tree() member, but only adds elements which are next to side...
Definition: elem.C:2167
virtual dof_id_type n_active_elem() const =0
A Node is like a Point, but with more information.
Definition: node.h:52
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 elemset variables, if any, into the provided data structures.
ExodusHeaderInfo read_header() const
Reads an ExodusII mesh file header, leaving this object&#39;s internal data structures unchanged...
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
std::map< int, std::string > id_to_ss_names
std::vector< std::string > get_complex_names(const std::vector< std::string > &names, bool write_complex_abs) const
This class is used as both an external data structure for passing around Exodus file header informati...
MPI_Info info
std::string get_block_name(int index)
Get the block name for the given block index if supplied in the mesh file.
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
std::map< int, std::string > id_to_edge_block_names
void read_sideset_info()
Reads information about all of the sidesets in the ExodusII mesh file.
void read_nodal_var_values(std::string nodal_var_name, int time_step)
Reads the nodal values for the variable &#39;nodal_var_name&#39; at the specified time into the &#39;nodal_var_va...
std::size_t n_edge_conds() const
ExodusII_IO_Helper(const ParallelObject &parent, bool v=false, bool run_only_on_proc0=true, bool single_precision=false)
Constructor.
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
void write_information_records(const std::vector< std::string > &records)
Writes the vector of information records.
virtual void read_var_names_impl(const char *var_type, int &count, std::vector< std::string > &result)
read_var_names() dispatches to this function.
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_as_dimension(unsigned dim)
Sets the value of _write_as_dimension.
std::map< dof_id_type, Real > nodal_var_values
const boundary_id_type side_id
std::string get_side_set_name(int index)
Get the side set name for the given side set index if supplied in the mesh file.
void sum(T &r) const
const std::map< boundary_id_type, std::string > & get_sideset_name_map() const
bool has_elem_integer(std::string_view name) const
Definition: mesh_base.C:638
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void get_sideset_data_indices(const MeshBase &mesh, 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...
const Parallel::Communicator & comm() const
virtual void write_nodesets(const MeshBase &mesh)
Writes the nodesets contained in "mesh".
void build_side_boundary_ids(std::vector< boundary_id_type > &b_ids) const
Builds the list of unique side boundary ids.
void read_bex_cv_blocks()
Reads the optional bex_cv_blocks from the ExodusII mesh file.
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 nodeset variables, if any, into the provided data structures.
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
std::vector< Real > time_steps
MappedOutputVector(const std::vector< Real > &vec_in, bool single_precision_in)
virtual void write_elements(const MeshBase &mesh, bool use_discontinuous=false)
Writes the elements contained in "mesh".
std::vector< std::set< subdomain_id_type > > get_complex_vars_active_subdomains(const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains, bool write_complex_abs) const
returns a "tripled" copy of vars_active_subdomains, which is necessary in the complex-valued case...
void write_elemsets(const MeshBase &mesh)
Write elemsets stored on the Mesh to the exo file.
const ExodusII_IO_Helper::Conversion & get_conversion(const ElemType type) const
void print_nodes(std::ostream &out_stream=libMesh::out)
Prints the nodal information, by default to libMesh::out.
The libMesh namespace provides an interface to certain functionality in the library.
std::vector< std::vector< char > > data_table
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void read_time_steps()
Reads and stores the timesteps in the &#39;time_steps&#39; array.
std::vector< int > node_sets_node_index
std::vector< char > & title
void write_nodal_values(int var_id, const std::vector< Real > &values, int timestep)
Writes the vector of values to a nodal variable.
bool _add_sides
Set to true iff we want to write separate "side" elements too.
void message(std::string_view msg)
Prints the message defined in msg.
void check_existing_vars(ExodusVarType type, std::vector< std::string > &names, std::vector< std::string > &names_from_file)
When appending: during initialization, check that variable names in the file match those you attempt ...
Real distance(const Point &p)
void read_nodes()
Reads the nodal data (x,y,z coordinates) from the ExodusII mesh file.
static const unsigned int type_to_n_nodes_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of nodes in the element...
Definition: elem.h:650
std::vector< std::string > nodeset_var_names
std::map< dof_id_type, dof_id_type > libmesh_node_num_to_exodus
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
void close() noexcept
Closes the ExodusII mesh file.
void write_element_values(const MeshBase &mesh, const std::vector< Real > &values, int timestep, const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains)
Writes the vector of values to the element variables.
void write_timestep(int timestep, Real time)
Writes the time for the timestep.
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
void get_elemsets(dof_id_type elemset_code, MeshBase::elemset_type &id_set_to_fill) const
Look up the element sets for a given elemset code and vice-versa.
Definition: mesh_base.C:429
void read_edge_blocks(MeshBase &mesh)
Read in edge blocks, storing information in the BoundaryInfo object.
std::tuple< dof_id_type, boundary_id_type > NodeBCTuple
As above, but the library creates and fills in a vector of (node-id, bc-id) pairs and returns it to t...
virtual void initialize_element_variables(std::vector< std::string > names, const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains)
Sets up the nodal variables.
void open(const char *filename, bool read_only)
Opens an ExodusII mesh file named filename.
unsigned int get_elem_integer_index(std::string_view name) const
Definition: mesh_base.C:626
processor_id_type n_processors() const
void libmesh_ignore(const Args &...)
void build_node_list(std::vector< dof_id_type > &node_id_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of nodes and ids for those nodes.
char * get_char_star(int i)
Provide access to the i&#39;th underlying char *.
const std::map< boundary_id_type, std::string > & get_nodeset_name_map() const
const std::string & get_edgeset_name(boundary_id_type id) const
int8_t boundary_id_type
Definition: id_types.h:51
void update()
Uses ex_update() to flush buffers to file.
void push_back_entry(const std::string &name)
Adds another name to the current data table.
dof_id_type id() const
Definition: dof_object.h:828
void read_all_nodesets()
New API that reads all nodesets simultaneously.
void read_global_values(std::vector< Real > &values, int timestep)
Reads the vector of global variables.
std::vector< std::string > global_var_names
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
MappedInputVector(std::vector< Real > &vec_in, bool single_precision_in)
void read_var_names(ExodusVarType type)
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)
Write nodeset data for the requested timestep.
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const =0
Similar to Elem::local_side_node(), but instead of a side id, takes an edge id and a node id on that ...
virtual dof_id_type max_elem_id() const =0
void subdomain_ids(std::set< subdomain_id_type > &ids, const bool global=true) const
Constructs a list of all subdomain identifiers in the local mesh if global == false, and in the global mesh if global == true (default).
Definition: mesh_base.C:996
void read_elemset_info()
Reads information about all of the elemsets in the ExodusII mesh file.
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
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)
Write elemset data for the requested timestep.
void write_element_values_element_major(const MeshBase &mesh, const std::vector< Real > &values, int timestep, const std::vector< std::set< subdomain_id_type >> &vars_active_subdomains, const std::vector< std::string > &derived_var_names, const std::map< subdomain_id_type, std::vector< std::string >> &subdomain_to_var_names)
Same as the function above, but assume the input &#39;values&#39; vector is in element-major order...
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
This class facilitates reading in vectors from Exodus file that may be of a different floating point ...
void read_nodeset_info()
Reads information about all of the nodesets in the ExodusII mesh file.
This is the ExodusII_IO_Helper class.
void set_coordinate_offset(Point p)
Allows you to set a vector that is added to the coordinates of all of the nodes.
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:1729
void write_global_values(const std::vector< Real > &values, int timestep)
Writes the vector of global variables.
std::vector< dof_id_type > _true_node_offsets
If we&#39;re adding "fake" sides to visualize SIDE_DISCONTINUOUS variables, we also need to know how many...
void build_shellface_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &shellface_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, shellfaces, and boundary ids for those shellfaces.
std::vector< char > title
An object whose state is distributed along a set of processors.
void read_elem_in_block(int block)
Reads all of the element connectivity for block block in the ExodusII mesh file.
const std::string & get_nodeset_name(boundary_id_type id) const
std::map< dof_id_type, dof_id_type > libmesh_elem_num_to_exodus
void initialize_global_variables(std::vector< std::string > names)
Sets up the global variables.
void read_sideset_data(const MeshBase &mesh, 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)
Read sideset variables, if any, into the provided data structures.
std::string enum_to_string(const T e)
std::vector< int > num_node_df_per_set
char ** get_char_star_star()
Provide access to the underlying C data table.
int get_node_set_id(int index)
Get the node set id for the given node set index.
virtual const Elem * elem_ptr(const dof_id_type i) const =0
void build_node_boundary_ids(std::vector< boundary_id_type > &b_ids) const
Builds the list of unique node boundary ids.
ExodusVarType
Wraps calls to exII::ex_get_var_names() and exII::ex_get_var_param().
std::map< std::string, ElemType > element_equivalence_map
Defines equivalence classes of Exodus element types that map to libmesh ElemTypes.
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
std::set< elemset_id_type > elemset_type
Typedef for the "set" container used to store elemset ids.
Definition: mesh_base.h:318
virtual void initialize(std::string title, const MeshBase &mesh, bool use_discontinuous=false)
Initializes the Exodus file.
virtual void create(std::string filename)
Opens an ExodusII mesh file named filename for writing.
std::vector< int > node_sets_node_list
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void read_num_time_steps()
Reads the number of timesteps currently stored in the Exodus file and stores it in the num_time_steps...
void read_block_info()
Reads information for all of the blocks in the ExodusII mesh file.
void read_elem_num_map()
Reads the optional node_num_map from the ExodusII mesh file.
void read_and_store_header_info()
Reads an ExodusII mesh file header, and stores required information on this object.
unsigned int spatial_dimension() const
Definition: mesh_base.C:543
void read_sideset(int id, int offset)
Reads information about sideset id and inserts it into the global sideset array at the position offse...
std::vector< dof_id_type > _added_side_node_offsets
If we&#39;re adding "fake" sides to visualize SIDE_DISCONTINUOUS variables, _added_side_node_offsets[p] g...
void build_edge_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &edge_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, edges, and boundary ids for those edges.
OStreamProxy out
NamesData(size_t n_strings, size_t string_length)
Constructor.
const std::set< boundary_id_type > & get_edge_boundary_ids() const
const std::string & get_sideset_name(boundary_id_type id) const
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
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)...
std::map< int, std::string > id_to_block_names
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
std::string & edgeset_name(boundary_id_type id)
int get_side_set_id(int index)
Get the side set id for the given side set index.
std::vector< std::string > elemset_var_names
void print_header()
Prints the ExodusII mesh file header, which includes the mesh title, the number of nodes...
std::map< int, std::string > id_to_ns_names
void initialize_nodal_variables(std::vector< std::string > names)
Sets up the nodal variables.
std::vector< int > num_nodes_per_set
std::vector< std::string > nodal_var_names
std::string get_node_set_name(int index)
Get the node set name for the given node set index if supplied in the mesh file.
std::map< int, std::string > id_to_elemset_names
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i)=0
virtual const Node * node_ptr(const dof_id_type i) const =0
std::vector< int > edge_block_ids
processor_id_type processor_id() const
void read_qa_records()
Reads the QA records from an ExodusII file.
This class is useful for managing anything that requires a char ** input/output in ExodusII file...
void read_elemental_var_values(std::string elemental_var_name, int time_step, std::map< dof_id_type, Real > &elem_var_value_map)
Reads elemental values for the variable &#39;elemental_var_name&#39; at the specified timestep into the &#39;elem...
static ElemType first_order_equivalent_type(const ElemType et)
Definition: elem.C:3066
std::vector< Real > node_sets_dist_fact
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
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...
int get_block_id(int index)
Get the block number for the given block index.
std::map< subdomain_id_type, std::vector< std::string > > get_complex_subdomain_to_var_names(const std::map< subdomain_id_type, std::vector< std::string >> &subdomain_to_var_names, bool write_complex_abs) const
Takes a map from subdomain id -> vector of active variable names as input and returns a corresponding...
std::vector< int > num_df_per_set
std::map< int, std::map< ElemType, ExodusII_IO_Helper::Conversion > > conversion_map
Associates libMesh ElemTypes with node/face/edge/etc.
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
std::vector< int > node_sets_dist_index
std::vector< int > elemset_id_list
void build_shellface_boundary_ids(std::vector< boundary_id_type > &b_ids) const
Builds the list of unique shellface boundary ids.
virtual void write_nodal_coordinates(const MeshBase &mesh, bool use_discontinuous=false)
Writes the nodal coordinates contained in "mesh".
std::vector< int > num_elem_df_per_set
void write_var_names_impl(const char *var_type, int &count, const std::vector< std::string > &names)
write_var_names() dispatches to this function.
std::vector< std::vector< long unsigned int > > bex_cv_conn
uint8_t dof_id_type
Definition: id_types.h:67
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...
static bool redundant_added_side(const Elem &elem, unsigned int side)
std::vector< std::vector< std::vector< Real > > > bex_dense_constraint_vecs