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