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