From d12e5123a19eaa9d7ce5a3bf374302f94dcef956 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Thu, 16 Jan 2025 07:53:26 +0000 Subject: [PATCH] More work towards supporting mixed-topology meshes (#3590) * Add very low level example code * Update to latest dolfinx * Add comments * Use vtk permutation * Call ffcx once * Remove unused * Fix some mypy moaning * Exit in parallel * Work on create_mesh [skip ci] * push back spans [skip ci] * add test * Adjust docstring * Enable no-partitioning mode * Use create_mesh * Update to API * Slight tidy up * More dolfinx integration * Simplify and add more docs * Fix for MPI * Start on generalising assembler * Used custom SP for now * Ranges * Check only one type of imap when Topology::index_map called * Fix default integration entities * Combine * Add sparsity build function * Simplify * Tidy * Format * Fix docs * Docs * Add FIXMEs * Add list of eles * Multiple dofmaps * Add new FS contructor * Add getters * Bind new constructor * Start on form * Add func to get cell types * Pass multiple ufcx_forms * Work on demo * Remove commented line * Start on adding cell_type to integral data * More work * More work * More hacking * More work * More work * Multiple kernels * Fix sparsity pattern creation * Add back rest of demo * Tidy * Tidy * Tidy * Tidy * Uncomment code * Tidy * Tidy * Add FIXME * Tidy * Tidy * Uncomment code * Simplify * Uncomment more code * Rename * Simplify * Tidy * Add comment * Get kernel properly * Comments * More work * Check number of elements/dofmaps * Add note * Tidy * Tidy * Tidy * Ruff fix * Add [] * Fix vector vs span * Ruff formatting * Ruff formatting * Doc fix * Add separate function * Add FIXME * mypy * Use std::size_t * Docs * Fix float type * Fix type * Fix test * Updates * Add check * Update cpp/dolfinx/fem/FunctionSpace.h * Style updates * Docs --------- Co-authored-by: Chris Richardson Co-authored-by: Garth N. Wells --- cpp/dolfinx/fem/Form.h | 112 ++++- cpp/dolfinx/fem/FunctionSpace.h | 120 +++-- cpp/dolfinx/fem/assemble_matrix_impl.h | 190 ++++---- cpp/dolfinx/fem/assembler.h | 27 +- cpp/dolfinx/fem/utils.h | 583 +++++++++++++------------ cpp/dolfinx/mesh/Geometry.h | 2 +- cpp/dolfinx/mesh/Topology.cpp | 15 +- cpp/dolfinx/mesh/Topology.h | 6 +- python/demo/demo_mixed-topology.py | 187 ++++++++ python/doc/source/demos.rst | 2 + python/dolfinx/fem/__init__.py | 21 + python/dolfinx/fem/forms.py | 93 +++- python/dolfinx/wrappers/assemble.cpp | 9 +- python/dolfinx/wrappers/fem.cpp | 22 +- python/dolfinx/wrappers/mesh.cpp | 1 + python/test/unit/fem/test_forms.py | 4 +- 16 files changed, 963 insertions(+), 431 deletions(-) create mode 100644 python/demo/demo_mixed-topology.py diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 58f3eef2ebd..c75a2af9e81 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -258,6 +258,35 @@ class Form return _function_spaces; } + /// @brief Get the kernel function for integral `i` on given domain + /// type. + /// @param[in] type Integral type. + /// @param[in] i The subdomain ID. + /// @param[in] kernel_idx Index of the kernel (we may have multiple kernels + /// for a given ID in mixed-topology meshes). + /// @return Function to call for `tabulate_tensor`. + std::function + kernel(IntegralType type, int i, int kernel_idx) const + { + const std::vector>& integrals + = _integrals[static_cast(type)]; + + // Get the range of integrals with a given ID + auto get_id = [](const auto& a) { return a.id; }; + auto start = std::ranges::lower_bound(integrals, i, std::less<>{}, get_id); + auto end = std::ranges::upper_bound(integrals, i, std::less<>{}, get_id); + + // Check that the kernel is valid and return it if so + if (start == integrals.end() or start->id != i + or std::distance(start, end) <= kernel_idx) + { + throw std::runtime_error("No kernel for requested domain index."); + } + + return std::next(start, kernel_idx)->kernel; + } + /// @brief Get the kernel function for integral `i` on given domain /// type. /// @param[in] type Integral type. @@ -267,13 +296,7 @@ class Form const geometry_type*, const int*, const uint8_t*)> kernel(IntegralType type, int i) const { - const auto& integrals = _integrals[static_cast(type)]; - auto it = std::ranges::lower_bound(integrals, i, std::less<>{}, - [](const auto& a) { return a.id; }); - if (it != integrals.end() and it->id == i) - return it->kernel; - else - throw std::runtime_error("No kernel for requested domain index."); + return kernel(type, i, 0); } /// @brief Get types of integrals in the form. @@ -310,8 +333,7 @@ class Form /// @param[in] i Index of the integral. std::vector active_coeffs(IntegralType type, std::size_t i) const { - assert(i < _integrals[static_cast(type)].size()); - return _integrals[static_cast(type)][i].coeffs; + return _integrals[static_cast(type)].at(i).coeffs; } /// @brief Get the IDs for integrals (kernels) for given integral type. @@ -324,9 +346,15 @@ class Form std::vector integral_ids(IntegralType type) const { std::vector ids; - const auto& integrals = _integrals[static_cast(type)]; + const std::vector>& integrals + = _integrals[static_cast(type)]; std::ranges::transform(integrals, std::back_inserter(ids), [](auto& integral) { return integral.id; }); + + // IDs may be repeated in mixed-topology meshes, so remove duplicates + std::sort(ids.begin(), ids.end()); + auto it = std::unique(ids.begin(), ids.end()); + ids.erase(it, ids.end()); return ids; } @@ -344,12 +372,14 @@ class Form /// local_facet_index_1)`. Data is flattened with row-major layout, /// `shape=(num_facets, 4)`. /// - /// @param[in] type Integral domain type. + /// @param[in] type Integral type. /// @param[in] i Integral ID, i.e. (sub)domain index. /// @return List of active entities for the given integral (kernel). std::span domain(IntegralType type, int i) const { - const auto& integrals = _integrals[static_cast(type)]; + // FIXME This should call domain with kernel_idx=0 + const std::vector>& integrals + = _integrals[static_cast(type)]; auto it = std::ranges::lower_bound(integrals, i, std::less<>{}, [](const auto& a) { return a.id; }); if (it != integrals.end() and it->id == i) @@ -358,24 +388,52 @@ class Form throw std::runtime_error("No mesh entities for requested domain index."); } + /// @brief Get the list of mesh entity indices for the ith integral + /// (kernel) of a given type. + /// @param[in] type Integral type. + /// @param[in] i Integral ID, i.e. (sub)domain index. + /// @param[in] kernel_idx Index of the kernel (we may have multiple kernels + /// for a given ID in mixed-topology meshes). + /// @return List of active entities for the given integral (kernel). + std::vector domain(IntegralType type, int i, + int kernel_idx) const + { + const std::vector>& integrals + = _integrals[static_cast(type)]; + auto get_id = [](const auto& a) { return a.id; }; + auto start = std::ranges::lower_bound(integrals, i, std::less<>{}, get_id); + auto end = std::ranges::upper_bound(integrals, i, std::less<>{}, get_id); + + // Check that the kernel is valid and return it if so + if (start == integrals.end() or start->id != i + or std::distance(start, end) <= kernel_idx) + { + throw std::runtime_error("No kernel for requested domain index."); + } + + return std::next(start, kernel_idx)->entities; + } + /// @brief Compute the list of entity indices in `mesh` for the ith /// integral (kernel) of a given type (i.e. cell, exterior facet, or /// interior facet). /// /// @param type Integral type. /// @param i Integral ID, i.e. the (sub)domain index. + /// @param kernel_idx Index of the kernel (we may have multiple + /// kernels for a given ID in mixed-topology meshes). /// @param mesh The mesh the entities are numbered with respect to. /// @return List of active entities in `mesh` for the given integral. - std::vector domain(IntegralType type, int i, + std::vector domain(IntegralType type, int i, int kernel_idx, const mesh::Mesh& mesh) const { // Hack to avoid passing shared pointer to this function std::shared_ptr> msh_ptr( &mesh, [](const mesh::Mesh*) {}); - std::span entities = domain(type, i); + std::vector entities = domain(type, i, kernel_idx); if (msh_ptr == _mesh) - return std::vector(entities.begin(), entities.end()); + return entities; else { std::span entity_map = _entity_maps.at(msh_ptr); @@ -416,6 +474,7 @@ class Form // Get the facet index const std::int32_t facet = c_to_f->links(entities[i])[entities[i + 1]]; + // Add cell and the local facet index mapped_entities.insert(mapped_entities.end(), {entity_map[facet], entities[i + 1]}); @@ -443,9 +502,9 @@ class Form } else if (codim == 1) { - // 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. + // 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. auto c_to_f = _mesh->topology()->connectivity(tdim, tdim - 1); assert(c_to_f); for (std::size_t i = 0; i < entities.size(); i += 2) @@ -453,6 +512,7 @@ class Form // Get the facet index const std::int32_t facet = c_to_f->links(entities[i])[entities[i + 1]]; + // Add cell and the local facet index mapped_entities.insert(mapped_entities.end(), {entity_map[facet], entities[i + 1]}); @@ -468,6 +528,20 @@ class Form } } + /// @brief Compute the list of entity indices in `mesh` for the ith + /// integral (kernel) of a given type (i.e. cell, exterior facet, or + /// interior facet). + /// + /// @param type Integral type. + /// @param i Integral ID, i.e. the (sub)domain index. + /// @param mesh The mesh the entities are numbered with respect to. + /// @return List of active entities in `mesh` for the given integral. + std::vector domain(IntegralType type, int i, + const mesh::Mesh& mesh) const + { + return domain(type, i, 0, mesh); + } + /// @brief Access coefficients. const std::vector< std::shared_ptr>>& @@ -531,5 +605,5 @@ class Form std::map>, std::vector> _entity_maps; -}; // namespace dolfinx::fem +}; } // namespace dolfinx::fem diff --git a/cpp/dolfinx/fem/FunctionSpace.h b/cpp/dolfinx/fem/FunctionSpace.h index e7b82f9c8a3..58518283281 100644 --- a/cpp/dolfinx/fem/FunctionSpace.h +++ b/cpp/dolfinx/fem/FunctionSpace.h @@ -44,12 +44,47 @@ class FunctionSpace FunctionSpace(std::shared_ptr> mesh, std::shared_ptr> element, std::shared_ptr dofmap) - : _mesh(mesh), _element(element), _dofmap(dofmap), + : _mesh(mesh), _dofmaps{dofmap}, _elements{element}, _id(boost::uuids::random_generator()()), _root_space_id(_id) { // Do nothing } + /// @brief Create function space for given mesh, elements and + /// degree-of-freedom maps. + /// @param[in] mesh Mesh that the space is defined on. + /// @param[in] elements Finite elements for the space, one for each cell type. + /// The elements must be ordered to be consistent with + /// mesh::topology::cell_types. + /// @param[in] dofmaps Degree-of-freedom maps for the space, one for each + /// element. The dofmaps must be ordered in the same way as the elements. + FunctionSpace( + std::shared_ptr> mesh, + std::vector>> elements, + std::vector> dofmaps) + : _mesh(mesh), _dofmaps(dofmaps), _elements(elements), + _id(boost::uuids::random_generator()()), _root_space_id(_id) + { + std::vector cell_types = mesh->topology()->cell_types(); + int num_cell_types = cell_types.size(); + if (elements.size() != num_cell_types) + { + throw std::runtime_error( + "Number of elements must match number of cell types"); + } + if (dofmaps.size() != num_cell_types) + { + throw std::runtime_error( + "Number of dofmaps must match number of cell types"); + } + for (std::size_t i = 0; i < num_cell_types; ++i) + { + if (elements.at(i)->cell_type() != cell_types.at(i)) + throw std::runtime_error( + "Element cell types must match mesh cell types"); + } + } + // Copy constructor (deleted) FunctionSpace(const FunctionSpace& V) = delete; @@ -76,19 +111,19 @@ class FunctionSpace FunctionSpace sub(const std::vector& component) const { assert(_mesh); - assert(_element); - assert(_dofmap); + assert(_elements.front()); + assert(_dofmaps.front()); // Check that component is valid if (component.empty()) throw std::runtime_error("Component must be non-empty"); // Extract sub-element - auto element = this->_element->extract_sub_element(component); + auto element = this->_elements.front()->extract_sub_element(component); // Extract sub dofmap - auto dofmap - = std::make_shared(_dofmap->extract_sub_dofmap(component)); + auto dofmap = std::make_shared( + _dofmaps.front()->extract_sub_dofmap(component)); // Create new sub space FunctionSpace sub_space(_mesh, element, dofmap); @@ -137,11 +172,11 @@ class FunctionSpace // Create collapsed DofMap auto [_collapsed_dofmap, collapsed_dofs] - = _dofmap->collapse(_mesh->comm(), *_mesh->topology()); + = _dofmaps.front()->collapse(_mesh->comm(), *_mesh->topology()); auto collapsed_dofmap = std::make_shared(std::move(_collapsed_dofmap)); - return {FunctionSpace(_mesh, _element, collapsed_dofmap), + return {FunctionSpace(_mesh, _elements.front(), collapsed_dofmap), std::move(collapsed_dofs)}; } @@ -154,8 +189,8 @@ class FunctionSpace /// 2-tensor. bool symmetric() const { - if (_element) - return _element->symmetric(); + if (_elements.front()) + return _elements.front()->symmetric(); return false; } @@ -178,8 +213,8 @@ class FunctionSpace "FunctionSpace that is a subspace."); } - assert(_element); - if (_element->is_mixed()) + assert(_elements.front()); + if (_elements.front()->is_mixed()) { throw std::runtime_error( "Cannot tabulate coordinates for a mixed FunctionSpace."); @@ -187,31 +222,32 @@ class FunctionSpace // Geometric dimension assert(_mesh); - assert(_element); + assert(_elements.front()); const std::size_t gdim = _mesh->geometry().dim(); const int tdim = _mesh->topology()->dim(); // Get dofmap local size - assert(_dofmap); - std::shared_ptr index_map = _dofmap->index_map; + assert(_dofmaps.front()); + std::shared_ptr index_map + = _dofmaps.front()->index_map; assert(index_map); - const int index_map_bs = _dofmap->index_map_bs(); - const int dofmap_bs = _dofmap->bs(); + const int index_map_bs = _dofmaps.front()->index_map_bs(); + const int dofmap_bs = _dofmaps.front()->bs(); - const int element_block_size = _element->block_size(); + const int element_block_size = _elements.front()->block_size(); const std::size_t scalar_dofs - = _element->space_dimension() / element_block_size; + = _elements.front()->space_dimension() / element_block_size; const std::int32_t num_dofs = index_map_bs * (index_map->size_local() + index_map->num_ghosts()) / dofmap_bs; // Get the dof coordinates on the reference element - if (!_element->interpolation_ident()) + if (!_elements.front()->interpolation_ident()) { throw std::runtime_error("Cannot evaluate dof coordinates - this element " "does not have pointwise evaluation."); } - const auto [X, Xshape] = _element->interpolation_points(); + const auto [X, Xshape] = _elements.front()->interpolation_points(); // Get coordinate map const CoordinateElement& cmap = _mesh->geometry().cmap(); @@ -245,7 +281,7 @@ class FunctionSpace const int num_cells = map->size_local() + map->num_ghosts(); std::span cell_info; - if (_element->needs_dof_transformations()) + if (_elements.front()->needs_dof_transformations()) { _mesh->topology_mutable()->create_entity_permutations(); cell_info = std::span(_mesh->topology()->get_cell_permutation_info()); @@ -264,7 +300,7 @@ class FunctionSpace // TODO: Check transform // Basis function reference-to-conforming transformation function auto apply_dof_transformation - = _element->template dof_transformation_fn( + = _elements.front()->template dof_transformation_fn( doftransform::standard); for (int c = 0; c < num_cells; ++c) @@ -282,7 +318,7 @@ class FunctionSpace x_b, std::span(cell_info.data(), cell_info.size()), c, x.extent(1)); // Get cell dofmap - auto dofs = _dofmap->cell_dofs(c); + auto dofs = _dofmaps.front()->cell_dofs(c); // Copy dof coordinates into vector if (!transpose) @@ -311,21 +347,49 @@ class FunctionSpace /// The finite element std::shared_ptr> element() const { - return _element; + if (_elements.size() > 1) + { + throw std::runtime_error( + "FunctionSpace has multiple elements, call `elements` instead."); + } + + return elements(0); + } + + /// The finite elements + std::shared_ptr> + elements(int cell_type_idx) const + { + return _elements.at(cell_type_idx); } /// The dofmap - std::shared_ptr dofmap() const { return _dofmap; } + std::shared_ptr dofmap() const + { + if (_dofmaps.size() > 1) + { + throw std::runtime_error( + "FunctionSpace has multiple dofmaps, call `dofmaps` instead."); + } + + return dofmaps(0); + } + + /// The dofmaps + std::shared_ptr dofmaps(int cell_type_idx) const + { + return _dofmaps.at(cell_type_idx); + } private: // The mesh std::shared_ptr> _mesh; // The finite element - std::shared_ptr> _element; + std::vector>> _elements; // The dofmap - std::shared_ptr _dofmap; + std::vector> _dofmaps; // The component w.r.t. to root space std::vector _component; diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 78582892f4f..b1abf6a6a80 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -492,7 +492,7 @@ void assemble_interior_facets( /// are applied. Matrix is not finalised. template void assemble_matrix( - la::MatSet auto mat_set, const Form& a, mdspan2_t x_dofmap, + la::MatSet auto mat_set, const Form& a, std::span> x, std::span constants, const std::map, std::pair, int>>& coefficients, @@ -510,90 +510,120 @@ void assemble_matrix( auto mesh1 = a.function_spaces().at(1)->mesh(); assert(mesh1); - // Get dofmap data - std::shared_ptr dofmap0 - = a.function_spaces().at(0)->dofmap(); - std::shared_ptr dofmap1 - = a.function_spaces().at(1)->dofmap(); - assert(dofmap0); - assert(dofmap1); - auto dofs0 = dofmap0->map(); - const int bs0 = dofmap0->bs(); - auto dofs1 = dofmap1->map(); - const int bs1 = dofmap1->bs(); - - auto element0 = a.function_spaces().at(0)->element(); - assert(element0); - auto element1 = a.function_spaces().at(1)->element(); - assert(element1); - fem::DofTransformKernel auto P0 - = element0->template dof_transformation_fn(doftransform::standard); - fem::DofTransformKernel auto P1T - = element1->template dof_transformation_right_fn( - doftransform::transpose); - - std::span cell_info0; - std::span cell_info1; - if (element0->needs_dof_transformations() - or element1->needs_dof_transformations() or a.needs_facet_permutations()) + // TODO: Mixed topology with exterior and interior facet integrals. + // + // NOTE: Can't just loop over cell types for interior facet integrals + // because we have a kernel per combination of comparable cell types, + // rather than one per cell type. Also, we need the dofmaps for two + // different cell types at the same time. + const int num_cell_types = mesh->topology()->cell_types().size(); + for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx) { - mesh0->topology_mutable()->create_entity_permutations(); - mesh1->topology_mutable()->create_entity_permutations(); - cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); - cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info()); - } + // Geometry dofmap and data + mdspan2_t x_dofmap = mesh->geometry().dofmap(cell_type_idx); + + // Get dofmap data + std::shared_ptr dofmap0 + = a.function_spaces().at(0)->dofmaps(cell_type_idx); + std::shared_ptr dofmap1 + = a.function_spaces().at(1)->dofmaps(cell_type_idx); + assert(dofmap0); + assert(dofmap1); + auto dofs0 = dofmap0->map(); + const int bs0 = dofmap0->bs(); + auto dofs1 = dofmap1->map(); + const int bs1 = dofmap1->bs(); + + auto element0 = a.function_spaces().at(0)->elements(cell_type_idx); + assert(element0); + auto element1 = a.function_spaces().at(1)->elements(cell_type_idx); + assert(element1); + fem::DofTransformKernel auto P0 + = element0->template dof_transformation_fn(doftransform::standard); + fem::DofTransformKernel auto P1T + = element1->template dof_transformation_right_fn( + doftransform::transpose); + + std::span cell_info0; + std::span cell_info1; + if (element0->needs_dof_transformations() + or element1->needs_dof_transformations() + or a.needs_facet_permutations()) + { + mesh0->topology_mutable()->create_entity_permutations(); + mesh1->topology_mutable()->create_entity_permutations(); + cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); + cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info()); + } - for (int i : a.integral_ids(IntegralType::cell)) - { - auto fn = a.kernel(IntegralType::cell, i); - assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); - impl::assemble_cells( - mat_set, x_dofmap, x, a.domain(IntegralType::cell, i), - {dofs0, bs0, a.domain(IntegralType::cell, i, *mesh0)}, P0, - {dofs1, bs1, a.domain(IntegralType::cell, i, *mesh1)}, P1T, bc0, bc1, - fn, coeffs, cstride, constants, cell_info0, cell_info1); - } + for (int i : a.integral_ids(IntegralType::cell)) + { + auto fn = a.kernel(IntegralType::cell, i, cell_type_idx); + assert(fn); + auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); + + impl::assemble_cells( + mat_set, x_dofmap, x, a.domain(IntegralType::cell, i, cell_type_idx), + {dofs0, bs0, a.domain(IntegralType::cell, i, cell_type_idx, *mesh0)}, + P0, + {dofs1, bs1, a.domain(IntegralType::cell, i, cell_type_idx, *mesh1)}, + P1T, bc0, bc1, fn, coeffs, cstride, constants, cell_info0, + cell_info1); + } - std::span perms; - if (a.needs_facet_permutations()) - { - mesh->topology_mutable()->create_entity_permutations(); - perms = std::span(mesh->topology()->get_facet_permutations()); - } + std::span perms; + if (a.needs_facet_permutations()) + { + mesh->topology_mutable()->create_entity_permutations(); + perms = std::span(mesh->topology()->get_facet_permutations()); + } - mesh::CellType cell_type = mesh->topology()->cell_type(); - int num_facets_per_cell - = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1); - for (int i : a.integral_ids(IntegralType::exterior_facet)) - { - auto fn = a.kernel(IntegralType::exterior_facet, i); - assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - impl::assemble_exterior_facets( - mat_set, x_dofmap, x, num_facets_per_cell, - a.domain(IntegralType::exterior_facet, i), - {dofs0, bs0, a.domain(IntegralType::exterior_facet, i, *mesh0)}, P0, - {dofs1, bs1, a.domain(IntegralType::exterior_facet, i, *mesh1)}, P1T, - bc0, bc1, fn, coeffs, cstride, constants, cell_info0, cell_info1, - perms); - } + mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx]; + int num_facets_per_cell + = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1); + for (int i : a.integral_ids(IntegralType::exterior_facet)) + { + if (num_cell_types > 1) + { + throw std::runtime_error("Exterior facet integrals with mixed " + "topology aren't supported yet"); + } - for (int i : a.integral_ids(IntegralType::interior_facet)) - { - const std::vector c_offsets = a.coefficient_offsets(); - auto fn = a.kernel(IntegralType::interior_facet, i); - assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::interior_facet, i}); - impl::assemble_interior_facets( - mat_set, x_dofmap, x, num_facets_per_cell, - a.domain(IntegralType::interior_facet, i), - {*dofmap0, bs0, a.domain(IntegralType::interior_facet, i, *mesh0)}, P0, - {*dofmap1, bs1, a.domain(IntegralType::interior_facet, i, *mesh1)}, P1T, - bc0, bc1, fn, coeffs, cstride, c_offsets, constants, cell_info0, - cell_info1, perms); + auto fn = a.kernel(IntegralType::exterior_facet, i); + assert(fn); + auto& [coeffs, cstride] + = coefficients.at({IntegralType::exterior_facet, i}); + impl::assemble_exterior_facets( + mat_set, x_dofmap, x, num_facets_per_cell, + a.domain(IntegralType::exterior_facet, i), + {dofs0, bs0, a.domain(IntegralType::exterior_facet, i, *mesh0)}, P0, + {dofs1, bs1, a.domain(IntegralType::exterior_facet, i, *mesh1)}, P1T, + bc0, bc1, fn, coeffs, cstride, constants, cell_info0, cell_info1, + perms); + } + + for (int i : a.integral_ids(IntegralType::interior_facet)) + { + if (num_cell_types > 1) + { + throw std::runtime_error("Interior facet integrals with mixed " + "topology aren't supported yet"); + } + + const std::vector c_offsets = a.coefficient_offsets(); + auto fn = a.kernel(IntegralType::interior_facet, i); + assert(fn); + auto& [coeffs, cstride] + = coefficients.at({IntegralType::interior_facet, i}); + impl::assemble_interior_facets( + mat_set, x_dofmap, x, num_facets_per_cell, + a.domain(IntegralType::interior_facet, i), + {*dofmap0, bs0, a.domain(IntegralType::interior_facet, i, *mesh0)}, + P0, + {*dofmap1, bs1, a.domain(IntegralType::interior_facet, i, *mesh1)}, + P1T, bc0, bc1, fn, coeffs, cstride, c_offsets, constants, cell_info0, + cell_info1, perms); + } } } diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index ab41da82792..4a15a3229e9 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -38,10 +38,10 @@ make_coefficients_span(const std::map, { using Key = typename std::remove_reference_t::key_type; std::map, int>> c; - std::ranges::transform( - coeffs, std::inserter(c, c.end()), - [](auto& e) -> typename decltype(c)::value_type - { return {e.first, {e.second.first, e.second.second}}; }); + std::ranges::transform(coeffs, std::inserter(c, c.end()), + [](auto& e) -> typename decltype(c)::value_type { + return {e.first, {e.second.first, e.second.second}}; + }); return c; } @@ -245,16 +245,15 @@ void assemble_matrix( assert(mesh); if constexpr (std::is_same_v>) { - impl::assemble_matrix(mat_add, a, mesh->geometry().dofmap(), - mesh->geometry().x(), constants, coefficients, - dof_marker0, dof_marker1); + impl::assemble_matrix(mat_add, a, mesh->geometry().x(), constants, + coefficients, dof_marker0, dof_marker1); } else { auto x = mesh->geometry().x(); std::vector> _x(x.begin(), x.end()); - impl::assemble_matrix(mat_add, a, mesh->geometry().dofmap(), _x, constants, - coefficients, dof_marker0, dof_marker1); + impl::assemble_matrix(mat_add, a, _x, constants, coefficients, dof_marker0, + dof_marker1); } } @@ -273,10 +272,12 @@ void assemble_matrix( const std::vector>>& bcs) { // Index maps for dof ranges - auto map0 = a.function_spaces().at(0)->dofmap()->index_map; - auto map1 = a.function_spaces().at(1)->dofmap()->index_map; - auto bs0 = a.function_spaces().at(0)->dofmap()->index_map_bs(); - auto bs1 = a.function_spaces().at(1)->dofmap()->index_map_bs(); + // NOTE: For mixed-topology meshes, there will be multiple DOF maps, + // but the index maps are the same. + auto map0 = a.function_spaces().at(0)->dofmaps(0)->index_map; + auto map1 = a.function_spaces().at(1)->dofmaps(0)->index_map; + auto bs0 = a.function_spaces().at(0)->dofmaps(0)->index_map_bs(); + auto bs1 = a.function_spaces().at(1)->dofmaps(0)->index_map_bs(); // Build dof markers std::vector dof_marker0, dof_marker1; diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index d0d5698744a..2dab25772ca 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -150,6 +150,34 @@ extract_function_spaces(const std::vector*>>& a) /// @return The corresponding sparsity pattern template la::SparsityPattern create_sparsity_pattern(const Form& a) +{ + std::shared_ptr mesh = a.mesh(); + assert(mesh); + + // Get index maps and block sizes from the DOF maps. Note that in + // mixed-topology meshes, despite there being multiple DOF maps, the + // index maps and block sizes are the same. + std::array, 2> dofmaps{ + *a.function_spaces().at(0)->dofmaps(0), + *a.function_spaces().at(1)->dofmaps(0)}; + + const std::array index_maps{dofmaps[0].get().index_map, + dofmaps[1].get().index_map}; + const std::array bs + = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()}; + + la::SparsityPattern pattern(mesh->comm(), index_maps, bs); + build_sparsity_pattern(pattern, a); + return pattern; +} + +/// @brief Build a sparsity pattern for a given form. +/// @note The pattern is not finalised, i.e. the caller is responsible +/// for calling SparsityPattern::assemble. +/// @param[in] pattern The sparsity pattern to add to +/// @param[in] a A bilinear form +template +void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) { if (a.rank() != 2) { @@ -157,13 +185,8 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) "Cannot create sparsity pattern. Form is not a bilinear."); } - // Get dof maps and mesh - std::array, 2> dofmaps{ - *a.function_spaces().at(0)->dofmap(), - *a.function_spaces().at(1)->dofmap()}; std::shared_ptr mesh = a.mesh(); assert(mesh); - std::shared_ptr mesh0 = a.function_spaces().at(0)->mesh(); assert(mesh0); std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh(); @@ -181,12 +204,6 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) common::Timer t0("Build sparsity"); - // Get common::IndexMaps for each dimension - const std::array index_maps{dofmaps[0].get().index_map, - dofmaps[1].get().index_map}; - const std::array bs - = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()}; - auto extract_cells = [](std::span facets) { assert(facets.size() % 2 == 0); @@ -197,48 +214,54 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) return cells; }; - // Create and build sparsity pattern - la::SparsityPattern pattern(mesh->comm(), index_maps, bs); - for (auto type : types) + const int num_cell_types = mesh->topology()->cell_types().size(); + for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx) { - std::vector ids = a.integral_ids(type); - switch (type) + std::array, 2> dofmaps{ + *a.function_spaces().at(0)->dofmaps(cell_type_idx), + *a.function_spaces().at(1)->dofmaps(cell_type_idx)}; + + // Create and build sparsity pattern + for (auto type : types) { - case IntegralType::cell: - for (int id : ids) - { - sparsitybuild::cells( - pattern, {a.domain(type, id, *mesh0), a.domain(type, id, *mesh1)}, - {{dofmaps[0], dofmaps[1]}}); - } - break; - case IntegralType::interior_facet: - for (int id : ids) + std::vector ids = a.integral_ids(type); + switch (type) { - sparsitybuild::interior_facets( - pattern, - {extract_cells(a.domain(type, id, *mesh0)), - extract_cells(a.domain(type, id, *mesh1))}, - {{dofmaps[0], dofmaps[1]}}); - } - break; - case IntegralType::exterior_facet: - for (int id : ids) - { - sparsitybuild::cells(pattern, - {extract_cells(a.domain(type, id, *mesh0)), - extract_cells(a.domain(type, id, *mesh1))}, - {{dofmaps[0], dofmaps[1]}}); + case IntegralType::cell: + for (int id : ids) + { + sparsitybuild::cells(pattern, + {a.domain(type, id, cell_type_idx, *mesh0), + a.domain(type, id, cell_type_idx, *mesh1)}, + {{dofmaps[0], dofmaps[1]}}); + } + break; + case IntegralType::interior_facet: + for (int id : ids) + { + sparsitybuild::interior_facets( + pattern, + {extract_cells(a.domain(type, id, *mesh0)), + extract_cells(a.domain(type, id, *mesh1))}, + {{dofmaps[0], dofmaps[1]}}); + } + break; + case IntegralType::exterior_facet: + for (int id : ids) + { + sparsitybuild::cells(pattern, + {extract_cells(a.domain(type, id, *mesh0)), + extract_cells(a.domain(type, id, *mesh1))}, + {{dofmaps[0], dofmaps[1]}}); + } + break; + default: + throw std::runtime_error("Unsupported integral type"); } - break; - default: - throw std::runtime_error("Unsupported integral type"); } } t0.stop(); - - return pattern; } /// Create an ElementDofLayout from a FiniteElement @@ -323,7 +346,7 @@ std::vector get_constant_names(const ufcx_form& ufcx_form); /// Use fem::create_form to create a fem::Form with coefficients and /// constants associated with the name/string. /// -/// @param[in] ufcx_form The UFCx form. +/// @param[in] ufcx_forms A list of UFCx forms, one for each cell type. /// @param[in] spaces Vector of function spaces. The number of spaces is /// equal to the rank of the form. /// @param[in] coefficients Coefficient fields in the form. @@ -337,7 +360,7 @@ std::vector get_constant_names(const ufcx_form& ufcx_form); /// @pre Each value in `subdomains` must be sorted by domain id. template > Form create_form_factory( - const ufcx_form& ufcx_form, + const std::vector>& ufcx_forms, const std::vector>>& spaces, const std::vector>>& coefficients, const std::vector>>& constants, @@ -349,29 +372,38 @@ Form create_form_factory( std::span>& entity_maps, std::shared_ptr> mesh = nullptr) { - if (ufcx_form.rank != (int)spaces.size()) - throw std::runtime_error("Wrong number of argument spaces for Form."); - if (ufcx_form.num_coefficients != (int)coefficients.size()) + for (const ufcx_form& ufcx_form : ufcx_forms) { - throw std::runtime_error( - "Mismatch between number of expected and provided Form coefficients."); - } - if (ufcx_form.num_constants != (int)constants.size()) - { - throw std::runtime_error( - "Mismatch between number of expected and provided Form constants."); + if (ufcx_form.rank != (int)spaces.size()) + throw std::runtime_error("Wrong number of argument spaces for Form."); + if (ufcx_form.num_coefficients != (int)coefficients.size()) + { + throw std::runtime_error("Mismatch between number of expected and " + "provided Form coefficients."); + } + if (ufcx_form.num_constants != (int)constants.size()) + { + throw std::runtime_error( + "Mismatch between number of expected and provided Form constants."); + } } // Check argument function spaces - for (std::size_t i = 0; i < spaces.size(); ++i) + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { - assert(spaces[i]->element()); - if (auto element_hash = ufcx_form.finite_element_hashes[i]; - element_hash != 0 - and element_hash != spaces[i]->element()->basix_element().hash()) + for (std::size_t i = 0; i < spaces.size(); ++i) { - throw std::runtime_error("Cannot create form. Elements are different to " - "those used to compile the form."); + assert(spaces[i]->elements(form_idx)); + if (auto element_hash + = ufcx_forms[form_idx].get().finite_element_hashes[i]; + element_hash != 0 + and element_hash + != spaces[i]->elements(form_idx)->basix_element().hash()) + { + throw std::runtime_error( + "Cannot create form. Elements are different to " + "those used to compile the form."); + } } } @@ -391,7 +423,10 @@ Form create_form_factory( assert(topology); const int tdim = topology->dim(); - const int* integral_offsets = ufcx_form.form_integral_offsets; + // NOTE: This assumes all forms in mixed-topology meshes have the same + // integral offsets. Since the UFL forms for each type of cell should be + // the same, I think this assumption is OK. + const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets; std::vector num_integrals_type(3); for (int i = 0; i < 3; ++i) num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i]; @@ -413,267 +448,277 @@ Form create_form_factory( // Attach cell kernels bool needs_facet_permutations = false; - std::vector default_cells; { - std::span ids(ufcx_form.form_integral_ids + std::vector default_cells; + std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[cell], num_integrals_type[cell]); auto itg = integrals.insert({IntegralType::cell, {}}); auto sd = subdomains.find(IntegralType::cell); - for (int i = 0; i < num_integrals_type[cell]; ++i) + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { - const int id = ids[i]; - ufcx_integral* integral - = ufcx_form.form_integrals[integral_offsets[cell] + i]; - assert(integral); - - // Build list of active coefficients - std::vector active_coeffs; - for (int j = 0; j < ufcx_form.num_coefficients; ++j) + const ufcx_form& ufcx_form = ufcx_forms[form_idx]; + for (int i = 0; i < num_integrals_type[cell]; ++i) { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); - } + const int id = ids[i]; + ufcx_integral* integral + = ufcx_form.form_integrals[integral_offsets[cell] + i]; + assert(integral); + + // Build list of active coefficients + std::vector active_coeffs; + for (int j = 0; j < ufcx_form.num_coefficients; ++j) + { + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); + } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; + kern_t k = nullptr; + if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float32; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex64); - } + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex64); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; + else if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float64; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex128); - } + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex128); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - if (!k) - { - throw std::runtime_error( - "UFCx kernel function is NULL. Check requested types."); - } + if (!k) + { + throw std::runtime_error( + "UFCx kernel function is NULL. Check requested types."); + } - // Build list of entities to assemble over - if (id == -1) - { - // Default kernel, operates on all (owned) cells - assert(topology->index_map(tdim)); - default_cells.resize(topology->index_map(tdim)->size_local(), 0); - std::iota(default_cells.begin(), default_cells.end(), 0); - itg.first->second.emplace_back(id, k, default_cells, active_coeffs); - } - else if (sd != subdomains.end()) - { - // NOTE: This requires that pairs are sorted - auto it - = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](const auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) - itg.first->second.emplace_back(id, k, it->second, active_coeffs); - } + // Build list of entities to assemble over + if (id == -1) + { + // Default kernel, operates on all (owned) cells + assert(topology->index_maps(tdim).at(form_idx)); + default_cells.resize( + topology->index_maps(tdim).at(form_idx)->size_local(), 0); + std::iota(default_cells.begin(), default_cells.end(), 0); + itg.first->second.emplace_back(id, k, default_cells, active_coeffs); + } + else if (sd != subdomains.end()) + { + // NOTE: This requires that pairs are sorted + auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + itg.first->second.emplace_back(id, k, it->second, active_coeffs); + } - if (integral->needs_facet_permutations) - needs_facet_permutations = true; + if (integral->needs_facet_permutations) + needs_facet_permutations = true; + } } } // Attach exterior facet kernels std::vector default_facets_ext; { - std::span ids(ufcx_form.form_integral_ids + std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[exterior_facet], num_integrals_type[exterior_facet]); auto itg = integrals.insert({IntegralType::exterior_facet, {}}); auto sd = subdomains.find(IntegralType::exterior_facet); - for (int i = 0; i < num_integrals_type[exterior_facet]; ++i) + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { - const int id = ids[i]; - ufcx_integral* integral - = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i]; - assert(integral); - std::vector active_coeffs; - for (int j = 0; j < ufcx_form.num_coefficients; ++j) + const ufcx_form& ufcx_form = ufcx_forms[form_idx]; + for (int i = 0; i < num_integrals_type[exterior_facet]; ++i) { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); - } + const int id = ids[i]; + ufcx_integral* integral + = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i]; + assert(integral); + std::vector active_coeffs; + for (int j = 0; j < ufcx_form.num_coefficients; ++j) + { + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); + } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; + kern_t k = nullptr; + if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float32; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex64); - } + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex64); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; + else if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float64; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex128); - } + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex128); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - assert(k); - - // Build list of entities to assembler over - const std::vector bfacets = mesh::exterior_facet_indices(*topology); - auto f_to_c = topology->connectivity(tdim - 1, tdim); - assert(f_to_c); - auto c_to_f = topology->connectivity(tdim, tdim - 1); - assert(c_to_f); - if (id == -1) - { - // Default kernel, operates on all (owned) exterior facets - default_facets_ext.reserve(2 * bfacets.size()); - for (std::int32_t f : bfacets) + assert(k); + + // Build list of entities to assembler over + const std::vector bfacets = mesh::exterior_facet_indices(*topology); + auto f_to_c = topology->connectivity(tdim - 1, tdim); + assert(f_to_c); + auto c_to_f = topology->connectivity(tdim, tdim - 1); + assert(c_to_f); + if (id == -1) + { + // Default kernel, operates on all (owned) exterior facets + default_facets_ext.reserve(2 * bfacets.size()); + for (std::int32_t f : bfacets) + { + // There will only be one pair for an exterior facet integral + auto pair + = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); + default_facets_ext.insert(default_facets_ext.end(), pair.begin(), + pair.end()); + } + itg.first->second.emplace_back(id, k, default_facets_ext, + active_coeffs); + } + else if (sd != subdomains.end()) { - // There will only be one pair for an exterior facet integral - auto pair - = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); - default_facets_ext.insert(default_facets_ext.end(), pair.begin(), - pair.end()); + // NOTE: This requires that pairs are sorted + auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + itg.first->second.emplace_back(id, k, it->second, active_coeffs); } - itg.first->second.emplace_back(id, k, default_facets_ext, - active_coeffs); - } - else if (sd != subdomains.end()) - { - // NOTE: This requires that pairs are sorted - auto it - = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](const auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) - itg.first->second.emplace_back(id, k, it->second, active_coeffs); - } - if (integral->needs_facet_permutations) - needs_facet_permutations = true; + if (integral->needs_facet_permutations) + needs_facet_permutations = true; + } } } // Attach interior facet kernels std::vector default_facets_int; { - std::span ids(ufcx_form.form_integral_ids + std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[interior_facet], num_integrals_type[interior_facet]); auto itg = integrals.insert({IntegralType::interior_facet, {}}); auto sd = subdomains.find(IntegralType::interior_facet); - - // Create indicator for interprocess facets - std::vector interprocess_marker; - if (num_integrals_type[interior_facet] > 0) + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { - assert(topology->index_map(tdim - 1)); - const std::vector& interprocess_facets - = topology->interprocess_facets(); - std::int32_t num_facets = topology->index_map(tdim - 1)->size_local() - + topology->index_map(tdim - 1)->num_ghosts(); - interprocess_marker.resize(num_facets, 0); - std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f) - { interprocess_marker[f] = 1; }); - } - - for (int i = 0; i < num_integrals_type[interior_facet]; ++i) - { - const int id = ids[i]; - ufcx_integral* integral - = ufcx_form.form_integrals[integral_offsets[interior_facet] + i]; - assert(integral); - std::vector active_coeffs; - for (int j = 0; j < ufcx_form.num_coefficients; ++j) + const ufcx_form& ufcx_form = ufcx_forms[form_idx]; + // Create indicator for interprocess facets + std::vector interprocess_marker; + if (num_integrals_type[interior_facet] > 0) { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); + assert(topology->index_map(tdim - 1)); + const std::vector& interprocess_facets + = topology->interprocess_facets(); + std::int32_t num_facets = topology->index_map(tdim - 1)->size_local() + + topology->index_map(tdim - 1)->num_ghosts(); + interprocess_marker.resize(num_facets, 0); + std::ranges::for_each(interprocess_facets, + [&interprocess_marker](auto f) + { interprocess_marker[f] = 1; }); } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) + for (int i = 0; i < num_integrals_type[interior_facet]; ++i) { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex64); - } + const int id = ids[i]; + ufcx_integral* integral + = ufcx_form.form_integrals[integral_offsets[interior_facet] + i]; + assert(integral); + std::vector active_coeffs; + for (int j = 0; j < ufcx_form.num_coefficients; ++j) + { + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); + } + + kern_t k = nullptr; + if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float32; +#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex64); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; + else if constexpr (std::is_same_v) + k = integral->tabulate_tensor_float64; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast::value_type*, const int*, - const unsigned char*)>(integral->tabulate_tensor_complex128); - } + else if constexpr (std::is_same_v>) + { + k = reinterpret_cast::value_type*, const int*, + const unsigned char*)>(integral->tabulate_tensor_complex128); + } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - assert(k); - - // Build list of entities to assembler over - auto f_to_c = topology->connectivity(tdim - 1, tdim); - assert(f_to_c); - auto c_to_f = topology->connectivity(tdim, tdim - 1); - assert(c_to_f); - if (id == -1) - { - // Default kernel, operates on all (owned) interior facets - assert(topology->index_map(tdim - 1)); - std::int32_t num_facets = topology->index_map(tdim - 1)->size_local(); - default_facets_int.reserve(4 * num_facets); - for (std::int32_t f = 0; f < num_facets; ++f) + assert(k); + + // Build list of entities to assembler over + auto f_to_c = topology->connectivity(tdim - 1, tdim); + assert(f_to_c); + auto c_to_f = topology->connectivity(tdim, tdim - 1); + assert(c_to_f); + if (id == -1) { - if (f_to_c->num_links(f) == 2) + // Default kernel, operates on all (owned) interior facets + assert(topology->index_map(tdim - 1)); + std::int32_t num_facets = topology->index_map(tdim - 1)->size_local(); + default_facets_int.reserve(4 * num_facets); + for (std::int32_t f = 0; f < num_facets; ++f) { - auto pairs - = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f); - default_facets_int.insert(default_facets_int.end(), pairs.begin(), - pairs.end()); - } - else if (interprocess_marker[f]) - { - throw std::runtime_error( - "Cannot compute interior facet integral over interprocess " - "facet. Please use ghost mode shared facet when creating the " - "mesh"); + if (f_to_c->num_links(f) == 2) + { + auto pairs + = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f); + default_facets_int.insert(default_facets_int.end(), pairs.begin(), + pairs.end()); + } + else if (interprocess_marker[f]) + { + throw std::runtime_error( + "Cannot compute interior facet integral over interprocess " + "facet. Please use ghost mode shared facet when creating the " + "mesh"); + } } + itg.first->second.emplace_back(id, k, default_facets_int, + active_coeffs); + } + else if (sd != subdomains.end()) + { + auto it = std::ranges::lower_bound(sd->second, id, std::less{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + itg.first->second.emplace_back(id, k, it->second, active_coeffs); } - itg.first->second.emplace_back(id, k, default_facets_int, - active_coeffs); - } - else if (sd != subdomains.end()) - { - auto it - = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](const auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) - itg.first->second.emplace_back(id, k, it->second, active_coeffs); - } - if (integral->needs_facet_permutations) - needs_facet_permutations = true; + if (integral->needs_facet_permutations) + needs_facet_permutations = true; + } } } @@ -744,7 +789,7 @@ Form create_form( throw std::runtime_error("Form constant \"" + name + "\" not provided."); } - return create_form_factory(ufcx_form, spaces, coeff_map, const_map, + return create_form_factory({ufcx_form}, spaces, coeff_map, const_map, subdomains, entity_maps, mesh); } diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index b5bffbc0dd5..ba628ad2fa3 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -269,7 +269,7 @@ create_geometry( // Build 'geometry' dofmap on the topology auto [_dof_index_map, bs, dofmaps] - = fem::build_dofmap_data(topology.index_map(topology.dim())->comm(), + = fem::build_dofmap_data(topology.index_maps(topology.dim())[0]->comm(), topology, dof_layouts, reorder_fn); auto dof_index_map = std::make_shared(std::move(_dof_index_map)); diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index fa7f27fefb8..d11378b71f8 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -768,7 +768,17 @@ const std::vector& Topology::entity_types(int dim) const //----------------------------------------------------------------------------- mesh::CellType Topology::cell_type() const { - return _entity_types.back().at(0); + std::vector cell_types = entity_types(dim()); + if (cell_types.size() > 1) + throw std::runtime_error( + "Multiple cell types of this dimension. Call cell_types " + "instead."); + return cell_types.front(); +} +//----------------------------------------------------------------------------- +std::vector Topology::cell_types() const +{ + return entity_types(dim()); } //----------------------------------------------------------------------------- std::vector> @@ -786,6 +796,9 @@ Topology::index_maps(int dim) const //----------------------------------------------------------------------------- std::shared_ptr Topology::index_map(int dim) const { + if (_entity_types[dim].size() > 1) + throw std::runtime_error( + "Multiple index maps of this dimension. Call index_maps instead."); return this->index_maps(dim).at(0); } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h index 7043bfa5e66..4966add5668 100644 --- a/cpp/dolfinx/mesh/Topology.h +++ b/cpp/dolfinx/mesh/Topology.h @@ -164,6 +164,10 @@ class Topology /// computed const std::vector& get_facet_permutations() const; + /// @brief Get the types of cells in the topology + /// @return The cell types + std::vector cell_types() const; + /// @brief List of inter-process facets of a given type. /// /// "Inter-process" facets are facets that are connected (1) to a cell @@ -194,7 +198,7 @@ class Topology /// existed. bool create_entities(int dim); - /// @brief Create connectivity between given pair of dimensions, `d0 + /// @brief Create connectivity between given pair of dimensions, `d0 /// -> d1`. /// @param[in] d0 Topological dimension. /// @param[in] d1 Topological dimension. diff --git a/python/demo/demo_mixed-topology.py b/python/demo/demo_mixed-topology.py new file mode 100644 index 00000000000..88b935bae36 --- /dev/null +++ b/python/demo/demo_mixed-topology.py @@ -0,0 +1,187 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.13.6 +# --- + +# # Poisson equation +# +# This demo illustrates how to solve a simple Helmholtz problem on a +# mixed-topology mesh. +# +# NOTE: Mixed-topology meshes are a work in progress and are not yet fully +# supported in DOLFINx. + +from mpi4py import MPI + +import numpy as np +from scipy.sparse.linalg import spsolve + +import basix +import dolfinx.cpp as _cpp +import ufl +from dolfinx.cpp.mesh import GhostMode, create_cell_partitioner, create_mesh +from dolfinx.fem import ( + FunctionSpace, + assemble_matrix, + coordinate_element, + mixed_topology_form, +) +from dolfinx.io.utils import cell_perm_vtk +from dolfinx.mesh import CellType, Mesh + +if MPI.COMM_WORLD.size > 1: + print("Not yet running in parallel") + exit(0) + + +# Create a mixed-topology mesh +nx = 16 +ny = 16 +nz = 16 +n_cells = nx * ny * nz + +cells: list = [[], []] +orig_idx: list = [[], []] +geom = [] + +if MPI.COMM_WORLD.rank == 0: + idx = 0 + for i in range(n_cells): + iz = i // (nx * ny) + j = i % (nx * ny) + iy = j // nx + ix = j % nx + + v0 = (iz * (ny + 1) + iy) * (nx + 1) + ix + v1 = v0 + 1 + v2 = v0 + (nx + 1) + v3 = v1 + (nx + 1) + v4 = v0 + (nx + 1) * (ny + 1) + v5 = v1 + (nx + 1) * (ny + 1) + v6 = v2 + (nx + 1) * (ny + 1) + v7 = v3 + (nx + 1) * (ny + 1) + if ix < nx / 2: + cells[0] += [v0, v1, v2, v3, v4, v5, v6, v7] + orig_idx[0] += [idx] + idx += 1 + else: + cells[1] += [v0, v1, v2, v4, v5, v6] + orig_idx[1] += [idx] + idx += 1 + cells[1] += [v1, v2, v3, v5, v6, v7] + orig_idx[1] += [idx] + idx += 1 + + n_points = (nx + 1) * (ny + 1) * (nz + 1) + sqxy = (nx + 1) * (ny + 1) + for v in range(n_points): + iz = v // sqxy + p = v % sqxy + iy = p // (nx + 1) + ix = p % (nx + 1) + geom += [[ix / nx, iy / ny, iz / nz]] + +cells_np = [np.array(c) for c in cells] +geomx = np.array(geom, dtype=np.float64) +hexahedron = coordinate_element(CellType.hexahedron, 1) +prism = coordinate_element(CellType.prism, 1) + +part = create_cell_partitioner(GhostMode.none) +mesh = create_mesh( + MPI.COMM_WORLD, cells_np, [hexahedron._cpp_object, prism._cpp_object], geomx, part +) + +# Create elements and dofmaps for each cell type +elements = [ + basix.create_element(basix.ElementFamily.P, basix.CellType.hexahedron, 1), + basix.create_element(basix.ElementFamily.P, basix.CellType.prism, 1), +] +elements_cpp = [_cpp.fem.FiniteElement_float64(e._e, None, True) for e in elements] +# NOTE: Both dofmaps have the same IndexMap, but different cell_dofs +dofmaps = _cpp.fem.create_dofmaps(mesh.comm, mesh.topology, elements_cpp) + +# Create C++ function space +V_cpp = _cpp.fem.FunctionSpace_float64(mesh, elements_cpp, dofmaps) + +# Create forms for each cell type. +# FIXME This hack is required at the moment because UFL does not yet know about +# mixed topology meshes. +a = [] +for i, cell_name in enumerate(["hexahedron", "prism"]): + print(f"Creating form for {cell_name}") + element = basix.ufl.wrap_element(elements[i]) + domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_name, 1, shape=(3,))) + V = FunctionSpace(Mesh(mesh, domain), element, V_cpp) + u, v = ufl.TrialFunction(V), ufl.TestFunction(V) + k = 12.0 + a += [(ufl.inner(ufl.grad(u), ufl.grad(v)) - k**2 * u * v) * ufl.dx] + +# Compile the form +# FIXME: For the time being, since UFL doesn't understand mixed topology meshes, +# we have to call mixed_topology_form instead of form. +a_form = mixed_topology_form(a, dtype=np.float64) + +# Assemble the matrix +A = assemble_matrix(a_form) + +# Solve +A_scipy = A.to_scipy() +b_scipy = np.ones(A_scipy.shape[1]) + +x = spsolve(A_scipy, b_scipy) +print(f"Solution vector norm {np.linalg.norm(x)}") + +# I/O +# Save to XDMF +xdmf = """ + + + + + +""" + +perm = [cell_perm_vtk(CellType.hexahedron, 8), cell_perm_vtk(CellType.prism, 6)] +topologies = ["Hexahedron", "Wedge"] + +for j in range(2): + vtk_topology = [] + geom_dm = mesh.geometry.dofmaps(j) + for c in geom_dm: + vtk_topology += list(c[perm[j]]) + topology_type = topologies[j] + + xdmf += f""" + + + + {" ".join(str(val) for val in vtk_topology)} + + + + + {" ".join(str(val) for val in mesh.geometry.x.flatten())} + + + + + {" ".join(str(val) for val in x)} + + + """ + +xdmf += """ + + + +""" + +fd = open("mixed-mesh.xdmf", "w") +fd.write(xdmf) +fd.close() diff --git a/python/doc/source/demos.rst b/python/doc/source/demos.rst index bbbcb889137..ec6b921beb5 100644 --- a/python/doc/source/demos.rst +++ b/python/doc/source/demos.rst @@ -37,6 +37,7 @@ PDEs (advanced) demos/demo_poisson_matrix_free.md demos/demo_pyamg.md demos/demo_hdg.md + demos/demo_mixed-topology.md Nonlinear problems @@ -114,3 +115,4 @@ List of all demos demos/demo_mixed-poisson.md demos/demo_pyamg.md demos/demo_hdg.md + demos/demo_mixed-topology.md diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index 8430eee2b9a..ef528c317f4 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -9,12 +9,14 @@ import numpy.typing as npt from dolfinx.cpp.fem import IntegralType, transpose_dofmap +from dolfinx.cpp.fem import build_sparsity_pattern as _build_sparsity_pattern from dolfinx.cpp.fem import compute_integration_domains as _compute_integration_domains from dolfinx.cpp.fem import create_interpolation_data as _create_interpolation_data from dolfinx.cpp.fem import create_sparsity_pattern as _create_sparsity_pattern from dolfinx.cpp.fem import discrete_curl as _discrete_curl from dolfinx.cpp.fem import discrete_gradient as _discrete_gradient from dolfinx.cpp.fem import interpolation_matrix as _interpolation_matrix +from dolfinx.cpp.la import SparsityPattern from dolfinx.cpp.mesh import Topology from dolfinx.fem.assemble import ( apply_lifting, @@ -41,6 +43,7 @@ extract_function_spaces, form, form_cpp_class, + mixed_topology_form, ) from dolfinx.fem.function import ( Constant, @@ -70,6 +73,23 @@ def create_sparsity_pattern(a: Form): return _create_sparsity_pattern(a._cpp_object) +def build_sparsity_pattern(pattern: SparsityPattern, a: Form): + """Build a sparsity pattern from a bilinear form. + + Args: + pattern: The sparsity pattern to add to + a: Bilinear form to build a sparsity pattern for. + + Returns: + Sparsity pattern for the form ``a``. + + Note: + The pattern is not finalised, i.e. the caller is responsible for + calling ``assemble`` on the sparsity pattern. + """ + return _build_sparsity_pattern(pattern, a._cpp_object) + + def create_interpolation_data( V_to: FunctionSpace, V_from: FunctionSpace, @@ -219,6 +239,7 @@ def compute_integration_domains( "functionspace", "locate_dofs_geometrical", "locate_dofs_topological", + "mixed_topology_form", "set_bc", "transpose_dofmap", ] diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 335ef212594..b3a0ee7e914 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -37,7 +37,7 @@ class Form: _cpp.fem.Form_float32, _cpp.fem.Form_float64, ] - _code: typing.Optional[str] + _code: typing.Optional[typing.Union[str, list[str]]] def __init__( self, @@ -48,8 +48,8 @@ def __init__( _cpp.fem.Form_float64, ], ufcx_form=None, - code: typing.Optional[str] = None, - module: typing.Optional[types.ModuleType] = None, + code: typing.Optional[typing.Union[str, list[str]]] = None, + module: typing.Optional[typing.Union[types.ModuleType, list[types.ModuleType]]] = None, ): """A finite element form. @@ -76,12 +76,12 @@ def ufcx_form(self): return self._ufcx_form @property - def code(self) -> typing.Union[str, None]: + def code(self) -> typing.Union[str, list[str], None]: """C code strings.""" return self._code @property - def module(self) -> typing.Union[types.ModuleType, None]: + def module(self) -> typing.Union[types.ModuleType, list[types.ModuleType], None]: """The CFFI module""" return self._module @@ -197,6 +197,87 @@ def form_cpp_class( } +def mixed_topology_form( + forms: typing.Iterable[ufl.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, +): + """ + Create a mixed-topology from from an array of Forms. + + # FIXME: This function is a temporary hack for mixed-topology meshes. It is needed + # because UFL does not know about mixed-topology meshes, so we need + # to pass a list of forms for each cell type. + + Args: + form: A list of UFL forms. Each form should be the same, just + defined on different cell types. + dtype: Scalar type to use for the compiled form. + form_compiler_options: See :func:`ffcx_jit ` + jit_options: See :func:`ffcx_jit `. + 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. + + Returns: + Compiled finite element Form. + """ + + if form_compiler_options is None: + form_compiler_options = dict() + + form_compiler_options["scalar_type"] = dtype + ftype = form_cpp_class(dtype) + + # Extract subdomain data from UFL form + sd = next(iter(forms)).subdomain_data() + (domain,) = list(sd.keys()) # Assuming single domain + + # Check that subdomain data for each integral type is the same + for data in sd.get(domain).values(): + assert all([d is data[0] for d in data if d is not None]) + + mesh = domain.ufl_cargo() + + ufcx_forms = [] + modules = [] + codes = [] + for form in forms: + ufcx_form, module, code = jit.ffcx_jit( + mesh.comm, + form, + form_compiler_options=form_compiler_options, + jit_options=jit_options, + ) + ufcx_forms.append(ufcx_form) + modules.append(module) + codes.append(code) + + # In a mixed-topology mesh, each form has the same C++ function space, + # so we can extract it from any of them + V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()] + + # TODO coeffs, constants, subdomains, entity_maps + f = ftype( + [module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)) for ufcx_form in ufcx_forms], + V, + [], + [], + {}, + {}, + mesh, + ) + return Form(f, ufcx_forms, codes, modules) + + def form( form: typing.Union[ufl.Form, typing.Iterable[ufl.Form]], dtype: npt.DTypeLike = default_scalar_type, @@ -303,7 +384,7 @@ def _form(form): _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()} f = ftype( - module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)), + [module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form))], V, coeffs, constants, diff --git a/python/dolfinx/wrappers/assemble.cpp b/python/dolfinx/wrappers/assemble.cpp index 107a5a24ea7..c897d383a77 100644 --- a/python/dolfinx/wrappers/assemble.cpp +++ b/python/dolfinx/wrappers/assemble.cpp @@ -209,8 +209,7 @@ void declare_assembly_functions(nb::module_& m) nb::arg("form"), "Pack coefficients for a Form."); m.def( "pack_constants", - [](const dolfinx::fem::Form& form) - { + [](const dolfinx::fem::Form& form) { return dolfinx_wrappers::as_nbarray(dolfinx::fem::pack_constants(form)); }, nb::arg("form"), "Pack constants for a Form."); @@ -272,9 +271,11 @@ void declare_assembly_functions(nb::module_& m) _bcs.push_back(*bc); } + // Get index map block size. Note that mixed-topology meshes + // will have multiple DOF maps, but the block sizes are the same. const std::array data_bs - = {a.function_spaces().at(0)->dofmap()->index_map_bs(), - a.function_spaces().at(1)->dofmap()->index_map_bs()}; + = {a.function_spaces().at(0)->dofmaps(0)->index_map_bs(), + a.function_spaces().at(1)->dofmaps(0)->index_map_bs()}; if (data_bs[0] != data_bs[1]) { diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index a80f0786410..b0c09ad05a5 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -74,6 +74,10 @@ void declare_function_space(nb::module_& m, std::string type) std::shared_ptr>, std::shared_ptr>(), nb::arg("mesh"), nb::arg("element"), nb::arg("dofmap")) + .def(nb::init>, + std::vector>>, + std::vector>>(), + nb::arg("mesh"), nb::arg("elements"), nb::arg("dofmaps")) .def("collapse", &dolfinx::fem::FunctionSpace::collapse) .def("component", &dolfinx::fem::FunctionSpace::component) .def("contains", &dolfinx::fem::FunctionSpace::contains, @@ -688,7 +692,7 @@ void declare_form(nb::module_& m, std::string type) nb::arg("entity_maps"), nb::arg("mesh").none()) .def( "__init__", - [](dolfinx::fem::Form* fp, std::uintptr_t form, + [](dolfinx::fem::Form* fp, std::vector forms, const std::vector< std::shared_ptr>>& spaces, const std::vector(form); + std::vector> ps; + for (auto form : forms) + ps.push_back(*(reinterpret_cast(form))); new (fp) dolfinx::fem::Form(dolfinx::fem::create_form_factory( - *p, spaces, coefficients, constants, sd, _entity_maps, + ps, spaces, coefficients, constants, sd, _entity_maps, mesh)); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), @@ -807,7 +813,7 @@ void declare_form(nb::module_& m, std::string type) } ufcx_form* p = reinterpret_cast(form); - return dolfinx::fem::create_form_factory(*p, spaces, coefficients, + return dolfinx::fem::create_form_factory({*p}, spaces, coefficients, constants, sd, {}, mesh); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), @@ -858,9 +864,11 @@ void declare_form(nb::module_& m, std::string type) nb::arg("constants"), nb::arg("subdomains"), nb::arg("entity_maps"), nb::arg("mesh"), "Create Form from a pointer to ufcx_form."); - m.def("create_sparsity_pattern", - &dolfinx::fem ::create_sparsity_pattern, nb::arg("a"), - "Create a sparsity pattern."); + m.def("create_sparsity_pattern", &dolfinx::fem::create_sparsity_pattern, + nb::arg("a"), "Create a sparsity pattern."); + + m.def("build_sparsity_pattern", &dolfinx::fem::build_sparsity_pattern, + nb::arg("pattern"), nb::arg("a"), "Build a sparsity pattern."); } template diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index f0198352adb..53199b6accb 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -650,6 +650,7 @@ void mesh(nb::module_& m) .def("index_map", &dolfinx::mesh::Topology::index_map, nb::arg("dim")) .def("index_maps", &dolfinx::mesh::Topology::index_maps, nb::arg("dim")) .def_prop_ro("cell_type", &dolfinx::mesh::Topology::cell_type) + .def_prop_ro("cell_types", &dolfinx::mesh::Topology::cell_types) .def_prop_ro( "entity_types", [](const dolfinx::mesh::Topology& self) diff --git a/python/test/unit/fem/test_forms.py b/python/test/unit/fem/test_forms.py index b1a2e5e09cf..f4ff8cccf33 100644 --- a/python/test/unit/fem/test_forms.py +++ b/python/test/unit/fem/test_forms.py @@ -89,7 +89,7 @@ def test_incorrect_element(): ) f = ftype( - module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)), + [module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form))], [space._cpp_object, space._cpp_object], [], [], @@ -101,7 +101,7 @@ def test_incorrect_element(): with pytest.raises(RuntimeError): f = ftype( - module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)), + [module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form))], [incorrect_space._cpp_object, incorrect_space._cpp_object], [], [],