diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 0601c18bd6..5cf7e4fe95 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -206,50 +206,87 @@ const std::vector prism_inverse_face_map = {4, 1, 2, 3, 5}; subdomain_id_type & subdomain_id_end, int & next_block_id) { - 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; - } + // Exodus requires every element block to contain a single element + // type. Group elements first by subdomain id and then by element + // type, so that a subdomain containing more than one element type + // (for example a mix of C0POLYHEDRON and HEX8 cells) is split across + // multiple blocks rather than written as a single, invalid block. + std::map>> elems_by_subdomain_type; + subdomain_id_type max_subdomain_id = 0; + bool have_elems = false; - // 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. + // don't put them in the map. That way the number of blocks + // should be correct. if (elem->infinite()) continue; - subdomain_map[ elem->subdomain_id() ].push_back(elem->id()); + const subdomain_id_type sbd_id = elem->subdomain_id(); + elems_by_subdomain_type[sbd_id][elem->type()].push_back(elem->id()); + max_subdomain_id = have_elems ? std::max(max_subdomain_id, sbd_id) : sbd_id; + have_elems = true; + } + + // Assign a block id to each (subdomain, element type) group. The + // first element type in each subdomain keeps the subdomain id as its + // block id, so single-type subdomains (the common case) are written + // exactly as before. Any additional element types in the same + // subdomain get synthesized block ids allocated above all existing + // subdomain ids. + std::map> subdomain_map; + subdomain_id_type next_synth_block_id = have_elems ? max_subdomain_id + 1 : 0; + + for (auto & [sbd_id, type_map] : elems_by_subdomain_type) + { + bool first_type = true; + for (auto & [elem_t, elem_ids] : type_map) + { + libmesh_ignore(elem_t); + const subdomain_id_type block_id = + first_type ? sbd_id : next_synth_block_id++; + first_type = false; + subdomain_map[block_id] = std::move(elem_ids); + } + } + + // Real element blocks occupy the ids below subdomain_id_end; any + // blocks synthesized for visualization sides are numbered after them. + if (!subdomain_map.empty()) + subdomain_id_end = subdomain_map.rbegin()->first + 1; + + // If we've been asked to add side elements, those 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 add an entry to keep track of + // the number of elements we'll add in each new block. + if (add_sides) + for (const auto & elem : mesh.active_element_ptr_range()) + { + if (elem->infinite()) + continue; - // 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)]; + const subdomain_id_type side_block_id = + cast_int(subdomain_id_end + elem->side_type(s)); + + // Guard the invariant above (also catches subdomain_id_type + // overflow wrapping a side block id back down onto an + // element block id). + libmesh_assert_greater_equal(side_block_id, subdomain_id_end); + + auto & marker = subdomain_map[side_block_id]; if (marker.empty()) marker.push_back(1); else ++marker[0]; } - } - - if (!add_sides && !subdomain_map.empty()) - subdomain_id_end = subdomain_map.rbegin()->first + 1; + } // Allocate optional block IDs after both mesh subdomains and any blocks // synthesized for visualization sides. @@ -3021,7 +3058,43 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti num_elem_this_blk_vec.push_back (cast_int(element_id_vec.size())); - std::string block_name = mesh.subdomain_name(subdomain_id); + // A block id normally *is* a subdomain id, but a subdomain that + // contains more than one element type is split across several + // blocks by build_subdomain_map(): only its first element type + // keeps the subdomain id, while every other type gets a + // synthesized block id allocated above all subdomain ids. We + // detect such synthesized blocks from the mismatch between the + // block id and the actual subdomain of the elements it holds. + const subdomain_id_type elem_subdomain_id = + mesh.elem_ref(element_id_vec[0]).subdomain_id(); + const bool is_synthesized_block = (subdomain_id != elem_subdomain_id); + + std::string block_name = mesh.subdomain_name(elem_subdomain_id); + + if (is_synthesized_block) + { + // Prefix the (original) subdomain's name with the element + // type, so the synthesized block is identifiable and its + // name stays distinct from the block that kept the + // subdomain id. This block becomes its own subdomain when + // the file is read back in. + const std::string type_suffix = Utility::enum_to_string(elem_t); + block_name = block_name.empty() + ? (std::to_string(elem_subdomain_id) + "_" + type_suffix) + : (block_name + "_" + type_suffix); + + // Informational: keep it to one rank so parallel (Nemesis) + // writes don't repeat it once per processor. + // NOTE might miss logs if only occurs on other ranks + if (this->processor_id() == 0) + libMesh::out << "ExodusII_IO: subdomain " << elem_subdomain_id + << " contains more than one element type; writing its " + << type_suffix << " elements to a new block \"" << block_name + << "\" (block id " << subdomain_id + << "), which will be read back as a separate subdomain." + << std::endl; + } + if (block_name.empty() && elem_t == C0POLYGON) block_name = "NSIDED_" + std::to_string(counter + 1); if (block_name.empty() && elem_t == C0POLYHEDRON) diff --git a/tests/mesh/exodus_test.C b/tests/mesh/exodus_test.C index 32700ea1e9..128a64b395 100644 --- a/tests/mesh/exodus_test.C +++ b/tests/mesh/exodus_test.C @@ -188,6 +188,7 @@ public: CPPUNIT_TEST(test_write_cube_header); CPPUNIT_TEST(test_write_hexagonal_prism_header); CPPUNIT_TEST(test_write_and_read_hexagonal_prism); + CPPUNIT_TEST(test_write_and_read_mixed_poly_hex); CPPUNIT_TEST_SUITE_END(); @@ -380,6 +381,91 @@ public: elem->node_id(side_nodes[n])); } } + + // A mesh with a C0POLYHEDRON and a HEX8 in the *same* subdomain. Exodus + // requires a single element type per block, so the writer has to split + // this subdomain into two blocks (one NFACED, one HEX8). Before that + // fix the whole subdomain was written as a single block of the first + // element's type, which errored (or, when the first element was the + // HEX8, crashed) instead of writing the polyhedra. + void test_write_and_read_mixed_poly_hex() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + + // A hexagonal-prism polyhedron, nodes 0..11, in subdomain 1. + const std::vector points = + { { 0, -2, 0}, {-1, -1, 0}, {-1, 1, 0}, { 0, 2, 0}, { 1, 1, 0}, { 1, -1, 0}, + { 0, -2, 1}, {-1, -1, 1}, {-1, 1, 1}, { 0, 2, 1}, { 1, 1, 1}, { 1, -1, 1} }; + for (auto p : index_range(points)) + mesh.add_point(points[p], p); + + const std::vector> nodes_on_side = + { {0, 1, 2, 3, 4, 5}, {0, 1, 7, 6}, {1, 2, 8, 7}, {2, 3, 9, 8}, + {3, 4, 10, 9}, {4, 5, 11, 10}, {5, 0, 6, 11}, {6, 7, 8, 9, 10, 11} }; + std::vector> sides(nodes_on_side.size()); + for (auto s : index_range(nodes_on_side)) + { + sides[s] = std::make_shared(nodes_on_side[s].size()); + for (auto i : index_range(nodes_on_side[s])) + sides[s]->set_node(i, mesh.node_ptr(nodes_on_side[s][i])); + } + std::unique_ptr mid_elem_node; + std::unique_ptr polyhedron = + std::make_unique(sides, mid_elem_node); + if (mid_elem_node) + mesh.add_node(std::move(mid_elem_node)); + polyhedron->set_id() = 0; + polyhedron->subdomain_id() = 1; + mesh.add_elem(std::move(polyhedron)); + + // A HEX8 in the *same* subdomain, nodes 12..19 (offset in x). + const std::vector hex_points = + { {5,0,0}, {6,0,0}, {6,1,0}, {5,1,0}, {5,0,1}, {6,0,1}, {6,1,1}, {5,1,1} }; + for (auto i : index_range(hex_points)) + mesh.add_point(hex_points[i], 12 + i); + std::unique_ptr hex = Elem::build(HEX8); + for (unsigned int i = 0; i != 8; ++i) + hex->set_node(i, mesh.node_ptr(12 + i)); + hex->set_id() = 1; + hex->subdomain_id() = 1; + mesh.add_elem(std::move(hex)); + + mesh.cache_elem_data(); + mesh.prepare_for_use(); + + { + ExodusII_IO exii(mesh); + exii.write("write_exodus_mixed_poly_hex.e"); + } + + TestCommWorld->barrier(); + + Mesh input_mesh(*TestCommWorld); + ExodusII_IO exii_input(input_mesh); + if (input_mesh.processor_id() == 0) + exii_input.read("write_exodus_mixed_poly_hex.e"); + + MeshCommunication().broadcast(input_mesh); + input_mesh.prepare_for_use(); + + CPPUNIT_ASSERT_EQUAL(cast_int(2), input_mesh.n_elem()); + + // Both element types must survive the round trip, in their own blocks. + unsigned int n_poly = 0, n_hex = 0; + for (const auto & elem : input_mesh.active_local_element_ptr_range()) + { + if (elem->type() == C0POLYHEDRON) + ++n_poly; + else if (elem->type() == HEX8) + ++n_hex; + } + input_mesh.comm().sum(n_poly); + input_mesh.comm().sum(n_hex); + CPPUNIT_ASSERT_EQUAL(1u, n_poly); + CPPUNIT_ASSERT_EQUAL(1u, n_hex); + } }; CPPUNIT_TEST_SUITE_REGISTRATION(ExodusC0PolyhedronTest);