diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp index 0760bc3aa68..d17b9353eae 100644 --- a/cpp/demo/codim_0_assembly/main.cpp +++ b/cpp/demo/codim_0_assembly/main.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -99,14 +100,11 @@ int main(int argc, char* argv[]) // `submesh`, we must provide a map from entities in `mesh` to // entities in `submesh`. This is simply the "inverse" of // `submesh_to_mesh`. - std::vector mesh_to_submesh(num_cells_local, -1); - for (std::size_t i = 0; i < submesh_to_mesh.size(); ++i) - mesh_to_submesh[submesh_to_mesh[i]] = i; - - std::map>, - std::span> - entity_maps = {{submesh, mesh_to_submesh}}; + const mesh::EntityMap entity_map(mesh->topology(), submesh->topology(), + tdim, submesh_to_mesh); + std::vector> entity_maps + = {std::make_shared(entity_map)}; // Next we compute the integration entities on the integration // domain `mesh` std::vector integration_entities diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index faa6339c654..f3fb368a126 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -1,4 +1,5 @@ -// Copyright (C) 2019-2025 Garth N. Wells and Chris Richardson +// Copyright (C) 2019-2025 Garth N. Wells, Chris Richardson, Joseph P. Dean and +// Jørgen S. Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -14,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -44,10 +46,17 @@ namespace impl { /// @brief Compute "`entity_map[entities[i]]`", with suitable handling /// for facets. -/// @param entities -/// @param entity_map -/// @param codim -/// @param c_to_f +/// @param entities List of integration entities. For cells this is a list of +/// local cell indices. For facets this is an mdspan of shape (num_facets, 2), +/// where the first column is the local cell index and the second column is +/// the local facet index (relative to the cell). +/// @param entity_map A map from an entity (local index relative to the process) +/// to a set of local entities in another domain. +/// @param codim The codimension of the mapped entities. If 0, the other domain +/// consists of cells of the input domain. If 1, the other domain consists of +/// facets of the input domain. +/// @param c_to_f If codimension is 1, this is the map from a cell in the input +/// domain to a facet of this domain. /// @return entity_map[entities[i]] std::vector compute_domain( auto entities, std::span entity_map, @@ -62,24 +71,31 @@ std::vector compute_domain( mapped_entities.reserve(entities.size()); if constexpr (entities.rank() == 1) { + // Map cells from integration mesh to argument/coefficient mesh std::span ents(entities.data_handle(), entities.size()); std::ranges::transform(ents, std::back_inserter(mapped_entities), [&entity_map](auto e) { return entity_map[e]; }); } else if (entities.rank() == 2) { + // Map facets from integration mesh to argument/coefficient mesh assert(codim.value() >= 0); if (codim.value() == 0) { + // Codim 0 case for (std::size_t i = 0; i < entities.extent(0); ++i) { - // Add cell and the local facet index + // Local facet index is preserved in submesh creation, so we just map + // the cell from the integration mesh to the argument/coefficient mesh + assert(static_cast(entities(i, 0)) < entity_map.size()); mapped_entities.insert(mapped_entities.end(), {entity_map[entities(i, 0)], entities(i, 1)}); } } else if (codim.value() == 1) { + // Codim 1 case + // In this case, the entity maps take facets in (`_mesh`) to cells // in `mesh`, so we need to get the facet number from the (cell, // local_facet pair) first. @@ -88,6 +104,7 @@ std::vector compute_domain( // Get the facet index, and add cell and the local facet index std::int32_t facet = c_to_f->get().links(entities(i, 0))[entities(i, 1)]; + assert(static_cast(facet) < entity_map.size()); mapped_entities.insert(mapped_entities.end(), {entity_map[facet], entities(i, 1)}); } @@ -97,7 +114,17 @@ std::vector compute_domain( } else throw std::runtime_error("Integral type not supported."); - + // Check that there are no negative values in the mapped entities + if (!mapped_entities.empty()) + { + auto min_val + = std::min_element(mapped_entities.begin(), mapped_entities.end()); + if (*min_val < 0) + { + throw std::runtime_error( + "Mapped entities contain negative values, check entity map."); + } + } return mapped_entities; } } // namespace impl @@ -221,8 +248,7 @@ class Form const std::vector>>& constants, bool needs_facet_permutations, - const std::map>, - std::span>& entity_maps) + const std::vector>& entity_maps) : _function_spaces(V), _integrals(std::forward(integrals)), _mesh(mesh), _coefficients(coefficients), _constants(constants), _needs_facet_permutations(needs_facet_permutations) @@ -230,32 +256,56 @@ class Form if (!_mesh) throw std::runtime_error("Form Mesh is null."); - // Check consistency of mesh(es) + // Check consistency of mesh(es) and find position of each mesh + // in the entity_map + std::vector argument_position(_function_spaces.size(), -1); + std::vector coefficient_position(_coefficients.size(), -1); { // Integration domain mesh is passed, so check that it is (1) // common for spaces and coefficients (2) or an entity_map is // available - for (auto& space : _function_spaces) + for (std::size_t i = 0; i < _function_spaces.size(); ++i) { - if (auto mesh0 = space->mesh(); - mesh0 != _mesh and !entity_maps.contains(mesh0)) + auto mesh0 = _function_spaces[i]->mesh(); + bool mesh_found = false; + for (std::size_t j = 0; j < entity_maps.size(); ++j) + { + assert(entity_maps[j]); + if (entity_maps[j]->contains(mesh0->topology())) + { + argument_position[i] = j; + mesh_found = true; + break; + } + } + if (mesh0 != _mesh and !mesh_found) { throw std::runtime_error( "Incompatible mesh. argument entity_maps must be provided."); } } - for (auto& c : coefficients) + for (std::size_t i = 0; i < _coefficients.size(); ++i) { - if (auto mesh0 = c->function_space()->mesh(); - mesh0 != _mesh and !entity_maps.contains(mesh0)) + auto mesh0 = _coefficients[i]->function_space()->mesh(); + bool mesh_found = false; + for (std::size_t j = 0; j < entity_maps.size(); ++j) + { + assert(entity_maps[j]); + if (entity_maps[j]->contains(mesh0->topology())) + { + coefficient_position[i] = j; + mesh_found = true; + break; + } + } + if (mesh0 != _mesh and !mesh_found) { throw std::runtime_error( "Incompatible mesh. coefficient entity_maps must be provided."); } } } - - for (auto& space : _function_spaces) + for (std::size_t i = 0; i < _function_spaces.size(); ++i) { // Working map: [integral type, domain ID, kernel_idx]->entities std::map, @@ -263,16 +313,31 @@ class Form std::span>> vdata; - if (auto mesh0 = space->mesh(); mesh0 == _mesh) + if (auto mesh0 = _function_spaces[i]->mesh(); mesh0 == _mesh) { for (auto& [key, integral] : _integrals) vdata.insert({key, std::span(integral.entities)}); } else { - auto it = entity_maps.find(mesh0); - assert(it != entity_maps.end()); - std::span entity_map = it->second; + std::int32_t apos = argument_position[i]; + assert(apos != -1); + auto emap = entity_maps[apos]; + assert(emap); + // Create an entity map from mesh0 to the integration domain + std::span entities0 + = emap->get_entities(mesh0->topology()); + std::span entities + = emap->get_entities(_mesh->topology()); + auto e_imap = _mesh->topology()->index_map(emap->dim()); + if (!e_imap) + throw std::runtime_error( + "No index map for entities, call `Topology::create_entities(" + + std::to_string(emap->dim()) + ")"); + std::vector entity_map( + e_imap->size_local() + e_imap->num_ghosts(), -1); + for (std::size_t i = 0; i < entities0.size(); ++i) + entity_map[entities[i]] = entities0[i]; for (auto& [key, itg] : _integrals) { auto [type, id, kernel_idx] = key; @@ -281,7 +346,7 @@ class Form { e = impl::compute_domain( md::mdspan(itg.entities.data(), itg.entities.size()), - entity_map); + std::span(entity_map)); } else if (type == IntegralType::exterior_facet or type == IntegralType::interior_facet) @@ -296,7 +361,7 @@ class Form md::mdspan>( itg.entities.data(), itg.entities.size() / 2, 2), - entity_map, codim, *c_to_f); + std::span(entity_map), codim, *c_to_f); } else throw std::runtime_error("Integral type not supported."); @@ -319,15 +384,31 @@ class Form } else { - auto it = entity_maps.find(mesh0); - assert(it != entity_maps.end()); - std::span entity_map = it->second; + std::int32_t cpos = coefficient_position[c]; + assert(cpos != -1); + // Create an entity map from mesh0 to the integration domain + auto emap = entity_maps[cpos]; + assert(emap); + std::span entities0 + = emap->get_entities(mesh0->topology()); + std::span entities + = emap->get_entities(_mesh->topology()); + auto e_imap = _mesh->topology()->index_map(emap->dim()); + if (!e_imap) + throw std::runtime_error( + "No index map for entities, call `Topology::create_entities(" + + std::to_string(entity_maps[cpos]->dim()) + ")"); + std::vector entity_map( + e_imap->size_local() + e_imap->num_ghosts(), -1); + for (std::size_t i = 0; i < entities0.size(); ++i) + entity_map[entities[i]] = entities0[i]; + std::vector e; if (type == IntegralType::cell) { e = impl::compute_domain( md::mdspan(integral.entities.data(), integral.entities.size()), - entity_map); + std::span(entity_map)); } else if (type == IntegralType::exterior_facet or type == IntegralType::interior_facet) @@ -342,7 +423,7 @@ class Form md::mdspan>( integral.entities.data(), integral.entities.size() / 2, 2), - entity_map, codim, *c_to_f); + std::span(entity_map), codim, *c_to_f); } else throw std::runtime_error("Integral type not supported."); diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 8ba38ef74c9..d43894539e5 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -370,8 +371,7 @@ Form create_form_factory( IntegralType, std::vector>>>& subdomains, - const std::map>, - std::span>& entity_maps, + const std::vector>& entity_maps, std::shared_ptr> mesh = nullptr) { for (const ufcx_form& ufcx_form : ufcx_forms) @@ -414,14 +414,6 @@ Form create_form_factory( mesh = spaces.front()->mesh(); if (!mesh) throw std::runtime_error("No mesh could be associated with the Form."); - for (auto& V : spaces) - { - if (mesh != V->mesh() and entity_maps.find(V->mesh()) == entity_maps.end()) - { - throw std::runtime_error( - "Incompatible mesh. entity_maps must be provided."); - } - } auto topology = mesh->topology(); assert(topology); @@ -790,8 +782,7 @@ Form create_form( IntegralType, std::vector>>>& subdomains, - const std::map>, - std::span>& entity_maps, + const std::vector>& entity_maps, std::shared_ptr> mesh = nullptr) { // Place coefficients in appropriate order @@ -850,8 +841,7 @@ Form create_form( IntegralType, std::vector>>>& subdomains, - const std::map>, - std::span>& entity_maps, + const std::vector>& entity_maps, std::shared_ptr> mesh = nullptr) { ufcx_form* form = fptr(); diff --git a/cpp/dolfinx/mesh/CMakeLists.txt b/cpp/dolfinx/mesh/CMakeLists.txt index 733b3202d4a..f15f6c436f5 100644 --- a/cpp/dolfinx/mesh/CMakeLists.txt +++ b/cpp/dolfinx/mesh/CMakeLists.txt @@ -10,6 +10,7 @@ set(HEADERS_mesh ${CMAKE_CURRENT_SOURCE_DIR}/permutationcomputation.h ${CMAKE_CURRENT_SOURCE_DIR}/topologycomputation.h ${CMAKE_CURRENT_SOURCE_DIR}/utils.h + ${CMAKE_CURRENT_SOURCE_DIR}/EntityMap.h PARENT_SCOPE ) diff --git a/cpp/dolfinx/mesh/EntityMap.h b/cpp/dolfinx/mesh/EntityMap.h new file mode 100644 index 00000000000..a8c86bc628f --- /dev/null +++ b/cpp/dolfinx/mesh/EntityMap.h @@ -0,0 +1,131 @@ +// Copyright (C) 2025 Jørgen S. Dokken +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once +#include "Topology.h" +#include +#include +#include +#include +namespace dolfinx::mesh +{ +/// @brief A map between entities of two meshes +class EntityMap +{ +public: + /// @brief Constructor of a map between a set of entities belonging to + /// two meshes. + /// @tparam U + /// @param topology0 The first topology in the mapping relation + /// @param topology1 + /// @param dim Topological dimension of the mapped entities + /// @param entities0 The entities belonging to the first mesh + /// @param entities1 The entities belonging to the second mesh + template + requires std::is_convertible_v, + std::vector> + EntityMap(std::shared_ptr topology0, + std::shared_ptr topology1, int dim, U&& entities0, + U&& entities1) + : _dim(dim), _topology0(topology0), + _entities0(std::forward(entities0)), _topology1(topology1), + _entities1(std::forward(entities1)) + { + if (_entities0.size() != _entities1.size()) + { + throw std::runtime_error("Entities must have the same size."); + } + } + + /// @brief Constructor of a map between a set of entities belonging to two + /// meshes. + /// + /// Entity `i` in mesh1 is assumed to map to `entities0[i]`. + /// + /// @tparam U + /// @param topology0 The first topology in the mapping relation + /// @param topology1 The second topology in the mapping relation + /// @param dim Topological dimension of the mapped entities + /// @param entities0 The entities belonging to the first mesh + template + requires std::is_convertible_v, + std::vector> + EntityMap(std::shared_ptr topology0, + std::shared_ptr topology1, int dim, U&& entities0) + : _dim(dim), _topology0(topology0), + _entities0(std::forward(entities0)), _topology1(topology1) + { + auto e_map = topology1->index_map(dim); + if (!e_map) + { + throw std::runtime_error( + "No index map for entities, call `Topology::create_entities(" + + std::to_string(dim) + ")"); + } + std::size_t num_ents + = static_cast(e_map->size_local() + e_map->num_ghosts()); + if (num_ents != _entities0.size()) + throw std::runtime_error("Size mismatch between entities and index map."); + _entities1.resize(_entities0.size()); + std::iota(_entities1.begin(), _entities1.end(), 0); + assert(_entities1.size() == _entities0.size()); + } + + /// Copy constructor + EntityMap(const EntityMap& map) = default; + + /// Move constructor + EntityMap(EntityMap&& map) = default; + + // Destructor + ~EntityMap() = default; + + /// TODO + bool contains(std::shared_ptr topology) const + { + return topology == _topology0 or topology == _topology1; + } + + /// TODO + std::uint8_t topology_index(std::shared_ptr topology) const + { + if (topology == _topology0) + return 0; + else if (topology == _topology1) + return 1; + else + throw std::runtime_error("Topology not in the map."); + } + + /// TODO + std::span + get_entities(std::shared_ptr topology) const + { + if (topology == _topology0) + return std::span(_entities0.data(), + _entities0.size()); + else if (topology == _topology1) + return std::span(_entities1.data(), + _entities1.size()); + else + throw std::runtime_error("Topology not in the map."); + } + + /// TODO + std::size_t dim() const { return _dim; } + +private: + std::size_t _dim; ///< Dimension of the entities + std::shared_ptr _topology0; ///< The first mesh + std::vector + _entities0; ///< Entities belonging to the first mesh + + std::shared_ptr _topology1; ///< The second mesh + std::vector + _entities1; ///< Entities belonging to the second mesh +}; + +} // namespace dolfinx::mesh \ No newline at end of file diff --git a/python/demo/demo_hdg.py b/python/demo/demo_hdg.py index 5fec78cc4e3..880dbbd9c63 100644 --- a/python/demo/demo_hdg.py +++ b/python/demo/demo_hdg.py @@ -41,7 +41,7 @@ import ufl from dolfinx import fem, mesh -from dolfinx.cpp.mesh import cell_num_entities +from dolfinx.cpp.mesh import EntityMap, cell_num_entities def par_print(comm, string): @@ -132,12 +132,16 @@ def u_e(x): # Create a cell integral measure over the facet mesh dx_f = ufl.Measure("dx", domain=facet_mesh) -# We write the mixed domain forms as integrals over msh. Hence, we must -# provide a map from facets in msh to cells in facet_mesh. This is the -# 'inverse' of facet_mesh_to_mesh, which we compute as follows: -mesh_to_facet_mesh = np.full(num_facets, -1) -mesh_to_facet_mesh[facet_mesh_to_mesh] = np.arange(len(facet_mesh_to_mesh)) -entity_maps = {facet_mesh: mesh_to_facet_mesh} +# We write the mixed domain forms as integrals over msh. +# We create the relation between the entities that we can pass to the form constructor +entity_maps = [ + EntityMap( + msh.topology._cpp_object, + facet_mesh.topology._cpp_object, + facet_mesh.topology.dim, + facet_mesh_to_mesh, + ) +] # Define forms h = ufl.CellDiameter(msh) @@ -170,6 +174,9 @@ def u_e(x): # Since the boundary condition is enforced in the facet space, we must # use the mesh_to_facet_mesh map to get the corresponding facets in # facet_mesh + +mesh_to_facet_mesh = np.full(num_facets, -1) +mesh_to_facet_mesh[facet_mesh_to_mesh] = np.arange(len(facet_mesh_to_mesh)) facet_mesh_boundary_facets = mesh_to_facet_mesh[msh_boundary_facets] # Get the dofs and apply the boundary condition diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py index 70dd5d52f36..219724fac81 100644 --- a/python/demo/demo_static-condensation.py +++ b/python/demo/demo_static-condensation.py @@ -196,7 +196,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, cu cells = np.arange(msh.topology.index_map(msh.topology.dim).size_local) integrals = {IntegralType.cell: [(-1, tabulate_A.address, cells, np.array([], dtype=np.int8))]} a_cond = Form( - formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, {}, mesh=msh._cpp_object) + formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, [], mesh=msh._cpp_object) ) A_cond = assemble_matrix(a_cond, bcs=[bc]) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 750eaacf245..1086acfcb55 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -202,7 +202,7 @@ def mixed_topology_form( dtype: npt.DTypeLike = default_scalar_type, form_compiler_options: typing.Optional[dict] = None, jit_options: typing.Optional[dict] = None, - entity_maps: typing.Optional[dict[Mesh, np.typing.NDArray[np.int32]]] = None, + entity_maps: typing.Optional[typing.Iterable[_cpp.mesh.EntityMap]] = None, ): """ Create a mixed-topology from from an array of Forms. @@ -220,12 +220,9 @@ def mixed_topology_form( entity_maps: If any trial functions, test functions, or coefficients in the form are not defined over the same mesh as the integration domain, `entity_maps` must be supplied. - For each key (a mesh, different to the integration domain - mesh) a map should be provided relating the entities in the - integration domain mesh to the entities in the key mesh e.g. - for a key-value pair (msh, emap) in `entity_maps`, `emap[i]` - is the entity in `msh` corresponding to entity `i` in the - integration domain mesh. + Each map is a relation between the a subset of entities for + two different mesh topologies. + Returns: Compiled finite element Form. @@ -272,7 +269,7 @@ def mixed_topology_form( [], [], {}, - {}, + [], mesh, ) return Form(f, ufcx_forms, codes, modules) @@ -283,7 +280,7 @@ def form( dtype: npt.DTypeLike = default_scalar_type, form_compiler_options: typing.Optional[dict] = None, jit_options: typing.Optional[dict] = None, - entity_maps: typing.Optional[dict[Mesh, np.typing.NDArray[np.int32]]] = None, + entity_maps: typing.Optional[typing.Iterable[_cpp.mesh.EntityMap]] = None, ): """Create a Form or an array of Forms. @@ -295,12 +292,8 @@ def form( entity_maps: If any trial functions, test functions, or coefficients in the form are not defined over the same mesh as the integration domain, ``entity_maps`` must be supplied. - For each key (a mesh, different to the integration domain - mesh) a map should be provided relating the entities in the - integration domain mesh to the entities in the key mesh e.g. - for a key-value pair ``(msh, emap)`` in ``entity_maps``, - ``emap[i]`` is the entity in ``msh`` corresponding to entity - ``i`` in the integration domain mesh. + Each map is a relation between the a subset of entities for two + different mesh topologies. Returns: Compiled finite element Form. @@ -379,9 +372,9 @@ def _form(form): } if entity_maps is None: - _entity_maps = dict() + _entity_maps = [] else: - _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()} + _entity_maps = entity_maps f = ftype( [module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form))], @@ -400,9 +393,9 @@ def _zero_form(form): assert len(V) > 0 msh = V[0].mesh if entity_maps is None: - _entity_maps = dict() + _entity_maps = [] else: - _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()} + _entity_maps = entity_maps f = ftype( spaces=V, integrals={}, @@ -571,7 +564,7 @@ def create_form( subdomains: dict[IntegralType, list[tuple[int, np.ndarray]]], coefficient_map: dict[ufl.Function, function.Function], constant_map: dict[ufl.Constant, function.Constant], - entity_maps: dict[Mesh, np.typing.NDArray[np.int32]] | None = None, + entity_maps: list[_cpp.mesh.EntityMap] | None = None, ) -> Form: """ Create a Form object from a data-independent compiled form. @@ -592,10 +585,7 @@ def create_form( an array of integers, where the i-th entry is the entity in the key mesh. """ - if entity_maps is None: - _entity_maps = {} - else: - _entity_maps = {m._cpp_object: emap for (m, emap) in entity_maps.items()} + _entity_maps = [] if entity_maps is None else entity_maps _subdomain_data = subdomains.copy() for _, idomain in _subdomain_data.items(): diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 594624f86b5..e00fcd095e1 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -644,9 +645,8 @@ void declare_form(nb::module_& m, std::string type) const std::vector< std::shared_ptr>>& constants, bool needs_permutation_data, - const std::map>, - nb::ndarray, - nb::c_contig>>& entity_maps, + const std::vector>& + entity_maps, std::shared_ptr> mesh) { std::map, @@ -669,15 +669,9 @@ void declare_form(nb::module_& m, std::string type) std::vector(c.data(), c.data() + c.size())}}); } } - - std::map>, - std::span> - _entity_maps; - for (auto& [msh, map] : entity_maps) - _entity_maps.emplace(msh, std::span(map.data(), map.size())); new (fp) dolfinx::fem::Form( spaces, std::move(_integrals), mesh, coefficients, constants, - needs_permutation_data, _entity_maps); + needs_permutation_data, entity_maps); }, nb::arg("spaces"), nb::arg("integrals"), nb::arg("coefficients"), nb::arg("constants"), nb::arg("need_permutation_data"), @@ -696,9 +690,8 @@ void declare_form(nb::module_& m, std::string type) std::vector, nb::c_contig>>>>& subdomains, - const std::map>, - nb::ndarray, - nb::c_contig>>& entity_maps, + const std::vector>& + entity_maps, std::shared_ptr> mesh) { std::map>, - std::span> - _entity_maps; - for (auto& [msh, map] : entity_maps) - _entity_maps.emplace(msh, std::span(map.data(), map.size())); std::vector> ps; for (auto form : forms) ps.push_back(*(reinterpret_cast(form))); new (fp) dolfinx::fem::Form(dolfinx::fem::create_form_factory( - ps, spaces, coefficients, constants, sd, _entity_maps, + ps, spaces, coefficients, constants, sd, entity_maps, mesh)); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), @@ -827,8 +814,7 @@ void declare_form(nb::module_& m, std::string type) std::vector>>>& subdomains, - const std::map>, - nb::ndarray>& + const std::vector>& entity_maps, std::shared_ptr> mesh = nullptr) { @@ -843,14 +829,9 @@ void declare_form(nb::module_& m, std::string type) x.emplace_back(id, std::span(idx.data(), idx.size())); sd.insert({itg, std::move(x)}); } - std::map>, - std::span> - _entity_maps; - for (auto& [msh, map] : entity_maps) - _entity_maps.emplace(msh, std::span(map.data(), map.size())); ufcx_form* p = reinterpret_cast(form); return dolfinx::fem::create_form( - *p, spaces, coefficients, constants, sd, _entity_maps, mesh); + *p, spaces, coefficients, constants, sd, entity_maps, mesh); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), nb::arg("constants"), nb::arg("subdomains"), nb::arg("entity_maps"), diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index a1707028680..1f158d21087 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -1,4 +1,5 @@ -// Copyright (C) 2017-2021 Chris N. Richardson and Garth N. Wells +// Copyright (C) 2017-2025 Chris N. Richardson, Garth N. Wells and Jørgen S. +// Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -12,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -123,6 +125,7 @@ void declare_meshtags(nb::module_& m, std::string type) template void declare_mesh(nb::module_& m, std::string type) { + std::string pyclass_geometry_name = std::string("Geometry_") + type; nb::class_>(m, pyclass_geometry_name.c_str(), "Geometry object") @@ -588,6 +591,43 @@ void mesh(nb::module_& m) m.def("compute_connectivity", &dolfinx::mesh::compute_connectivity, nb::arg("topology"), nb::arg("d0"), nb::arg("d1")); + // dolfinx::mesh::EntityMap class + nb::class_(m, "EntityMap", "EntityMap object") + .def( + "__init__", + [](dolfinx::mesh::EntityMap* self, + std::shared_ptr topology0, + std::shared_ptr topology1, int dim, + nb::ndarray, nb::c_contig> + entities0, + nb::ndarray, nb::c_contig> + entities1) + { + new (self) dolfinx::mesh::EntityMap( + topology0, topology1, dim, + std::vector(entities0.data(), + entities0.data() + entities0.size()), + std::vector(entities1.data(), + entities1.data() + entities1.size())); + }, + nb::arg("topology0"), nb::arg("topology1"), nb::arg("dim"), + nb::arg("entities0"), nb::arg("entities1")) + .def( + "__init__", + [](dolfinx::mesh::EntityMap* self, + std::shared_ptr topology0, + std::shared_ptr topology1, int dim, + nb::ndarray, nb::c_contig> + entities0) + { + new (self) dolfinx::mesh::EntityMap( + topology0, topology1, dim, + std::vector(entities0.data(), + entities0.data() + entities0.size())); + }, + nb::arg("topology0"), nb::arg("topology1"), nb::arg("dim"), + nb::arg("entities0")); + // dolfinx::mesh::Topology class nb::class_(m, "Topology", nb::dynamic_attr(), "Topology object") diff --git a/python/test/unit/fem/test_assemble_mesh_independent_form.py b/python/test/unit/fem/test_assemble_mesh_independent_form.py index 0c7ed671ca3..3a36de9c30c 100644 --- a/python/test/unit/fem/test_assemble_mesh_independent_form.py +++ b/python/test/unit/fem/test_assemble_mesh_independent_form.py @@ -113,10 +113,6 @@ def create_and_integrate(N, compiled_form): submesh, sub_to_parent, _, _ = dolfinx.mesh.create_submesh( mesh, mesh.topology.dim - 1, facets ) - imap = mesh.topology.index_map(mesh.topology.dim - 1) - num_facets = imap.size_local + imap.num_ghosts - parent_to_sub = np.full(num_facets, -1, dtype=np.int32) - parent_to_sub[sub_to_parent] = np.arange(sub_to_parent.size, dtype=np.int32) def g(x): return -3 * x[1] ** 3 + x[0] @@ -131,8 +127,14 @@ def g(x): dolfinx.fem.IntegralType.exterior_facet, mesh.topology, sub_to_parent ) subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(subdomain_id, facet_entities)]} + entity_map = dolfinx.cpp.mesh.EntityMap( + mesh.topology._cpp_object, + submesh.topology._cpp_object, + submesh.topology.dim, + sub_to_parent, + ) form = dolfinx.fem.create_form( - compiled_form, [Vh], mesh, subdomains, {w: wh}, {}, {submesh: parent_to_sub} + compiled_form, [Vh], mesh, subdomains, {w: wh}, {}, [entity_map] ) # Compute exact solution diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index c2cffad73d8..45a332ce7a4 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -13,6 +13,7 @@ import ufl from dolfinx import default_scalar_type, fem, la +from dolfinx.cpp.mesh import EntityMap from dolfinx.fem import compute_integration_domains from dolfinx.mesh import ( CellType, @@ -250,25 +251,21 @@ def coeff_expr(x): c = msh.comm.allreduce(fem.assemble_scalar(M), op=MPI.SUM) # Assemble a mixed-domain form using msh as integration domain. - # Entity maps must map cells in msh (the integration domain mesh, - # defined by the integration measure) to cells in smsh. - cell_imap = msh.topology.index_map(tdim) - num_cells = cell_imap.size_local + cell_imap.num_ghosts - msh_to_smsh = np.full(num_cells, -1) - msh_to_smsh[smsh_to_msh] = np.arange(len(smsh_to_msh)) - entity_maps = {smsh: np.array(msh_to_smsh, dtype=np.int32)} - a0 = fem.form(a_ufl(u, q, f, g, measure_msh), entity_maps=entity_maps) + # We create an entity map that relates the entities in the submesh to those of the parent mesh + entity_map = EntityMap( + msh.topology._cpp_object, smsh.topology._cpp_object, smsh.topology.dim, smsh_to_msh + ) + a0 = fem.form(a_ufl(u, q, f, g, measure_msh), entity_maps=[entity_map]) A0 = fem.assemble_matrix(a0, bcs=[bc]) A0.scatter_reverse() assert np.isclose(A0.squared_norm(), A.squared_norm()) - L0 = fem.form(L_ufl(q, f, g, measure_msh), entity_maps=entity_maps) + L0 = fem.form(L_ufl(q, f, g, measure_msh), entity_maps=[entity_map]) b0 = fem.assemble_vector(L0) fem.apply_lifting(b0.array, [a0], bcs=[[bc]]) b0.scatter_reverse(la.InsertMode.add) assert np.isclose(la.norm(b0), la.norm(b)) - - M0 = fem.form(M_ufl(f, g, measure_msh), entity_maps=entity_maps) + M0 = fem.form(M_ufl(f, g, measure_msh), entity_maps=[entity_map]) c0 = msh.comm.allreduce(fem.assemble_scalar(M0), op=MPI.SUM) assert np.isclose(c0, c) @@ -278,21 +275,18 @@ def coeff_expr(x): # Create the measure (this time defined over the submesh) measure_smsh = create_measure(smsh, integral_type) - # Entity maps must map cells in smsh (the integration domain mesh) to - # cells in msh - entity_maps = {msh: np.array(smsh_to_msh, dtype=np.int32)} - a1 = fem.form(a_ufl(u, q, f, g, measure_smsh), entity_maps=entity_maps) + a1 = fem.form(a_ufl(u, q, f, g, measure_smsh), entity_maps=[entity_map]) A1 = fem.assemble_matrix(a1, bcs=[bc]) A1.scatter_reverse() assert np.isclose(A1.squared_norm(), A.squared_norm()) - L1 = fem.form(L_ufl(q, f, g, measure_smsh), entity_maps=entity_maps) + L1 = fem.form(L_ufl(q, f, g, measure_smsh), entity_maps=[entity_map]) b1 = fem.assemble_vector(L1) fem.apply_lifting(b1.array, [a1], bcs=[[bc]]) b1.scatter_reverse(la.InsertMode.add) assert np.isclose(la.norm(b1), la.norm(b)) - M1 = fem.form(M_ufl(f, g, measure_smsh), entity_maps=entity_maps) + M1 = fem.form(M_ufl(f, g, measure_smsh), entity_maps=[entity_map]) c1 = msh.comm.allreduce(fem.assemble_scalar(M1), op=MPI.SUM) assert np.isclose(c1, c) @@ -360,29 +354,26 @@ def coeff_expr(x): M = fem.form(M_ufl(f, f, ds)) c = msh.comm.allreduce(fem.assemble_scalar(M), op=MPI.SUM) - # Since msh is the integration domain, we must pass entity maps taking - # facets in msh to cells in smsh. This is simply the inverse of smsh_to_msh. - facet_imap = msh.topology.index_map(tdim - 1) - num_facets = facet_imap.size_local + facet_imap.num_ghosts - msh_to_smsh = np.full(num_facets, -1) - msh_to_smsh[smsh_to_msh] = np.arange(len(smsh_to_msh)) - entity_maps = {smsh: msh_to_smsh} + # We create the realation between the submesh and the parent mesh + entity_map = EntityMap( + msh.topology._cpp_object, smsh.topology._cpp_object, smsh.topology.dim, smsh_to_msh + ) # Create forms and compare - a1 = fem.form(a_ufl(u, vbar, f, g, ds), entity_maps=entity_maps) + a1 = fem.form(a_ufl(u, vbar, f, g, ds), entity_maps=[entity_map]) A1 = fem.assemble_matrix(a1, bcs=[bc]) A1.scatter_reverse() assert np.isclose(A.squared_norm(), A1.squared_norm()) - L1 = fem.form(L_ufl(vbar, f, g, ds), entity_maps=entity_maps) + L1 = fem.form(L_ufl(vbar, f, g, ds), entity_maps=[entity_map]) b1 = fem.assemble_vector(L1) fem.apply_lifting(b1.array, [a1], bcs=[[bc]]) b1.scatter_reverse(la.InsertMode.add) assert np.isclose(la.norm(b), la.norm(b1)) - M1 = fem.form(M_ufl(f, g, ds), entity_maps=entity_maps) + M1 = fem.form(M_ufl(f, g, ds), entity_maps=[entity_map]) c1 = msh.comm.allreduce(fem.assemble_scalar(M1), op=MPI.SUM) assert np.isclose(c, c1) @@ -505,23 +496,57 @@ def f_right(x): ) # We create an entity map from the parent mesh to the submesh, where - # the cell on either side of the interface is mapped to the same cell + # the cell on either side of the interface is mapped to the same cell. mesh.topology.create_connectivity(tdim - 1, tdim) f_to_c = mesh.topology.connectivity(tdim - 1, tdim) parent_to_left = np.full(num_cells_local, -1, dtype=np.int32) parent_to_right = np.full(num_cells_local, -1, dtype=np.int32) parent_to_left[left_to_parent] = np.arange(len(left_to_parent)) parent_to_right[right_to_parent] = np.arange(len(right_to_parent)) + left_cells = [] + right_cells = [] + parent_left_cells = [] + parent_right_cells = [] + # FIXME: this would be way easier if we transfer the meshtag to the relevant submesh(es) for tag in [4, 5]: + # Loop through facets of the parent mesh for facet in facet_tag.find(tag): + # Find cells on the interface cells = f_to_c.links(facet) assert len(cells) == 2 + + # Map cells to submesh(es) + # Not every submesh has a cell on the interface left_map = parent_to_left[cells] + if (left_val := max(left_map)) > -1: + left_cells.extend([left_val for _ in left_map]) + parent_left_cells.extend(cells) right_map = parent_to_right[cells] - parent_to_left[cells] = max(left_map) - parent_to_right[cells] = max(right_map) + if (right_val := max(right_map)) > -1: + right_cells.extend([right_val for _ in right_map]) + parent_right_cells.extend(cells) + + # Create entity maps + parent_left_cells = np.asarray(parent_left_cells, dtype=np.int32) + parent_right_cells = np.asarray(parent_right_cells, dtype=np.int32) + left_cells = np.asarray(left_cells, dtype=np.int32) + right_cells = np.asarray(right_cells, dtype=np.int32) + entity_map_left = EntityMap( + mesh.topology._cpp_object, + left_mesh.topology._cpp_object, + tdim, + parent_left_cells, + left_cells, + ) + entity_map_right = EntityMap( + mesh.topology._cpp_object, + right_mesh.topology._cpp_object, + tdim, + parent_right_cells, + right_cells, + ) + entity_maps = [entity_map_left, entity_map_right] - entity_maps = {left_mesh: parent_to_left, right_mesh: parent_to_right} J_compiled = fem.form(J, entity_maps=entity_maps) J_local = fem.assemble_scalar(J_compiled) J_sum = mesh.comm.allreduce(J_local, op=MPI.SUM) @@ -586,17 +611,15 @@ def test_mixed_measures(): v = ufl.TestFunction(V) q = ufl.TestFunction(Q) + entity_maps = [ + EntityMap(msh.topology._cpp_object, smsh.topology._cpp_object, tdim, smsh_to_msh) + ] # First, assemble a block vector using both dx_msh and dx_smsh L = [fem.form(ufl.inner(2.3, v) * dx_msh), fem.form(ufl.inner(1.3, q) * dx_smsh)] b0 = assemble_vector(L, kind=PETSc.Vec.Type.MPI) b0.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # Now, assemble the same vector using only dx_msh - cell_imap = msh.topology.index_map(tdim) - num_cells = cell_imap.size_local + cell_imap.num_ghosts - msh_to_smsh = np.full(num_cells, -1) - msh_to_smsh[smsh_to_msh] = np.arange(len(smsh_to_msh)) - entity_maps = {smsh: msh_to_smsh} L = [ fem.form(ufl.inner(2.3, v) * dx_msh), fem.form(ufl.inner(1.3, q) * dx_msh(1), entity_maps=entity_maps), @@ -641,7 +664,6 @@ def test_interior_facet_codim_1(msh): fdim = tdim - 1 msh.topology.create_connectivity(fdim, tdim) facet_imap = msh.topology.index_map(fdim) - num_facets = facet_imap.size_local + facet_imap.num_ghosts # Mark all local and owned interior facets and "unmark" exterior facets facet_vector = la.vector(facet_imap, 1, dtype=np.int32) @@ -655,9 +677,9 @@ def test_interior_facet_codim_1(msh): submesh, sub_to_parent, _, _ = create_submesh(msh, fdim, interior_facets) # Create inverse map - mesh_to_submesh = np.full(num_facets, -1) - mesh_to_submesh[sub_to_parent] = np.arange(len(sub_to_parent)) - entity_maps = {submesh: mesh_to_submesh} + entity_maps = [ + EntityMap(msh.topology._cpp_object, submesh.topology._cpp_object, fdim, sub_to_parent) + ] def assemble_interior_facet_formulation(formulation, entity_maps): F = fem.form(formulation, entity_maps=entity_maps) diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py index 3ea1b8dbd85..9fec7e6706f 100644 --- a/python/test/unit/fem/test_custom_jit_kernels.py +++ b/python/test/unit/fem/test_custom_jit_kernels.py @@ -102,11 +102,11 @@ def test_numba_assembly(dtype): formtype = form_cpp_class(dtype) a = Form( formtype( - [V._cpp_object, V._cpp_object], integrals, [], [], False, {}, mesh=mesh._cpp_object + [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object ) ) integrals = {IntegralType.cell: [(-1, k1.address, cells, active_coeffs)]} - L = Form(formtype([V._cpp_object], integrals, [], [], False, {}, mesh=mesh._cpp_object)) + L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object)) A = dolfinx.fem.assemble_matrix(a) A.scatter_reverse() @@ -141,7 +141,7 @@ def test_coefficient(dtype): formtype = form_cpp_class(dtype) L = Form( formtype( - [V._cpp_object], integrals, [vals._cpp_object], [], False, {}, mesh=mesh._cpp_object + [V._cpp_object], integrals, [vals._cpp_object], [], False, [], mesh=mesh._cpp_object ) ) @@ -279,14 +279,14 @@ def test_cffi_assembly(): integrals = {IntegralType.cell: [(-1, ptrA, cells, active_coeffs)]} a = Form( _cpp.fem.Form_float64( - [V._cpp_object, V._cpp_object], integrals, [], [], False, {}, mesh=mesh._cpp_object + [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object ) ) ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL")) integrals = {IntegralType.cell: [(-1, ptrL, cells, active_coeffs)]} L = Form( - _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, {}, mesh=mesh._cpp_object) + _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object) ) A = fem.assemble_matrix(a) diff --git a/python/test/unit/fem/test_forms.py b/python/test/unit/fem/test_forms.py index 05f2408dd05..f6b2865984e 100644 --- a/python/test/unit/fem/test_forms.py +++ b/python/test/unit/fem/test_forms.py @@ -94,7 +94,7 @@ def test_incorrect_element(): [], [], {IntegralType.cell: []}, - {}, + [], mesh._cpp_object, ) dolfinx.fem.Form(f, ufcx_form, code) @@ -106,7 +106,7 @@ def test_incorrect_element(): [], [], {IntegralType.cell: []}, - {}, + [], mesh._cpp_object, ) dolfinx.fem.Form(f, ufcx_form, code)