Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow specifying attribute name when reading meshtag #3257

Open
wants to merge 34 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 27 commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
9c22fac
Allow to specify attribute name when reading meshtag
mleoni-pf Jun 6, 2024
4efecce
Add python interface
mleoni-pf Jun 6, 2024
361a202
Disambiguate functions
mleoni-pf Jun 7, 2024
983af85
Fix Python wrapper
mleoni-pf Jun 7, 2024
9c31b90
Applying review changes
mleoni-pf Jun 7, 2024
6950dd0
Linting
mleoni-pf Jun 7, 2024
b7fa6c4
Only one function in the python wrapper and interface
mleoni-pf Jun 7, 2024
cf379a4
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Jun 7, 2024
ca55e25
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Jun 20, 2024
fc18895
Applied code review changes
mleoni-pf Jun 20, 2024
6075ef4
Added unit test
mleoni-pf Jun 20, 2024
fa3a840
Fix comments
mleoni-pf Jun 20, 2024
2501a9e
Removed unused variables
mleoni-pf Jun 20, 2024
dbf155a
Merge branch 'main' into mleoni/readNamedMeshtags
garth-wells Jun 26, 2024
26aed10
Merge branch 'main' into mleoni/readNamedMeshtags
garth-wells Jun 26, 2024
7719809
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Jul 1, 2024
b4a5187
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Jul 8, 2024
47e9fdc
Merge branch 'main' into mleoni/readNamedMeshtags
garth-wells Jul 17, 2024
60829f0
Using ASCII file to circumvent CI issue
mleoni-pf Jul 17, 2024
af42617
Merge branch 'mleoni/readNamedMeshtags' of github.com:mleoni-pf/dolfi…
mleoni-pf Jul 17, 2024
9df7574
Specify ASCII encoding
mleoni-pf Jul 17, 2024
90d5b48
Fix a bug and create mesh within test
mleoni-pf Jul 17, 2024
b4ed95a
Dodge issue 3316
mleoni-pf Jul 17, 2024
3ed4bf7
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Jul 31, 2024
edd09a1
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Aug 26, 2024
5b13c3b
Merge branch 'FEniCS:main' into mleoni/readNamedMeshtags
mleoni-pf Aug 29, 2024
c4f87c6
Merge remote-tracking branch 'origin/main' into mleoni/readNamedMeshtags
garth-wells Sep 12, 2024
80fe197
Merge branch 'main' into mleoni/readNamedMeshtags
jhale Sep 22, 2024
ae2d766
Merge branch 'main' into mleoni/readNamedMeshtags
garth-wells Sep 24, 2024
4457147
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Sep 26, 2024
3bbeb7e
Applied review suggestions
mleoni-pf Sep 26, 2024
7ee33ec
Fix missing renaming
mleoni-pf Sep 27, 2024
48e0835
Fix python wrapper
mleoni-pf Sep 27, 2024
9c4c8cd
Merge branch 'main' into mleoni/readNamedMeshtags
mleoni-pf Sep 27, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 32 additions & 4 deletions cpp/dolfinx/io/XDMFFile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -322,8 +322,9 @@ template void XDMFFile::write_meshtags(const mesh::MeshTags<std::int32_t>&,
/// @endcond
//-----------------------------------------------------------------------------
mesh::MeshTags<std::int32_t>
XDMFFile::read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
std::string xpath)
XDMFFile::read_meshtags_by_label(const mesh::Mesh<double>& 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();
Expand All @@ -336,8 +337,28 @@ XDMFFile::read_meshtags(const mesh::Mesh<double>& 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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this just iterating through the attributes? If yes, can it be made simpler using pugixml functionality, see https://pugixml.org/docs/manual.html#access.basic.

The code is hard to follow - at a minimum it needs a comment on what it is doing.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I remember looking at this. The problem was, if I recall correctly, that I am searching for an attribute by the name of "Name" and the "basic functionality" functions do not allow to search for a name. If I wanted to loop using next_attribute I would then need to check if the attribute's name is "Name" and, if not, continue.

I will add documentation but the code doesn't look to exotic to me [because I wrote it, of course XD].

{
pugi::xml_attribute hint;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this for?

Copy link
Contributor Author

@mleoni-pf mleoni-pf Jun 20, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's something that pugixml requires, like a hint to try to speed-up the search of a node.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The hint is not set or re-used - does the code compile without it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it can because hint is also an output parameter of the function.

pugi::xml_attribute name = attribute_node.attribute("Name", hint);
if (std::string(name.value()) == attribute_label)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Creating the string is a bit clumsy, maybe use https://en.cppreference.com/w/cpp/string/char_traits/compare.

Copy link
Contributor Author

@mleoni-pf mleoni-pf Jun 20, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using the function you suggested requires providing a number of characters to compare, so it would have to be something like

if (auto len = std::max(std::strlen(name.value()), attribute_label.size());
    std::char_traits<char>::compare(name.value(), attribute_label.c_str(),
                                    len) == 0)

right?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think name.value() == attribute_label might work, see https://en.cppreference.com/w/cpp/string/basic_string/operator_cmp.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done!

{
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<std::int32_t>(
_comm.comm(), values_data_node, _h5_id);

Expand Down Expand Up @@ -376,6 +397,13 @@ XDMFFile::read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
return meshtags;
}
//-----------------------------------------------------------------------------
mesh::MeshTags<std::int32_t>
XDMFFile::read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
std::string xpath)
{
return read_meshtags_by_label(mesh, name, std::string(), xpath);
}
//-----------------------------------------------------------------------------
std::pair<mesh::CellType, int> XDMFFile::read_cell_type(std::string grid_name,
std::string xpath)
{
Expand Down
14 changes: 13 additions & 1 deletion cpp/dolfinx/io/XDMFFile.h
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,21 @@ class XDMFFile
const mesh::Geometry<T>& x, std::string geometry_xpath,
std::string xpath = "/Xdmf/Domain");

/// Read MeshTags by name
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a bit confusing - the docstring says "Read MeshTags by name" but the function name is "read_meshtags_by_label". Is it by name or by label?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The issue comes from the following fact: mesh tags are identified by "a word", so it is natural that this word be called the "name" of the mesh tag. The problem is that there is already another name parameter which refers to the name of the mesh node in the same file. I picked the word "label" to distinguish the two internally, in the code implementation, but a user should not be concerned with this implementation detail: a user saved a mesh tag giving it a "name" and wants to read it back by the same "name".
The fact remains that the user still needs to call a function whose name includes the word "label", not "name", which is indeed a bit confusing. I will change everything to be called more consistently, the code shouldn't become too unintelligible on the development side.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the way, re-reading this conversation I noticed that in the first review, four months ago, you asked me to change the method name to avoid the word "name" in favour of "label". I had forgotten about that and, based on your comment, I reverted the change. Which one should I keep? To a user, a mesh tag's "label" means presumably very little...

/// @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<std::int32_t>
read_meshtags_by_label(const mesh::Mesh<double>& 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<std::int32_t>
read_meshtags(const mesh::Mesh<double>& mesh, std::string name,
Expand Down
1 change: 1 addition & 0 deletions cpp/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
refinement/interval.cpp
Expand Down
91 changes: 91 additions & 0 deletions cpp/test/mesh/read_named_meshtags.cpp
Original file line number Diff line number Diff line change
@@ -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 <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);
std::ranges::fill(domain_values, domain_value);

std::vector<std::int32_t> material_values(n_cells);
std::ranges::fill(material_values, 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"));

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());
}
11 changes: 9 additions & 2 deletions python/dolfinx/io/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down
21 changes: 13 additions & 8 deletions python/dolfinx/wrappers/io.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<T>& self) { self.close(); })
.def(
"write", [](dolfinx::io::VTXWriter<T>& self, double t)
{ self.write(t); }, nb::arg("t"));
"write",
[](dolfinx::io::VTXWriter<T>& self, double t) { self.write(t); },
nb::arg("t"));
}

{
Expand Down Expand Up @@ -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<T>& self) { self.close(); })
.def(
"write", [](dolfinx::io::FidesWriter<T>& self, double t)
{ self.write(t); }, nb::arg("t"));
"write",
[](dolfinx::io::FidesWriter<T>& self, double t) { self.write(t); },
nb::arg("t"));
}
#endif
}
Expand Down Expand Up @@ -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<float>(xdmf_file);
xdmf_real_fn<double>(xdmf_file);
Expand Down
Loading