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