Allow specifying attribute name when reading meshtag

cpp/dolfinx/io/XDMFFile.cpp

@@ -335,6 +335,7 @@ template void XDMFFile::write_meshtags(const mesh::MeshTags<std::int32_t>&,
 //-----------------------------------------------------------------------------
 mesh::MeshTags<std::int32_t>
 XDMFFile::read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
+                        std::optional<std::string> attribute_name,
                         std::string xpath)
 {
   spdlog::info("XDMF read meshtags ({})", name);
@@ -350,6 +351,20 @@ XDMFFile::read_meshtags(const mesh::Mesh<double>& 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<std::int32_t>(
       _comm.comm(), values_data_node, _h5_id);
 

cpp/dolfinx/io/XDMFFile.h

@@ -12,7 +12,9 @@
 #include <dolfinx/mesh/cell_types.h>
 #include <filesystem>
 #include <memory>
+#include <optional>
 #include <string>
+#include <utility>
 #include <variant>
 
 namespace pugi
@@ -43,16 +45,16 @@ 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
 /// (http://www.xdmf.org) format. It creates an XML file that describes
 /// 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<double>& 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<std::int32_t>
   read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
+                std::optional<std::string> attribute_name,
                 std::string xpath = "/Xdmf/Domain");
 
   /// Write Information

cpp/test/CMakeLists.txt

@@ -37,13 +37,15 @@ add_executable(
   vector.cpp
   matrix.cpp
   io.cpp
+  common/CIFailure.cpp
   common/sub_systems_manager.cpp
   common/index_map.cpp
   common/sort.cpp
+  common/CIFailure.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

cpp/test/mesh/read_named_meshtags.cpp

@@ -0,0 +1,85 @@
+// 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 <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/io/XDMFFile.h>
+#include <dolfinx/la/Vector.h>
+#include <dolfinx/mesh/MeshTags.h>
+#include <dolfinx/mesh/generation.h>
+#include <dolfinx/mesh/utils.h>
+#include <numeric>
+#include <string>
+
+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::Mesh<double>>(
+      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<std::int32_t> indices(n_cells);
+  std::iota(std::begin(indices), std::end(indices), 0);
+
+  std::vector<std::int32_t> domain_values(n_cells, domain_value);
+  std::vector<std::int32_t> material_values(n_cells, material_value);
+
+  mesh::MeshTags<std::int32_t> mt_domains(mesh->topology(), 2, indices,
+                                          domain_values);
+  mt_domains.name = "domain";
+
+  mesh::MeshTags<std::int32_t> 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<mesh::Mesh<double>>(meshFile.read_mesh(
+      fem::CoordinateElement<double>(mesh::CellType::triangle, 1),
+      mesh::GhostMode::none, "mesh"));
+
+  mesh::MeshTags<std::int32_t> mt_first
+      = meshFile.read_meshtags(*mesh, "material", {});
+  CHECK(mt_first.values().front() == material_value);
+
+  mesh::MeshTags<std::int32_t> mt_domain
+      = meshFile.read_meshtags(*mesh, "domain", "domain", "/Xdmf/Domain");
+  CHECK(mt_domain.values().front() == domain_value);
+
+  mesh::MeshTags<std::int32_t> mt_material
+      = meshFile.read_meshtags(*mesh, "material", "material", "/Xdmf/Domain");
+  CHECK(mt_material.values().front() == material_value);
+
+  CHECK_THROWS(meshFile.read_meshtags(*mesh, "mesh", "missing"));
+}
+} // namespace
+
+TEST_CASE("Read meshtag by name", "[read_meshtag_by_name]")
+{
+  test_read_named_meshtags();
+}

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,

python/dolfinx/wrappers/io.cpp

@@ -24,6 +24,7 @@
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/complex.h>
 #include <nanobind/stl/filesystem.h>
+#include <nanobind/stl/optional.h>
 #include <nanobind/stl/pair.h>
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/string.h>
@@ -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,