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