diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp index 25226a0166..08ece3eae9 100644 --- a/cpp/dolfinx/io/XDMFFile.cpp +++ b/cpp/dolfinx/io/XDMFFile.cpp @@ -334,8 +334,9 @@ template void XDMFFile::write_meshtags(const mesh::MeshTags&, /// @endcond //----------------------------------------------------------------------------- mesh::MeshTags -XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, - std::string xpath) +XDMFFile::read_meshtags_by_label(const mesh::Mesh& mesh, + std::string name, std::string attribute_label, + std::string xpath) { spdlog::info("XDMF read meshtags ({})", name); pugi::xml_node node = _xml_doc->select_node(xpath.c_str()).node(); @@ -348,8 +349,28 @@ XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, const auto [entities, eshape] = read_topology_data(name, xpath); - pugi::xml_node values_data_node - = grid_node.child("Attribute").child("DataItem"); + pugi::xml_node attribute_node = grid_node.child("Attribute"); + pugi::xml_node values_data_node = attribute_node.child("DataItem"); + if (!attribute_label.empty()) + { + bool found = false; + while (!found and attribute_node) + { + pugi::xml_attribute hint; + pugi::xml_attribute name = attribute_node.attribute("Name", hint); + if (std::string(name.value()) == attribute_label) + { + values_data_node = attribute_node.child("DataItem"); + found = true; + } + attribute_node = attribute_node.next_sibling("Attribute"); + } + if (!found) + { + throw std::runtime_error("Attribute with name '" + attribute_label + + "' not found."); + } + } const std::vector values = xdmf_utils::get_dataset( _comm.comm(), values_data_node, _h5_id); @@ -388,6 +409,13 @@ XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, return meshtags; } //----------------------------------------------------------------------------- +mesh::MeshTags +XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, + std::string xpath) +{ + return read_meshtags_by_label(mesh, name, std::string(), xpath); +} +//----------------------------------------------------------------------------- std::pair XDMFFile::read_cell_type(std::string grid_name, std::string xpath) { diff --git a/cpp/dolfinx/io/XDMFFile.h b/cpp/dolfinx/io/XDMFFile.h index 35104bc6a7..140258f0f6 100644 --- a/cpp/dolfinx/io/XDMFFile.h +++ b/cpp/dolfinx/io/XDMFFile.h @@ -164,9 +164,21 @@ class XDMFFile const mesh::Geometry& x, std::string geometry_xpath, std::string xpath = "/Xdmf/Domain"); + /// Read MeshTags by name + /// @param[in] mesh The Mesh that the data is defined on + /// @param[in] name Name of the grid node in the xml file. E.g. "Material" in + /// Grid Name="Material" GridType="Uniform" + /// @param[in] attribute_label The name of the attribute to read + /// @param[in] xpath XPath where MeshTags Grid is stored in file + mesh::MeshTags + read_meshtags_by_label(const mesh::Mesh& mesh, std::string name, + std::string attribute_label, + std::string xpath = "/Xdmf/Domain"); + /// Read MeshTags /// @param[in] mesh The Mesh that the data is defined on - /// @param[in] name + /// @param[in] name Name of the grid node in the xml file. E.g. "Material" in + /// Grid Name="Material" GridType="Uniform" /// @param[in] xpath XPath where MeshTags Grid is stored in file mesh::MeshTags read_meshtags(const mesh::Mesh& mesh, std::string name, diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index 4a3dd8f381..901bd29247 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -42,6 +42,7 @@ add_executable( common/sort.cpp fem/functionspace.cpp mesh/distributed_mesh.cpp + mesh/read_named_meshtags.cpp mesh/generation.cpp common/CIFailure.cpp mesh/refinement/interval.cpp diff --git a/cpp/test/mesh/read_named_meshtags.cpp b/cpp/test/mesh/read_named_meshtags.cpp new file mode 100644 index 0000000000..805ed29e6d --- /dev/null +++ b/cpp/test/mesh/read_named_meshtags.cpp @@ -0,0 +1,91 @@ +// Copyright (C) 2024 Massimiliano Leoni +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later +// +// Unit tests for reading meshtags by name + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dolfinx; + +namespace +{ + +void test_read_named_meshtags() +{ + const std::string mesh_file = "Domain.xdmf"; + + constexpr std::int32_t domain_value = 1; + constexpr std::int32_t material_value = 2; + + // Create mesh + auto part = mesh::create_cell_partitioner(mesh::GhostMode::none); + auto mesh = std::make_shared>( + mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {3, 3}, + mesh::CellType::triangle, part)); + + const std::int32_t n_cells = mesh->topology()->index_map(2)->size_local(); + std::vector indices(n_cells); + std::iota(std::begin(indices), std::end(indices), 0); + + std::vector domain_values(n_cells); + std::ranges::fill(domain_values, domain_value); + + std::vector material_values(n_cells); + std::ranges::fill(material_values, material_value); + + mesh::MeshTags mt_domains(mesh->topology(), 2, indices, + domain_values); + mt_domains.name = "domain"; + + mesh::MeshTags mt_materials(mesh->topology(), 2, indices, + material_values); + mt_materials.name = "material"; + + io::XDMFFile file(mesh->comm(), mesh_file, "w", io::XDMFFile::Encoding::HDF5); + file.write_mesh(*mesh); + file.write_meshtags(mt_domains, mesh->geometry(), + "/Xdmf/Domain/mesh/Geometry"); + file.write_meshtags(mt_materials, mesh->geometry(), + "/Xdmf/Domain/Grid/Geometry"); + + file.close(); + + io::XDMFFile meshFile(MPI_COMM_WORLD, mesh_file, "r", + io::XDMFFile::Encoding::HDF5); + mesh = std::make_shared>(meshFile.read_mesh( + fem::CoordinateElement(mesh::CellType::triangle, 1), + mesh::GhostMode::none, "mesh")); + + const auto mt_domain = meshFile.read_meshtags_by_label( + *mesh, "domain", "domain", "/Xdmf/Domain"); + + CHECK(mt_domain.values().front() == domain_value); + + const auto mt_material = meshFile.read_meshtags_by_label( + *mesh, "material", "material", "/Xdmf/Domain"); + + CHECK(mt_material.values().front() == material_value); + + CHECK_THROWS(meshFile.read_meshtags_by_label(*mesh, "mesh", "missing")); +} + +} // namespace + +TEST_CASE("Read meshtag by name", "[read_meshtag_by_name]") +{ + CHECK_NOTHROW(test_read_named_meshtags()); +} diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index 005e28385a..95d2035316 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -276,8 +276,15 @@ def read_mesh( ) return Mesh(msh, domain) - def read_meshtags(self, mesh, name, xpath="/Xdmf/Domain"): - mt = super().read_meshtags(mesh._cpp_object, name, xpath) + def read_meshtags(self, mesh, name, attribute_label=None, xpath="/Xdmf/Domain"): + """Read meshtags with a specific name as specified in the xdmf file. If + `attribute_label` is `None`, read the first attribute in the file. + If `attribute_label` is not `None` but no attributes have the provided name, + throw an error. If multiple attributes have the provided name, read the first + one found.""" + mt = super().read_meshtags( + mesh._cpp_object, name, attribute_label if attribute_label is not None else "", xpath + ) return MeshTags(mt) diff --git a/python/dolfinx/wrappers/io.cpp b/python/dolfinx/wrappers/io.cpp index ed7ae36816..6956fc6bdd 100644 --- a/python/dolfinx/wrappers/io.cpp +++ b/python/dolfinx/wrappers/io.cpp @@ -142,8 +142,9 @@ void declare_vtx_writer(nb::module_& m, std::string type) nb::arg("policy") = dolfinx::io::VTXMeshPolicy::update) .def("close", [](dolfinx::io::VTXWriter& self) { self.close(); }) .def( - "write", [](dolfinx::io::VTXWriter& self, double t) - { self.write(t); }, nb::arg("t")); + "write", + [](dolfinx::io::VTXWriter& self, double t) { self.write(t); }, + nb::arg("t")); } { @@ -182,8 +183,9 @@ void declare_vtx_writer(nb::module_& m, std::string type) nb::arg("policy") = dolfinx::io::FidesMeshPolicy::update) .def("close", [](dolfinx::io::FidesWriter& self) { self.close(); }) .def( - "write", [](dolfinx::io::FidesWriter& self, double t) - { self.write(t); }, nb::arg("t")); + "write", + [](dolfinx::io::FidesWriter& self, double t) { self.write(t); }, + nb::arg("t")); } #endif } @@ -308,15 +310,18 @@ void io(nb::module_& m) nb::arg("name") = "mesh", nb::arg("xpath") = "/Xdmf/Domain") .def("read_cell_type", &dolfinx::io::XDMFFile::read_cell_type, nb::arg("name") = "mesh", nb::arg("xpath") = "/Xdmf/Domain") - .def("read_meshtags", &dolfinx::io::XDMFFile::read_meshtags, - nb::arg("mesh"), nb::arg("name"), nb::arg("xpath") = "/Xdmf/Domain") + .def("read_meshtags", &dolfinx::io::XDMFFile::read_meshtags_by_label, + nb::arg("mesh"), nb::arg("name"), nb::arg("attribute_label"), + nb::arg("xpath")) .def("write_information", &dolfinx::io::XDMFFile::write_information, nb::arg("name"), nb::arg("value"), nb::arg("xpath") = "/Xdmf/Domain") .def("read_information", &dolfinx::io::XDMFFile::read_information, nb::arg("name"), nb::arg("xpath") = "/Xdmf/Domain") .def_prop_ro( - "comm", [](dolfinx::io::XDMFFile& self) - { return MPICommWrapper(self.comm()); }, nb::keep_alive<0, 1>()); + "comm", + [](dolfinx::io::XDMFFile& self) + { return MPICommWrapper(self.comm()); }, + nb::keep_alive<0, 1>()); xdmf_real_fn(xdmf_file); xdmf_real_fn(xdmf_file);