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