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