Allow specifying attribute name when reading meshtag (#3257)

* 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 <[email protected]> * Update cpp/dolfinx/io/XDMFFile.h Co-authored-by: Jørgen Schartum Dokken <[email protected]> * 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 <[email protected]> * 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 <[email protected]> Co-authored-by: Garth N. Wells <[email protected]>
FEniCS · Nov 22, 2024 · 13bc37f · 13bc37f
1 parent c0e48a1
commit 13bc37f
Show file tree

Hide file tree

Showing 7 changed files with 189 additions and 20 deletions.
diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/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);
 

diff --git a/cpp/dolfinx/io/XDMFFile.h b/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

diff --git 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
@@ -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 <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_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::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_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::Mesh<double>>(mesh_file.read_mesh(
+      fem::CoordinateElement<double>(mesh::CellType::triangle, 1),
+      mesh::GhostMode::none, "mesh"));
+
+  mesh::MeshTags<std::int32_t> mt_first
+      = mesh_file.read_meshtags(*mesh, "material", {});
+  CHECK(mt_first.values().front() == material_value);
+
+  mesh::MeshTags<std::int32_t> mt_domain
+      = mesh_file.read_meshtags(*mesh, "domain", "domain", "/Xdmf/Domain");
+  CHECK(mt_domain.values().front() == domain_value);
+
+  mesh::MeshTags<std::int32_t> 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
@@ -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
@@ -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,

diff --git 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)