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