diff --git a/include/mesh/exodusII_io.h b/include/mesh/exodusII_io.h index 5296c9bd393..e2e6c409703 100644 --- a/include/mesh/exodusII_io.h +++ b/include/mesh/exodusII_io.h @@ -140,6 +140,45 @@ class ExodusII_IO : public MeshInput, */ void write_complex_magnitude (bool val); + /** + * If true, this flag enforces the following behaviors: + * + * .) When reading an Exodus file, instead of setting the + * elem_num_map and node_num_map ids directly on Nodes and Elems + * read in from the Exodus file, set their unique_ids + * instead. Normally libmesh chooses the unique_ids automatically, + * but in this case we override that choice. As a consequence, + * the libMesh Elems/Nodes will be numbered based on their + * implied ordering in the exo file (sequentially by block and + * without any gaps in the numbering). + * .) When writing an Exodus file, populate the elem_num_map and + * node_num_map with the unique_ids of the Elems/Nodes being + * written. + * + * There are a few reasons why one might want to employ this + * behavior. One case is if you are using a ReplicatedMesh and the + * {node,elem}_num_map contains "large" (much larger than num_nodes) + * ids. This could happen either because the {node,elem}_num_map + * starts with a large id, or because there is a large gap in the + * {node,elem}_num_map numbering. Because the ReplicatedMesh always + * constructs std::vectors of size "largest_id" rather than size + * "num_nodes" or "num_elems", this mesh will be stored + * inefficiently, allocating nullptrs in the std::vector for the + * "missing" Nodes/Elems. Setting this flag to true will ensure that + * only vectors of size "num_nodes" and "num_elems" are allocated. + * + * Another case is when you delete Nodes/Elems during a simulation + * and/or renumber them, but don't want the original Nodes/Elems to + * be renumbered when the Mesh is written back out to an Exodus + * file. Setting this flag to true causes the original Node/Elem + * ids to be "saved" in the unique_id fields, so that they can later + * be written back to disk with their original numbering. + * + * This capability is relatively new, and should be considered + * experimental. + */ + void set_unique_ids_from_maps (bool val); + /** * By default, we only write out the elements physically stored in * the mesh. If we have any SIDE_DISCONTINUOUS variables, however, @@ -670,6 +709,14 @@ class ExodusII_IO : public MeshInput, */ bool _write_complex_abs; + /** + * Set Elem/Node unique_ids based on the elem_num_map and node_num_map + * contents during reading, populate those maps with unique_ids during + * writing. + * See also: docs for set_unique_ids_from_maps(bool) function. + */ + bool _set_unique_ids_from_maps; + /** * Set to true (false is the default) to generate independent nodes * for every Bezier Extraction element. diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index 8a74d34506a..79fc1d40809 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -61,6 +61,7 @@ namespace libMesh // Forward declarations class MeshBase; +class DofObject; /** * This is the \p ExodusII_IO_Helper class. This class hides the @@ -305,6 +306,44 @@ class ExodusII_IO_Helper : public ParallelObject int time_step, std::map & elem_var_value_map); + /** + * Helper function that takes a (1-based) Exodus node/elem id and + * determines the corresponding libMesh Node/Elem id. Takes into account + * whether the user has chosen to set the Node/Elem unique ids based on + * the {node,elem}_num_map or to let libMesh set them. + */ + dof_id_type get_libmesh_node_id(int exodus_node_id); + dof_id_type get_libmesh_elem_id(int exodus_elem_id); + + /** + * Helper function that conditionally sets the unique_id of the + * passed-in Node/Elem. Calling this function does nothing if + * _set_unique_ids_from_maps == false, otherwise it sets the + * unique_id based on the entries of the {node,elem_num_map}. The + * input index is assumed to be a zero-based index into the + * {node,elem}_num_map array. + */ + void conditionally_set_node_unique_id( + MeshBase & mesh, Node * node, int zero_based_node_num_map_index); + void conditionally_set_elem_unique_id( + MeshBase & mesh, Elem * elem, int zero_based_elem_num_map_index); + +private: + + /** + * Internal implementation for the two sets of functions above. + */ + dof_id_type get_libmesh_id( + int exodus_id, + const std::vector & num_map); + + void set_dof_object_unique_id( + MeshBase & mesh, + DofObject * dof_object, + int exodus_mapped_id); + +public: + /** * Opens an \p ExodusII mesh file named \p filename for writing. */ @@ -850,6 +889,10 @@ class ExodusII_IO_Helper : public ParallelObject // On/Off message flag bool verbose; + // Same as the ExodusII_IO flag by the same name. This flag is + // also set whenever ExodusII_IO::set_unique_ids_from_maps() is called. + bool set_unique_ids_from_maps; + // This flag gets set after the Exodus file has been successfully opened for writing. // Both the create() and open() (if called with EX_WRITE) functions may set this flag. bool opened_for_writing; diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index c539616a510..7e9c8a8edab 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -458,7 +458,7 @@ class MeshBase : public ParallelObject /** * \returns The next unique id to be used. */ - unique_id_type next_unique_id() { return _next_unique_id; } + unique_id_type next_unique_id() const { return _next_unique_id; } /** * Sets the next available unique id to be used. On a diff --git a/src/mesh/exodusII_io.C b/src/mesh/exodusII_io.C index b62db85cb87..3cb657da79f 100644 --- a/src/mesh/exodusII_io.C +++ b/src/mesh/exodusII_io.C @@ -134,6 +134,7 @@ ExodusII_IO::ExodusII_IO (MeshBase & mesh, #endif _allow_empty_variables(false), _write_complex_abs(true), + _set_unique_ids_from_maps(false), _disc_bex(false) { // if !LIBMESH_HAVE_EXODUS_API, we didn't use this @@ -157,6 +158,7 @@ ExodusII_IO::ExodusII_IO (const MeshBase & mesh, #endif _allow_empty_variables(false), _write_complex_abs(true), + _set_unique_ids_from_maps(false), _disc_bex(false) { // if !LIBMESH_HAVE_EXODUS_API, we didn't use this @@ -325,23 +327,30 @@ void ExodusII_IO::read (const std::string & fname) std::unordered_map spline_nodeelem_ptrs; // Loop over the nodes, create Nodes with local processor_id 0. - for (int i=0; inum_nodes; i++) + for (auto i : make_range(exio_helper->num_nodes)) { - // Use the node_num_map to get the correct ID for Exodus - int exodus_id = exio_helper->node_num_map[i]; + // Determine the libmesh node id implied by "i". The + // get_libmesh_node_id() helper function expects a 1-based + // Exodus node id, so we construct the "implied" Exodus node id + // from "i" by adding 1. + auto libmesh_node_id = exio_helper->get_libmesh_node_id(/*exodus_node_id=*/i+1); // Catch the node that was added to the mesh - Node * added_node = mesh.add_point (Point(exio_helper->x[i], exio_helper->y[i], exio_helper->z[i]), exodus_id-1); + Node * added_node = mesh.add_point (Point(exio_helper->x[i], exio_helper->y[i], exio_helper->z[i]), libmesh_node_id); - // If the Mesh assigned an ID different from what is in the - // Exodus file, we should probably error. - libmesh_error_msg_if(added_node->id() != static_cast(exodus_id-1), + // Sanity check: throw an error if the Mesh assigned an ID to + // the Node which does not match the libmesh_id we just determined. + libmesh_error_msg_if(added_node->id() != static_cast(libmesh_node_id), "Error! Mesh assigned node ID " << added_node->id() << " which is different from the (zero-based) Exodus ID " - << exodus_id-1 + << libmesh_node_id << "!"); + // If the _set_unique_ids_from_maps flag is true, set the + // unique_id for "node", otherwise do nothing. + exio_helper->conditionally_set_node_unique_id(mesh, added_node, i); + // If we have a set of spline weights, these nodes are going to // be used as control points for Bezier elements, and we need // to attach a NodeElem to each to make sure it doesn't get @@ -361,6 +370,13 @@ void ExodusII_IO::read (const std::string & fname) // Give the NodeElem ids at the end, so we can match any // existing ids in the file for other elements + // + // We don't set the unique_id for this NodeElem here even if + // the user has set the _set_unique_ids_from_maps flag + // because these NodeElems don't have entries in the + // elem_num_map. Therefore, we just let the Mesh assign + // whatever unique_id is "next" as the Elem is added to the + // Mesh. elem->set_id() = exio_helper->end_elem_id() + i; elem->set_node(0, added_node); @@ -449,19 +465,13 @@ void ExodusII_IO::read (const std::string & fname) // Assign the current subdomain to this Elem uelem->subdomain_id() = static_cast(subdomain_id); - // Use the elem_num_map to obtain the ID of this element in - // the Exodus file. Make sure we aren't reading garbage if - // the file is corrupt. - libmesh_error_msg_if(std::size_t(j) >= exio_helper->elem_num_map.size(), - "Error: Trying to read Exodus file with more elements than elem_num_map entries.\n"); - int exodus_id = exio_helper->elem_num_map[j]; + // Determine the libmesh elem id implied by "j". The + // ExodusII_IO_Helper::get_libmesh_elem_id() helper function + // expects a 1-based Exodus elem id, so we construct the + // "implied" Exodus elem id from "j" by adding 1. + auto libmesh_elem_id = exio_helper->get_libmesh_elem_id(/*exodus_elem_id=*/j+1); - // Assign this element the same ID it had in the Exodus - // file, but make it zero-based by subtracting 1. Note: - // some day we could use 1-based numbering in libmesh and - // thus match the Exodus numbering exactly, but at the - // moment libmesh is zero-based. - uelem->set_id(exodus_id-1); + uelem->set_id(libmesh_elem_id); // Record that we have seen an element of dimension uelem->dim() elems_of_dimension[uelem->dim()] = true; @@ -469,13 +479,17 @@ void ExodusII_IO::read (const std::string & fname) // Catch the Elem pointer that the Mesh throws back Elem * elem = mesh.add_elem(std::move(uelem)); - // If the Mesh assigned an ID different from what is in the - // Exodus file, we should probably error. - libmesh_error_msg_if(elem->id() != static_cast(exodus_id-1), + // If the _set_unique_ids_from_maps flag is true, set the + // unique_id for "elem", otherwise do nothing. + exio_helper->conditionally_set_elem_unique_id(mesh, elem, j); + + // If the Mesh assigned an ID different from the one we + // tried to give it, we should probably error. + libmesh_error_msg_if(elem->id() != static_cast(libmesh_elem_id), "Error! Mesh assigned ID " << elem->id() << " which is different from the (zero-based) Exodus ID " - << exodus_id-1 + << libmesh_elem_id << "!"); // Assign extra integer IDs @@ -528,17 +542,14 @@ void ExodusII_IO::read (const std::string & fname) { for (int k=0; knum_nodes_per_elem; k++) { - // global index - int gi = (elem_num)*exio_helper->num_nodes_per_elem + conv.get_node_map(k); - - // The entries in 'connect' are actually (1-based) - // indices into the node_num_map, so to get the right - // node ID we: - // 1.) Subtract 1 from connect[gi] - // 2.) Pass it through node_num_map to get the corresponding Exodus ID - // 3.) Subtract 1 from that, since libmesh node numbering is "zero"-based, - // even when the Exodus node numbering doesn't start with 1. - int libmesh_node_id = exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1; + // Get index into this block's connectivity array + int gi = elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k); + + // Get the 1-based Exodus node id from the "connect" array + auto exodus_node_id = exio_helper->connect[gi]; + + // Convert this index to a libMesh Node id + auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id); // Set the node pointer in the Elem elem->set_node(k, mesh.node_ptr(libmesh_node_id)); @@ -609,10 +620,13 @@ void ExodusII_IO::read (const std::string & fname) libmesh_vector_at(coef_vec, elem_node_index); // Get the libMesh node corresponding to that row - const int gi = (elem_num)*exio_helper->bex_num_elem_cvs + - spline_node_index; - const dof_id_type libmesh_node_id = - exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1; + const int gi = elem_num * exio_helper->bex_num_elem_cvs + spline_node_index; + + // Get the 1-based Exodus node id from the "connect" array + auto exodus_node_id = exio_helper->connect[gi]; + + // Convert this index to a libMesh Node id + auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id); if (coef != 0) // Ignore irrelevant spline nodes key.emplace_back(libmesh_node_id, coef); @@ -703,15 +717,10 @@ void ExodusII_IO::read (const std::string & fname) for (auto e : index_range(exio_helper->elem_list)) { - // The numbers in the Exodus file sidesets should be thought - // of as (1-based) indices into the elem_num_map array. So, - // to get the right element ID we have to: - // 1.) Subtract 1 from elem_list[e] (to get a zero-based index) - // 2.) Pass it through elem_num_map (to get the corresponding Exodus ID) - // 3.) Subtract 1 from that, since libmesh is "zero"-based, - // even when the Exodus numbering doesn't start with 1. + // Call helper function to get the libmesh Elem id for the + // e'th entry in the current elem_list. dof_id_type libmesh_elem_id = - cast_int(exio_helper->elem_num_map[exio_helper->elem_list[e] - 1] - 1); + exio_helper->get_libmesh_elem_id(exio_helper->elem_list[e]); // Set any relevant node/edge maps for this element Elem & elem = mesh.elem_ref(libmesh_elem_id); @@ -804,14 +813,10 @@ void ExodusII_IO::read (const std::string & fname) std::map elem_to_elemsets; for (auto e : index_range(exio_helper->elemset_list)) { - // Follow standard (see sideset case above) approach for - // converting the ids stored in the elemset_list to - // libmesh Elem ids. - // - // TODO: this should be moved to a helper function so we - // don't duplicate the code. - dof_id_type libmesh_elem_id = - cast_int(exio_helper->elem_num_map[exio_helper->elemset_list[e] - 1] - 1); + // Call helper function to get the libmesh Elem id for the + // e'th entry in the current elemset_list. + dof_id_type libmesh_elem_id = + exio_helper->get_libmesh_elem_id(exio_helper->elemset_list[e]); // Get a pointer to this Elem Elem * elem = mesh.elem_ptr(libmesh_elem_id); @@ -887,22 +892,9 @@ void ExodusII_IO::read (const std::string & fname) for (int i=0; inum_nodes_per_set[nodeset]; ++i) { - int exodus_id = exio_helper->node_sets_node_list[i + offset]; - - // It's possible for nodesets to have invalid ids in them - // by accident. Instead of possibly accessing past the - // end of node_num_map, let's make sure we have that many - // entries. - libmesh_error_msg_if(static_cast(exodus_id - 1) >= exio_helper->node_num_map.size(), - "Invalid Exodus node id " << exodus_id - << " found in nodeset " << nodeset_id); - - // As before, the entries in 'node_list' are 1-based - // indices into the node_num_map array, so we have to map - // them. See comment above. - int libmesh_node_id = exio_helper->node_num_map[exodus_id - 1] - 1; - mesh.get_boundary_info().add_node(cast_int(libmesh_node_id), - nodeset_id); + int exodus_node_id = exio_helper->node_sets_node_list[i + offset]; + auto libmesh_node_id = exio_helper->get_libmesh_node_id(exodus_node_id); + mesh.get_boundary_info().add_node(libmesh_node_id, nodeset_id); } } } @@ -971,7 +963,15 @@ void ExodusII_IO::write_complex_magnitude (bool val) _write_complex_abs = val; } +void ExodusII_IO::set_unique_ids_from_maps (bool val) +{ + _set_unique_ids_from_maps = val; + // Set this flag on the helper object as well. The helper needs to know about this + // flag, since it sometimes needs to construct libmesh Node ids from nodal connectivity + // arrays (see e.g. ExodusII_IO_Helper::read_edge_blocks()). + exio_helper->set_unique_ids_from_maps = val; +} void ExodusII_IO::use_mesh_dimension_instead_of_spatial_dimension(bool val) { diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 9f3b6078c82..bc49f9f462d 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -286,6 +286,7 @@ ExodusII_IO_Helper::ExodusII_IO_Helper(const ParallelObject & parent, num_nodal_vars(0), num_elem_vars(0), verbose(v), + set_unique_ids_from_maps(false), opened_for_writing(false), opened_for_reading(false), _run_only_on_proc0(run_only_on_proc0), @@ -1408,15 +1409,13 @@ void ExodusII_IO_Helper::read_edge_blocks(MeshBase & mesh) // Loop over indices in connectivity array, build edge elements, // look them up in the edge_map. - for (unsigned int i=0, sz=connect.size(); iconnect[i+n]; + dof_id_type libmesh_node_id = this->get_libmesh_node_id(exodus_node_id); edge->set_node(n, mesh.node_ptr(libmesh_node_id)); } @@ -1487,6 +1486,12 @@ void ExodusII_IO_Helper::read_elem_num_map () if (num_elem) { + // The elem_num_map may contain ids larger than num_elem. In + // other words, the elem_num_map is not necessarily just a + // permutation of the "trivial" 1,2,3,... mapping, it can + // contain effectively "any" numbers. Therefore, to get + // "_end_elem_id", we need to check what the max entry in the + // elem_num_map is. auto it = std::max_element(elem_num_map.begin(), elem_num_map.end()); _end_elem_id = *it; } @@ -1863,7 +1868,7 @@ void ExodusII_IO_Helper::read_nodal_var_values(std::string nodal_var_name, int t } // Clear out any previously read nodal variable values - nodal_var_values.clear(); + this->nodal_var_values.clear(); std::vector unmapped_nodal_var_values(num_nodes); @@ -1878,18 +1883,21 @@ void ExodusII_IO_Helper::read_nodal_var_values(std::string nodal_var_name, int t MappedInputVector(unmapped_nodal_var_values, _single_precision).data()); EX_CHECK_ERR(ex_err, "Error reading nodal variable values!"); - for (unsigned i=0; i(num_nodes); i++) + for (auto i : make_range(num_nodes)) { - libmesh_assert_less(i, this->node_num_map.size()); - - // Use the node_num_map to obtain the ID of this node in the Exodus file, - // and remember to subtract 1 since libmesh is zero-based and Exodus is 1-based. - const unsigned mapped_node_id = this->node_num_map[i] - 1; - - libmesh_assert_less(i, unmapped_nodal_var_values.size()); + // Determine the libmesh node id implied by "i". The + // get_libmesh_node_id() helper function expects a 1-based + // Exodus node id, so we construct the "implied" Exodus node id + // from "i" by adding 1. + // + // If the user has set the "set_unique_ids_from_maps" flag to + // true, then calling get_libmesh_node_id(i+1) will just return + // i, otherwise it will determine the value (with error + // checking) using this->node_num_map. + auto libmesh_node_id = this->get_libmesh_node_id(/*exodus_node_id=*/i+1); // Store the nodal value in the map. - nodal_var_values[mapped_node_id] = unmapped_nodal_var_values[i]; + this->nodal_var_values[libmesh_node_id] = unmapped_nodal_var_values[i]; } } @@ -2127,12 +2135,14 @@ void ExodusII_IO_Helper::read_elemental_var_values(std::string elemental_var_nam for (unsigned j=0; j(num_elem_this_blk); j++) { - // Use the elem_num_map to obtain the ID of this element in the Exodus file, - // and remember to subtract 1 since libmesh is zero-based and Exodus is 1-based. - unsigned mapped_elem_id = this->elem_num_map[ex_el_num] - 1; + // Determine the libmesh id of the element with zero-based + // index "ex_el_num". This function expects a one-based + // index, so we add 1 to ex_el_num when we pass it in. + auto libmesh_elem_id = + this->get_libmesh_elem_id(ex_el_num + 1); // Store the elemental value in the map. - elem_var_value_map[mapped_elem_id] = block_elem_var_values[j]; + elem_var_value_map[libmesh_elem_id] = block_elem_var_values[j]; // Go to the next sequential element ID. ex_el_num++; @@ -2141,6 +2151,100 @@ void ExodusII_IO_Helper::read_elemental_var_values(std::string elemental_var_nam } + +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); +} + +dof_id_type ExodusII_IO_Helper::get_libmesh_elem_id(int exodus_elem_id) +{ + return this->get_libmesh_id(exodus_elem_id, this->elem_num_map); +} + +dof_id_type +ExodusII_IO_Helper::get_libmesh_id(int exodus_id, + const std::vector & num_map) +{ + // The input exodus_id is assumed to be a (1-based) index into + // the {node,elem}_num_map, so in order to use exodus_id as an index + // in C++, we need to first make it zero-based. + auto exodus_id_zero_based = + cast_int(exodus_id - 1); + + // Throw an informative error message rather than accessing past the + // end of the provided num_map. If we are setting Elem unique_ids + // based on the num_map, we don't need to do this check. + if (!this->set_unique_ids_from_maps) + libmesh_error_msg_if(exodus_id_zero_based >= num_map.size(), + "Cannot get LibMesh id for Exodus id: " << exodus_id); + + // If the user set the flag which stores Exodus node/elem ids as + // unique_ids instead of regular ids, then the libmesh id we are + // looking for is actually just "exodus_id_zero_based". Otherwise, + // we need to look up the Node/Elem's id in the provided num_map, + // *and* then subtract 1 from that because the entries in the + // num_map are also 1-based. + dof_id_type libmesh_id = + this->set_unique_ids_from_maps ? + cast_int(exodus_id_zero_based) : + cast_int(num_map[exodus_id_zero_based] - 1); + + return libmesh_id; +} + + + +void +ExodusII_IO_Helper:: +conditionally_set_node_unique_id(MeshBase & mesh, Node * node, int zero_based_node_num_map_index) +{ + this->set_dof_object_unique_id(mesh, node, libmesh_vector_at(this->node_num_map, zero_based_node_num_map_index)); +} + +void +ExodusII_IO_Helper:: +conditionally_set_elem_unique_id(MeshBase & mesh, Elem * elem, int zero_based_elem_num_map_index) +{ + this->set_dof_object_unique_id(mesh, elem, libmesh_vector_at(this->elem_num_map, zero_based_elem_num_map_index)); +} + +void +ExodusII_IO_Helper::set_dof_object_unique_id( + MeshBase & mesh, + DofObject * dof_object, + int exodus_mapped_id) +{ + if (this->set_unique_ids_from_maps) + { + // Exodus ids are always 1-based while libmesh ids are always + // 0-based, so to make a libmesh unique_id here, we subtract 1 + // from the exodus_mapped_id to make it 0-based. + auto exodus_mapped_id_zero_based = + cast_int(exodus_mapped_id - 1); + + // Set added_node's unique_id to "exodus_mapped_id_zero_based". + dof_object->set_unique_id(cast_int(exodus_mapped_id_zero_based)); + + // Normally the Mesh is responsible for setting the unique_ids + // of Nodes/Elems in a consistent manner, so when we set the unique_id + // of a Node/Elem manually based on the {node,elem}_num_map, we need to + // make sure that the "next" unique id assigned by the Mesh + // will still be valid. We do this by making sure that the + // next_unique_id is greater than the one we set manually. The + // APIs for doing this are only defined when unique ids are + // enabled. +#ifdef LIBMESH_ENABLE_UNIQUE_ID + unique_id_type next_unique_id = mesh.next_unique_id(); + mesh.set_next_unique_id(std::max(next_unique_id, static_cast(exodus_mapped_id_zero_based + 1))); +#else + // Avoid compiler warnings about the unused variable + libmesh_ignore(mesh); +#endif + } +} + + // For Writing Solutions void ExodusII_IO_Helper::create(std::string filename) @@ -2433,12 +2537,22 @@ void ExodusII_IO_Helper::write_nodal_coordinates(const MeshBase & mesh, bool use #endif }; - // And in the node_num_map - since the nodes aren't organized in - // blocks, libmesh will always write out the identity map - // here... unless there has been some refinement and coarsening, or - // node deletion without a corresponding call to contract(). You - // need to write this any time there could be 'holes' in the node - // numbering, so we write it every time. + // And in the node_num_map. If the user has set the + // _set_unique_ids_from_maps flag, then we will write the Node + // unique_ids to the node_num_map, otherwise we will just write a + // trivial node_num_map, since in that we don't write the unique_id + // information to the Exodus file. In other words, set the + // _set_unique_ids_from_maps flag to true on both the reading and + // writing ExodusII_IO objects if you want to preserve the + // node_num_map information without actually renumbering the Nodes + // in libmesh according to the node_num_map. + // + // One reason why you might not want to actually renumber the Nodes + // in libmesh according to the node_num_map is that it can introduce + // undesirable large "gaps" in the numbering, e.g. Nodes numbered + // [0, 1, 1000, 10001] which is not ideal for the ReplicatedMesh + // _nodes data structure, which stores the Nodes in a contiguous + // array based on Node id. // Let's skip the node_num_map in the discontinuous and add_sides // cases, since we're effectively duplicating nodes for the sake of @@ -2462,17 +2576,27 @@ void ExodusII_IO_Helper::write_nodal_coordinates(const MeshBase & mesh, bool use // Fill in node_num_map entry with the proper (1-based) node // id, unless we're not going to be able to keep the map up - // later. + // later. If the user has chosen to _set_unique_ids_from_maps, + // then we fill up the node_num_map with (1-based) unique + // ids rather than node ids. if (!_add_sides) - node_num_map.push_back(node.id() + 1); + { + if (this->set_unique_ids_from_maps) + node_num_map.push_back(node.unique_id() + 1); + else + node_num_map.push_back(node.id() + 1); + } - // Also map the zero-based libmesh node id to the 1-based - // Exodus ID it will be assigned (this is equivalent to the - // current size of the x vector). + // Also map the zero-based libmesh node id to the (1-based) + // index in the node_num_map it corresponds to + // (this is equivalent to the current size of the "x" vector, + // so we just use x.size()). This map is used to look up + // an Exodus Node id given a libMesh Node id, so it does + // involve unique_ids. libmesh_node_num_to_exodus[ cast_int(node.id()) ] = cast_int(x.size()); - } + } // end for (node_ptr) } - else + else // use_discontinuous { for (const auto & elem : mesh.active_element_ptr_range()) for (const Node & node : elem->node_ref_range()) @@ -2616,8 +2740,9 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti ++counter; } - elem_num_map.resize(num_elem); - std::vector::iterator curr_elem_map_end = elem_num_map.begin(); + // Here we reserve() space so that we can push_back() onto the + // elem_num_map in the loops below. + this->elem_num_map.reserve(num_elem); // In the case of discontinuous plotting we initialize a map from // (element, node) pairs to the corresponding discontinuous node index. @@ -2820,7 +2945,16 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // We need these later if we're adding fake sides, but we don't need // to recalculate it. auto num_elem_this_blk_it = num_elem_this_blk_vec.begin(); + + // We write "fake" ids to the elem_num_map when adding fake sides. + // I don't think it's too important exactly what fake ids are used, + // as long as they don't conflict with any other ids that are + // already in the elem_num_map. auto next_fake_id = mesh.max_elem_id() + 1; // 1-based numbering in Exodus +#ifdef LIBMESH_ENABLE_UNIQUE_ID + if (this->set_unique_ids_from_maps) + next_fake_id = mesh.next_unique_id(); +#endif for (auto & [subdomain_id, element_id_vec] : subdomain_map) { @@ -2841,9 +2975,9 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // element ids to add. if (subdomain_id < subdomain_id_end) { - connect.resize(element_id_vec.size()*num_nodes_per_elem); + connect.resize(element_id_vec.size()*num_nodes_per_elem); - for (auto i : index_range(element_id_vec)) + for (auto i : index_range(element_id_vec)) { unsigned int elem_id = element_id_vec[i]; libmesh_elem_num_to_exodus[elem_id] = ++libmesh_elem_num_to_exodus_counter; // 1-based indexing for Exodus @@ -2891,24 +3025,22 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti libmesh_map_find(discontinuous_node_indices, std::make_pair(elem_id, elem_node_index)); } - } - } + } // end for(j) + + // push_back() either elem_id+1 or the current Elem's + // unique_id+1 into the elem_num_map, depending on the value + // of the set_unique_ids_from_maps flag. + if (this->set_unique_ids_from_maps) + this->elem_num_map.push_back(elem.unique_id() + 1); + else + this->elem_num_map.push_back(elem_id + 1); - // This transform command stores its result in a range that - // begins at the third argument, so this command is adding - // values to the elem_num_map vector starting from - // curr_elem_map_end. Here we add 1 to each id to make a - // 1-based exodus file. - curr_elem_map_end = std::transform - (element_id_vec.begin(), - element_id_vec.end(), - curr_elem_map_end, - [](dof_id_type id){return id+1;}); + } // end for(i) } - // If this is a "fake" block of added sides, we build those as - // we go. - else + else // subdomain_id >= subdomain_id_end { + // If this is a "fake" block of added sides, we build those as + // we go. libmesh_assert(_add_sides); libmesh_assert(num_elem_this_blk_it != num_elem_this_blk_vec.end()); @@ -2944,12 +3076,11 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti } } - auto old_curr_map_end = curr_elem_map_end; - curr_elem_map_end += num_elem_this_blk; - - std::generate - (old_curr_map_end, curr_elem_map_end, - [&next_fake_id](){return next_fake_id++;}); + // Store num_elem_this_blk "fake" ids into the + // elem_num_map. Use a traditional for-loop to avoid unused + // variable warnings about the loop counter. + for (int i=0; ielem_num_map.push_back(next_fake_id++); } ++num_elem_this_blk_it; @@ -2962,7 +3093,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) // write out the element number map that we created ex_err = exII::ex_put_elem_num_map(ex_id, elem_num_map.data()); diff --git a/tests/Makefile.am b/tests/Makefile.am index acf0ad8239d..b7c9f423971 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -217,6 +217,8 @@ data = matrices/geom_1_extraction_op.h5 \ meshes/non_manifold_junction1.exo \ meshes/non_manifold_junction2.exo \ meshes/non_manifold_junction3.exo \ + meshes/nontrivial_elem_num_map.exo \ + meshes/nontrivial_node_num_map.exo \ meshes/Cluster_34.stl \ meshes/engraving.stl \ meshes/quad4_tri3_smoothed.xda.gz \ diff --git a/tests/Makefile.in b/tests/Makefile.in index d2ab60962dd..37f2b705cb6 100644 --- a/tests/Makefile.in +++ b/tests/Makefile.in @@ -2386,6 +2386,8 @@ data = matrices/geom_1_extraction_op.h5 \ meshes/non_manifold_junction1.exo \ meshes/non_manifold_junction2.exo \ meshes/non_manifold_junction3.exo \ + meshes/nontrivial_elem_num_map.exo \ + meshes/nontrivial_node_num_map.exo \ meshes/Cluster_34.stl \ meshes/engraving.stl \ meshes/quad4_tri3_smoothed.xda.gz \ diff --git a/tests/mesh/mesh_input.C b/tests/mesh/mesh_input.C index 4dd3a2ff749..c08c37e6663 100644 --- a/tests/mesh/mesh_input.C +++ b/tests/mesh/mesh_input.C @@ -112,6 +112,8 @@ public: CPPUNIT_TEST( testExodusCopyNodalSolutionReplicated ); CPPUNIT_TEST( testExodusCopyElementSolutionReplicated ); CPPUNIT_TEST( testExodusReadHeader ); + CPPUNIT_TEST( testExodusSetNodeUniqueIdsFromMaps ); + CPPUNIT_TEST( testExodusSetElemUniqueIdsFromMaps ); #if LIBMESH_DIM > 2 CPPUNIT_TEST( testExodusIGASidesets ); CPPUNIT_TEST( testLowOrderEdgeBlocks ); @@ -333,6 +335,130 @@ public: #ifdef LIBMESH_HAVE_EXODUS_API + + void testExodusSetNodeUniqueIdsFromMaps_implementation( + bool set_unique_ids, + const std::vector & expected_unique_ids) + { + // This test requires that libmesh is compiled with unique_ids enabled +#ifdef LIBMESH_ENABLE_UNIQUE_ID + { + ReplicatedMesh mesh(*TestCommWorld); + ExodusII_IO exii(mesh); + + // Set Node/Elem unique ids based on the node/elem_num_map + exii.set_unique_ids_from_maps(set_unique_ids); + + // Read the mesh + exii.read("meshes/nontrivial_node_num_map.exo"); + + // Verify the results. + for (auto i : index_range(expected_unique_ids)) + { + // Debugging: + // libMesh::out << "unique_id for node " << i + // << " = " << mesh.node_ptr(i)->unique_id() + // << std::endl; + + CPPUNIT_ASSERT_EQUAL(mesh.node_ptr(i)->unique_id(), expected_unique_ids[i]); + } + } +#else + // Prevent compiler warnings about unused variables when + // unique_ids are not enabled. + libmesh_ignore(set_unique_ids, expected_unique_ids); +#endif // LIBMESH_ENABLE_UNIQUE_ID + } // end testExodusSetNodeUniqueIdsFromMaps_implementation() + + + + void testExodusSetNodeUniqueIdsFromMaps() + { + LOG_UNIT_TEST; + + // node_num_map = 1, 3, 9, 8, 2, 5, 7, 6, 4 + // The input is a (zero-based) version of the node_num_map as it + // exists in the file. + this->testExodusSetNodeUniqueIdsFromMaps_implementation( + /*set_unique_ids=*/true, + /*expected_unique_ids=*/{0, 2, 8, 7, 1, 4, 6, 5, 3}); + + // The input is the (zero-based) _position_ of the (one-based) id + // "i" in the node_num_map, as it exists in the file. For example, + // zero-based Node id 3 corresponds to one-based Node id 4, which + // appears in (zero-based) position 8 in the node_num_map, hence: + // 3 -> 8, etc. + this->testExodusSetNodeUniqueIdsFromMaps_implementation( + /*set_unique_ids=*/false, + /*expected_unique_ids=*/{0, 4, 1, 8, 5, 7, 6, 3, 2}); + } + + void testExodusSetElemUniqueIdsFromMaps_implementation( + bool set_unique_ids, + const std::vector & expected_unique_ids) + { + // This test requires that libmesh is compiled with unique_ids enabled +#ifdef LIBMESH_ENABLE_UNIQUE_ID + { + ReplicatedMesh mesh(*TestCommWorld); + ExodusII_IO exii(mesh); + + // Set Node/Elem unique ids based on the node/elem_num_map + exii.set_unique_ids_from_maps(set_unique_ids); + + // Read the mesh + exii.read("meshes/nontrivial_elem_num_map.exo"); + + // Verify the results. + auto expected_it = expected_unique_ids.begin(); + for (const auto & elem : mesh.element_ptr_range()) + { + // Debugging: + // libMesh::out << "unique_id for Elem " << elem->id() + // << " = " << elem->unique_id() + // << std::endl; + + CPPUNIT_ASSERT_EQUAL(elem->unique_id(), *expected_it++); + } + } +#else + // Prevent compiler warnings about unused variables when + // unique_ids are not enabled. + libmesh_ignore(set_unique_ids, expected_unique_ids); +#endif // LIBMESH_ENABLE_UNIQUE_ID + } + + void testExodusSetElemUniqueIdsFromMaps() + { + LOG_UNIT_TEST; + + // The mesh used in this test has a non-trivial elem_num_map with + // a non-contiguous numbering that has a "gap" at the beginning + // and is also missing "16": + // elem_num_map = 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, + // 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, + // 45, 46 ; + { + // When we assign (zero-based) elem_num_map entries as unique_ids, the + // unique_ids are just a zero-based version of the entries above. + std::vector expected_unique_ids = { + 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45}; + this->testExodusSetElemUniqueIdsFromMaps_implementation( + /*set_unique_ids=*/true, expected_unique_ids); + } + + { + // When libmesh assigns unique_ids, all the Nodes (50) are + // numbered first, so the first Elem is assigned a unique_id of + // 50 and the rest are sequential from there. + std::vector expected_unique_ids(/*size=*/35); + std::iota(expected_unique_ids.begin(), expected_unique_ids.end(), /*start=*/50); + this->testExodusSetElemUniqueIdsFromMaps_implementation( + /*set_unique_ids=*/false, expected_unique_ids); + } + } + void testExodusReadHeader () { LOG_UNIT_TEST; diff --git a/tests/meshes/nontrivial_elem_num_map.exo b/tests/meshes/nontrivial_elem_num_map.exo new file mode 100644 index 00000000000..39b7778c398 Binary files /dev/null and b/tests/meshes/nontrivial_elem_num_map.exo differ diff --git a/tests/meshes/nontrivial_node_num_map.exo b/tests/meshes/nontrivial_node_num_map.exo new file mode 100644 index 00000000000..cfec214e9c7 Binary files /dev/null and b/tests/meshes/nontrivial_node_num_map.exo differ