From 13bc37f8db1ed162aac185ec8679e31c7f2277af Mon Sep 17 00:00:00 2001 From: mleoni-pf <160495781+mleoni-pf@users.noreply.github.com> Date: Fri, 22 Nov 2024 14:39:51 +0100 Subject: [PATCH] Allow specifying attribute name when reading meshtag (#3257) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Safe default value for unbound variables in dolfinx.conf * Allow to specify attribute name when reading meshtag * Add python interface * Disambiguate functions * Fix Python wrapper * Applying review changes * Linting * Only one function in the python wrapper and interface * Applied code review changes * Added unit test * Fix comments * Removed unused variables * Using ASCII file to circumvent CI issue * Specify ASCII encoding * Fix a bug and create mesh within test * Dodge issue 3316 * Applied review suggestions * Fix missing renaming * Fix python wrapper * Simplify code * Update cpp/dolfinx/io/XDMFFile.h Co-authored-by: Jørgen Schartum Dokken * Update cpp/dolfinx/io/XDMFFile.h Co-authored-by: Jørgen Schartum Dokken * Sprinkling some const * Apply revision * Fix doc * Fix doc * Ruff format * Const and ref fixes * Ruff fixes * Update python/dolfinx/io/utils.py Co-authored-by: Jørgen Schartum Dokken * Add type hints * Fix dispatch * Fixes * Add missing include * Tidy cmake order * Minor fixes * Added test * Ruff fixes * Ruff fixes --------- Co-authored-by: Jørgen Schartum Dokken Co-authored-by: Garth N. Wells --- cpp/dolfinx/io/XDMFFile.cpp | 15 ++++ cpp/dolfinx/io/XDMFFile.h | 25 ++++--- cpp/test/CMakeLists.txt | 3 +- cpp/test/mesh/read_named_meshtags.cpp | 86 +++++++++++++++++++++++ python/dolfinx/io/utils.py | 37 ++++++++-- python/dolfinx/wrappers/io.cpp | 4 +- python/test/unit/io/test_xdmf_meshtags.py | 39 ++++++++++ 7 files changed, 189 insertions(+), 20 deletions(-) create mode 100644 cpp/test/mesh/read_named_meshtags.cpp diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp index 25226a01669..aee81c258a0 100644 --- a/cpp/dolfinx/io/XDMFFile.cpp +++ b/cpp/dolfinx/io/XDMFFile.cpp @@ -335,6 +335,7 @@ template void XDMFFile::write_meshtags(const mesh::MeshTags&, //----------------------------------------------------------------------------- mesh::MeshTags XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, + std::optional attribute_name, std::string xpath) { spdlog::info("XDMF read meshtags ({})", name); @@ -350,6 +351,20 @@ XDMFFile::read_meshtags(const mesh::Mesh& mesh, std::string name, pugi::xml_node values_data_node = grid_node.child("Attribute").child("DataItem"); + if (attribute_name) + { + // Search for a child that contains an attribute with the requested name + pugi::xml_node attribute_node = grid_node.find_child( + [attr_name = *attribute_name](auto n) + { return n.attribute("Name").value() == attr_name; }); + if (!attribute_node) + { + throw std::runtime_error("Attribute with name '" + *attribute_name + + "' not found."); + } + else + values_data_node = attribute_node.child("DataItem"); + } const std::vector values = xdmf_utils::get_dataset( _comm.comm(), values_data_node, _h5_id); diff --git a/cpp/dolfinx/io/XDMFFile.h b/cpp/dolfinx/io/XDMFFile.h index 35104bc6a7c..18d510867bd 100644 --- a/cpp/dolfinx/io/XDMFFile.h +++ b/cpp/dolfinx/io/XDMFFile.h @@ -12,7 +12,9 @@ #include #include #include +#include #include +#include #include namespace pugi @@ -43,7 +45,7 @@ class MeshTags; namespace dolfinx::io { -/// Read and write mesh::Mesh, fem::Function and other objects in +/// @brief Read and write mesh::Mesh, fem::Function and other objects in /// XDMF. /// /// This class supports the output of meshes and functions in XDMF @@ -51,8 +53,8 @@ namespace dolfinx::io /// the data and points to a HDF5 file that stores the actual problem /// data. Output of data in parallel is supported. /// -/// XDMF is not suitable for higher order geometries, as their currently -/// only supports 1st and 2nd order geometries. +/// XDMF is not suitable for higher-order geometries; it only supports +/// 1st- and 2nd-order geometries. class XDMFFile { public: @@ -63,12 +65,9 @@ class XDMFFile ASCII }; - /// Default encoding type - static const Encoding default_encoding = Encoding::HDF5; - /// Constructor XDMFFile(MPI_Comm comm, const std::filesystem::path& filename, - std::string file_mode, Encoding encoding = default_encoding); + std::string file_mode, Encoding encoding = Encoding::HDF5); /// Move constructor XDMFFile(XDMFFile&&) = default; @@ -97,10 +96,10 @@ class XDMFFile void write_geometry(const mesh::Geometry& geometry, std::string name, std::string xpath = "/Xdmf/Domain"); - /// Read in Mesh + /// Read Mesh /// @param[in] element Element that describes the geometry of a cell /// @param[in] mode The type of ghosting/halo to use for the mesh when - /// distributed in parallel + /// distributed in parallel /// @param[in] name /// @param[in] xpath XPath where Mesh Grid is located /// @return A Mesh distributed on the same communicator as the @@ -165,11 +164,15 @@ class XDMFFile std::string xpath = "/Xdmf/Domain"); /// Read MeshTags - /// @param[in] mesh The Mesh that the data is defined on - /// @param[in] name + /// @param[in] mesh Mesh that the input data is defined on + /// @param[in] name Name of the grid node in the xml-scheme of the + /// XDMF-file. E.g. "Material" in Grid Name="Material" + /// GridType="Uniform" + /// @param[in] attribute_name Name of the attribute to read /// @param[in] xpath XPath where MeshTags Grid is stored in file mesh::MeshTags read_meshtags(const mesh::Mesh& mesh, std::string name, + std::optional attribute_name, std::string xpath = "/Xdmf/Domain"); /// Write Information diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index 4a3dd8f3817..ada545792a1 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -37,13 +37,14 @@ add_executable( vector.cpp matrix.cpp io.cpp + common/CIFailure.cpp common/sub_systems_manager.cpp common/index_map.cpp common/sort.cpp fem/functionspace.cpp mesh/distributed_mesh.cpp mesh/generation.cpp - common/CIFailure.cpp + mesh/read_named_meshtags.cpp mesh/refinement/interval.cpp mesh/refinement/option.cpp mesh/refinement/rectangle.cpp diff --git a/cpp/test/mesh/read_named_meshtags.cpp b/cpp/test/mesh/read_named_meshtags.cpp new file mode 100644 index 00000000000..cf3830aa1ae --- /dev/null +++ b/cpp/test/mesh/read_named_meshtags.cpp @@ -0,0 +1,86 @@ +// 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_name = "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, domain_value); + std::vector material_values(n_cells, 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_name, "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 mesh_file(MPI_COMM_WORLD, mesh_file_name, "r", + io::XDMFFile::Encoding::HDF5); + mesh = std::make_shared>(mesh_file.read_mesh( + fem::CoordinateElement(mesh::CellType::triangle, 1), + mesh::GhostMode::none, "mesh")); + + mesh::MeshTags mt_first + = mesh_file.read_meshtags(*mesh, "material", {}); + CHECK(mt_first.values().front() == material_value); + + mesh::MeshTags mt_domain + = mesh_file.read_meshtags(*mesh, "domain", "domain", "/Xdmf/Domain"); + CHECK(mt_domain.values().front() == domain_value); + + mesh::MeshTags mt_material + = mesh_file.read_meshtags(*mesh, "material", "material", "/Xdmf/Domain"); + CHECK(mt_material.values().front() == material_value); + + CHECK_THROWS(mesh_file.read_meshtags(*mesh, "mesh", "missing")); + mesh_file.close(); +} +} // namespace + +TEST_CASE("Read meshtag by name", "[read_meshtag_by_name]") +{ + test_read_named_meshtags(); +} diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index 2f5a2b90f26..f08eb362499 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -122,7 +122,6 @@ class VTKFile(_cpp.io.VTKFile): VTK supports arbitrary order Lagrange finite elements for the geometry description. XDMF is the preferred format for geometry order <= 2. - """ def __enter__(self): @@ -205,8 +204,32 @@ 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: Mesh, + name: str, + attribute_name: typing.Optional[str] = None, + xpath: str = "/Xdmf/Domain", + ) -> MeshTags: + """Read MeshTags with a specific name as specified in the XMDF file. + + Args: + mesh: Mesh that the input data is defined on. + name: Name of the grid node in the xml-scheme of the + XDMF-file. + attribute_name: The name of the attribute to read. If + ``attribute_name`` is empty, reads the first attribute in + the file. If ``attribute_name`` is not empty but no + attributes have the provided name, throws an error. If + multiple attributes have the provided name, reads the + first one found. + xpath: XPath where MeshTags Grid is stored in file. + + Returns: + A MeshTags object containing the requested data read from + file. + """ + mt = super().read_meshtags(mesh._cpp_object, name, attribute_name, xpath) return MeshTags(mt) @@ -215,12 +238,12 @@ def distribute_entity_data( ) -> tuple[npt.NDArray[np.int64], np.ndarray]: """Given a set of mesh entities and values, distribute them to the process that owns the entity. - The entities are described by the global vertex indices of the mesh. These entity indices are - using the original input ordering. + The entities are described by the global vertex indices of the mesh. + These entity indices are using the original input ordering. Returns: - Entities owned by the process (and their local entity-to-vertex indices) and the - corresponding values. + Entities owned by the process (and their local entity-to-vertex + indices) and the corresponding values. """ return _cpp.io.distribute_entity_data( mesh.topology._cpp_object, diff --git a/python/dolfinx/wrappers/io.cpp b/python/dolfinx/wrappers/io.cpp index 2ce136ea732..35ab47e2378 100644 --- a/python/dolfinx/wrappers/io.cpp +++ b/python/dolfinx/wrappers/io.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -271,7 +272,8 @@ void io(nb::module_& m) .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") + nb::arg("mesh"), nb::arg("name"), nb::arg("attribute_name").none(), + 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, diff --git a/python/test/unit/io/test_xdmf_meshtags.py b/python/test/unit/io/test_xdmf_meshtags.py index d8d18f48da5..37f8e4411fa 100644 --- a/python/test/unit/io/test_xdmf_meshtags.py +++ b/python/test/unit/io/test_xdmf_meshtags.py @@ -99,3 +99,42 @@ def test_3d(tempdir, cell_type, encoding): num_facets = int(tree.findall(".//Grid[@Name='facets']/Topology")[0].get("NumberOfElements")) assert num_lines == lines_local assert num_facets == facets_local + + +@pytest.mark.parametrize("cell_type", celltypes_3D) +@pytest.mark.parametrize("encoding", encodings) +def test_read_named_meshtags(tempdir, cell_type, encoding): + domain_value = 1 + material_value = 2 + + filename = Path(tempdir, "named_meshtags.xdmf") + comm = MPI.COMM_WORLD + mesh = create_unit_cube(comm, 4, 4, 4, cell_type) + + indices = np.arange(mesh.topology.index_map(3).size_local) + domain_values = np.full(indices.shape, domain_value, dtype=np.int32) + mt_domains = meshtags(mesh, 3, indices, domain_values) + mt_domains.name = "domain" + + material_values = np.full(indices.shape, material_value, dtype=np.int32) + mt_materials = meshtags(mesh, 3, indices, material_values) + mt_materials.name = "material" + + with XDMFFile(comm, filename, "w", encoding=encoding) as file: + file.write_mesh(mesh) + file.write_meshtags(mt_domains, mesh.geometry) + file.write_meshtags(mt_materials, mesh.geometry) + + with XDMFFile(comm, filename, "r", encoding=encoding) as file: + mesh_in = file.read_mesh() + tdim = mesh_in.topology.dim + mesh_in.topology.create_connectivity(1, tdim) + + mt_first_in = file.read_meshtags(mesh_in, "material") + assert all(v == material_value for v in mt_first_in.values) + + mt_domains_in = file.read_meshtags(mesh_in, "domain") + assert all(v == domain_value for v in mt_domains_in.values) + + mt_materials_in = file.read_meshtags(mesh_in, "material") + assert all(v == material_value for v in mt_materials_in.values)