From dda406ba2d0cbf9c49c9b4026233090e8fdf74be Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 14:01:59 +0000 Subject: [PATCH 1/7] Add transfer tags to submesh --- cpp/dolfinx/fem/Function.h | 14 +- cpp/dolfinx/fem/interpolate.h | 34 ++-- cpp/dolfinx/mesh/EntityMap.cpp | 50 ------ cpp/dolfinx/mesh/EntityMap.h | 55 ++++++- cpp/dolfinx/mesh/Topology.h | 7 + cpp/dolfinx/mesh/utils.h | 149 +++++++++++++++++- python/dolfinx/mesh.py | 51 ++++++ .../dolfinx/wrappers/dolfinx_wrappers/mesh.h | 14 ++ python/dolfinx/wrappers/mesh.cpp | 2 +- 9 files changed, 294 insertions(+), 82 deletions(-) diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h index cfdec86e843..f53e775b727 100644 --- a/cpp/dolfinx/fem/Function.h +++ b/cpp/dolfinx/fem/Function.h @@ -164,7 +164,7 @@ class Function std::pair, std::vector>( md::mdspan>)>& f, - CellRange auto&& cells) + mesh::CellRange auto&& cells) { assert(_function_space); assert(_function_space->element()); @@ -249,7 +249,7 @@ class Function /// /// @pre `cells0` and `cells1` must have the same length. void interpolate(const Function& u0, - CellRange auto&& cells0, CellRange auto&& cells1) + mesh::CellRange auto&& cells0, mesh::CellRange auto&& cells1) { fem::interpolate(*this, cells1, u0, cells0); } @@ -261,7 +261,7 @@ class Function /// @param[in] u Function to be interpolated. /// @param[in] cells Cells to interpolate from. void interpolate(const Function& u, - CellRange auto&& cells) + mesh::CellRange auto&& cells) { fem::interpolate(*this, u, cells); } @@ -297,7 +297,7 @@ class Function /// /// @pre `cells0` `cells1` must have the same length. void interpolate(const Expression& e0, - CellRange auto&& cells0, CellRange auto&& cells1) + mesh::CellRange auto&& cells0, mesh::CellRange auto&& cells1) { // Extract mesh const mesh::Mesh* mesh0 = nullptr; @@ -398,7 +398,7 @@ class Function /// interpolate from if `e0` has Function coefficients. If no mesh can /// be associated with `e0` then the mesh associated with `this` is used. void interpolate(const Expression& e0, - CellRange auto&& cells) + mesh::CellRange auto&& cells) { interpolate(e0, cells, cells); } @@ -432,7 +432,7 @@ class Function /// interpolation points of `this` with cells in `u`. Can be computed /// with `fem::create_interpolation_data`. void interpolate(const Function& u, - CellRange auto&& cells, double tol, int maxit, + mesh::CellRange auto&& cells, double tol, int maxit, const geometry::PointOwnershipData& interpolation_data) { fem::interpolate(*this, u, cells, tol, maxit, interpolation_data); @@ -455,7 +455,7 @@ class Function /// @param[in] maxit Maximum number of Newton iterations in non-affine /// pull-back. If the mesh geometry is affine this argument is ignored. void eval(std::span x, std::array xshape, - CellRange auto&& cells, std::span u, + mesh::CellRange auto&& cells, std::span u, std::array ushape, double tol, int maxit) const { if (cells.empty()) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index 731c475fdac..f87771add6b 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -36,11 +36,7 @@ concept MDSpan = requires(T x, std::size_t idx) { { x.extent(1) } -> std::integral; }; -/// Requirement on range of cell indices. -template -concept CellRange = std::ranges::input_range and std::ranges::sized_range - and std::is_integral_v< - std::remove_const_t>>; + /// @brief Compute the evaluation points in the physical space at which /// an expression should be computed to interpolate it in a finite @@ -55,7 +51,7 @@ concept CellRange = std::ranges::input_range and std::ranges::sized_range template std::vector interpolation_coords(const fem::FiniteElement& element, const mesh::Geometry& geometry, - CellRange auto&& cells) + mesh::CellRange auto&& cells) { // Find CoordinateElement appropriate to element auto cmap_index = [&geometry](mesh::CellType cell_type) @@ -139,7 +135,7 @@ std::vector interpolation_coords(const fem::FiniteElement& element, /// calling `interpolation_coords`. template void interpolate(Function& u, std::span f, - std::array fshape, CellRange auto&& cells); + std::array fshape, mesh::CellRange auto&& cells); namespace impl { @@ -368,8 +364,8 @@ void interpolation_apply(U&& Pi, V&& data, std::span coeffs, int bs) /// elements must share the same basis function map. Neither is checked /// by the function. template -void interpolate_same_map(Function& u1, CellRange auto&& cells1, - const Function& u0, CellRange auto&& cells0) +void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, + const Function& u0, mesh::CellRange auto&& cells0) { auto V0 = u0.function_space(); assert(V0); @@ -469,9 +465,9 @@ void interpolate_same_map(Function& u1, CellRange auto&& cells1, /// @pre Function%s `u1` and `u0` must share the same mesh. This is not /// checked by the function. template -void interpolate_nonmatching_maps(Function& u1, CellRange auto&& cells1, +void interpolate_nonmatching_maps(Function& u1, mesh::CellRange auto&& cells1, const Function& u0, - CellRange auto&& cells0) + mesh::CellRange auto&& cells0) { // Get mesh auto V0 = u0.function_space(); @@ -719,7 +715,7 @@ void interpolate_nonmatching_maps(Function& u1, CellRange auto&& cells1, /// @param [out] coeffs Output Function coefficients. template void point_evaluation(const FiniteElement& element, bool symmetric, - const DofMap& dofmap, CellRange auto&& cells, + const DofMap& dofmap, mesh::CellRange auto&& cells, std::span cell_info, std::span f, std::array fshape, std::span coeffs) @@ -821,7 +817,7 @@ void point_evaluation(const FiniteElement& element, bool symmetric, /// @param [out] coeffs Output Function coefficients. template void identity_mapped_evaluation(const FiniteElement& element, bool symmetric, - const DofMap& dofmap, CellRange auto&& cells, + const DofMap& dofmap, mesh::CellRange auto&& cells, std::span cell_info, std::span f, std::array fshape, @@ -897,7 +893,7 @@ void identity_mapped_evaluation(const FiniteElement& element, bool symmetric, /// @param [out] coeffs Output Function coefficients. template void piola_mapped_evaluation(const FiniteElement& element, bool symmetric, - const DofMap& dofmap, CellRange auto&& cells, + const DofMap& dofmap, mesh::CellRange auto&& cells, std::span cell_info, std::span f, std::array fshape, @@ -1077,7 +1073,7 @@ void piola_mapped_evaluation(const FiniteElement& element, bool symmetric, template geometry::PointOwnershipData create_interpolation_data( const mesh::Geometry& geometry0, const FiniteElement& element0, - const mesh::Mesh& mesh1, CellRange auto&& cells, T padding) + const mesh::Mesh& mesh1, mesh::CellRange auto&& cells, T padding) { // Collect all the points at which values are needed to define the // interpolating function @@ -1097,7 +1093,7 @@ geometry::PointOwnershipData create_interpolation_data( template void interpolate(Function& u, std::span f, - std::array fshape, CellRange auto&& cells) + std::array fshape, mesh::CellRange auto&& cells) { // TODO: Index for mixed-topology, zero for now const int index = 0; @@ -1180,7 +1176,7 @@ void interpolate(Function& u, std::span f, /// fem::create_interpolation_data. template void interpolate(Function& u1, const Function& u0, - CellRange auto&& cells, double tol, int maxit, + mesh::CellRange auto&& cells, double tol, int maxit, const geometry::PointOwnershipData& interpolation_data) { auto mesh1 = u1.function_space()->mesh(); @@ -1257,8 +1253,8 @@ void interpolate(Function& u1, const Function& u0, /// /// @pre `cells0` and `cells1` have the same size. template -void interpolate(Function& u1, CellRange auto&& cells1, - const Function& u0, CellRange auto&& cells0) +void interpolate(Function& u1, mesh::CellRange auto&& cells1, + const Function& u0, mesh::CellRange auto&& cells0) { if (cells0.size() != cells1.size()) throw std::runtime_error("Length of cell lists do not match."); diff --git a/cpp/dolfinx/mesh/EntityMap.cpp b/cpp/dolfinx/mesh/EntityMap.cpp index ac563df08c3..916278a3aff 100644 --- a/cpp/dolfinx/mesh/EntityMap.cpp +++ b/cpp/dolfinx/mesh/EntityMap.cpp @@ -6,7 +6,6 @@ #include "EntityMap.h" #include "Topology.h" -#include #include #include #include @@ -25,55 +24,6 @@ std::shared_ptr EntityMap::sub_topology() const { return _sub_topology; } -//----------------------------------------------------------------------------- -std::vector -EntityMap::sub_topology_to_topology(std::span entities, - bool inverse) const -{ - if (!inverse) - { - // In this case, we want to map from entity indices in - // `_sub_topology` to corresponding entities in `_topology`. Hence, - // for each index in `entities`, we get the corresponding index in - // `_topology` using `_sub_topology_to_topology` - auto mapped - = entities - | std::views::transform([this](std::int32_t i) - { return _sub_topology_to_topology[i]; }); - return std::vector(mapped.begin(), mapped.end()); - } - else - { - // In this case, we are mapping from entity indices in `_topology` - // to entity indices in `_sub_topology`. Hence, we first need to - // construct the "inverse" of `_sub_topology_to_topology` - std::unordered_map topology_to_sub_topology; - topology_to_sub_topology.reserve(_sub_topology_to_topology.size()); - for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i) - { - topology_to_sub_topology.insert( - {_sub_topology_to_topology[i], static_cast(i)}); - } - // For each entity index in `entities` (which are indices in - // `_topology`), get the corresponding entity in `_sub_topology`. - // Since `_sub_topology` consists of a subset of entities in - // `_topology`, there are entities in topology that may not exist in - // `_sub_topology`. If this is the case, mark those entities with - // -1. - - auto mapped - = entities - | std::views::transform( - [&topology_to_sub_topology](std::int32_t i) - { - // Map the entity if it exists. If it doesn't, mark with - // -1. - auto it = topology_to_sub_topology.find(i); - return (it != topology_to_sub_topology.end()) ? it->second : -1; - }); - return std::vector(mapped.begin(), mapped.end()); - } -} //----------------------------------------------------------------------------- } // namespace dolfinx::mesh diff --git a/cpp/dolfinx/mesh/EntityMap.h b/cpp/dolfinx/mesh/EntityMap.h index 07aad6299ba..9711519f2e1 100644 --- a/cpp/dolfinx/mesh/EntityMap.h +++ b/cpp/dolfinx/mesh/EntityMap.h @@ -1,4 +1,4 @@ -// Copyright (C) 2025 Jørgen S. Dokken and Joseph P. Dean +// Copyright (C) 2025-2026 Jørgen S. Dokken and Joseph P. Dean // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -10,6 +10,7 @@ #include #include #include +#include #include namespace dolfinx::mesh @@ -99,9 +100,55 @@ class EntityMap /// `this->sub_topology()`. /// @return A list of mapped entity indices. Entities that do not /// exist in the target topology are marked as -1. - std::vector - sub_topology_to_topology(std::span entities, - bool inverse) const; +std::vector +sub_topology_to_topology(CellRange auto&& entities, + bool inverse) const +{ + if (!inverse) + { + // In this case, we want to map from entity indices in + // `_sub_topology` to corresponding entities in `_topology`. Hence, + // for each index in `entities`, we get the corresponding index in + // `_topology` using `_sub_topology_to_topology` + auto mapped + = entities + | std::views::transform([this](std::int32_t i) + { return _sub_topology_to_topology[i]; }); + return std::vector(mapped.begin(), mapped.end()); + } + else + { + // In this case, we are mapping from entity indices in `_topology` + // to entity indices in `_sub_topology`. Hence, we first need to + // construct the "inverse" of `_sub_topology_to_topology` + std::unordered_map topology_to_sub_topology; + topology_to_sub_topology.reserve(_sub_topology_to_topology.size()); + for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i) + { + topology_to_sub_topology.insert( + {_sub_topology_to_topology[i], static_cast(i)}); + } + + // For each entity index in `entities` (which are indices in + // `_topology`), get the corresponding entity in `_sub_topology`. + // Since `_sub_topology` consists of a subset of entities in + // `_topology`, there are entities in topology that may not exist in + // `_sub_topology`. If this is the case, mark those entities with + // -1. + + auto mapped + = entities + | std::views::transform( + [&topology_to_sub_topology](std::int32_t i) + { + // Map the entity if it exists. If it doesn't, mark with + // -1. + auto it = topology_to_sub_topology.find(i); + return (it != topology_to_sub_topology.end()) ? it->second : -1; + }); + return std::vector(mapped.begin(), mapped.end()); + } +} private: // Dimension of the entities diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h index 11e55d0b785..814919adf5e 100644 --- a/cpp/dolfinx/mesh/Topology.h +++ b/cpp/dolfinx/mesh/Topology.h @@ -18,6 +18,7 @@ #include #include #include +#include namespace dolfinx::common { @@ -26,6 +27,12 @@ class IndexMap; namespace dolfinx::mesh { + /// Requirement on range of cell indices. +template +concept CellRange = std::ranges::input_range and std::ranges::sized_range + and std::is_integral_v< + std::remove_const_t>>; + enum class CellType : std::int8_t; /// @brief Topology stores the topology of a mesh, consisting of mesh diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index b67c40aa756..3381b815ff7 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1,4 +1,4 @@ -// Copyright (C) 2019-2024 Garth N. Wells +// Copyright (C) 2019-2026 Garth N. Wells and Jørgen S. Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -8,6 +8,7 @@ #include "EntityMap.h" #include "Mesh.h" +#include "MeshTags.h" #include "Topology.h" #include "graphbuild.h" #include @@ -21,6 +22,7 @@ #include #include #include +#include #include #include @@ -1446,4 +1448,149 @@ create_submesh(const Mesh& mesh, int dim, std::move(subx_to_x_dofmap)}; } +/// @brief Transfer a meshtags object from a parent to a submesh. +/// +/// @note All entities not tagged in the parent mesh is tagged with the +/// default value in the submesh, which is `std::numeric_limits::max()`. +/// +/// @param[in] tags The meshtags object on the parent mesh. +/// @param[in] submesh_topology The topology of the submesh. +/// @param[in] vertex_map Map from submesh vertex to parent mesh vertex. +/// @param[in] cell_map Map from submesh cell to parent mesh entity of the same +/// dimension as the meshtags. +/// @return A meshtags object on the submesh, and the bi-directional entity map +/// between the tagged entities in the submesh and the parent mesh. +template +std::tuple, EntityMap> transfer_meshtags_to_submesh( + const MeshTags& tags, + std::shared_ptr submesh_topology, + const EntityMap& vertex_map, const EntityMap& cell_map) +{ + + int tag_dim = tags.dim(); + int submesh_tdim = submesh_topology->dim(); + auto topology = tags.topology(); + if (tag_dim > submesh_tdim) + { + throw std::runtime_error("Tag dimension must be less than or equal to " + "submesh dimension"); + } + + std::shared_ptr sub_cell_imap + = submesh_topology->index_map(submesh_tdim); + if (!sub_cell_imap) + { + throw std::runtime_error( + std::format("Entities of dimension {} does not exist in mesh topology.", + submesh_tdim)); + } + + // Create a map from parent entity to submesh cell + std::int32_t submesh_num_cells + = sub_cell_imap->size_local() + sub_cell_imap->num_ghosts(); + auto sub_cells = std::ranges::views::iota(0, submesh_num_cells); + std::vector sub_cell_to_parent_entity + = cell_map.sub_topology_to_topology(sub_cells, false); + + // Create a full lookup for all cells on the parent mesh, as the tag can have + // entities that are not in the submesh + auto e_to_v = topology->connectivity(tag_dim, 0); + if (!e_to_v) + { + throw std::runtime_error("Missing connectivity between " + + std::to_string(tag_dim) + " and 0"); + } + + std::vector parent_entity_to_sub_cell(e_to_v->num_nodes(), -1); + for (std::int32_t i = 0; i < sub_cell_to_parent_entity.size(); ++i) + parent_entity_to_sub_cell[sub_cell_to_parent_entity[i]] = i; + + // Get map from submesh vertex to parent vertex + std::vector sub_to_parent_vertex; + { + auto sub_vertex_map = submesh_topology->index_map(0); + std::int32_t num_sub_vertices + = sub_vertex_map->size_local() + sub_vertex_map->num_ghosts(); + auto sub_vertices = std::ranges::views::iota(0, num_sub_vertices); + + sub_to_parent_vertex + = vertex_map.sub_topology_to_topology(sub_vertices, false); + } + // Access various connectivity maps + auto sub_e_to_v = submesh_topology->connectivity(tag_dim, 0); + auto sub_c_to_e = submesh_topology->connectivity(submesh_tdim, tag_dim); + auto e_to_sub_cell = topology->connectivity(tag_dim, submesh_tdim); + if (!e_to_sub_cell) + { + throw std::runtime_error("Missing connectivity between " + + std::to_string(tag_dim) + " and " + + std::to_string(submesh_tdim)); + } + + auto sub_entity_imap = submesh_topology->index_map(tag_dim); + + // Prepare sub entity to parent map + std::size_t num_sum_entities + = sub_entity_imap->size_local() + sub_entity_imap->num_ghosts(); + std::vector sub_entity_to_parent(num_sum_entities, -1); + std::vector submesh_values(num_sum_entities, + std::numeric_limits::max()); + std::vector submesh_indices(num_sum_entities); + std::iota(submesh_indices.begin(), submesh_indices.end(), 0); + + std::span tagged_entities = tags.indices(); + std::span tagged_values = tags.values(); + // For each entity in the tag, find all cells of the submesh connected to this + // entity + for (std::size_t i = 0; i < tagged_entities.size(); ++i) + { + auto entity = tagged_entities[i]; + auto entity_vertices = e_to_v->links(entity); + bool entity_found = false; + for (auto parent_cell : e_to_sub_cell->links(entity)) + { + if (entity_found) + break; + std::int32_t sub_cell = parent_entity_to_sub_cell[parent_cell]; + if (sub_cell != -1) + { + // For a cell in the sub mesh find all attached entities, and + // define them by their vertices in the sub mesh + for (auto sub_entity : sub_c_to_e->links(sub_cell)) + { + if (entity_found) + break; + auto sub_vertices = sub_e_to_v->links(sub_entity); + bool entity_matches = true; + for (auto sub_vertex : sub_vertices) + { + if (std::find(entity_vertices.begin(), entity_vertices.end(), + sub_to_parent_vertex[sub_vertex]) + == entity_vertices.end()) + { + entity_matches = false; + break; + } + } + if (entity_matches) + { + // Found entity in submesh with the same vertices as in the parent + // mesh + submesh_values[sub_entity] = tagged_values[i]; + entity_found = true; + sub_entity_to_parent[sub_entity] = entity; + } + } + } + } + } + + MeshTags new_meshtag(submesh_topology, tag_dim, submesh_indices, + submesh_values); + new_meshtag.name = tags.name; + EntityMap sub_entity_map(tags.topology(), submesh_topology, tag_dim, + submesh_indices); + return std::make_tuple(new_meshtag, sub_entity_map); +} + } // namespace dolfinx::mesh diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 186c45015d0..14f7d313f69 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -78,6 +78,7 @@ "to_string", "to_type", "transfer_meshtag", + "transfer_meshtags_to_submesh" "uniform_refine", ] @@ -1268,3 +1269,53 @@ def create_geometry( if (dtype := np.dtype(element.dtype)) != x.dtype: raise ValueError(f"Mismatch in x dtype ({x.dtype}) and coordinate element ({dtype})") return Geometry(ftype(index_map, dofmap, element._cpp_object, x, input_global_indices)) + + +def transfer_meshtags_to_submesh( + entity_tag: MeshTags, + submesh: Mesh, + vertex_to_parent: EntityMap, + cell_to_parent: EntityMap, +) -> tuple[MeshTags, EntityMap]: + """ + Transfer a ``entity_tag`` from a parent mesh to a ``submesh``. + + Note: + All entities not tagged in the parent mesh will be tagged with the numerical + maximum of the data type of ``entity_tag`` in the submesh. + + Args: + entity_tag: Tag to transfer + submesh: Submesh to transfer tag to + vertex_to_parent: Mapping from submesh vertices to parent mesh vertices + cell_to_parent: Mapping from submesh cells to parent entities + Returns: + A tuple (submesh_tag, sub_to_parent_entity_map) where: ``submesh_tag`` is the tag on the + submesh and ``sub_to_parent_entity_map`` is a mapping from submesh entities in the tag to + the corresponding entities in the parent. + """ + dim = entity_tag.dim + sub_tdim = submesh.topology.dim + if dim > sub_tdim: + raise RuntimeError( + f"Cannot transfer meshtags of dimension {dim} to submesh with topological dimension" + ) + + submesh.topology.create_connectivity(sub_tdim, sub_tdim) + submesh.topology.create_connectivity(entity_tag.dim, 0) + submesh.topology.create_connectivity(sub_tdim, entity_tag.dim) + entity_tag.topology.create_connectivity(dim, 0) + entity_tag.topology.create_connectivity(dim, sub_tdim) + dtype = entity_tag.values.dtype + if dtype == np.int32: + ftype = _cpp.mesh.transfer_meshtags_to_submesh_int32 + elif dtype == np.int64: + ftype = _cpp.mesh.transfer_meshtags_to_submesh_int64 + elif dtype == np.float64: + ftype = _cpp.mesh.transfer_meshtags_to_submesh_float64 + else: + raise NotImplementedError(f"MeshTags with dtype {dtype} not supported for transfer to submesh.") + cpp_tag, sub_to_parent_entity_map = ftype( + entity_tag._cpp_object, submesh.topology._cpp_object, vertex_to_parent._cpp_object, cell_to_parent._cpp_object + ) + return MeshTags(cpp_tag), sub_to_parent_entity_map \ No newline at end of file diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h index 24f22d14247..86298c5c1b4 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h @@ -131,6 +131,20 @@ void declare_meshtags(nb::module_& m, const std::string& type) return dolfinx::mesh::create_meshtags( topology, dim, entities, std::span(values.data(), values.size())); }); + std::string pyfunc_name = "transfer_meshtags_to_submesh_" + type; + m.def( + pyfunc_name.c_str(), + [](const dolfinx::mesh::MeshTags& tags, + std::shared_ptr submesh_topology, + const dolfinx::mesh::EntityMap& vertex_map, + const dolfinx::mesh::EntityMap& cell_map) + { + return dolfinx::mesh::transfer_meshtags_to_submesh( + tags, submesh_topology, + vertex_map, cell_map); + }, + nanobind::arg("tags"), nanobind::arg("submesh_topology"), + nanobind::arg("vertex_map"), nanobind::arg("cell_map")); } template diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index d86035de9bd..cbb6d0aec05 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -181,7 +181,7 @@ void mesh(nb::module_& m) { std::vector mapped_entities = self.sub_topology_to_topology( - std::span(entities.data(), entities.size()), inverse); + std::span(entities.data(), entities.size()), inverse); return as_nbarray(std::move(mapped_entities)); }, nb::arg("entities"), nb::arg("inverse")) From d04822e54cfacd10aaf12d5215fa613468649380 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 14:04:22 +0000 Subject: [PATCH 2/7] Wrap as python object --- python/dolfinx/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 14f7d313f69..7f4d60f59dd 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -1318,4 +1318,4 @@ def transfer_meshtags_to_submesh( cpp_tag, sub_to_parent_entity_map = ftype( entity_tag._cpp_object, submesh.topology._cpp_object, vertex_to_parent._cpp_object, cell_to_parent._cpp_object ) - return MeshTags(cpp_tag), sub_to_parent_entity_map \ No newline at end of file + return MeshTags(cpp_tag), EntityMap(sub_to_parent_entity_map) \ No newline at end of file From 67a625016515ad36c3449e767f7b9b6aed80cb96 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 14:17:42 +0000 Subject: [PATCH 3/7] Fix uint --- cpp/dolfinx/mesh/utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 3381b815ff7..13d11317a3a 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1502,7 +1502,7 @@ std::tuple, EntityMap> transfer_meshtags_to_submesh( } std::vector parent_entity_to_sub_cell(e_to_v->num_nodes(), -1); - for (std::int32_t i = 0; i < sub_cell_to_parent_entity.size(); ++i) + for (std::size_t i = 0; i < sub_cell_to_parent_entity.size(); ++i) parent_entity_to_sub_cell[sub_cell_to_parent_entity[i]] = i; // Get map from submesh vertex to parent vertex From 45ed1fb23e10bfccd301d7ae4fc8df054adf57d8 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 15:32:20 +0000 Subject: [PATCH 4/7] Add early exist for edim == submesh tdim --- cpp/dolfinx/mesh/utils.h | 167 ++++++++++++++++++----------- python/dolfinx/mesh.py | 33 +++--- python/test/unit/mesh/test_mesh.py | 37 +++++++ 3 files changed, 160 insertions(+), 77 deletions(-) diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 13d11317a3a..ca94b7db5fd 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -1450,23 +1450,18 @@ create_submesh(const Mesh& mesh, int dim, /// @brief Transfer a meshtags object from a parent to a submesh. /// -/// @note All entities not tagged in the parent mesh is tagged with the -/// default value in the submesh, which is `std::numeric_limits::max()`. -/// /// @param[in] tags The meshtags object on the parent mesh. /// @param[in] submesh_topology The topology of the submesh. /// @param[in] vertex_map Map from submesh vertex to parent mesh vertex. /// @param[in] cell_map Map from submesh cell to parent mesh entity of the same /// dimension as the meshtags. -/// @return A meshtags object on the submesh, and the bi-directional entity map -/// between the tagged entities in the submesh and the parent mesh. +/// @return A meshtags object on the submesh. template -std::tuple, EntityMap> transfer_meshtags_to_submesh( +MeshTags transfer_meshtags_to_submesh( const MeshTags& tags, std::shared_ptr submesh_topology, const EntityMap& vertex_map, const EntityMap& cell_map) { - int tag_dim = tags.dim(); int submesh_tdim = submesh_topology->dim(); auto topology = tags.topology(); @@ -1475,7 +1470,6 @@ std::tuple, EntityMap> transfer_meshtags_to_submesh( throw std::runtime_error("Tag dimension must be less than or equal to " "submesh dimension"); } - std::shared_ptr sub_cell_imap = submesh_topology->index_map(submesh_tdim); if (!sub_cell_imap) @@ -1494,14 +1488,16 @@ std::tuple, EntityMap> transfer_meshtags_to_submesh( // Create a full lookup for all cells on the parent mesh, as the tag can have // entities that are not in the submesh - auto e_to_v = topology->connectivity(tag_dim, 0); - if (!e_to_v) + auto parent_entity_imap = topology->index_map(submesh_tdim); + if (!parent_entity_imap) { - throw std::runtime_error("Missing connectivity between " - + std::to_string(tag_dim) + " and 0"); + throw std::runtime_error(std::format( + "Entities of dimension {} does not exist in parent mesh topology.", + submesh_tdim)); } - - std::vector parent_entity_to_sub_cell(e_to_v->num_nodes(), -1); + std::size_t num_parent_entities + = parent_entity_imap->size_local() + parent_entity_imap->num_ghosts(); + std::vector parent_entity_to_sub_cell(num_parent_entities, -1); for (std::size_t i = 0; i < sub_cell_to_parent_entity.size(); ++i) parent_entity_to_sub_cell[sub_cell_to_parent_entity[i]] = i; @@ -1519,23 +1515,50 @@ std::tuple, EntityMap> transfer_meshtags_to_submesh( // Access various connectivity maps auto sub_e_to_v = submesh_topology->connectivity(tag_dim, 0); auto sub_c_to_e = submesh_topology->connectivity(submesh_tdim, tag_dim); - auto e_to_sub_cell = topology->connectivity(tag_dim, submesh_tdim); - if (!e_to_sub_cell) + auto sub_entity_imap = submesh_topology->index_map(tag_dim); + auto e_to_v = topology->connectivity(tag_dim, 0); + std::shared_ptr> + e_to_sub_cell = nullptr; + if (tag_dim != submesh_tdim) { - throw std::runtime_error("Missing connectivity between " - + std::to_string(tag_dim) + " and " - + std::to_string(submesh_tdim)); + e_to_sub_cell = topology->connectivity(tag_dim, submesh_tdim); + if (!e_to_sub_cell) + { + throw std::runtime_error( + std::format("Missing connectivity between {} and {} in parent mesh", + tag_dim, submesh_tdim)); + } } - auto sub_entity_imap = submesh_topology->index_map(tag_dim); + if (!sub_e_to_v) + { + throw std::runtime_error(std::format( + "Missing connectivity between {} and {} in submesh", tag_dim, 0)); + } + if (!sub_c_to_e) + { + throw std::runtime_error( + std::format("Missing connectivity between {} and {} in submesh", + submesh_tdim, tag_dim)); + } + if (!sub_entity_imap) + { + throw std::runtime_error(std::format( + "Entities of dimension {} does not exist in submesh topology.", + tag_dim)); + } + if (!e_to_v) + { + throw std::runtime_error( + std::format("Missing connectivity between {} and 0", tag_dim)); + } // Prepare sub entity to parent map - std::size_t num_sum_entities + std::size_t num_sub_entities = sub_entity_imap->size_local() + sub_entity_imap->num_ghosts(); - std::vector sub_entity_to_parent(num_sum_entities, -1); - std::vector submesh_values(num_sum_entities, - std::numeric_limits::max()); - std::vector submesh_indices(num_sum_entities); + constexpr T max_val = std::numeric_limits::max(); + std::vector submesh_values(num_sub_entities, max_val); + std::vector submesh_indices(num_sub_entities); std::iota(submesh_indices.begin(), submesh_indices.end(), 0); std::span tagged_entities = tags.indices(); @@ -1545,52 +1568,76 @@ std::tuple, EntityMap> transfer_meshtags_to_submesh( for (std::size_t i = 0; i < tagged_entities.size(); ++i) { auto entity = tagged_entities[i]; - auto entity_vertices = e_to_v->links(entity); - bool entity_found = false; - for (auto parent_cell : e_to_sub_cell->links(entity)) + + auto find_and_map_sub_entity + = [tag_dim, submesh_tdim, &e_to_v, &parent_entity_to_sub_cell, + &sub_to_parent_vertex, &sub_e_to_v, &sub_c_to_e, + &e_to_sub_cell](std::size_t i, std::int32_t entity) { - if (entity_found) - break; - std::int32_t sub_cell = parent_entity_to_sub_cell[parent_cell]; - if (sub_cell != -1) + // Fast exit if the tag dimension is the same as the submesh dimension, + // as we can directly map the parent entity to the submesh cell + if (tag_dim == submesh_tdim) + return parent_entity_to_sub_cell[entity]; + + auto entity_vertices = e_to_v->links(entity); + auto parent_sub_cells = e_to_sub_cell->links(entity); + // Check if any of the sub cells are in the submesh, and if so check if + // any of the entities in the submesh connected to the same cell share all + // vertices with the parent entity. + for (auto parent_cell : parent_sub_cells) { - // For a cell in the sub mesh find all attached entities, and - // define them by their vertices in the sub mesh + std::int32_t sub_cell = parent_entity_to_sub_cell[parent_cell]; + if (sub_cell == -1) + continue; + for (auto sub_entity : sub_c_to_e->links(sub_cell)) { - if (entity_found) - break; auto sub_vertices = sub_e_to_v->links(sub_entity); - bool entity_matches = true; - for (auto sub_vertex : sub_vertices) - { - if (std::find(entity_vertices.begin(), entity_vertices.end(), - sub_to_parent_vertex[sub_vertex]) - == entity_vertices.end()) - { - entity_matches = false; - break; - } - } + + // Replace the innermost loop with std::all_of + bool entity_matches + = std::all_of(sub_vertices.begin(), sub_vertices.end(), + [&](auto sub_vertex) + { + auto parent_v = sub_to_parent_vertex[sub_vertex]; + return std::find(entity_vertices.begin(), + entity_vertices.end(), parent_v) + != entity_vertices.end(); + }); + + // If a match is found, apply values and exit the lambda immediately if (entity_matches) - { - // Found entity in submesh with the same vertices as in the parent - // mesh - submesh_values[sub_entity] = tagged_values[i]; - entity_found = true; - sub_entity_to_parent[sub_entity] = entity; - } + return sub_entity; } } - } + return -1; + }; + + // Execute the search for the current entity + std::int32_t sub_entity = find_and_map_sub_entity(i, entity); + if (sub_entity != -1) + submesh_values[sub_entity] = tagged_values[i]; } - MeshTags new_meshtag(submesh_topology, tag_dim, submesh_indices, - submesh_values); + // Filter out the entities that were never mapped (values still equal max) + std::vector filtered_indices; + std::vector filtered_values; + filtered_indices.reserve(num_sub_entities); + filtered_values.reserve(num_sub_entities); + for (std::size_t i = 0; i < submesh_values.size(); ++i) + { + if (submesh_values[i] != max_val) + { + filtered_indices.push_back(submesh_indices[i]); + filtered_values.push_back(submesh_values[i]); + } + } + filtered_indices.shrink_to_fit(); + filtered_values.shrink_to_fit(); + MeshTags new_meshtag(submesh_topology, tag_dim, filtered_indices, + filtered_values); new_meshtag.name = tags.name; - EntityMap sub_entity_map(tags.topology(), submesh_topology, tag_dim, - submesh_indices); - return std::make_tuple(new_meshtag, sub_entity_map); + return new_meshtag; } } // namespace dolfinx::mesh diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 7f4d60f59dd..464ce24d6f9 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -78,7 +78,7 @@ "to_string", "to_type", "transfer_meshtag", - "transfer_meshtags_to_submesh" + "transfer_meshtags_to_submesh", "uniform_refine", ] @@ -1276,23 +1276,17 @@ def transfer_meshtags_to_submesh( submesh: Mesh, vertex_to_parent: EntityMap, cell_to_parent: EntityMap, -) -> tuple[MeshTags, EntityMap]: - """ - Transfer a ``entity_tag`` from a parent mesh to a ``submesh``. - - Note: - All entities not tagged in the parent mesh will be tagged with the numerical - maximum of the data type of ``entity_tag`` in the submesh. +) -> MeshTags: + """Transfer a ``entity_tag`` from a parent mesh to a ``submesh``. Args: entity_tag: Tag to transfer submesh: Submesh to transfer tag to - vertex_to_parent: Mapping from submesh vertices to parent mesh vertices + vertex_to_parent: Mapping from submesh vertices to parent + mesh vertices cell_to_parent: Mapping from submesh cells to parent entities Returns: - A tuple (submesh_tag, sub_to_parent_entity_map) where: ``submesh_tag`` is the tag on the - submesh and ``sub_to_parent_entity_map`` is a mapping from submesh entities in the tag to - the corresponding entities in the parent. + The transferred meshtags object on the submesh. """ dim = entity_tag.dim sub_tdim = submesh.topology.dim @@ -1306,7 +1300,7 @@ def transfer_meshtags_to_submesh( submesh.topology.create_connectivity(sub_tdim, entity_tag.dim) entity_tag.topology.create_connectivity(dim, 0) entity_tag.topology.create_connectivity(dim, sub_tdim) - dtype = entity_tag.values.dtype + dtype = entity_tag.values.dtype if dtype == np.int32: ftype = _cpp.mesh.transfer_meshtags_to_submesh_int32 elif dtype == np.int64: @@ -1314,8 +1308,13 @@ def transfer_meshtags_to_submesh( elif dtype == np.float64: ftype = _cpp.mesh.transfer_meshtags_to_submesh_float64 else: - raise NotImplementedError(f"MeshTags with dtype {dtype} not supported for transfer to submesh.") - cpp_tag, sub_to_parent_entity_map = ftype( - entity_tag._cpp_object, submesh.topology._cpp_object, vertex_to_parent._cpp_object, cell_to_parent._cpp_object + raise NotImplementedError( + f"MeshTags with dtype {dtype} not supported for transfer to submesh." + ) + cpp_tag = ftype( + entity_tag._cpp_object, + submesh.topology._cpp_object, + vertex_to_parent._cpp_object, + cell_to_parent._cpp_object, ) - return MeshTags(cpp_tag), EntityMap(sub_to_parent_entity_map) \ No newline at end of file + return MeshTags(cpp_tag) diff --git a/python/test/unit/mesh/test_mesh.py b/python/test/unit/mesh/test_mesh.py index 08804d65ede..253c3ed4aa0 100644 --- a/python/test/unit/mesh/test_mesh.py +++ b/python/test/unit/mesh/test_mesh.py @@ -36,6 +36,7 @@ exterior_facet_indices, locate_entities, locate_entities_boundary, + transfer_meshtags_to_submesh, ) @@ -782,3 +783,39 @@ def test_mesh_single_process_distribution(partitioner): adj = mesh.topology.connectivity(*conn) for i in range(adj.num_nodes): assert adj.links(i).size == 2 + + +@pytest.mark.parametrize("codim", [0, 1, 2, 3]) +def test_transfer_to_submesh(codim): + mesh = create_unit_cube(MPI.COMM_WORLD, 8, 4, 5) + tdim = mesh.topology.dim + assert tdim - codim >= 0 + entities = locate_entities(mesh, tdim - codim, lambda x: x[0] >= 0.5) + submesh, entity_map, vertex_map, _node_map = create_submesh(mesh, tdim - codim, entities) + + def marker1(x): + return x[1] <= 0.5 + + def marker2(x): + return x[2] < 0.5 + + for i in range(submesh.topology.dim + 1): + mesh.topology.create_entities(i) + pe_map = mesh.topology.index_map(i) + num_parent_entities = pe_map.size_local + pe_map.num_ghosts + values = np.zeros(num_parent_entities, dtype=np.int32) + values[locate_entities(mesh, i, marker1)] = 1 + values[locate_entities(mesh, i, marker2)] = 2 + et_indices = np.flatnonzero(values) + et_values = values[et_indices] + et = dolfinx.mesh.meshtags(mesh, i, et_indices, et_values) + + sub_et = transfer_meshtags_to_submesh(et, submesh, vertex_map, entity_map) + ref_one = locate_entities(submesh, i, marker1) + ref_two = np.sort(locate_entities(submesh, i, marker2)) + ref_one = np.setdiff1d(ref_one, ref_two, assume_unique=True) + + marked1 = sub_et.find(1) + marked2 = sub_et.find(2) + np.testing.assert_allclose(marked1, ref_one) + np.testing.assert_allclose(marked2, ref_two) From 7ab9e2a5cef2fb1280da4cf78edec84591babdb9 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 15:38:02 +0000 Subject: [PATCH 5/7] Clang formatting --- cpp/dolfinx/fem/interpolate.h | 17 ++++--- cpp/dolfinx/mesh/EntityMap.h | 92 +++++++++++++++++------------------ cpp/dolfinx/mesh/Topology.h | 4 +- 3 files changed, 58 insertions(+), 55 deletions(-) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index f87771add6b..5f24cc79d03 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -36,8 +36,6 @@ concept MDSpan = requires(T x, std::size_t idx) { { x.extent(1) } -> std::integral; }; - - /// @brief Compute the evaluation points in the physical space at which /// an expression should be computed to interpolate it in a finite /// element space. @@ -135,7 +133,8 @@ std::vector interpolation_coords(const fem::FiniteElement& element, /// calling `interpolation_coords`. template void interpolate(Function& u, std::span f, - std::array fshape, mesh::CellRange auto&& cells); + std::array fshape, + mesh::CellRange auto&& cells); namespace impl { @@ -365,7 +364,8 @@ void interpolation_apply(U&& Pi, V&& data, std::span coeffs, int bs) /// by the function. template void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, - const Function& u0, mesh::CellRange auto&& cells0) + const Function& u0, + mesh::CellRange auto&& cells0) { auto V0 = u0.function_space(); assert(V0); @@ -465,7 +465,8 @@ void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, /// @pre Function%s `u1` and `u0` must share the same mesh. This is not /// checked by the function. template -void interpolate_nonmatching_maps(Function& u1, mesh::CellRange auto&& cells1, +void interpolate_nonmatching_maps(Function& u1, + mesh::CellRange auto&& cells1, const Function& u0, mesh::CellRange auto&& cells0) { @@ -817,7 +818,8 @@ void point_evaluation(const FiniteElement& element, bool symmetric, /// @param [out] coeffs Output Function coefficients. template void identity_mapped_evaluation(const FiniteElement& element, bool symmetric, - const DofMap& dofmap, mesh::CellRange auto&& cells, + const DofMap& dofmap, + mesh::CellRange auto&& cells, std::span cell_info, std::span f, std::array fshape, @@ -1093,7 +1095,8 @@ geometry::PointOwnershipData create_interpolation_data( template void interpolate(Function& u, std::span f, - std::array fshape, mesh::CellRange auto&& cells) + std::array fshape, + mesh::CellRange auto&& cells) { // TODO: Index for mixed-topology, zero for now const int index = 0; diff --git a/cpp/dolfinx/mesh/EntityMap.h b/cpp/dolfinx/mesh/EntityMap.h index 9711519f2e1..8efedd679a6 100644 --- a/cpp/dolfinx/mesh/EntityMap.h +++ b/cpp/dolfinx/mesh/EntityMap.h @@ -9,8 +9,8 @@ #include "Topology.h" #include #include -#include #include +#include #include namespace dolfinx::mesh @@ -100,55 +100,55 @@ class EntityMap /// `this->sub_topology()`. /// @return A list of mapped entity indices. Entities that do not /// exist in the target topology are marked as -1. -std::vector -sub_topology_to_topology(CellRange auto&& entities, - bool inverse) const -{ - if (!inverse) - { - // In this case, we want to map from entity indices in - // `_sub_topology` to corresponding entities in `_topology`. Hence, - // for each index in `entities`, we get the corresponding index in - // `_topology` using `_sub_topology_to_topology` - auto mapped - = entities - | std::views::transform([this](std::int32_t i) - { return _sub_topology_to_topology[i]; }); - return std::vector(mapped.begin(), mapped.end()); - } - else + std::vector sub_topology_to_topology(CellRange auto&& entities, + bool inverse) const { - // In this case, we are mapping from entity indices in `_topology` - // to entity indices in `_sub_topology`. Hence, we first need to - // construct the "inverse" of `_sub_topology_to_topology` - std::unordered_map topology_to_sub_topology; - topology_to_sub_topology.reserve(_sub_topology_to_topology.size()); - for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i) + if (!inverse) { - topology_to_sub_topology.insert( - {_sub_topology_to_topology[i], static_cast(i)}); + // In this case, we want to map from entity indices in + // `_sub_topology` to corresponding entities in `_topology`. Hence, + // for each index in `entities`, we get the corresponding index in + // `_topology` using `_sub_topology_to_topology` + auto mapped + = entities + | std::views::transform([this](std::int32_t i) + { return _sub_topology_to_topology[i]; }); + return std::vector(mapped.begin(), mapped.end()); + } + else + { + // In this case, we are mapping from entity indices in `_topology` + // to entity indices in `_sub_topology`. Hence, we first need to + // construct the "inverse" of `_sub_topology_to_topology` + std::unordered_map topology_to_sub_topology; + topology_to_sub_topology.reserve(_sub_topology_to_topology.size()); + for (std::size_t i = 0; i < _sub_topology_to_topology.size(); ++i) + { + topology_to_sub_topology.insert( + {_sub_topology_to_topology[i], static_cast(i)}); + } + + // For each entity index in `entities` (which are indices in + // `_topology`), get the corresponding entity in `_sub_topology`. + // Since `_sub_topology` consists of a subset of entities in + // `_topology`, there are entities in topology that may not exist in + // `_sub_topology`. If this is the case, mark those entities with + // -1. + + auto mapped = entities + | std::views::transform( + [&topology_to_sub_topology](std::int32_t i) + { + // Map the entity if it exists. If it doesn't, mark + // with -1. + auto it = topology_to_sub_topology.find(i); + return (it != topology_to_sub_topology.end()) + ? it->second + : -1; + }); + return std::vector(mapped.begin(), mapped.end()); } - - // For each entity index in `entities` (which are indices in - // `_topology`), get the corresponding entity in `_sub_topology`. - // Since `_sub_topology` consists of a subset of entities in - // `_topology`, there are entities in topology that may not exist in - // `_sub_topology`. If this is the case, mark those entities with - // -1. - - auto mapped - = entities - | std::views::transform( - [&topology_to_sub_topology](std::int32_t i) - { - // Map the entity if it exists. If it doesn't, mark with - // -1. - auto it = topology_to_sub_topology.find(i); - return (it != topology_to_sub_topology.end()) ? it->second : -1; - }); - return std::vector(mapped.begin(), mapped.end()); } -} private: // Dimension of the entities diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h index 814919adf5e..c3b4f1ce18c 100644 --- a/cpp/dolfinx/mesh/Topology.h +++ b/cpp/dolfinx/mesh/Topology.h @@ -7,6 +7,7 @@ #pragma once #include +#include #include #include #include @@ -18,7 +19,6 @@ #include #include #include -#include namespace dolfinx::common { @@ -27,7 +27,7 @@ class IndexMap; namespace dolfinx::mesh { - /// Requirement on range of cell indices. +/// Requirement on range of cell indices. template concept CellRange = std::ranges::input_range and std::ranges::sized_range and std::is_integral_v< From 654addc76e0a94a6dd43d262e7ea54f7d71423c4 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 22 Apr 2026 15:40:11 +0000 Subject: [PATCH 6/7] Clang format bindings --- python/dolfinx/wrappers/dolfinx_wrappers/mesh.h | 3 +-- python/dolfinx/wrappers/mesh.cpp | 4 +++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h index 86298c5c1b4..ca71c264b03 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h @@ -140,8 +140,7 @@ void declare_meshtags(nb::module_& m, const std::string& type) const dolfinx::mesh::EntityMap& cell_map) { return dolfinx::mesh::transfer_meshtags_to_submesh( - tags, submesh_topology, - vertex_map, cell_map); + tags, submesh_topology, vertex_map, cell_map); }, nanobind::arg("tags"), nanobind::arg("submesh_topology"), nanobind::arg("vertex_map"), nanobind::arg("cell_map")); diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index cbb6d0aec05..f68c1bb68a5 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -181,7 +181,9 @@ void mesh(nb::module_& m) { std::vector mapped_entities = self.sub_topology_to_topology( - std::span(entities.data(), entities.size()), inverse); + std::span(entities.data(), + entities.size()), + inverse); return as_nbarray(std::move(mapped_entities)); }, nb::arg("entities"), nb::arg("inverse")) From 35bdf80b6935dda8cbf71797f8bb68c4d9a30d21 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 22 Apr 2026 20:24:14 +0000 Subject: [PATCH 7/7] Use C++20 to clarify logic --- cpp/dolfinx/mesh/EntityMap.h | 4 +-- cpp/dolfinx/mesh/utils.h | 52 +++++++++++++++++++----------------- 2 files changed, 29 insertions(+), 27 deletions(-) diff --git a/cpp/dolfinx/mesh/EntityMap.h b/cpp/dolfinx/mesh/EntityMap.h index 8efedd679a6..eeb64c67df8 100644 --- a/cpp/dolfinx/mesh/EntityMap.h +++ b/cpp/dolfinx/mesh/EntityMap.h @@ -110,7 +110,7 @@ class EntityMap // for each index in `entities`, we get the corresponding index in // `_topology` using `_sub_topology_to_topology` auto mapped - = entities + = std::forward(entities) | std::views::transform([this](std::int32_t i) { return _sub_topology_to_topology[i]; }); return std::vector(mapped.begin(), mapped.end()); @@ -135,7 +135,7 @@ class EntityMap // `_sub_topology`. If this is the case, mark those entities with // -1. - auto mapped = entities + auto mapped = std::forward(entities) | std::views::transform( [&topology_to_sub_topology](std::int32_t i) { diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index ca94b7db5fd..c4a59b1fe78 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -147,7 +147,7 @@ compute_vertex_coords_boundary(const mesh::Mesh& mesh, int dim, // Get first cell and find position const std::int32_t c = v_to_c->links(v).front(); auto cell_vertices = c_to_v->links(c); - auto it = std::find(cell_vertices.begin(), cell_vertices.end(), v); + auto it = std::ranges::find(cell_vertices, v); assert(it != cell_vertices.end()); const std::size_t local_pos = std::distance(cell_vertices.begin(), it); @@ -1567,43 +1567,45 @@ MeshTags transfer_meshtags_to_submesh( // entity for (std::size_t i = 0; i < tagged_entities.size(); ++i) { - auto entity = tagged_entities[i]; - auto find_and_map_sub_entity = [tag_dim, submesh_tdim, &e_to_v, &parent_entity_to_sub_cell, &sub_to_parent_vertex, &sub_e_to_v, &sub_c_to_e, - &e_to_sub_cell](std::size_t i, std::int32_t entity) + &e_to_sub_cell](std::int32_t entity) { // Fast exit if the tag dimension is the same as the submesh dimension, // as we can directly map the parent entity to the submesh cell if (tag_dim == submesh_tdim) return parent_entity_to_sub_cell[entity]; + // Given an entity in the parent meshtag, find all submesh-cells that are + // entities in parent mesh that contain this entity. auto entity_vertices = e_to_v->links(entity); auto parent_sub_cells = e_to_sub_cell->links(entity); - // Check if any of the sub cells are in the submesh, and if so check if - // any of the entities in the submesh connected to the same cell share all - // vertices with the parent entity. - for (auto parent_cell : parent_sub_cells) + auto submesh_cells + = parent_sub_cells + | std::views::transform([&parent_entity_to_sub_cell](auto c) + { return parent_entity_to_sub_cell[c]; }) + | std::views::filter([](auto sub_cell) { return sub_cell != -1; }); + for (auto sub_cell : submesh_cells) { - std::int32_t sub_cell = parent_entity_to_sub_cell[parent_cell]; - if (sub_cell == -1) - continue; - for (auto sub_entity : sub_c_to_e->links(sub_cell)) { - auto sub_vertices = sub_e_to_v->links(sub_entity); - - // Replace the innermost loop with std::all_of - bool entity_matches - = std::all_of(sub_vertices.begin(), sub_vertices.end(), - [&](auto sub_vertex) - { - auto parent_v = sub_to_parent_vertex[sub_vertex]; - return std::find(entity_vertices.begin(), - entity_vertices.end(), parent_v) - != entity_vertices.end(); - }); + // Convert submesh entity vertices to parent vertices + auto parent_vertices + = sub_e_to_v->links(sub_entity) + | std::views::transform([&](auto v) + { return sub_to_parent_vertex[v]; }); + + // Check if all parent vertices of the submesh entity are in the + // parent entity + bool entity_matches = std::ranges::all_of( + parent_vertices, + [&](auto p_v) + { + // With C++23 this can use std::ranges::contains + return std::ranges::find(entity_vertices, p_v) + != std::ranges::end(entity_vertices); + }); // If a match is found, apply values and exit the lambda immediately if (entity_matches) @@ -1614,7 +1616,7 @@ MeshTags transfer_meshtags_to_submesh( }; // Execute the search for the current entity - std::int32_t sub_entity = find_and_map_sub_entity(i, entity); + std::int32_t sub_entity = find_and_map_sub_entity(tagged_entities[i]); if (sub_entity != -1) submesh_values[sub_entity] = tagged_values[i]; }