From 7dc50cf7066fe1720b634a2797fec5282a8464e3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 12 Feb 2026 15:27:08 -0700 Subject: [PATCH 1/6] Factor out ExodusII_IO_Helper::build_subdomain_map With the "local" option added it works for Nemesis too. Hopefully this will let us get add_sides support into Nemesis eventually, but also let me get multi-block-subdomain support into both at once. --- include/mesh/exodusII_io_helper.h | 19 +++- include/mesh/nemesis_io_helper.h | 5 -- src/mesh/exodusII_io_helper.C | 145 ++++++++++++++---------------- src/mesh/nemesis_io_helper.C | 27 ++---- 4 files changed, 91 insertions(+), 105 deletions(-) diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index fbf57f5dc3..b0aa519011 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -1104,13 +1104,30 @@ class ExodusII_IO_Helper : public ParallelObject int & count, std::vector & result); -private: +protected: /** * Set to true iff we want to write separate "side" elements too. */ bool _add_sides = false; + /** + * Map of subdomains to element numbers. + */ + std::map> _subdomain_map; + + /** + * One beyond the last "real" subdomain id, in cases where we have + * "fake" subdomain ids for added sides + */ + subdomain_id_type _subdomain_id_end; + + /** + * Method for constructing _subdomain_map + */ + void build_subdomain_map(const MeshBase & mesh, bool local); + +private: /** * write_var_names() dispatches to this function. */ diff --git a/include/mesh/nemesis_io_helper.h b/include/mesh/nemesis_io_helper.h index f88a78627f..962f4260df 100644 --- a/include/mesh/nemesis_io_helper.h +++ b/include/mesh/nemesis_io_helper.h @@ -435,11 +435,6 @@ class Nemesis_IO_Helper : public ExodusII_IO_Helper */ std::set nodes_attached_to_local_elems; - /** - * Map of subdomains to element numbers. - */ - std::map> subdomain_map; - /** * This is the block connectivity, i.e. for each subdomain (block) there * is an element connectivity list. This map associates the block ID to that vector. diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 3dcbeaa6dc..668da15cbe 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -198,57 +198,6 @@ const std::vector prism_inverse_face_map = {4, 1, 2, 3, 5}; #endif } - - std::map> - build_subdomain_map(const MeshBase & mesh, bool add_sides, subdomain_id_type & subdomain_id_end) - { - std::map> subdomain_map; - - // If we've been asked to add side elements, those will go in - // their own blocks. - if (add_sides) - { - std::set sbd_ids; - mesh.subdomain_ids(sbd_ids); - if (!sbd_ids.empty()) - subdomain_id_end = *sbd_ids.rbegin()+1; - } - - // Loop through element and map between block and element vector. - for (const auto & elem : mesh.active_element_ptr_range()) - { - // We skip writing infinite elements to the Exodus file, so - // don't put them in the subdomain_map. That way the number of - // blocks should be correct. - if (elem->infinite()) - continue; - - subdomain_map[ elem->subdomain_id() ].push_back(elem->id()); - - // If we've been asked to add side elements, those will go in their own - // blocks. We don't have any ids to list for elements that don't - // explicitly exist in the mesh, but we do an entry to keep - // track of the number of elements we'll add in each new block. - if (add_sides) - for (auto s : elem->side_index_range()) - { - if (EquationSystems::redundant_added_side(*elem,s)) - continue; - - auto & marker = - subdomain_map[subdomain_id_end + elem->side_type(s)]; - if (marker.empty()) - marker.push_back(1); - else - ++marker[0]; - } - } - - if (!add_sides && !subdomain_map.empty()) - subdomain_id_end = subdomain_map.rbegin()->first + 1; - - return subdomain_map; - } } // end anonymous namespace @@ -2356,8 +2305,7 @@ void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh num_edge = bi.n_edge_conds(); // We need to know about all processors' subdomains - subdomain_id_type subdomain_id_end = 0; - auto subdomain_map = build_subdomain_map(mesh, _add_sides, subdomain_id_end); + build_subdomain_map(mesh, false); num_elem = n_active_elem; num_nodes = 0; @@ -2478,7 +2426,7 @@ void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh num_side_sets = cast_int(unique_side_boundaries.size()); num_node_sets = cast_int(unique_node_boundaries.size()); - num_elem_blk = cast_int(subdomain_map.size()); + num_elem_blk = cast_int(this->_subdomain_map.size()); if (str_title.size() > MAX_LINE_LENGTH) { @@ -2702,16 +2650,11 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti { LOG_SCOPE("write_elements()", "ExodusII_IO_Helper"); - // Map from block ID to a vector of element IDs in that block. Element - // IDs are now of type dof_id_type, subdomain IDs are of type subdomain_id_type. - subdomain_id_type subdomain_id_end = 0; - auto subdomain_map = build_subdomain_map(mesh, _add_sides, subdomain_id_end); - if ((_run_only_on_proc0) && (this->processor_id() != 0)) return; // element map vector - num_elem_blk = cast_int(subdomain_map.size()); + num_elem_blk = cast_int(this->_subdomain_map.size()); block_ids.resize(num_elem_blk); std::vector elem_blk_id; @@ -2732,15 +2675,15 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // counter indexes into the block_ids vector unsigned int counter = 0; - for (auto & [subdomain_id, element_id_vec] : subdomain_map) + for (auto & [subdomain_id, element_id_vec] : this->_subdomain_map) { block_ids[counter] = subdomain_id; - const ElemType elem_t = (subdomain_id >= subdomain_id_end) ? - ElemType(subdomain_id - subdomain_id_end) : + const ElemType elem_t = (subdomain_id >= this->_subdomain_id_end) ? + ElemType(subdomain_id - this->_subdomain_id_end) : mesh.elem_ref(element_id_vec[0]).type(); - if (subdomain_id >= subdomain_id_end) + if (subdomain_id >= this->_subdomain_id_end) { libmesh_assert(_add_sides); libmesh_assert(element_id_vec.size() == 1); @@ -2993,14 +2936,14 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti next_fake_id = mesh.next_unique_id(); #endif - for (auto & [subdomain_id, element_id_vec] : subdomain_map) + for (auto & [subdomain_id, element_id_vec] : this->_subdomain_map) { // Use the first element in the block to get representative // information for a "real" block. Note that Exodus assumes all // elements in a block are of the same type! We are using that // same assumption here! - const ElemType elem_t = (subdomain_id >= subdomain_id_end) ? - ElemType(subdomain_id - subdomain_id_end) : + const ElemType elem_t = (subdomain_id >= this->_subdomain_id_end) ? + ElemType(subdomain_id - this->_subdomain_id_end) : mesh.elem_ref(element_id_vec[0]).type(); const auto & conv = get_conversion(elem_t); @@ -3010,7 +2953,7 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // If this is a *real* block, we just loop over vectors of // element ids to add. - if (subdomain_id < subdomain_id_end) + if (subdomain_id < this->_subdomain_id_end) { connect.resize(element_id_vec.size()*num_nodes_per_elem); @@ -3130,7 +3073,7 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti nullptr, // elem_edge_conn (unused) nullptr); // elem_face_conn (unused) EX_CHECK_ERR(ex_err, "Error writing element connectivities"); - } // end for (auto & [subdomain_id, element_id_vec] : subdomain_map) + } // end for (auto & [subdomain_id, element_id_vec] : this->_subdomain_map) // write out the element number map that we created ex_err = exII::ex_put_elem_num_map(ex_id, elem_num_map.data()); @@ -4509,12 +4452,7 @@ void ExodusII_IO_Helper::write_element_values ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_BLOCK, &num_elem_vars); EX_CHECK_ERR(ex_err, "Error reading number of elemental variables."); - // We will eventually loop over the element blocks (subdomains) and - // write the data one block at a time. Build a data structure that - // maps each subdomain to a list of element ids it contains. - std::map> subdomain_map; - for (const auto & elem : mesh.active_element_ptr_range()) - subdomain_map[elem->subdomain_id()].push_back(elem->id()); + this->build_subdomain_map(mesh, false); // Use mesh.n_elem() to access into the values vector rather than // the number of elements the Exodus writer thinks the mesh has, @@ -4533,13 +4471,13 @@ void ExodusII_IO_Helper::write_element_values for (unsigned int var_id=0; var_id(num_elem_vars); ++var_id) { // The size of the subdomain map is the number of blocks. - auto it = subdomain_map.begin(); + auto it = this->_subdomain_map.begin(); // Reference to the set of active subdomains for the current variable. const auto & active_subdomains = vars_active_subdomains[var_id]; - for (unsigned int j=0; it!=subdomain_map.end(); ++it, ++j) + for (unsigned int j=0; it!=this->_subdomain_map.end(); ++it, ++j) { // Skip any variable/subdomain pairs that are inactive. // Note that if active_subdomains is empty, it is interpreted @@ -4957,6 +4895,59 @@ get_complex_subdomain_to_var_names +void ExodusII_IO_Helper::build_subdomain_map(const MeshBase & mesh, bool local) +{ + // Start from scratch + this->_subdomain_map.clear(); + this->_subdomain_id_end = 0; + + // If we've been asked to add side elements, those will go in + // their own blocks. + if (this->_add_sides) + { + std::set sbd_ids; + mesh.subdomain_ids(sbd_ids); + if (!sbd_ids.empty()) + this->_subdomain_id_end = *sbd_ids.rbegin()+1; + } + + // Loop through element and map between block and element vector. + const auto range = local ? mesh.active_local_element_ptr_range() : + mesh.active_element_ptr_range(); + for (const auto & elem : range) + { + // We skip writing infinite elements to the Exodus file, so + // don't put them in the subdomain_map. That way the number of + // blocks should be correct. + if (elem->infinite()) + continue; + + this->_subdomain_map[ elem->subdomain_id() ].push_back(elem->id()); + + // If we've been asked to add side elements, those will go in their own + // blocks. We don't have any ids to list for elements that don't + // explicitly exist in the mesh, but we do an entry to keep + // track of the number of elements we'll add in each new block. + if (this->_add_sides) + for (auto s : elem->side_index_range()) + { + if (EquationSystems::redundant_added_side(*elem,s)) + continue; + + auto & marker = + this->_subdomain_map[this->_subdomain_id_end + elem->side_type(s)]; + if (marker.empty()) + marker.push_back(1); + else + ++marker[0]; + } + } + + if (!this->_add_sides && !this->_subdomain_map.empty()) + this->_subdomain_id_end = this->_subdomain_map.rbegin()->first + 1; +} + + int ExodusII_IO_Helper::Conversion::get_node_map(int i) const { if (!node_map) diff --git a/src/mesh/nemesis_io_helper.C b/src/mesh/nemesis_io_helper.C index 3f2811f09b..5dbe677333 100644 --- a/src/mesh/nemesis_io_helper.C +++ b/src/mesh/nemesis_io_helper.C @@ -1666,20 +1666,7 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh) // Make sure there is no leftover information in the subdomain_map, and reserve // enough space to store the elements we need. - this->subdomain_map.clear(); - for (const auto & [sbd_id, cnt] : local_subdomain_counts) - { - if (verbose) - { - libMesh::out << "[" << this->processor_id() << "] " - << "local_subdomain_counts [" << static_cast(sbd_id) << "]= " - << cnt - << std::endl; - } - - this->subdomain_map[sbd_id].reserve(cnt); - } - + this->build_subdomain_map(pmesh, true); // First loop over the elements to figure out which elements are in which subdomain for (const auto & elem : pmesh.active_local_element_ptr_range()) @@ -1687,10 +1674,6 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh) // Grab the nodes while we're here. for (auto n : elem->node_index_range()) this->nodes_attached_to_local_elems.insert( elem->node_id(n) ); - - subdomain_id_type cur_subdomain = elem->subdomain_id(); - - this->subdomain_map[cur_subdomain].push_back(elem->id()); } // Set num_nodes which is used by exodusII_io_helper @@ -1736,7 +1719,7 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh) this->libmesh_elem_num_to_exodus.clear(); // Now loop over each subdomain and get a unique numbering for the elements - for (auto & [block_id, elem_ids_this_subdomain] : subdomain_map) + for (auto & [block_id, elem_ids_this_subdomain] : this->_subdomain_map) { block_ids.push_back(block_id); @@ -2308,7 +2291,7 @@ void Nemesis_IO_Helper::write_elements(const MeshBase & mesh, bool /*use_discont subdomain_id_type block = cast_int(it->first); const std::vector & this_block_connectivity = it->second; - std::vector & elements_in_this_block = subdomain_map[block]; + std::vector & elements_in_this_block = this->_subdomain_map[block]; // Use the first element in this block to get representative information. // Note that Exodus assumes all elements in a block are of the same type! @@ -2662,10 +2645,10 @@ Nemesis_IO_Helper::write_element_values(const MeshBase & mesh, { const subdomain_id_type sbd_id = cast_int(sbd_id_int); - auto it = subdomain_map.find(sbd_id); + auto it = this->_subdomain_map.find(sbd_id); const std::vector empty_vec; const std::vector & elem_ids = - (it == subdomain_map.end()) ? empty_vec : it->second; + (it == this->_subdomain_map.end()) ? empty_vec : it->second; // Possibly skip this (variable, subdomain) combination. Also, check that there is at // least one element on the subdomain... Indeed, it is possible to have zero elements, From 6e0c302fbb000ec20037c1d27aedfbeff41d1255 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 13 Feb 2026 08:22:32 -0700 Subject: [PATCH 2/6] Move Exodus helper function implementations to .C These are long enough that they might not be inlined regardless, and this helper class is complicated enough with the declarations alone --- include/mesh/exodusII_io_helper.h | 24 ++---------------------- src/mesh/exodusII_io_helper.C | 28 ++++++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 22 deletions(-) diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index b0aa519011..7cd0d3b2ac 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -942,34 +942,14 @@ class ExodusII_IO_Helper : public ParallelObject * each processor. This plus a node id gives a valid nodal solution * vector index. */ - dof_id_type node_id_to_vec_id(dof_id_type n) const - { - if (_added_side_node_offsets.empty()) - return n; - - // Find the processor id that has node_id in the parallel vec - const auto lb = std::upper_bound(_true_node_offsets.begin(), - _true_node_offsets.end(), n); - libmesh_assert(lb != _true_node_offsets.end()); - const processor_id_type p = lb - _true_node_offsets.begin(); - - return n + (p ? _added_side_node_offsets[p-1] : 0); - } + dof_id_type node_id_to_vec_id(dof_id_type n) const; /* * Returns the sum of both added node "offsets" on processors 0 * through p-1 and real nodes added on processors 0 to p. * This is the starting index for added nodes' data. */ - dof_id_type added_node_offset_on(processor_id_type p) const - { - libmesh_assert (p < _true_node_offsets.size()); - const dof_id_type added_node_offsets = - (_added_side_node_offsets.empty() || !p) ? 0 : - _added_side_node_offsets[p-1]; - return _true_node_offsets[p] + added_node_offsets; - } - + dof_id_type added_node_offset_on(processor_id_type p) const; protected: /** diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 668da15cbe..05c77e9c92 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -4948,6 +4948,34 @@ void ExodusII_IO_Helper::build_subdomain_map(const MeshBase & mesh, bool local) } + +dof_id_type ExodusII_IO_Helper::node_id_to_vec_id(dof_id_type n) const +{ + if (_added_side_node_offsets.empty()) + return n; + + // Find the processor id that has node_id in the parallel vec + const auto lb = std::upper_bound(_true_node_offsets.begin(), + _true_node_offsets.end(), n); + libmesh_assert(lb != _true_node_offsets.end()); + const processor_id_type p = lb - _true_node_offsets.begin(); + + return n + (p ? _added_side_node_offsets[p-1] : 0); +} + + + +dof_id_type ExodusII_IO_Helper::added_node_offset_on(processor_id_type p) const +{ + libmesh_assert (p < _true_node_offsets.size()); + const dof_id_type added_node_offsets = + (_added_side_node_offsets.empty() || !p) ? 0 : + _added_side_node_offsets[p-1]; + return _true_node_offsets[p] + added_node_offsets; +} + + + int ExodusII_IO_Helper::Conversion::get_node_map(int i) const { if (!node_map) From 25f5d87a8baa0b0c25250a5c005a8ff12bcb0582 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 13 Feb 2026 08:39:39 -0700 Subject: [PATCH 3/6] Factor out added_side node offset calculations --- include/mesh/exodusII_io_helper.h | 6 ++ src/mesh/exodusII_io_helper.C | 118 ++++++++++++++++-------------- 2 files changed, 70 insertions(+), 54 deletions(-) diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index 7cd0d3b2ac..3990ee1435 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -952,6 +952,12 @@ class ExodusII_IO_Helper : public ParallelObject dof_id_type added_node_offset_on(processor_id_type p) const; protected: + /** + * Calculate _added_side_node_offsets needed to add "fake" side + * elements to the given mesh + */ + void calculate_added_side_node_offsets(const MeshBase & mesh); + /** * When appending: during initialization, check that variable names * in the file match those you attempt to initialize with. diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 05c77e9c92..754ce900bf 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -2310,60 +2310,9 @@ void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh num_elem = n_active_elem; num_nodes = 0; - // If we're adding face elements they'll need copies of their nodes. - // We also have to count of how many nodes (and gaps between nodes!) - // are on each processor, to calculate offsets for any nodal data - // writing later. - _added_side_node_offsets.clear(); - if (_add_sides) - { - dof_id_type num_side_elem = 0; - dof_id_type num_local_side_nodes = 0; - - for (const auto & elem : mesh.active_local_element_ptr_range()) - { - for (auto s : elem->side_index_range()) - { - if (EquationSystems::redundant_added_side(*elem,s)) - continue; - - num_side_elem++; - num_local_side_nodes += elem->nodes_on_side(s).size(); - } - } - - mesh.comm().sum(num_side_elem); - num_elem += num_side_elem; - - mesh.comm().allgather(num_local_side_nodes, _added_side_node_offsets); - const processor_id_type n_proc = mesh.n_processors(); - libmesh_assert_equal_to(n_proc, _added_side_node_offsets.size()); - - for (auto p : make_range(n_proc-1)) - _added_side_node_offsets[p+1] += _added_side_node_offsets[p]; - - num_nodes = _added_side_node_offsets[n_proc-1]; - - dof_id_type n_local_nodes = cast_int - (std::distance(mesh.local_nodes_begin(), - mesh.local_nodes_end())); - dof_id_type n_total_nodes = n_local_nodes; - mesh.comm().sum(n_total_nodes); - - const dof_id_type max_nn = mesh.max_node_id(); - const dof_id_type n_gaps = max_nn - n_total_nodes; - const dof_id_type gaps_per_processor = n_gaps / n_proc; - const dof_id_type remainder_gaps = n_gaps % n_proc; - - n_local_nodes = n_local_nodes + // Actual nodes - gaps_per_processor + // Our even share of gaps - (mesh.processor_id() < remainder_gaps); // Leftovers - - mesh.comm().allgather(n_local_nodes, _true_node_offsets); - for (auto p : make_range(n_proc-1)) - _true_node_offsets[p+1] += _true_node_offsets[p]; - libmesh_assert_equal_to(_true_node_offsets[n_proc-1], mesh.max_node_id()); - } + // If we're adding face elements they'll need copies of their nodes, + // and we'll need to manage the extra copies. + this->calculate_added_side_node_offsets(mesh); // If _write_as_dimension is nonzero, use it to set num_dim in the Exodus file. if (_write_as_dimension) @@ -4949,6 +4898,67 @@ void ExodusII_IO_Helper::build_subdomain_map(const MeshBase & mesh, bool local) +void ExodusII_IO_Helper::calculate_added_side_node_offsets(const MeshBase & mesh) +{ + + // If we're adding face elements they'll need copies of their nodes. + // We also have to count of how many nodes (and gaps between nodes!) + // are on each processor, to calculate offsets for any nodal data + // writing later. + _added_side_node_offsets.clear(); + if (!_add_sides) + return; + + dof_id_type num_side_elem = 0; + dof_id_type num_local_side_nodes = 0; + + for (const auto & elem : mesh.active_local_element_ptr_range()) + { + for (auto s : elem->side_index_range()) + { + if (EquationSystems::redundant_added_side(*elem,s)) + continue; + + num_side_elem++; + num_local_side_nodes += elem->nodes_on_side(s).size(); + } + } + + mesh.comm().sum(num_side_elem); + num_elem += num_side_elem; + + mesh.comm().allgather(num_local_side_nodes, _added_side_node_offsets); + const processor_id_type n_proc = mesh.n_processors(); + libmesh_assert_equal_to(n_proc, _added_side_node_offsets.size()); + + for (auto p : make_range(n_proc-1)) + _added_side_node_offsets[p+1] += _added_side_node_offsets[p]; + + num_nodes = _added_side_node_offsets[n_proc-1]; + + dof_id_type n_local_nodes = cast_int + (std::distance(mesh.local_nodes_begin(), + mesh.local_nodes_end())); + dof_id_type n_total_nodes = n_local_nodes; + mesh.comm().sum(n_total_nodes); + + const dof_id_type max_nn = mesh.max_node_id(); + const dof_id_type n_gaps = max_nn - n_total_nodes; + const dof_id_type gaps_per_processor = n_gaps / n_proc; + const dof_id_type remainder_gaps = n_gaps % n_proc; + + n_local_nodes = n_local_nodes + // Actual nodes + gaps_per_processor + // Our even share of gaps + (mesh.processor_id() < remainder_gaps); // Leftovers + + mesh.comm().allgather(n_local_nodes, _true_node_offsets); + for (auto p : make_range(n_proc-1)) + _true_node_offsets[p+1] += _true_node_offsets[p]; + libmesh_assert_equal_to(_true_node_offsets[n_proc-1], mesh.max_node_id()); +} + + + dof_id_type ExodusII_IO_Helper::node_id_to_vec_id(dof_id_type n) const { if (_added_side_node_offsets.empty()) From 452b82514bb7ea22369a3a1fd0252ad01ec8d9be Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 13 Feb 2026 10:04:44 -0700 Subject: [PATCH 4/6] Remove redundant specifier --- include/mesh/exodusII_io_helper.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index 3990ee1435..008fe7098b 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -1090,8 +1090,6 @@ class ExodusII_IO_Helper : public ParallelObject int & count, std::vector & result); -protected: - /** * Set to true iff we want to write separate "side" elements too. */ From d67790761162c39a0ccc2530b01baf6b02f6c82e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 13 Feb 2026 10:05:23 -0700 Subject: [PATCH 5/6] Factor out read_extra_integers --- include/mesh/exodusII_io_helper.h | 12 ++++++ src/mesh/exodusII_io.C | 49 +++------------------ src/mesh/exodusII_io_helper.C | 71 +++++++++++++++++++++++++++++++ 3 files changed, 88 insertions(+), 44 deletions(-) diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index 008fe7098b..97942358f9 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -306,6 +306,18 @@ class ExodusII_IO_Helper : public ParallelObject int time_step, std::map & elem_var_value_map); + /** + * Reads element "extra integer" values, using the variable names + * specified in \p extra_integer_var_names, and converting the + * stored values to integer format. The returned value is a vector + * (indexed in the same way as extra_integer_vars) of maps from + * element id to extra integer value. + * + * Currently values from the last time step are used. + */ + std::vector> read_extra_integers( + const std::vector & extra_integer_var_names); + /** * Helper function that takes a (1-based) Exodus node/elem id and * determines the corresponding libMesh Node/Elem id. Takes into account diff --git a/src/mesh/exodusII_io.C b/src/mesh/exodusII_io.C index ef2b60d290..3a986575fe 100644 --- a/src/mesh/exodusII_io.C +++ b/src/mesh/exodusII_io.C @@ -406,13 +406,11 @@ void ExodusII_IO::read (const std::string & fname) // have already been attached to spline control nodes. mesh.reserve_elem(exio_helper->w.size() + exio_helper->num_elem); - // Read variables for extra integer IDs - std::vector> elem_ids(extra_ids.size()); - // We use the last time step to load the IDs exio_helper->read_num_time_steps(); - unsigned int last_step = exio_helper->num_time_steps; - for (auto i : index_range(extra_ids)) - exio_helper->read_elemental_var_values(_extra_integer_vars[i], last_step, elem_ids[i]); + + // Read variables for extra integer IDs + auto extra_integer_values = + exio_helper->read_extra_integers(_extra_integer_vars); // Read in the element connectivity for each block. int nelem_last_block = 0; @@ -504,44 +502,7 @@ void ExodusII_IO::read (const std::string & fname) // Assign extra integer IDs for (auto & id : extra_ids) - { - const Real v = elem_ids[id][elem->id()]; - - if (v == Real(-1)) - { - elem->set_extra_integer(id, DofObject::invalid_id); - continue; - } - - // Ignore FE_INVALID here even if we've enabled FPEs; a - // thrown exception is preferred over an FPE signal. - FPEDisabler disable_fpes; - const long long iv = std::llround(v); - - // Check if the real number is outside of the range we can - // convert exactly - - long long max_representation = 1; - max_representation = (max_representation << std::min(std::numeric_limits::digits, - std::numeric_limits::digits)); - libmesh_error_msg_if(iv > max_representation, - "Error! An element integer value higher than " - << max_representation - << " was found! Exodus uses real numbers for storing element " - " integers, which can only represent integers from 0 to " - << max_representation - << "."); - - libmesh_error_msg_if(iv < 0, - "Error! An element integer value less than -1" - << " was found! Exodus uses real numbers for storing element " - " integers, which can only represent integers from 0 to " - << max_representation - << "."); - - - elem->set_extra_integer(id, cast_int(iv)); - } + elem->set_extra_integer(id, extra_integer_values[id][elem->id()]); // Set all the nodes for this element // diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 754ce900bf..000c2e257c 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -2131,6 +2131,77 @@ void ExodusII_IO_Helper::read_elemental_var_values(std::string elemental_var_nam +std::vector> +ExodusII_IO_Helper::read_extra_integers + (const std::vector & extra_integer_var_names) +{ + std::vector> + extra_integer_values(extra_integer_var_names.size()); + + if (extra_integer_var_names.empty()) + return extra_integer_values; + + // Make sure we have an up-to-date count of the number of time steps in the file. + this->read_num_time_steps(); + + // Prepare to check if each real number is outside of the + // range we can convert exactly + + const int exodus_digits = _single_precision ? + std::numeric_limits::digits : + std::numeric_limits::digits; + + const int shift = std::min(std::numeric_limits::digits, + exodus_digits); + + const long long max_representation = 1LL << shift; + + for (auto i : index_range(extra_integer_var_names)) + { + std::map raw_vals; + + // Read element extra "integers" as doubles from the last time step + this->read_elemental_var_values(extra_integer_var_names[i], this->num_time_steps, raw_vals); + + // Convert doubles to actual integers, at least within the + // largest convex subset of doubles where this is a bijection. + auto & values = extra_integer_values[i]; + + for (auto [elem_id, extra_val] : raw_vals) + { + if (extra_val == Real(-1)) + { + values[elem_id] = DofObject::invalid_id; + continue; + } + + // Ignore FE_INVALID here even if we've enabled FPEs; a + // thrown exception is preferred over an FPE signal. + FPEDisabler disable_fpes; + const long long int_val = std::llround(extra_val); + + libmesh_error_msg_if(int_val > max_representation, + "Error! An element integer value higher than " + << max_representation + << " was found! Exodus uses real numbers for storing element " + " integers, which can only represent integers from 0 to " + << max_representation << "."); + + libmesh_error_msg_if(int_val < 0, + "Error! An element integer value less than -1" + << " was found! Exodus uses real numbers for storing element " + " integers, which can only represent integers from 0 to " + << max_representation << "."); + + values[elem_id] = cast_int(int_val); + } + } + + return extra_integer_values; +} + + + dof_id_type ExodusII_IO_Helper::get_libmesh_node_id(int exodus_node_id) { return this->get_libmesh_id(exodus_node_id, this->node_num_map); From 7a54db01c1f86e5e4f317aa9d709161abe23f3d1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 19 Feb 2026 10:08:24 -0600 Subject: [PATCH 6/6] Be more efficient when we can get away with it --- src/mesh/exodusII_io_helper.C | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 000c2e257c..499c9b9f18 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -4472,7 +4472,11 @@ void ExodusII_IO_Helper::write_element_values ex_err = exII::ex_get_variable_param(ex_id, exII::EX_ELEM_BLOCK, &num_elem_vars); EX_CHECK_ERR(ex_err, "Error reading number of elemental variables."); - this->build_subdomain_map(mesh, false); + // We might be appending values to an existing file, in which case + // we haven't done an initialize() and we need to build the + // subdomain map here. + if (this->_subdomain_map.empty()) + this->build_subdomain_map(mesh, false); // Use mesh.n_elem() to access into the values vector rather than // the number of elements the Exodus writer thinks the mesh has,