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 29 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 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
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 @@ -334,8 +334,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 @@ -348,8 +349,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.

{
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?

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.

{
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 @@ -388,6 +409,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?

/// @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
mesh/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