From 01dfc5b57517f4c8d4c981b17050e537c1039132 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 8 Jul 2026 09:58:54 +0200 Subject: [PATCH 01/20] Use geometry type consistently --- cpp/dolfinx/fem/assemble_expression_impl.h | 6 +-- cpp/dolfinx/fem/assemble_matrix_impl.h | 40 ++++++--------- cpp/dolfinx/fem/assemble_scalar_impl.h | 44 +++++++--------- cpp/dolfinx/fem/assemble_vector_impl.h | 50 ++++++------------- cpp/dolfinx/fem/assembler.h | 38 +++----------- cpp/dolfinx/fem/interpolate.h | 8 +-- cpp/dolfinx/fem/traits.h | 7 ++- cpp/dolfinx/fem/utils.h | 18 ++++--- .../dolfinx/wrappers/dolfinx_wrappers/fem.h | 4 +- 9 files changed, 81 insertions(+), 134 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_expression_impl.h b/cpp/dolfinx/fem/assemble_expression_impl.h index 7f3dcd9ddd2..80a3729428f 100644 --- a/cpp/dolfinx/fem/assemble_expression_impl.h +++ b/cpp/dolfinx/fem/assemble_expression_impl.h @@ -61,11 +61,11 @@ namespace dolfinx::fem::impl /// @param[in] perms Entity permutation information for use in `fn`. template void tabulate_expression( - std::span values, fem::FEkernel auto fn, + std::span values, fem::FEkernel auto fn, std::array Xshape, std::size_t value_size, std::size_t num_argument_dofs, md::mdspan> x_dofmap, - std::span> x, + std::span x, md::mdspan> coeffs, std::span constants, fem::MDSpan2 auto entities, std::span cell_info, @@ -153,7 +153,7 @@ void tabulate_expression( /// expression values at the evaluation points. template void tabulate_expression( - std::span values, fem::FEkernel auto fn, + std::span values, fem::FEkernel auto fn, std::array Xshape, std::size_t value_size, md::mdspan> coeffs, std::span constants, const mesh::Mesh& mesh, diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 410e8b5d5f2..cd20de18d7a 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -62,18 +62,16 @@ using mdspan2_t = md::mdspan>; /// function mesh. /// @param cell_info1 Cell permutation information for the trial /// function mesh. -template +template void assemble_cells_matrix( la::MatSet auto mat_set, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, std::span cells, std::tuple> dofmap0, fem::DofTransformKernel auto P0, std::tuple> dofmap1, fem::DofTransformKernel auto P1T, std::span bc0, - std::span bc1, FEkernel auto kernel, + std::span bc1, FEkernel auto kernel, md::mdspan> coeffs, std::span constants, std::span cell_info0, std::span cell_info1) @@ -90,7 +88,7 @@ void assemble_cells_matrix( const int ndim0 = bs0 * num_dofs0; const int ndim1 = bs1 * num_dofs1; std::vector Ae(ndim0 * ndim1); - std::vector> cdofs(3 * x_dofmap.extent(1)); + std::vector cdofs(3 * x_dofmap.extent(1)); // Iterate over active cells assert(cells0.size() == cells.size()); @@ -226,12 +224,10 @@ void assemble_cells_matrix( /// function mesh. /// @param[in] perms Entity permutation integer. Empty if entity /// permutations are not required. -template +template void assemble_entities( la::MatSet auto mat_set, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> entities, @@ -245,7 +241,7 @@ void assemble_entities( std::extents>> dofmap1, fem::DofTransformKernel auto P1T, std::span bc0, - std::span bc1, FEkernel auto kernel, + std::span bc1, FEkernel auto kernel, md::mdspan> coeffs, std::span constants, std::span cell_info0, std::span cell_info1, @@ -258,7 +254,7 @@ void assemble_entities( const auto [dmap1, bs1, entities1] = dofmap1; // Data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); + std::vector cdofs(3 * x_dofmap.extent(1)); const int num_dofs0 = dmap0.extent(1); const int num_dofs1 = dmap1.extent(1); const int ndim0 = bs0 * num_dofs0; @@ -393,12 +389,10 @@ void assemble_entities( /// function mesh. /// @param[in] perms Facet permutation integer. Empty if facet /// permutations are not required. -template +template void assemble_interior_facets( la::MatSet auto mat_set, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> facets, @@ -412,7 +406,7 @@ void assemble_interior_facets( std::extents>> dofmap1, fem::DofTransformKernel auto P1T, std::span bc0, - std::span bc1, FEkernel auto kernel, + std::span bc1, FEkernel auto kernel, md::mdspan> coeffs, @@ -427,7 +421,7 @@ void assemble_interior_facets( const auto [dmap1, bs1, facets1] = dofmap1; // Data structures used in assembly - using X = scalar_value_t; + using X = U; std::vector cdofs(2 * x_dofmap.extent(1) * 3); std::span cdofs0(cdofs.data(), x_dofmap.extent(1) * 3); std::span cdofs1(cdofs.data() + x_dofmap.extent(1) * 3, @@ -612,9 +606,7 @@ void assemble_interior_facets( template void assemble_matrix( la::MatSet auto mat_set, const Form& a, - md::mdspan, - md::extents> - x, + md::mdspan> x, std::span constants, const std::map, std::pair, int>>& coefficients, @@ -687,7 +679,7 @@ void assemble_matrix( std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx); auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); assert(cells.size() * cstride == coeffs.size()); - impl::assemble_cells_matrix( + impl::assemble_cells_matrix( mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0}, P0, {dofs1, bs1, cells1}, P1T, bc0, bc1, fn, md::mdspan(coeffs.data(), cells.size(), cstride), constants, @@ -732,7 +724,7 @@ void assemble_matrix( std::span facets0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0); std::span facets1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0); assert((facets.size() / 4) * 2 * cstride == coeffs.size()); - impl::assemble_interior_facets( + impl::assemble_interior_facets( mat_set, x_dofmap, x, mdspanx22_t(facets.data(), facets.size() / 4, 2, 2), {*dofmap0, bs0, @@ -777,7 +769,7 @@ void assemble_matrix( std::span e1 = a.domain_arg(itg_type, 1, i, 0); mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); assert((entities.size() / 2) * cstride == coeffs.size()); - impl::assemble_entities( + impl::assemble_entities( mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0, {dofs1, bs1, entities1}, P1T, bc0, bc1, fn, md::mdspan(coeffs.data(), entities.extent(0), cstride), constants, diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 23070de9f96..d2c50dbffbe 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -17,20 +17,20 @@ #include #include #include +#include #include namespace dolfinx::fem::impl { /// Assemble functional over cells -template -T assemble_cells(mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, - std::span cells, FEkernel auto fn, - std::span constants, - md::mdspan> coeffs, - std::span> cdofs_b) +template +T assemble_cells( + mdspan2_t x_dofmap, + md::mdspan> x, + std::span cells, FEkernel auto fn, + std::span constants, + md::mdspan> coeffs, + std::span> cdofs_b) { T value(0); if (cells.empty()) @@ -65,19 +65,17 @@ T assemble_cells(mdspan2_t x_dofmap, /// However, entities may be attached to more than one cell. This function /// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen /// from cell used to define the entity. -template +template T assemble_entities( mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> entities, - FEkernel auto fn, std::span constants, + FEkernel auto fn, std::span constants, md::mdspan> coeffs, md::mdspan> perms, - std::span> cdofs_b) + std::span> cdofs_b) { T value(0); if (entities.empty()) @@ -106,21 +104,19 @@ T assemble_entities( } /// Assemble functional over interior facets -template +template T assemble_interior_facets( mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> facets, - FEkernel auto fn, std::span constants, + FEkernel auto fn, std::span constants, md::mdspan> coeffs, md::mdspan> perms, - std::span> cdofs_b) + std::span> cdofs_b) { T value(0); if (facets.empty()) @@ -160,9 +156,7 @@ T assemble_interior_facets( template T assemble_scalar( const fem::Form& M, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, std::span constants, const std::map, std::pair, int>>& coefficients, @@ -171,7 +165,7 @@ T assemble_scalar( std::shared_ptr> mesh = M.mesh(); assert(mesh); - std::vector> cdofs_b(2 * 3 * x_dofmap.extent(1)); + std::vector cdofs_b(2 * 3 * x_dofmap.extent(1)); T value = 0; for (int i = 0; i < M.num_integrals(IntegralType::cell, cell_type_idx); ++i) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index fe61b3a6959..9dfa55ebb4e 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -60,18 +60,16 @@ using mdspan2_t = md::mdspan>; /// coefficient for cell `i`. /// @param[in] cell_info0 Cell permutation information for the test /// function mesh. -template ::value_type> requires std::is_same_v::value_type, T> void assemble_cells( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, std::span cells, std::tuple> dofmap, - FEkernel auto kernel, std::span constants, + FEkernel auto kernel, std::span constants, md::mdspan> coeffs, std::span cell_info0) { @@ -82,7 +80,7 @@ void assemble_cells( assert(_bs < 0 or _bs == bs); // Create data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); + std::vector cdofs(3 * x_dofmap.extent(1)); std::vector be(bs * dmap.extent(1)); // Iterate over active cells @@ -153,14 +151,12 @@ void assemble_cells( /// function mesh. /// @param[in] perms Entity permutation integer. Empty if entity /// permutations are not required. -template ::value_type> requires std::is_same_v::value_type, T> void assemble_entities( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> entities, @@ -168,7 +164,7 @@ void assemble_entities( md::mdspan>> dofmap, - FEkernel auto kernel, std::span constants, + FEkernel auto kernel, std::span constants, md::mdspan> coeffs, std::span cell_info0, md::mdspan> perms) @@ -181,7 +177,7 @@ void assemble_entities( // Create data structures used in assembly const int num_dofs = dmap.extent(1); - std::vector> cdofs(3 * x_dofmap.extent(1)); + std::vector cdofs(3 * x_dofmap.extent(1)); std::vector be(bs * num_dofs); assert(entities0.size() == entities.size()); for (std::size_t f = 0; f < entities.extent(0); ++f) @@ -248,14 +244,12 @@ void assemble_entities( /// function mesh. /// @param[in] perms Facet permutation integer. Empty if facet /// permutations are not required. -template ::value_type> requires std::is_same_v::value_type, T> void assemble_interior_facets( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, + md::mdspan> x, md::mdspan> facets, @@ -263,14 +257,14 @@ void assemble_interior_facets( md::mdspan>> dofmap, - FEkernel auto kernel, std::span constants, + FEkernel auto kernel, std::span constants, md::mdspan> coeffs, std::span cell_info0, md::mdspan> perms) { - using X = scalar_value_t; + using X = U; if (facets.empty()) return; @@ -564,9 +558,7 @@ template ::value_type, T> void assemble_vector( V&& b, const Form& L, - md::mdspan, - md::extents> - x, + md::mdspan> x, std::span constants, const std::map, std::pair, int>>& coefficients) @@ -762,22 +754,12 @@ void assemble_vector( std::pair, int>>& coefficients) { using mdspanx3_t - = md::mdspan, - md::extents>; + = md::mdspan>; std::shared_ptr> mesh = L.mesh(); assert(mesh); auto x = mesh->geometry().x(); - if constexpr (std::is_same_v>) - { - impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3), - constants, coefficients); - } - else - { - std::vector> _x(x.begin(), x.end()); - impl::assemble_vector(b, L, mdspanx3_t(_x.data(), _x.size() / 3, 3), - constants, coefficients); - } + impl::assemble_vector(b, L, mdspanx3_t(x.data(), x.size() / 3, 3), constants, + coefficients); } } // namespace dolfinx::fem::impl diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 0654e108e0d..b9259b8da40 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -176,8 +176,7 @@ T assemble_scalar( std::pair, int>>& coefficients) { using mdspanx3_t - = md::mdspan, - md::extents>; + = md::mdspan>; std::shared_ptr> mesh = M.mesh(); assert(mesh); @@ -191,19 +190,9 @@ T assemble_scalar( // Geometry dofmap and data md::mdspan> x_dofmap = mesh->geometry().dofmaps().at(cell_type_idx); - if constexpr (std::is_same_v>) - { - val += impl::assemble_scalar(M, x_dofmap, - mdspanx3_t(x.data(), x.size() / 3, 3), - constants, coefficients, cell_type_idx); - } - else - { - std::vector> _x(x.begin(), x.end()); - val += impl::assemble_scalar(M, x_dofmap, - mdspanx3_t(_x.data(), _x.size() / 3, 3), - constants, coefficients, cell_type_idx); - } + val += impl::assemble_scalar(M, x_dofmap, + mdspanx3_t(x.data(), x.size() / 3, 3), + constants, coefficients, cell_type_idx); } return val; } @@ -469,25 +458,14 @@ void assemble_matrix( { common::Timer t_assm("[Assemble Matrix]"); using mdspanx3_t - = md::mdspan, - md::extents>; + = md::mdspan>; std::shared_ptr> mesh = a.mesh(); assert(mesh); std::span x = mesh->geometry().x(); - if constexpr (std::is_same_v>) - { - impl::assemble_matrix( - mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3), constants, - coefficients, dof_marker0, dof_marker1); - } - else - { - std::vector> _x(x.begin(), x.end()); - impl::assemble_matrix( - mat_add, a, mdspanx3_t(_x.data(), _x.size() / 3, 3), constants, - coefficients, dof_marker0, dof_marker1); - } + impl::assemble_matrix( + mat_add, a, mdspanx3_t(x.data(), x.size() / 3, 3), constants, + coefficients, dof_marker0, dof_marker1); } /// @brief Assemble bilinear form into a matrix diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index a3f6c2f4a8e..a1fdbc65a34 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -310,7 +310,9 @@ void scatter_values(MPI_Comm comm, std::span src_ranks, template void interpolation_apply(U&& Pi, V&& data, std::span coeffs, int bs) { - using X = typename dolfinx::scalar_value_t; + // Geometry (real) scalar type, taken from the interpolation operator Pi + // rather than scalar_value_t so it is independent of the value scalar T. + using X = typename std::remove_cvref_t::value_type; // Compute coefficients = Pi * x (matrix-vector multiply) if (bs == 1) @@ -421,7 +423,7 @@ void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, auto [i_m, im_shape] = element1->create_interpolation_operator(*element0); // Iterate over mesh and interpolate on each cell - using X = typename dolfinx::scalar_value_t; + using X = U; assert(cells0.size() == cells1.size()); for (auto cell0_it = cells0.begin(), cell1_it = cells1.begin(); cell0_it != cells0.end() and cell1_it != cells1.end(); @@ -663,7 +665,7 @@ void interpolate_nonmatching_maps(Function& u1, coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k]; // Evaluate v at the interpolation points (physical space values) - using X = typename dolfinx::scalar_value_t; + using X = U; for (std::size_t p = 0; p < Xshape[0]; ++p) { for (int k = 0; k < bs0; ++k) diff --git a/cpp/dolfinx/fem/traits.h b/cpp/dolfinx/fem/traits.h index 3a8d521e1de..d6915de30ae 100644 --- a/cpp/dolfinx/fem/traits.h +++ b/cpp/dolfinx/fem/traits.h @@ -25,10 +25,9 @@ concept DofTransformKernel /// /// Kernel functions that can be passed to an assembler for execution /// must satisfy this concept. -template -concept FEkernel - = std::is_invocable_v*, - const int*, const std::uint8_t*, void*>; +template > +concept FEkernel = std::is_invocable_v; /// @brief Concept for mdspan of rank 1 or 2. template diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 721b79eadb0..25493241b74 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -970,17 +970,18 @@ Expression create_expression( static_cast(e.entity_dimension)}; std::vector value_shape(e.value_shape, e.value_shape + e.num_components); - std::function*, - const int*, const std::uint8_t*, void*)> + std::function tabulate_tensor = nullptr; if constexpr (std::is_same_v) tabulate_tensor = e.tabulate_tensor_float32; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS else if constexpr (std::is_same_v>) { - tabulate_tensor = reinterpret_cast*, const int*, - const unsigned char*, void*)>(e.tabulate_tensor_complex64); + tabulate_tensor + = reinterpret_cast( + e.tabulate_tensor_complex64); } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS else if constexpr (std::is_same_v) @@ -988,9 +989,10 @@ Expression create_expression( #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS else if constexpr (std::is_same_v>) { - tabulate_tensor = reinterpret_cast*, const int*, - const unsigned char*, void*)>(e.tabulate_tensor_complex128); + tabulate_tensor + = reinterpret_cast( + e.tabulate_tensor_complex128); } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS else diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 530635651bf..6eb30dfb7bb 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -333,11 +333,9 @@ void declare_function_space(nb::module_& m, std::string type) } // Declare DirichletBC objects for type T -template +template > void declare_objects(nb::module_& m, std::string type) { - using U = typename dolfinx::scalar_value_t; - // dolfinx::fem::DirichletBC std::string pyclass_name = std::string("DirichletBC_") + type; nb::class_> dirichletbc( From 4dea4827d7f1084c0485bff8ed2bd92ce82824a0 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 8 Jul 2026 10:39:13 +0200 Subject: [PATCH 02/20] One more geom type --- python/dolfinx/wrappers/dolfinx_wrappers/fem.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 6eb30dfb7bb..4a1433f2a59 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -687,7 +687,7 @@ void declare_objects(nb::module_& m, std::string type) { auto tabulate_expression_ptr = (void (*)(T*, const T*, const T*, - const typename geom_type::value_type*, + const U*, const int*, const std::uint8_t*, void*))fn_addr; new (ex) dolfinx::fem::Expression( coefficients, constants, std::span(X.data(), X.size()), @@ -780,7 +780,7 @@ void declare_form(nb::module_& m, std::string type) { auto kn_ptr = (void (*)(T*, const T*, const T*, - const typename geom_type::value_type*, + const U*, const int*, const std::uint8_t*, void*))ptr; _integrals.insert( {{type, id, 0}, From 6be352f5d97f9fb0e8c4662694be6496e941bcca Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 8 Jul 2026 10:46:40 +0200 Subject: [PATCH 03/20] More in extract kernel --- cpp/dolfinx/fem/utils.h | 6 +++--- python/dolfinx/wrappers/dolfinx_wrappers/fem.h | 11 ----------- 2 files changed, 3 insertions(+), 14 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 25493241b74..8714a20edcc 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -568,7 +568,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - impl::kernel_t k = impl::extract_kernel(integral); + impl::kernel_t k = impl::extract_kernel(integral); if (!k) { throw std::runtime_error( @@ -648,7 +648,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - impl::kernel_t k = impl::extract_kernel(integral); + impl::kernel_t k = impl::extract_kernel(integral); assert(k); // Build list of entities to assembler over @@ -774,7 +774,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - impl::kernel_t k = impl::extract_kernel(integral); + impl::kernel_t k = impl::extract_kernel(integral); // Build list of entities to assembler over auto e_to_c = topology->connectivity(dim, tdim); diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 4a1433f2a59..9eb9d719a0b 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -58,17 +58,6 @@ namespace dolfinx_wrappers namespace nb = nanobind; namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE; -template -struct geom_type -{ - typedef T value_type; -}; -template -struct geom_type> -{ - typedef typename T::value_type value_type; -}; - // Copy a container of pointers to a // std::vector> auto ptr_to_ref_wrapper_vec(auto& x) From ff0ea507d266219282dc9d11671f6ae16eb20c46 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 8 Jul 2026 10:57:35 +0200 Subject: [PATCH 04/20] Clang format --- python/dolfinx/wrappers/dolfinx_wrappers/fem.h | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 9eb9d719a0b..0c6b012d448 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -675,9 +675,8 @@ void declare_objects(nb::module_& m, std::string type) argument_space) { auto tabulate_expression_ptr - = (void (*)(T*, const T*, const T*, - const U*, - const int*, const std::uint8_t*, void*))fn_addr; + = (void (*)(T*, const T*, const T*, const U*, const int*, + const std::uint8_t*, void*))fn_addr; new (ex) dolfinx::fem::Expression( coefficients, constants, std::span(X.data(), X.size()), {X.shape(0), X.shape(1)}, tabulate_expression_ptr, value_shape, @@ -768,9 +767,8 @@ void declare_form(nb::module_& m, std::string type) for (auto& [id, ptr, e, c] : kernels) { auto kn_ptr - = (void (*)(T*, const T*, const T*, - const U*, - const int*, const std::uint8_t*, void*))ptr; + = (void (*)(T*, const T*, const T*, const U*, const int*, + const std::uint8_t*, void*))ptr; _integrals.insert( {{type, id, 0}, {kn_ptr, From 67bb87535ab2618acf5ce8901fab7cba26ae0abe Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 8 Jul 2026 13:21:15 +0200 Subject: [PATCH 05/20] Leftover templates --- python/dolfinx/wrappers/dolfinx_wrappers/fem.h | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 0c6b012d448..11f2ba48ef3 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -728,11 +728,9 @@ void declare_objects(nb::module_& m, std::string type) "Create Expression from a pointer to ufc_form."); } -template +template > void declare_form(nb::module_& m, std::string type) { - using U = typename dolfinx::scalar_value_t; - // dolfinx::fem::Form std::string pyclass_name_form = std::string("Form_") + type; nb::class_>(m, pyclass_name_form.c_str(), @@ -758,7 +756,7 @@ void declare_form(nb::module_& m, std::string type) std::shared_ptr> mesh) { std::map, - dolfinx::fem::integral_data> + dolfinx::fem::integral_data> _integrals; // Loop over kernel for each entity type From eab634ab765df68dec0484e65267821b6443b5bd Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 06:34:57 +0200 Subject: [PATCH 06/20] Use common type registry --- cpp/dolfinx/fem/interpolate.h | 4 +- cpp/dolfinx/la/MatrixCSR.h | 2 +- python/dolfinx/fem/__init__.py | 3 ++ python/dolfinx/fem/bcs.py | 21 ++++---- python/dolfinx/fem/forms.py | 12 ++++- python/dolfinx/fem/function.py | 51 +++++++------------ python/dolfinx/fem/typemap.py | 91 ++++++++++++++++++++++++++++++++++ 7 files changed, 137 insertions(+), 47 deletions(-) create mode 100644 python/dolfinx/fem/typemap.py diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index a1fdbc65a34..cbe7a5ad0ab 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -423,7 +423,7 @@ void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, auto [i_m, im_shape] = element1->create_interpolation_operator(*element0); // Iterate over mesh and interpolate on each cell - using X = U; + using X = scalar_value_t; assert(cells0.size() == cells1.size()); for (auto cell0_it = cells0.begin(), cell1_it = cells1.begin(); cell0_it != cells0.end() and cell1_it != cells1.end(); @@ -665,7 +665,7 @@ void interpolate_nonmatching_maps(Function& u1, coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k]; // Evaluate v at the interpolation points (physical space values) - using X = U; + using X = scalar_value_t; for (std::size_t p = 0; p < Xshape[0]; ++p) { for (int k = 0; k < bs0; ++k) diff --git a/cpp/dolfinx/la/MatrixCSR.h b/cpp/dolfinx/la/MatrixCSR.h index 46b117f4683..56a526c2070 100644 --- a/cpp/dolfinx/la/MatrixCSR.h +++ b/cpp/dolfinx/la/MatrixCSR.h @@ -347,7 +347,7 @@ class MatrixCSR { const std::size_t nrows = num_all_rows(); const std::size_t ncols = _index_maps[1]->size_global(); - std::vector A(nrows * ncols * _bs[0] * _bs[1], 0.0); + std::vector A(nrows * ncols * _bs[0] * _bs[1], value_type(0)); for (std::size_t r = 0; r < nrows; ++r) { for (std::int32_t j = _row_ptr[r]; j < _row_ptr[r + 1]; ++j) diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index 7b45c04735b..65541f8b9cd 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -48,6 +48,7 @@ FunctionSpace, functionspace, ) +from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.fem.utils import ( build_sparsity_pattern, compute_integration_domains, @@ -99,8 +100,10 @@ "locate_dofs_geometrical", "locate_dofs_topological", "mixed_topology_form", + "get_cpp_type", "pack_coefficients", "pack_constants", + "register_cpp_type", "set_bc", "transpose_dofmap", ] diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index b8bb1be1aed..84128e59d94 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -21,6 +21,7 @@ import dolfinx from dolfinx import cpp as _cpp from dolfinx.fem.function import Constant, Function, FunctionSpace +from dolfinx.fem.typemap import get_cpp_type from dolfinx.typing import Scalar @@ -200,19 +201,19 @@ def dirichletbc( try: dtype = value.dtype - if np.issubdtype(dtype, np.float32): - bctype = _cpp.fem.DirichletBC_float32 - elif np.issubdtype(dtype, np.float64): - bctype = _cpp.fem.DirichletBC_float64 - elif np.issubdtype(dtype, np.complex64): - bctype = _cpp.fem.DirichletBC_complex64 - elif np.issubdtype(dtype, np.complex128): - bctype = _cpp.fem.DirichletBC_complex128 - else: - raise NotImplementedError(f"Type {value.dtype} not supported.") except AttributeError: raise AttributeError("Boundary condition value must have a dtype attribute.") + # Geometry type is fixed by the function space (or the value's space), + # defaulting to matched precision when neither carries one. + if V is not None: + geometry_dtype = np.dtype(V.element.dtype) + elif isinstance(value, Function): + geometry_dtype = np.dtype(value.function_space.element.dtype) + else: + geometry_dtype = np.dtype(dtype).type(0).real.dtype + bctype = get_cpp_type("DirichletBC", dtype, geometry_dtype) + # Unwrap value object, if required if isinstance(value, np.ndarray): _value = value diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 96513d9cc36..76774f2f77e 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -25,6 +25,7 @@ from dolfinx import default_scalar_type, jit from dolfinx.fem import IntegralType from dolfinx.fem.function import Constant, Function, FunctionSpace +from dolfinx.fem.typemap import get_cpp_type from dolfinx.typing import Scalar if typing.TYPE_CHECKING: @@ -280,7 +281,6 @@ def mixed_topology_form( 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() @@ -295,6 +295,9 @@ def mixed_topology_form( raise RuntimeError("Expecting to find a Mesh in the form.") comm = mesh.comm if jit_comm is None else jit_comm + # Geometry type is fixed by the mesh. + ftype = get_cpp_type("Form", dtype, mesh.geometry.x.dtype) + ufcx_forms = [] modules = [] codes = [] @@ -364,7 +367,6 @@ def form( form_compiler_options = dict() form_compiler_options["scalar_type"] = dtype - ftype = form_cpp_class(dtype) def _form(form): """Compile a single UFL form.""" @@ -381,6 +383,9 @@ def _form(form): raise RuntimeError("Expecting to find a Mesh in the form.") comm = msh.comm if jit_comm is None else jit_comm + # Geometry type is fixed by the mesh. + ftype = get_cpp_type("Form", dtype, msh.geometry.x.dtype) + ufcx_form, module, code = jit.ffcx_jit( comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options ) @@ -443,6 +448,9 @@ def _zero_form(form): assert len(V) > 0 msh = V[0].mesh + # Geometry type is fixed by the mesh. + ftype = get_cpp_type("Form", dtype, msh.geometry.x.dtype) + f = ftype( spaces=V, integrals={}, diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 0d26d9d1f12..f35ddfccdc0 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -22,6 +22,7 @@ from dolfinx import default_scalar_type, jit, la from dolfinx.fem.dofmap import DofMap from dolfinx.fem.element import FiniteElement, finiteelement +from dolfinx.fem.typemap import get_cpp_type from dolfinx.geometry import PointOwnershipData from dolfinx.typing import Real, Scalar @@ -205,17 +206,13 @@ def __init__( else: raise RuntimeError("Expressions with more that one Argument not allowed.") - def _create_expression(dtype): - if np.issubdtype(dtype, np.float32): - return _cpp.fem.create_expression_float32 - elif np.issubdtype(dtype, np.float64): - return _cpp.fem.create_expression_float64 - elif np.issubdtype(dtype, np.complex64): - return _cpp.fem.create_expression_complex64 - elif np.issubdtype(dtype, np.complex128): - return _cpp.fem.create_expression_complex128 - else: - raise NotImplementedError(f"Type {dtype} not supported.") + # Geometry type is fixed by the expression's mesh. + expr_domains = ufl.domain.extract_domains(e) + if len(expr_domains) > 0: + geometry_dtype = expr_domains[0].ufl_cargo().geometry.x.dtype + else: + geometry_dtype = np.dtype(dtype).type(0).real.dtype + create_expression = get_cpp_type("Expression", dtype, geometry_dtype) _entity_maps = ( [entity_map._cpp_object for entity_map in entity_maps] @@ -223,7 +220,7 @@ def _create_expression(dtype): else [] ) ffi = module.ffi - self._cpp_object = _create_expression(dtype)( + self._cpp_object = create_expression( ffi.cast("uintptr_t", ffi.addressof(self._ufcx_expression)), coeffs, constants, @@ -376,27 +373,14 @@ def __init__( if dtype is None: dtype = default_scalar_type - assert np.issubdtype(V.element.dtype, np.dtype(dtype).type(0).real.dtype), ( - "Incompatible FunctionSpace dtype and requested dtype." - ) - - # Create cpp Function - def functiontype(dtype): - if np.issubdtype(dtype, np.float32): - return _cpp.fem.Function_float32 - elif np.issubdtype(dtype, np.float64): - return _cpp.fem.Function_float64 - elif np.issubdtype(dtype, np.complex64): - return _cpp.fem.Function_complex64 - elif np.issubdtype(dtype, np.complex128): - return _cpp.fem.Function_complex128 - else: - raise NotImplementedError(f"Type {dtype} not supported.") - + # Scalar type (dtype) is independent of the geometry type, which + # is fixed by the function space. + geometry_dtype = V.element.dtype + cpp_type = get_cpp_type("Function", dtype, geometry_dtype) if x is not None: - self._cpp_object = functiontype(dtype)(V._cpp_object, x._cpp_object) # type: ignore + self._cpp_object = cpp_type(V._cpp_object, x._cpp_object) # type: ignore else: - self._cpp_object = functiontype(dtype)(V._cpp_object) # type: ignore + self._cpp_object = cpp_type(V._cpp_object) # type: ignore # Initialize the ufl.FunctionSpace super().__init__(V.ufl_function_space()) @@ -651,18 +635,21 @@ def functionspace( | tuple[str, int, tuple] | tuple[str, int, tuple, bool] ), + dtype: np.dtype = None, ) -> FunctionSpace: """Create a finite element function space. Args: mesh: Mesh that space is defined on. element: Finite element description. + dtype: Scalar type of the element. Returns: A function space. """ # Create UFL element - dtype = mesh.geometry.x.dtype + if dtype is None: + dtype = mesh.geometry.x.dtype try: e = ElementMetaData(*element) # type: ignore ufl_e = basix.ufl.element( diff --git a/python/dolfinx/fem/typemap.py b/python/dolfinx/fem/typemap.py new file mode 100644 index 00000000000..85aec2b4512 --- /dev/null +++ b/python/dolfinx/fem/typemap.py @@ -0,0 +1,91 @@ +# Copyright (C) 2026 Michal Habera +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +"""Registry from (object kind, scalar dtype, geometry dtype) to compiled type. + +``Function``, ``Expression``, ``Form`` and ``DirichletBC`` are templated on +both a scalar and a mesh geometry type. DOLFINx registers the matched-precision +pairs; downstream libraries register further pairs (e.g. fp16 dofs on an fp64 +mesh) via :func:`register_cpp_type`. +""" + +from __future__ import annotations + +import typing + +import numpy as np +import numpy.typing as npt + +from dolfinx import cpp as _cpp + +__all__ = ["get_cpp_type", "register_cpp_type"] + +# Object kinds templated on (scalar, geometry). +_KINDS = ("Function", "Expression", "Form", "DirichletBC") + +# (kind, scalar dtype, geometry dtype) -> compiled C++ class or factory. +_registry: dict[tuple[str, np.dtype, np.dtype], typing.Any] = {} + + +def register_cpp_type( + kind: str, + scalar_dtype: npt.DTypeLike, + geometry_dtype: npt.DTypeLike, + cpp_type: typing.Any, +) -> None: + """Register the compiled C++ type for an object kind and dtype pair. + + Args: + kind: Object kind, one of ``"Function"``, ``"Expression"``, + ``"Form"`` or ``"DirichletBC"``. + scalar_dtype: Scalar (degree-of-freedom) type. + geometry_dtype: Mesh geometry type. + cpp_type: Compiled C++ class or factory that the wrapper builds + with for this pair. + """ + if kind not in _KINDS: + raise ValueError(f"Unknown object kind '{kind}'. Expected one of {_KINDS}.") + _registry[(kind, np.dtype(scalar_dtype), np.dtype(geometry_dtype))] = cpp_type + + +def get_cpp_type( + kind: str, scalar_dtype: npt.DTypeLike, geometry_dtype: npt.DTypeLike +) -> typing.Any: + """Look up the compiled C++ type for an object kind and dtype pair. + + Args: + kind: Object kind (see :func:`register_cpp_type`). + scalar_dtype: Scalar (degree-of-freedom) type. + geometry_dtype: Mesh geometry type. + + Returns: + The registered C++ class or factory. + """ + scalar, geometry = np.dtype(scalar_dtype), np.dtype(geometry_dtype) + try: + return _registry[(kind, scalar, geometry)] + except KeyError: + raise NotImplementedError( + f"No {kind} type registered for scalar dtype '{scalar}' on geometry " + f"dtype '{geometry}'. Downstream libraries can register one with " + f"dolfinx.fem.register_cpp_type()." + ) from None + + +# Register the matched-precision built-ins (geometry == real part of scalar). +# ``Expression`` is built through a factory function, the others through their +# class constructor. +_MATCHED: tuple[tuple[npt.DTypeLike, npt.DTypeLike, str], ...] = ( + (np.float32, np.float32, "float32"), + (np.float64, np.float64, "float64"), + (np.complex64, np.float32, "complex64"), + (np.complex128, np.float64, "complex128"), +) + +for _scalar, _geom, _name in _MATCHED: + register_cpp_type("Function", _scalar, _geom, getattr(_cpp.fem, f"Function_{_name}")) + register_cpp_type("Expression", _scalar, _geom, getattr(_cpp.fem, f"create_expression_{_name}")) + register_cpp_type("Form", _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) + register_cpp_type("DirichletBC", _scalar, _geom, getattr(_cpp.fem, f"DirichletBC_{_name}")) From 0c3646b56bf61cff1e3740fe679456f2737aaacc Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 06:46:56 +0200 Subject: [PATCH 07/20] Lint fix --- python/dolfinx/fem/__init__.py | 2 +- python/dolfinx/fem/typemap.py | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index 65541f8b9cd..9e9ec025548 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -95,12 +95,12 @@ "form", "form_cpp_class", "functionspace", + "get_cpp_type", "interpolate_geometry", "interpolation_matrix", "locate_dofs_geometrical", "locate_dofs_topological", "mixed_topology_form", - "get_cpp_type", "pack_coefficients", "pack_constants", "register_cpp_type", diff --git a/python/dolfinx/fem/typemap.py b/python/dolfinx/fem/typemap.py index 85aec2b4512..81d5636ea01 100644 --- a/python/dolfinx/fem/typemap.py +++ b/python/dolfinx/fem/typemap.py @@ -3,12 +3,12 @@ # This file is part of DOLFINx (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -"""Registry from (object kind, scalar dtype, geometry dtype) to compiled type. +"""Registry from (kind, scalar dtype, geometry dtype) to compiled type. -``Function``, ``Expression``, ``Form`` and ``DirichletBC`` are templated on -both a scalar and a mesh geometry type. DOLFINx registers the matched-precision -pairs; downstream libraries register further pairs (e.g. fp16 dofs on an fp64 -mesh) via :func:`register_cpp_type`. +``Function``, ``Expression``, ``Form`` and ``DirichletBC`` are templated +on both a scalar and a mesh geometry type. DOLFINx registers the +matched-precision pairs; downstream libraries register further pairs +(e.g. fp16 dofs on an fp64 mesh) via :func:`register_cpp_type`. """ from __future__ import annotations @@ -74,9 +74,9 @@ def get_cpp_type( ) from None -# Register the matched-precision built-ins (geometry == real part of scalar). -# ``Expression`` is built through a factory function, the others through their -# class constructor. +# Register the matched-precision built-ins (geometry == real part of +# scalar). ``Expression`` is built through a factory, the others through +# their class constructor. _MATCHED: tuple[tuple[npt.DTypeLike, npt.DTypeLike, str], ...] = ( (np.float32, np.float32, "float32"), (np.float64, np.float64, "float64"), From ece6cfd3c24ec1bea2362fea0657cbbe3a195dab Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 07:08:58 +0200 Subject: [PATCH 08/20] Extract dtype from geometry array --- python/dolfinx/fem/bcs.py | 8 ++++---- python/dolfinx/fem/function.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index 84128e59d94..2be2085590c 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -204,12 +204,12 @@ def dirichletbc( except AttributeError: raise AttributeError("Boundary condition value must have a dtype attribute.") - # Geometry type is fixed by the function space (or the value's space), - # defaulting to matched precision when neither carries one. + # Geometry type is the mesh geometry type of the function space (or the + # value's space), defaulting to matched precision when neither has one. if V is not None: - geometry_dtype = np.dtype(V.element.dtype) + geometry_dtype = np.dtype(V.mesh.geometry.x.dtype) elif isinstance(value, Function): - geometry_dtype = np.dtype(value.function_space.element.dtype) + geometry_dtype = np.dtype(value.function_space.mesh.geometry.x.dtype) else: geometry_dtype = np.dtype(dtype).type(0).real.dtype bctype = get_cpp_type("DirichletBC", dtype, geometry_dtype) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index f35ddfccdc0..5234760a87f 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -374,8 +374,8 @@ def __init__( dtype = default_scalar_type # Scalar type (dtype) is independent of the geometry type, which - # is fixed by the function space. - geometry_dtype = V.element.dtype + # is fixed by the mesh. + geometry_dtype = V.mesh.geometry.x.dtype cpp_type = get_cpp_type("Function", dtype, geometry_dtype) if x is not None: self._cpp_object = cpp_type(V._cpp_object, x._cpp_object) # type: ignore From f89680a35db2d63fc005a3a3f61fde7aa4a971b8 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 08:58:02 +0200 Subject: [PATCH 09/20] Fix dispatch for default floats --- cpp/dolfinx/fem/utils.h | 8 ++++---- python/dolfinx/fem/function.py | 4 +--- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 8714a20edcc..10d69b21ed9 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -973,10 +973,10 @@ Expression create_expression( std::function tabulate_tensor = nullptr; - if constexpr (std::is_same_v) + if constexpr (std::is_same_v && std::is_same_v) tabulate_tensor = e.tabulate_tensor_float32; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v> && std::is_same_v) { tabulate_tensor = reinterpret_cast create_expression( e.tabulate_tensor_complex64); } #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) + else if constexpr (std::is_same_v && std::is_same_v) tabulate_tensor = e.tabulate_tensor_float64; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v> && std::is_same_v) { tabulate_tensor = reinterpret_cast FunctionSpace: """Create a finite element function space. @@ -648,8 +647,7 @@ def functionspace( A function space. """ # Create UFL element - if dtype is None: - dtype = mesh.geometry.x.dtype + dtype = mesh.geometry.x.dtype try: e = ElementMetaData(*element) # type: ignore ufl_e = basix.ufl.element( From 2431f18cd2a05b8a8899a29507c8212856101af0 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 20:37:14 +0200 Subject: [PATCH 10/20] Split Constant --- .../dolfinx/wrappers/dolfinx_wrappers/fem.h | 58 ++++++++++--------- python/dolfinx/wrappers/fem.cpp | 5 ++ 2 files changed, 36 insertions(+), 27 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 11f2ba48ef3..226a7c29edf 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -321,6 +321,37 @@ void declare_function_space(nb::module_& m, std::string type) } } +template +void declare_constant(nb::module_& m, std::string type) +{ + // dolfinx::fem::Constant + std::string pyclass_name_constant = std::string("Constant_") + type; + nb::class_>( + m, pyclass_name_constant.c_str(), + "Value constant with respect to integration domain") + .def( + "__init__", + [](dolfinx::fem::Constant* cp, + nb::ndarray c) + { + std::vector shape(c.shape_ptr(), + c.shape_ptr() + c.ndim()); + new (cp) + dolfinx::fem::Constant(std::span(c.data(), c.size()), shape); + }, + nb::arg("c").noconvert(), "Create a constant from a value array") + .def_prop_ro("dtype", [](const dolfinx::fem::Constant&) + { return dolfinx_wrappers::numpy_dtype_v; }) + .def_prop_ro( + "value", + [](dolfinx::fem::Constant& self) + { + return nb::ndarray( + self.value.data(), self.shape.size(), self.shape.data()); + }, + nb::rv_policy::reference_internal); +} + // Declare DirichletBC objects for type T template > void declare_objects(nb::module_& m, std::string type) @@ -628,33 +659,6 @@ void declare_objects(nb::module_& m, std::string type) .def_prop_ro("function_space", &dolfinx::fem::Function::function_space); - // dolfinx::fem::Constant - std::string pyclass_name_constant = std::string("Constant_") + type; - nb::class_>( - m, pyclass_name_constant.c_str(), - "Value constant with respect to integration domain") - .def( - "__init__", - [](dolfinx::fem::Constant* cp, - nb::ndarray c) - { - std::vector shape(c.shape_ptr(), - c.shape_ptr() + c.ndim()); - new (cp) - dolfinx::fem::Constant(std::span(c.data(), c.size()), shape); - }, - nb::arg("c").noconvert(), "Create a constant from a value array") - .def_prop_ro("dtype", [](const dolfinx::fem::Constant&) - { return dolfinx_wrappers::numpy_dtype_v; }) - .def_prop_ro( - "value", - [](dolfinx::fem::Constant& self) - { - return nb::ndarray( - self.value.data(), self.shape.size(), self.shape.data()); - }, - nb::rv_policy::reference_internal); - // dolfinx::fem::Expression std::string pyclass_name_expr = std::string("Expression_") + type; nb::class_>(m, pyclass_name_expr.c_str(), diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 423d136743e..b1b43b3d76c 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -162,6 +162,11 @@ void fem(nb::module_& m) .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral") .value("ridge", dolfinx::fem::IntegralType::ridge, "ridge integral"); + declare_constant(m, "float32"); + declare_constant(m, "float64"); + declare_constant>(m, "complex64"); + declare_constant>(m, "complex128"); + declare_objects(m, "float32"); declare_objects(m, "float64"); declare_objects>(m, "complex64"); From 00730b422faabddf7e77e34e8af101302a1af598 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 9 Jul 2026 20:41:49 +0200 Subject: [PATCH 11/20] Lint --- cpp/dolfinx/fem/utils.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 10d69b21ed9..ebea0ee4540 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -976,7 +976,8 @@ Expression create_expression( if constexpr (std::is_same_v && std::is_same_v) tabulate_tensor = e.tabulate_tensor_float32; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v> && std::is_same_v) + else if constexpr (std::is_same_v> + && std::is_same_v) { tabulate_tensor = reinterpret_cast create_expression( else if constexpr (std::is_same_v && std::is_same_v) tabulate_tensor = e.tabulate_tensor_float64; #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v> && std::is_same_v) + else if constexpr (std::is_same_v> + && std::is_same_v) { tabulate_tensor = reinterpret_cast Date: Sun, 12 Jul 2026 13:37:37 +0200 Subject: [PATCH 12/20] Use class-local registries --- python/dolfinx/fem/bcs.py | 12 +++-- python/dolfinx/fem/forms.py | 14 ++++-- python/dolfinx/fem/function.py | 15 ++++-- python/dolfinx/fem/typemap.py | 88 ++++++++-------------------------- 4 files changed, 51 insertions(+), 78 deletions(-) diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index 2be2085590c..733c5c569eb 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -13,7 +13,7 @@ from __future__ import annotations from collections.abc import Callable, Iterable -from typing import Generic +from typing import ClassVar, Generic import numpy as np import numpy.typing as npt @@ -21,7 +21,7 @@ import dolfinx from dolfinx import cpp as _cpp from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import get_cpp_type +from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type from dolfinx.typing import Scalar @@ -99,6 +99,7 @@ class DirichletBC(Generic[Scalar]): The conditions are imposed on a linear system. """ + _cpp_registry: ClassVar[dict] = {} _cpp_object: ( _cpp.fem.DirichletBC_complex64 | _cpp.fem.DirichletBC_complex128 @@ -212,7 +213,7 @@ def dirichletbc( geometry_dtype = np.dtype(value.function_space.mesh.geometry.x.dtype) else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - bctype = get_cpp_type("DirichletBC", dtype, geometry_dtype) + bctype = get_cpp_type(DirichletBC, dtype, geometry_dtype) # Unwrap value object, if required if isinstance(value, np.ndarray): @@ -250,3 +251,8 @@ def _bc_space(V, bcs): return [bc for bc in bcs if V.contains(bc.function_space)] return [_bc_space(V, bcs) if V is not None else [] for V in spaces] + + +# Register the matched-precision built-ins. +for _scalar, _geom, _name in MATCHED_PRECISIONS: + register_cpp_type(DirichletBC, _scalar, _geom, getattr(_cpp.fem, f"DirichletBC_{_name}")) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 76774f2f77e..31b4f43ec0c 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -25,7 +25,7 @@ from dolfinx import default_scalar_type, jit from dolfinx.fem import IntegralType from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import get_cpp_type +from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type from dolfinx.typing import Scalar if typing.TYPE_CHECKING: @@ -38,6 +38,7 @@ class Form(typing.Generic[Scalar]): """A finite element form.""" + _cpp_registry: typing.ClassVar[dict] = {} _cpp_object: ( _cpp.fem.Form_complex64 | _cpp.fem.Form_complex128 @@ -296,7 +297,7 @@ def mixed_topology_form( comm = mesh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = get_cpp_type("Form", dtype, mesh.geometry.x.dtype) + ftype = get_cpp_type(Form, dtype, mesh.geometry.x.dtype) ufcx_forms = [] modules = [] @@ -384,7 +385,7 @@ def _form(form): comm = msh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = get_cpp_type("Form", dtype, msh.geometry.x.dtype) + ftype = get_cpp_type(Form, dtype, msh.geometry.x.dtype) ufcx_form, module, code = jit.ffcx_jit( comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options @@ -449,7 +450,7 @@ def _zero_form(form): msh = V[0].mesh # Geometry type is fixed by the mesh. - ftype = get_cpp_type("Form", dtype, msh.geometry.x.dtype) + ftype = get_cpp_type(Form, dtype, msh.geometry.x.dtype) f = ftype( spaces=V, @@ -814,3 +815,8 @@ def derivative_block( "F must be either a UFL form (with rank zero or one), or a sequence of " "rank-one UFL forms." ) + + +# Register the matched-precision built-ins. +for _scalar, _geom, _name in MATCHED_PRECISIONS: + register_cpp_type(Form, _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index ca0a47ddbbc..23159c6b7a9 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -22,7 +22,7 @@ from dolfinx import default_scalar_type, jit, la from dolfinx.fem.dofmap import DofMap from dolfinx.fem.element import FiniteElement, finiteelement -from dolfinx.fem.typemap import get_cpp_type +from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type from dolfinx.geometry import PointOwnershipData from dolfinx.typing import Real, Scalar @@ -113,6 +113,7 @@ class Expression(Generic[Scalar]): """ + _cpp_registry: typing.ClassVar[dict] = {} _ufl_expression: ufl.core.expr.Expr _argument_space: FunctionSpace | None _cpp_object: ( @@ -212,7 +213,7 @@ def __init__( geometry_dtype = expr_domains[0].ufl_cargo().geometry.x.dtype else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - create_expression = get_cpp_type("Expression", dtype, geometry_dtype) + create_expression = get_cpp_type(Expression, dtype, geometry_dtype) _entity_maps = ( [entity_map._cpp_object for entity_map in entity_maps] @@ -338,6 +339,7 @@ class Function(ufl.Coefficient, Generic[Scalar]): """ + _cpp_registry: typing.ClassVar[dict] = {} _cpp_object: ( _cpp.fem.Function_complex64 | _cpp.fem.Function_complex128 @@ -376,7 +378,7 @@ def __init__( # Scalar type (dtype) is independent of the geometry type, which # is fixed by the mesh. geometry_dtype = V.mesh.geometry.x.dtype - cpp_type = get_cpp_type("Function", dtype, geometry_dtype) + cpp_type = get_cpp_type(Function, dtype, geometry_dtype) if x is not None: self._cpp_object = cpp_type(V._cpp_object, x._cpp_object) # type: ignore else: @@ -848,3 +850,10 @@ def tabulate_dof_coordinates(self) -> npt.NDArray[Real]: degrees-of-freedom. """ return self._cpp_object.tabulate_dof_coordinates() + + +# Register the matched-precision built-ins. ``Expression`` is built +# through a factory, ``Function`` through its class constructor. +for _scalar, _geom, _name in MATCHED_PRECISIONS: + register_cpp_type(Function, _scalar, _geom, getattr(_cpp.fem, f"Function_{_name}")) + register_cpp_type(Expression, _scalar, _geom, getattr(_cpp.fem, f"create_expression_{_name}")) diff --git a/python/dolfinx/fem/typemap.py b/python/dolfinx/fem/typemap.py index 81d5636ea01..89969720f2d 100644 --- a/python/dolfinx/fem/typemap.py +++ b/python/dolfinx/fem/typemap.py @@ -3,89 +3,41 @@ # This file is part of DOLFINx (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -"""Registry from (kind, scalar dtype, geometry dtype) to compiled type. - -``Function``, ``Expression``, ``Form`` and ``DirichletBC`` are templated -on both a scalar and a mesh geometry type. DOLFINx registers the -matched-precision pairs; downstream libraries register further pairs -(e.g. fp16 dofs on an fp64 mesh) via :func:`register_cpp_type`. +"""Registry from (scalar dtype, geometry dtype) to compiled C++ type. + +``Function``, ``Expression``, ``Form`` and ``DirichletBC`` each own a +static ``_cpp_registry`` dict (bound to the class) mapping a +``(scalar, geometry)`` dtype pair to the compiled C++ type or factory. +DOLFINx registers the matched-precision pairs; downstream libraries add +further pairs (e.g. fp16 dofs on an fp64 mesh) with +:func:`register_cpp_type`. """ -from __future__ import annotations - -import typing - import numpy as np -import numpy.typing as npt - -from dolfinx import cpp as _cpp - -__all__ = ["get_cpp_type", "register_cpp_type"] - -# Object kinds templated on (scalar, geometry). -_KINDS = ("Function", "Expression", "Form", "DirichletBC") -# (kind, scalar dtype, geometry dtype) -> compiled C++ class or factory. -_registry: dict[tuple[str, np.dtype, np.dtype], typing.Any] = {} +def register_cpp_type(cls, scalar_dtype, geometry_dtype, cpp_type): + """Register ``cpp_type`` on ``cls`` for a (scalar, geometry) pair.""" + cls._cpp_registry[np.dtype(scalar_dtype), np.dtype(geometry_dtype)] = cpp_type -def register_cpp_type( - kind: str, - scalar_dtype: npt.DTypeLike, - geometry_dtype: npt.DTypeLike, - cpp_type: typing.Any, -) -> None: - """Register the compiled C++ type for an object kind and dtype pair. - Args: - kind: Object kind, one of ``"Function"``, ``"Expression"``, - ``"Form"`` or ``"DirichletBC"``. - scalar_dtype: Scalar (degree-of-freedom) type. - geometry_dtype: Mesh geometry type. - cpp_type: Compiled C++ class or factory that the wrapper builds - with for this pair. - """ - if kind not in _KINDS: - raise ValueError(f"Unknown object kind '{kind}'. Expected one of {_KINDS}.") - _registry[(kind, np.dtype(scalar_dtype), np.dtype(geometry_dtype))] = cpp_type - - -def get_cpp_type( - kind: str, scalar_dtype: npt.DTypeLike, geometry_dtype: npt.DTypeLike -) -> typing.Any: - """Look up the compiled C++ type for an object kind and dtype pair. - - Args: - kind: Object kind (see :func:`register_cpp_type`). - scalar_dtype: Scalar (degree-of-freedom) type. - geometry_dtype: Mesh geometry type. - - Returns: - The registered C++ class or factory. - """ - scalar, geometry = np.dtype(scalar_dtype), np.dtype(geometry_dtype) +def get_cpp_type(cls, scalar_dtype, geometry_dtype): + """Compiled C++ type registered on ``cls`` for a dtype pair.""" + key = np.dtype(scalar_dtype), np.dtype(geometry_dtype) try: - return _registry[(kind, scalar, geometry)] + return cls._cpp_registry[key] except KeyError: raise NotImplementedError( - f"No {kind} type registered for scalar dtype '{scalar}' on geometry " - f"dtype '{geometry}'. Downstream libraries can register one with " - f"dolfinx.fem.register_cpp_type()." + f"No compiled {cls.__name__} for scalar '{key[0]}' on geometry " + f"'{key[1]}'. Register one with dolfinx.fem.register_cpp_type()." ) from None -# Register the matched-precision built-ins (geometry == real part of -# scalar). ``Expression`` is built through a factory, the others through -# their class constructor. -_MATCHED: tuple[tuple[npt.DTypeLike, npt.DTypeLike, str], ...] = ( +# Matched-precision built-ins (geometry == real part of scalar) as +# (scalar, geometry, C++ suffix). +MATCHED_PRECISIONS = ( (np.float32, np.float32, "float32"), (np.float64, np.float64, "float64"), (np.complex64, np.float32, "complex64"), (np.complex128, np.float64, "complex128"), ) - -for _scalar, _geom, _name in _MATCHED: - register_cpp_type("Function", _scalar, _geom, getattr(_cpp.fem, f"Function_{_name}")) - register_cpp_type("Expression", _scalar, _geom, getattr(_cpp.fem, f"create_expression_{_name}")) - register_cpp_type("Form", _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) - register_cpp_type("DirichletBC", _scalar, _geom, getattr(_cpp.fem, f"DirichletBC_{_name}")) From a39eeae74b822fde30cd801bdc57ba3e85ecbebd Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Sun, 12 Jul 2026 20:01:05 +0200 Subject: [PATCH 13/20] Simplify --- python/dolfinx/fem/bcs.py | 11 ++++++++--- python/dolfinx/fem/forms.py | 11 ++++++++--- python/dolfinx/fem/function.py | 18 +++++++++++++----- python/dolfinx/fem/typemap.py | 10 ---------- 4 files changed, 29 insertions(+), 21 deletions(-) diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index 733c5c569eb..9dd47472a33 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -21,7 +21,7 @@ import dolfinx from dolfinx import cpp as _cpp from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type +from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.typing import Scalar @@ -253,6 +253,11 @@ def _bc_space(V, bcs): return [_bc_space(V, bcs) if V is not None else [] for V in spaces] -# Register the matched-precision built-ins. -for _scalar, _geom, _name in MATCHED_PRECISIONS: +# Register the matched-precision built-ins (geometry == real scalar part). +for _scalar, _geom, _name in ( + (np.float32, np.float32, "float32"), + (np.float64, np.float64, "float64"), + (np.complex64, np.float32, "complex64"), + (np.complex128, np.float64, "complex128"), +): register_cpp_type(DirichletBC, _scalar, _geom, getattr(_cpp.fem, f"DirichletBC_{_name}")) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 31b4f43ec0c..60dfd026c91 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -25,7 +25,7 @@ from dolfinx import default_scalar_type, jit from dolfinx.fem import IntegralType from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type +from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.typing import Scalar if typing.TYPE_CHECKING: @@ -817,6 +817,11 @@ def derivative_block( ) -# Register the matched-precision built-ins. -for _scalar, _geom, _name in MATCHED_PRECISIONS: +# Register the matched-precision built-ins (geometry == real scalar part). +for _scalar, _geom, _name in ( + (np.float32, np.float32, "float32"), + (np.float64, np.float64, "float64"), + (np.complex64, np.float32, "complex64"), + (np.complex128, np.float64, "complex128"), +): register_cpp_type(Form, _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 23159c6b7a9..ad3658f753a 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -22,7 +22,7 @@ from dolfinx import default_scalar_type, jit, la from dolfinx.fem.dofmap import DofMap from dolfinx.fem.element import FiniteElement, finiteelement -from dolfinx.fem.typemap import MATCHED_PRECISIONS, get_cpp_type, register_cpp_type +from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.geometry import PointOwnershipData from dolfinx.typing import Real, Scalar @@ -210,7 +210,9 @@ def __init__( # Geometry type is fixed by the expression's mesh. expr_domains = ufl.domain.extract_domains(e) if len(expr_domains) > 0: - geometry_dtype = expr_domains[0].ufl_cargo().geometry.x.dtype + expr_domain = expr_domains[0] + assert isinstance(expr_domain, ufl.Mesh) + geometry_dtype = expr_domain.ufl_cargo().geometry.x.dtype else: geometry_dtype = np.dtype(dtype).type(0).real.dtype create_expression = get_cpp_type(Expression, dtype, geometry_dtype) @@ -852,8 +854,14 @@ def tabulate_dof_coordinates(self) -> npt.NDArray[Real]: return self._cpp_object.tabulate_dof_coordinates() -# Register the matched-precision built-ins. ``Expression`` is built -# through a factory, ``Function`` through its class constructor. -for _scalar, _geom, _name in MATCHED_PRECISIONS: +# Register the matched-precision built-ins (geometry == real part of +# scalar). ``Expression`` is built through a factory, ``Function`` +# through its class constructor. +for _scalar, _geom, _name in ( + (np.float32, np.float32, "float32"), + (np.float64, np.float64, "float64"), + (np.complex64, np.float32, "complex64"), + (np.complex128, np.float64, "complex128"), +): register_cpp_type(Function, _scalar, _geom, getattr(_cpp.fem, f"Function_{_name}")) register_cpp_type(Expression, _scalar, _geom, getattr(_cpp.fem, f"create_expression_{_name}")) diff --git a/python/dolfinx/fem/typemap.py b/python/dolfinx/fem/typemap.py index 89969720f2d..3830507330d 100644 --- a/python/dolfinx/fem/typemap.py +++ b/python/dolfinx/fem/typemap.py @@ -31,13 +31,3 @@ def get_cpp_type(cls, scalar_dtype, geometry_dtype): f"No compiled {cls.__name__} for scalar '{key[0]}' on geometry " f"'{key[1]}'. Register one with dolfinx.fem.register_cpp_type()." ) from None - - -# Matched-precision built-ins (geometry == real part of scalar) as -# (scalar, geometry, C++ suffix). -MATCHED_PRECISIONS = ( - (np.float32, np.float32, "float32"), - (np.float64, np.float64, "float64"), - (np.complex64, np.float32, "complex64"), - (np.complex128, np.float64, "complex128"), -) From 9778283d295cbf10f6d31bb2e3ca402e232d5c17 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 14 Jul 2026 16:06:22 +0200 Subject: [PATCH 14/20] Assert when extracting from ufcx structs --- cpp/dolfinx/fem/interpolate.h | 4 ++-- cpp/dolfinx/fem/kernel.h | 3 +++ cpp/dolfinx/fem/utils.h | 36 ++++++++++++++--------------------- 3 files changed, 19 insertions(+), 24 deletions(-) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index cbe7a5ad0ab..85ec5cc3190 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -423,7 +423,7 @@ void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, auto [i_m, im_shape] = element1->create_interpolation_operator(*element0); // Iterate over mesh and interpolate on each cell - using X = scalar_value_t; + using X = U; // geometry (real) type, independent of the value scalar T assert(cells0.size() == cells1.size()); for (auto cell0_it = cells0.begin(), cell1_it = cells1.begin(); cell0_it != cells0.end() and cell1_it != cells1.end(); @@ -665,7 +665,7 @@ void interpolate_nonmatching_maps(Function& u1, coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k]; // Evaluate v at the interpolation points (physical space values) - using X = scalar_value_t; + using X = U; // geometry (real) type, independent of the value scalar T for (std::size_t p = 0; p < Xshape[0]; ++p) { for (int k = 0; k < bs0; ++k) diff --git a/cpp/dolfinx/fem/kernel.h b/cpp/dolfinx/fem/kernel.h index 3c949864ca1..a889043db0f 100644 --- a/cpp/dolfinx/fem/kernel.h +++ b/cpp/dolfinx/fem/kernel.h @@ -34,6 +34,9 @@ using kernel_t = std::function> constexpr kernel_t extract_kernel(const ufcx_integral* integral) { + static_assert(std::is_same_v>, + "UFCx kernels require geometry type U == scalar_value_t."); + if constexpr (std::is_same_v) return integral->tabulate_tensor_float32; else if constexpr (std::is_same_v) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index cdac727a217..292d9f07a6e 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -970,32 +970,24 @@ Expression create_expression( static_cast(e.entity_dimension)}; std::vector value_shape(e.value_shape, e.value_shape + e.num_components); + + static_assert(std::is_same_v>, + "UFCx kernels require geometry type U == scalar_value_t."); + + using kptr_t = void (*)(T*, const T*, const T*, const U*, const int*, + const std::uint8_t*, void*); std::function tabulate_tensor = nullptr; - if constexpr (std::is_same_v && std::is_same_v) - tabulate_tensor = e.tabulate_tensor_float32; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v> - && std::is_same_v) - { - tabulate_tensor - = reinterpret_cast( - e.tabulate_tensor_complex64); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v && std::is_same_v) - tabulate_tensor = e.tabulate_tensor_float64; + if constexpr (std::is_same_v) + tabulate_tensor = reinterpret_cast(e.tabulate_tensor_float32); + else if constexpr (std::is_same_v) + tabulate_tensor = reinterpret_cast(e.tabulate_tensor_float64); #ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v> - && std::is_same_v) - { - tabulate_tensor - = reinterpret_cast( - e.tabulate_tensor_complex128); - } + else if constexpr (std::is_same_v>) + tabulate_tensor = reinterpret_cast(e.tabulate_tensor_complex64); + else if constexpr (std::is_same_v>) + tabulate_tensor = reinterpret_cast(e.tabulate_tensor_complex128); #endif // DOLFINX_NO_STDC_COMPLEX_KERNELS else throw std::runtime_error("Type not supported."); From 569b0b855935fa9c6aaa5581738e454b8f7348ff Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 14 Jul 2026 16:11:17 +0200 Subject: [PATCH 15/20] Lint --- cpp/dolfinx/fem/interpolate.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index 85ec5cc3190..87434c4ff3f 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -423,7 +423,7 @@ void interpolate_same_map(Function& u1, mesh::CellRange auto&& cells1, auto [i_m, im_shape] = element1->create_interpolation_operator(*element0); // Iterate over mesh and interpolate on each cell - using X = U; // geometry (real) type, independent of the value scalar T + using X = U; // geometry (real) type, independent of the value scalar T assert(cells0.size() == cells1.size()); for (auto cell0_it = cells0.begin(), cell1_it = cells1.begin(); cell0_it != cells0.end() and cell1_it != cells1.end(); @@ -665,7 +665,7 @@ void interpolate_nonmatching_maps(Function& u1, coeffs0[dof_bs0 * i + k] = array0[dof_bs0 * dofs0[i] + k]; // Evaluate v at the interpolation points (physical space values) - using X = U; // geometry (real) type, independent of the value scalar T + using X = U; // geometry (real) type, independent of the value scalar T for (std::size_t p = 0; p < Xshape[0]; ++p) { for (int k = 0; k < bs0; ++k) From 3f63447660ac8720b259c04a220bfd76bc24a945 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 14 Jul 2026 16:28:14 +0200 Subject: [PATCH 16/20] Guard ufcx stricts in nanobind wrappers --- .../dolfinx/wrappers/dolfinx_wrappers/fem.h | 64 +++++++++++++------ 1 file changed, 44 insertions(+), 20 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 226a7c29edf..7d48bd2d369 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -721,11 +721,17 @@ void declare_objects(nb::module_& m, std::string type) const std::vector& entity_maps, std::shared_ptr> argument_space) { - const ufcx_expression* p - = reinterpret_cast(expression); - return dolfinx::fem::create_expression( - *p, coefficients, constants, ptr_to_ref_wrapper_vec(entity_maps), - argument_space); + if constexpr (std::is_same_v>) + { + const ufcx_expression* p + = reinterpret_cast(expression); + return dolfinx::fem::create_expression( + *p, coefficients, constants, ptr_to_ref_wrapper_vec(entity_maps), + argument_space); + } + else + throw std::runtime_error( + "create_expression requires geometry type scalar_value_t."); }, nb::arg("expression"), nb::arg("coefficients"), nb::arg("constants"), nb::arg("entity_maps"), nb::arg("argument_space").none(), @@ -818,14 +824,20 @@ void declare_form(nb::module_& m, std::string type) sd.insert({itg, std::move(x)}); } - std::vector> ps; - ps.reserve(forms.size()); - 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, - ptr_to_ref_wrapper_vec(entity_maps), mesh)); + if constexpr (std::is_same_v>) + { + std::vector> ps; + ps.reserve(forms.size()); + 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, + ptr_to_ref_wrapper_vec(entity_maps), mesh)); + } + else + throw std::runtime_error( + "create_form requires geometry type scalar_value_t."); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), nb::arg("constants"), nb::arg("subdomains"), nb::arg("entity_maps"), @@ -899,9 +911,15 @@ void declare_form(nb::module_& m, std::string type) sd.insert({itg, std::move(x)}); } - ufcx_form* p = reinterpret_cast(form); - return dolfinx::fem::create_form_factory( - {*p}, spaces, coefficients, constants, sd, {}, std::move(mesh)); + if constexpr (std::is_same_v>) + { + ufcx_form* p = reinterpret_cast(form); + return dolfinx::fem::create_form_factory( + {*p}, spaces, coefficients, constants, sd, {}, std::move(mesh)); + } + else + throw std::runtime_error( + "create_form requires geometry type scalar_value_t."); }, nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), nb::arg("constants"), nb::arg("subdomains"), nb::arg("mesh"), @@ -937,10 +955,16 @@ void declare_form(nb::module_& m, std::string type) sd.insert({itg, std::move(x)}); } - ufcx_form* p = reinterpret_cast(form); - return dolfinx::fem::create_form( - *p, spaces, coefficients, constants, sd, - ptr_to_ref_wrapper_vec(entity_maps), std::move(mesh)); + if constexpr (std::is_same_v>) + { + ufcx_form* p = reinterpret_cast(form); + return dolfinx::fem::create_form( + *p, spaces, coefficients, constants, sd, + ptr_to_ref_wrapper_vec(entity_maps), std::move(mesh)); + } + else + throw std::runtime_error( + "create_form requires geometry type scalar_value_t."); }, // NOLINTEND(performance-no-int-to-ptr) nb::arg("form"), nb::arg("spaces"), nb::arg("coefficients"), From 4028fe30c969fe3d30bf7051d3477e408467d16e Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 16 Jul 2026 16:49:25 +0200 Subject: [PATCH 17/20] Re-route dtype dispatches through registry only --- python/dolfinx/fem/__init__.py | 3 - python/dolfinx/fem/bcs.py | 37 ++++++--- python/dolfinx/fem/element.py | 65 +++++++++++---- python/dolfinx/fem/forms.py | 54 +++++++------ python/dolfinx/fem/function.py | 143 +++++++++++++++++++++++---------- python/dolfinx/fem/typemap.py | 33 -------- 6 files changed, 202 insertions(+), 133 deletions(-) delete mode 100644 python/dolfinx/fem/typemap.py diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index 9e9ec025548..7b45c04735b 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -48,7 +48,6 @@ FunctionSpace, functionspace, ) -from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.fem.utils import ( build_sparsity_pattern, compute_integration_domains, @@ -95,7 +94,6 @@ "form", "form_cpp_class", "functionspace", - "get_cpp_type", "interpolate_geometry", "interpolation_matrix", "locate_dofs_geometrical", @@ -103,7 +101,6 @@ "mixed_topology_form", "pack_coefficients", "pack_constants", - "register_cpp_type", "set_bc", "transpose_dofmap", ] diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index 92452e70a1f..d3534e3f24d 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -21,7 +21,6 @@ import dolfinx from dolfinx import cpp as _cpp from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.typing import Scalar @@ -99,7 +98,13 @@ class DirichletBC(Generic[Scalar]): The conditions are imposed on a linear system. """ - _cpp_registry: ClassVar[dict] = {} + # Matched-precision built-ins (geometry == real scalar part). + _cpp_registry: ClassVar[dict] = { + (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.DirichletBC_float32, + (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.DirichletBC_float64, + (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.DirichletBC_complex64, + (np.dtype(np.complex128), np.dtype(np.float64)): _cpp.fem.DirichletBC_complex128, + } _cpp_object: ( _cpp.fem.DirichletBC_complex64 | _cpp.fem.DirichletBC_complex128 @@ -107,6 +112,22 @@ class DirichletBC(Generic[Scalar]): | _cpp.fem.DirichletBC_float64 ) + @classmethod + def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): + """Register a compiled C++ type for a (scalar, geometry) pair.""" + cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, scalar_dtype, geometry_dtype): + """Compiled C++ type registered for a (scalar, geometry) pair.""" + try: + return cls._cpp_registry[scalar_dtype, geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " + f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__(self, bc): """Initialise a Dirichlet boundary condition. @@ -213,7 +234,7 @@ def dirichletbc( geometry_dtype = np.dtype(value.function_space.mesh.geometry.x.dtype) else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - bctype = get_cpp_type(DirichletBC, dtype, geometry_dtype) + bctype = DirichletBC.get_cpp_type(dtype, geometry_dtype) # Unwrap value object, if required if isinstance(value, np.ndarray): @@ -251,13 +272,3 @@ def _bc_space(V, bcs): return [bc for bc in bcs if V.contains(bc.function_space)] return [_bc_space(V, bcs) if V is not None else [] for V in spaces] - - -# Register the matched-precision built-ins (geometry == real scalar part). -for _scalar, _geom, _name in ( - (np.float32, np.float32, "float32"), - (np.float64, np.float64, "float64"), - (np.complex64, np.float32, "complex64"), - (np.complex128, np.float64, "complex128"), -): - register_cpp_type(DirichletBC, _scalar, _geom, getattr(_cpp.fem, f"DirichletBC_{_name}")) diff --git a/python/dolfinx/fem/element.py b/python/dolfinx/fem/element.py index 1d82d341aff..181c78aee42 100644 --- a/python/dolfinx/fem/element.py +++ b/python/dolfinx/fem/element.py @@ -6,7 +6,7 @@ """Finite elements.""" from functools import singledispatch -from typing import Generic +from typing import ClassVar, Generic import numpy as np import numpy.typing as npt @@ -20,8 +20,29 @@ class CoordinateElement(Generic[Real]): """Coordinate element describing the geometry map for mesh cells.""" + # Built-in geometry types. + _cpp_registry: ClassVar[dict] = { + np.dtype(np.float32): _cpp.fem.CoordinateElement_float32, + np.dtype(np.float64): _cpp.fem.CoordinateElement_float64, + } _cpp_object: _cpp.fem.CoordinateElement_float32 | _cpp.fem.CoordinateElement_float64 + @classmethod + def register_cpp_type(cls, cpp_type, geometry_dtype): + """Register a compiled C++ type for a geometry dtype.""" + cls._cpp_registry[geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, geometry_dtype): + """Compiled C++ type registered for a geometry dtype.""" + try: + return cls._cpp_registry[geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " + f"Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, cmap: _cpp.fem.CoordinateElement_float32 | _cpp.fem.CoordinateElement_float64 ): @@ -144,12 +165,8 @@ def coordinate_element( Returns: A coordinate element. """ - if np.issubdtype(dtype, np.float32): - return CoordinateElement(_cpp.fem.CoordinateElement_float32(celltype, degree, variant)) - elif np.issubdtype(dtype, np.float64): - return CoordinateElement(_cpp.fem.CoordinateElement_float64(celltype, degree, variant)) - else: - raise RuntimeError("Unsupported dtype.") + cpp_type = CoordinateElement.get_cpp_type(np.dtype(dtype)) + return CoordinateElement(cpp_type(celltype, degree, variant)) @coordinate_element.register(basix.finite_element.FiniteElement) @@ -164,17 +181,36 @@ def _(e: basix.finite_element.FiniteElement) -> CoordinateElement: Returns: A coordinate element. """ - try: - return CoordinateElement(_cpp.fem.CoordinateElement_float32(e._e)) - except TypeError: - return CoordinateElement(_cpp.fem.CoordinateElement_float64(e._e)) + cpp_type = CoordinateElement.get_cpp_type(e.dtype) + return CoordinateElement(cpp_type(e._e)) class FiniteElement(Generic[Real]): """A finite element.""" + # Built-in geometry types. + _cpp_registry: ClassVar[dict] = { + np.dtype(np.float32): _cpp.fem.FiniteElement_float32, + np.dtype(np.float64): _cpp.fem.FiniteElement_float64, + } _cpp_object: _cpp.fem.FiniteElement_float32 | _cpp.fem.FiniteElement_float64 + @classmethod + def register_cpp_type(cls, cpp_type, geometry_dtype): + """Register a compiled C++ type for a geometry dtype.""" + cls._cpp_registry[geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, geometry_dtype): + """Compiled C++ type registered for a geometry dtype.""" + try: + return cls._cpp_registry[geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " + f"Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, cpp_object: _cpp.fem.FiniteElement_float32 | _cpp.fem.FiniteElement_float64, @@ -346,12 +382,7 @@ def finiteelement( the selected element. FiniteElement_dtype: Geometry type of the element. """ - if np.issubdtype(FiniteElement_dtype, np.float32): - CppElement = _cpp.fem.FiniteElement_float32 - elif np.issubdtype(FiniteElement_dtype, np.float64): - CppElement = _cpp.fem.FiniteElement_float64 - else: - raise ValueError(f"Unsupported dtype: {FiniteElement_dtype}") + CppElement = FiniteElement.get_cpp_type(np.dtype(FiniteElement_dtype)) if ufl_e.is_mixed: elements = [ diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 69cc2effc48..e5c8845a917 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -25,7 +25,6 @@ from dolfinx import default_scalar_type, jit from dolfinx.fem import IntegralType from dolfinx.fem.function import Constant, Function, FunctionSpace -from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.typing import Scalar if typing.TYPE_CHECKING: @@ -38,7 +37,13 @@ class Form(typing.Generic[Scalar]): """A finite element form.""" - _cpp_registry: typing.ClassVar[dict] = {} + # Matched-precision built-ins (geometry == real scalar part). + _cpp_registry: typing.ClassVar[dict] = { + (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.Form_float32, + (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.Form_float64, + (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.Form_complex64, + (np.dtype(np.complex128), np.dtype(np.float64)): _cpp.fem.Form_complex128, + } _cpp_object: ( _cpp.fem.Form_complex64 | _cpp.fem.Form_complex128 @@ -47,6 +52,22 @@ class Form(typing.Generic[Scalar]): ) _code: str | list[str] | None + @classmethod + def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): + """Register a compiled C++ type for a (scalar, geometry) pair.""" + cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, scalar_dtype, geometry_dtype): + """Compiled C++ type registered for a (scalar, geometry) pair.""" + try: + return cls._cpp_registry[scalar_dtype, geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " + f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, form: _cpp.fem.Form_complex64 @@ -224,16 +245,9 @@ def form_cpp_class( This function is for advanced usage, typically when writing custom kernels using Numba or C. """ - if np.issubdtype(dtype, np.float32): - return _cpp.fem.Form_float32 - elif np.issubdtype(dtype, np.float64): - return _cpp.fem.Form_float64 - elif np.issubdtype(dtype, np.complex64): - return _cpp.fem.Form_complex64 - elif np.issubdtype(dtype, np.complex128): - return _cpp.fem.Form_complex128 - else: - raise NotImplementedError(f"Type {dtype} not supported.") + # Geometry is the real part of the scalar type (matched precision). + scalar_dtype = np.dtype(dtype) + return Form.get_cpp_type(scalar_dtype, scalar_dtype.type(0).real.dtype) _ufl_to_dolfinx_domain = { @@ -297,7 +311,7 @@ def mixed_topology_form( comm = mesh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = get_cpp_type(Form, dtype, mesh.geometry.x.dtype) + ftype = Form.get_cpp_type(np.dtype(dtype), mesh.geometry.x.dtype) ufcx_forms = [] modules = [] @@ -385,7 +399,7 @@ def _form(form): comm = msh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = get_cpp_type(Form, dtype, msh.geometry.x.dtype) + ftype = Form.get_cpp_type(np.dtype(dtype), msh.geometry.x.dtype) ufcx_form, module, code = jit.ffcx_jit( comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options @@ -450,7 +464,7 @@ def _zero_form(form): msh = V[0].mesh # Geometry type is fixed by the mesh. - ftype = get_cpp_type(Form, dtype, msh.geometry.x.dtype) + ftype = Form.get_cpp_type(np.dtype(dtype), msh.geometry.x.dtype) f = ftype( spaces=V, @@ -815,13 +829,3 @@ def derivative_block( "F must be either a UFL form (with rank zero or one), or a sequence of " "rank-one UFL forms." ) - - -# Register the matched-precision built-ins (geometry == real scalar part). -for _scalar, _geom, _name in ( - (np.float32, np.float32, "float32"), - (np.float64, np.float64, "float64"), - (np.complex64, np.float32, "complex64"), - (np.complex128, np.float64, "complex128"), -): - register_cpp_type(Form, _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 11c967b944c..759528df915 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -22,7 +22,6 @@ from dolfinx import default_scalar_type, jit, la from dolfinx.fem.dofmap import DofMap from dolfinx.fem.element import FiniteElement, finiteelement -from dolfinx.fem.typemap import get_cpp_type, register_cpp_type from dolfinx.geometry import PointOwnershipData from dolfinx.typing import Real, Scalar @@ -36,6 +35,13 @@ class Constant(ufl.Constant, Generic[Scalar]): """A constant with respect to a domain.""" + # Built-in scalar types. + _cpp_registry: typing.ClassVar[dict] = { + np.dtype(np.float32): _cpp.fem.Constant_float32, + np.dtype(np.float64): _cpp.fem.Constant_float64, + np.dtype(np.complex64): _cpp.fem.Constant_complex64, + np.dtype(np.complex128): _cpp.fem.Constant_complex128, + } _cpp_object: ( _cpp.fem.Constant_complex64 | _cpp.fem.Constant_complex128 @@ -43,6 +49,22 @@ class Constant(ufl.Constant, Generic[Scalar]): | _cpp.fem.Constant_float64 ) + @classmethod + def register_cpp_type(cls, cpp_type, scalar_dtype): + """Register a compiled C++ type for a scalar dtype.""" + cls._cpp_registry[scalar_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, scalar_dtype): + """Compiled C++ type registered for a scalar dtype.""" + try: + return cls._cpp_registry[scalar_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for scalar '{scalar_dtype}'. " + f"Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, domain, @@ -57,18 +79,10 @@ def __init__( c = np.asarray(c) super().__init__(domain, c.shape) try: - if np.issubdtype(c.dtype, np.complex64): - self._cpp_object = _cpp.fem.Constant_complex64(c) - elif np.issubdtype(c.dtype, np.complex128): - self._cpp_object = _cpp.fem.Constant_complex128(c) - elif np.issubdtype(c.dtype, np.float32): - self._cpp_object = _cpp.fem.Constant_float32(c) - elif np.issubdtype(c.dtype, np.float64): - self._cpp_object = _cpp.fem.Constant_float64(c) - else: - raise RuntimeError("Unsupported dtype") + cpp_type = Constant.get_cpp_type(c.dtype) except AttributeError as err: raise AttributeError("Constant value must have a dtype attribute.") from err + self._cpp_object = cpp_type(c) @property def value(self): @@ -113,7 +127,14 @@ class Expression(Generic[Scalar]): """ - _cpp_registry: typing.ClassVar[dict] = {} + # Matched-precision built-ins (geometry == real scalar part). An + # ``Expression`` is built through a factory, not a class constructor. + _cpp_registry: typing.ClassVar[dict] = { + (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.create_expression_float32, + (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.create_expression_float64, + (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.create_expression_complex64, + (np.dtype(np.complex128), np.dtype(np.float64)): _cpp.fem.create_expression_complex128, + } _ufl_expression: ufl.core.expr.Expr _argument_space: FunctionSpace | None _cpp_object: ( @@ -124,6 +145,22 @@ class Expression(Generic[Scalar]): ) _code: str + @classmethod + def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): + """Register a C++ factory for a (scalar, geometry) pair.""" + cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, scalar_dtype, geometry_dtype): + """C++ factory registered for a (scalar, geometry) pair.""" + try: + return cls._cpp_registry[scalar_dtype, geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " + f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, e: ufl.core.expr.Expr, @@ -215,7 +252,7 @@ def __init__( geometry_dtype = expr_domain.ufl_cargo().geometry.x.dtype else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - create_expression = get_cpp_type(Expression, dtype, geometry_dtype) + create_expression = Expression.get_cpp_type(np.dtype(dtype), geometry_dtype) _entity_maps = ( [entity_map._cpp_object for entity_map in entity_maps] @@ -341,7 +378,13 @@ class Function(ufl.Coefficient, Generic[Scalar]): """ - _cpp_registry: typing.ClassVar[dict] = {} + # Matched-precision built-ins (geometry == real scalar part). + _cpp_registry: typing.ClassVar[dict] = { + (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.Function_float32, + (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.Function_float64, + (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.Function_complex64, + (np.dtype(np.complex128), np.dtype(np.float64)): _cpp.fem.Function_complex128, + } _cpp_object: ( _cpp.fem.Function_complex64 | _cpp.fem.Function_complex128 @@ -351,6 +394,22 @@ class Function(ufl.Coefficient, Generic[Scalar]): _x: la.Vector[Scalar] + @classmethod + def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): + """Register a compiled C++ type for a (scalar, geometry) pair.""" + cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, scalar_dtype, geometry_dtype): + """Compiled C++ type registered for a (scalar, geometry) pair.""" + try: + return cls._cpp_registry[scalar_dtype, geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " + f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, V: FunctionSpace, @@ -380,7 +439,7 @@ def __init__( # Scalar type (dtype) is independent of the geometry type, which # is fixed by the mesh. geometry_dtype = V.mesh.geometry.x.dtype - cpp_type = get_cpp_type(Function, dtype, geometry_dtype) + cpp_type = Function.get_cpp_type(np.dtype(dtype), geometry_dtype) if x is not None: self._cpp_object = cpp_type(V._cpp_object, x._cpp_object) # type: ignore else: @@ -645,7 +704,6 @@ def functionspace( Args: mesh: Mesh that space is defined on. element: Finite element description. - dtype: Scalar type of the element. Returns: A function space. @@ -688,11 +746,9 @@ def functionspace( "Mesh and element dtype are not compatible." ) - # Initialize the cpp.FunctionSpace - try: - cppV = _cpp.fem.FunctionSpace_float64(mesh._cpp_object, element._cpp_object, cpp_dofmap) # type: ignore - except TypeError: - cppV = _cpp.fem.FunctionSpace_float32(mesh._cpp_object, element._cpp_object, cpp_dofmap) # type: ignore + # Initialize the cpp.FunctionSpace. Geometry type is fixed by the mesh. + cpp_type = FunctionSpace.get_cpp_type(dtype) + cppV = cpp_type(mesh._cpp_object, element._cpp_object, cpp_dofmap) # type: ignore return FunctionSpace(mesh, ufl_e, cppV) @@ -700,9 +756,30 @@ def functionspace( class FunctionSpace(ufl.FunctionSpace, Generic[Real]): """A space on which Functions (fields) can be defined.""" + # Built-in geometry types. + _cpp_registry: typing.ClassVar[dict] = { + np.dtype(np.float32): _cpp.fem.FunctionSpace_float32, + np.dtype(np.float64): _cpp.fem.FunctionSpace_float64, + } _cpp_object: _cpp.fem.FunctionSpace_float32 | _cpp.fem.FunctionSpace_float64 _mesh: Mesh[Real] + @classmethod + def register_cpp_type(cls, cpp_type, geometry_dtype): + """Register a compiled C++ type for a geometry dtype.""" + cls._cpp_registry[geometry_dtype] = cpp_type + + @classmethod + def get_cpp_type(cls, geometry_dtype): + """Compiled C++ type registered for a geometry dtype.""" + try: + return cls._cpp_registry[geometry_dtype] + except KeyError: + raise NotImplementedError( + f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " + f"Register one with {cls.__name__}.register_cpp_type()." + ) from None + def __init__( self, mesh: Mesh[Real], @@ -748,14 +825,9 @@ def clone(self) -> FunctionSpace[Real]: Returns: A new function space that shares data """ # noqa: D301 - try: - Vcpp = _cpp.fem.FunctionSpace_float64( - self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap - ) # type: ignore - except TypeError: - Vcpp = _cpp.fem.FunctionSpace_float32( - self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap - ) # type: ignore + # Geometry type is fixed by the mesh. + cpp_type = FunctionSpace.get_cpp_type(self._mesh.geometry.x.dtype) + Vcpp = cpp_type(self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap) # type: ignore return FunctionSpace(self._mesh, self.ufl_element(), Vcpp) @property @@ -852,16 +924,3 @@ def tabulate_dof_coordinates(self) -> npt.NDArray[Real]: degrees-of-freedom. """ return self._cpp_object.tabulate_dof_coordinates() - - -# Register the matched-precision built-ins (geometry == real part of -# scalar). ``Expression`` is built through a factory, ``Function`` -# through its class constructor. -for _scalar, _geom, _name in ( - (np.float32, np.float32, "float32"), - (np.float64, np.float64, "float64"), - (np.complex64, np.float32, "complex64"), - (np.complex128, np.float64, "complex128"), -): - register_cpp_type(Function, _scalar, _geom, getattr(_cpp.fem, f"Function_{_name}")) - register_cpp_type(Expression, _scalar, _geom, getattr(_cpp.fem, f"create_expression_{_name}")) diff --git a/python/dolfinx/fem/typemap.py b/python/dolfinx/fem/typemap.py deleted file mode 100644 index 3830507330d..00000000000 --- a/python/dolfinx/fem/typemap.py +++ /dev/null @@ -1,33 +0,0 @@ -# Copyright (C) 2026 Michal Habera -# -# This file is part of DOLFINx (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -"""Registry from (scalar dtype, geometry dtype) to compiled C++ type. - -``Function``, ``Expression``, ``Form`` and ``DirichletBC`` each own a -static ``_cpp_registry`` dict (bound to the class) mapping a -``(scalar, geometry)`` dtype pair to the compiled C++ type or factory. -DOLFINx registers the matched-precision pairs; downstream libraries add -further pairs (e.g. fp16 dofs on an fp64 mesh) with -:func:`register_cpp_type`. -""" - -import numpy as np - - -def register_cpp_type(cls, scalar_dtype, geometry_dtype, cpp_type): - """Register ``cpp_type`` on ``cls`` for a (scalar, geometry) pair.""" - cls._cpp_registry[np.dtype(scalar_dtype), np.dtype(geometry_dtype)] = cpp_type - - -def get_cpp_type(cls, scalar_dtype, geometry_dtype): - """Compiled C++ type registered on ``cls`` for a dtype pair.""" - key = np.dtype(scalar_dtype), np.dtype(geometry_dtype) - try: - return cls._cpp_registry[key] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{key[0]}' on geometry " - f"'{key[1]}'. Register one with dolfinx.fem.register_cpp_type()." - ) from None From afb3038e60e720baf3bfed104076bad9815e1de8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 16 Jul 2026 17:05:14 +0000 Subject: [PATCH 18/20] Add libaec overlay port to avoid gitlab.dkrz.de rate limiting (HTTP 429) --- cpp/.vcpkg-overlay/README.md | 21 ++++++++--- cpp/.vcpkg-overlay/libaec/portfile.cmake | 44 ++++++++++++++++++++++++ cpp/.vcpkg-overlay/libaec/usage | 7 ++++ cpp/.vcpkg-overlay/libaec/vcpkg.json | 17 +++++++++ 4 files changed, 84 insertions(+), 5 deletions(-) create mode 100644 cpp/.vcpkg-overlay/libaec/portfile.cmake create mode 100644 cpp/.vcpkg-overlay/libaec/usage create mode 100644 cpp/.vcpkg-overlay/libaec/vcpkg.json diff --git a/cpp/.vcpkg-overlay/README.md b/cpp/.vcpkg-overlay/README.md index e8e98953a3b..17672e51f8b 100644 --- a/cpp/.vcpkg-overlay/README.md +++ b/cpp/.vcpkg-overlay/README.md @@ -1,14 +1,25 @@ -# vcpkg overlay port for Intel MPI +# vcpkg overlay ports + +## Intel MPI This vcpkg overlay port contains scripts for installing Intel MPI on Windows (only). MSMPI, which is used by default with vcpkg, does not support the MPI3 standard. Using this port requires that Intel OneAPI binaries are already installed. On Unix systems the built-in OpenMPI or MPICH ports can be used. -From the root of this repository it can be activated by e.g.: - - cmake -DCMAKE_TOOLCHAIN_FILE=%VCPKG_ROOT%/scripts/buildsystems/vcpkg.cmake -DVCPKG_OVERLAY_PORTS="cpp/.vcpkg-overlay" -B build-dir -S cpp/ - This overlay port was adapted from the original at: https://github.com/arcaneframework/framework-ci + +## libaec + +The upstream vcpkg port for `libaec` downloads from `https://gitlab.dkrz.de`, +which can intermittently fail with HTTP 429 (rate limiting) errors in CI. +This overlay port redirects the download to the GitHub mirror at +https://github.com/MathisRosenhauer/libaec. + +## Usage + +From the root of this repository it can be activated by e.g.: + + cmake -DCMAKE_TOOLCHAIN_FILE=%VCPKG_ROOT%/scripts/buildsystems/vcpkg.cmake -DVCPKG_OVERLAY_PORTS="cpp/.vcpkg-overlay" -B build-dir -S cpp/ diff --git a/cpp/.vcpkg-overlay/libaec/portfile.cmake b/cpp/.vcpkg-overlay/libaec/portfile.cmake new file mode 100644 index 00000000000..65c34cc787a --- /dev/null +++ b/cpp/.vcpkg-overlay/libaec/portfile.cmake @@ -0,0 +1,44 @@ +vcpkg_from_github( + OUT_SOURCE_PATH SOURCE_PATH + REPO MathisRosenhauer/libaec + REF "v${VERSION}" + SHA512 55bd605590015e0f903a231268265051336c172c18935fdfecead5630454a99811f3f196117556ad1b62f3052d54894fd8e257a9ea9f45cc03f59c93cde00cda + HEAD_REF master +) + +string(COMPARE EQUAL "${VCPKG_LIBRARY_LINKAGE}" "static" BUILD_STATIC) + +vcpkg_cmake_configure( + SOURCE_PATH "${SOURCE_PATH}" + OPTIONS + -DBUILD_STATIC_LIBS=${BUILD_STATIC} + -Dlibaec_INSTALL_CMAKEDIR=share/${PORT} +) +vcpkg_cmake_install() +vcpkg_copy_pdbs() +vcpkg_cmake_config_fixup() +vcpkg_replace_string("${CURRENT_PACKAGES_DIR}/share/libaec/libaec-config.cmake" + "if(libaec_USE_STATIC_LIBS)" + "if(\"${BUILD_STATIC}\") # forced by vcpkg" +) + +# Compatibility with user's CMake < 3.18 (vcpkg claims support for >= 3.16): +# Make imported targets global so that libaec-config.cmake can create ALIAS targets. +set(_target_file "libaec_shared-targets") +if(BUILD_STATIC) + set(_target_file "libaec_static-targets") +endif() +file(READ "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" libaec_targets) +string(REGEX REPLACE " (SHARED|STATIC) IMPORTED" " \\1 IMPORTED \${libaec_maybe_global}" libaec_targets "${libaec_targets}") +file(WRITE "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" "set(libaec_maybe_global \"\") +if(CMAKE_VERSION VERSION_LESS 3.18) + set(libaec_maybe_global \"GLOBAL\") +endif() +${libaec_targets} +" +) + +file(REMOVE_RECURSE "${CURRENT_PACKAGES_DIR}/debug/include") + +file(INSTALL "${CURRENT_PORT_DIR}/usage" DESTINATION "${CURRENT_PACKAGES_DIR}/share/${PORT}") +vcpkg_install_copyright(FILE_LIST "${SOURCE_PATH}/LICENSE.txt") diff --git a/cpp/.vcpkg-overlay/libaec/usage b/cpp/.vcpkg-overlay/libaec/usage new file mode 100644 index 00000000000..a6aeb5f6a63 --- /dev/null +++ b/cpp/.vcpkg-overlay/libaec/usage @@ -0,0 +1,7 @@ +libaec provides CMake targets: + + find_package(libaec CONFIG REQUIRED) + # libaec API + target_link_libraries(main PRIVATE libaec::aec) + # szip compatible API + target_link_libraries(main PRIVATE libaec::sz) diff --git a/cpp/.vcpkg-overlay/libaec/vcpkg.json b/cpp/.vcpkg-overlay/libaec/vcpkg.json new file mode 100644 index 00000000000..69c73b3559d --- /dev/null +++ b/cpp/.vcpkg-overlay/libaec/vcpkg.json @@ -0,0 +1,17 @@ +{ + "name": "libaec", + "version": "1.1.7", + "description": "Adaptive Entropy Coding library", + "homepage": "https://github.com/MathisRosenhauer/libaec", + "license": "BSD-2-Clause", + "dependencies": [ + { + "name": "vcpkg-cmake", + "host": true + }, + { + "name": "vcpkg-cmake-config", + "host": true + } + ] +} From 01a3ffede59a369521f96c3a8038c8992351b708 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Thu, 16 Jul 2026 19:09:56 +0200 Subject: [PATCH 19/20] Revert "Add libaec overlay port to avoid gitlab.dkrz.de rate limiting (HTTP 429)" This reverts commit afb3038e60e720baf3bfed104076bad9815e1de8. --- cpp/.vcpkg-overlay/README.md | 21 +++-------- cpp/.vcpkg-overlay/libaec/portfile.cmake | 44 ------------------------ cpp/.vcpkg-overlay/libaec/usage | 7 ---- cpp/.vcpkg-overlay/libaec/vcpkg.json | 17 --------- 4 files changed, 5 insertions(+), 84 deletions(-) delete mode 100644 cpp/.vcpkg-overlay/libaec/portfile.cmake delete mode 100644 cpp/.vcpkg-overlay/libaec/usage delete mode 100644 cpp/.vcpkg-overlay/libaec/vcpkg.json diff --git a/cpp/.vcpkg-overlay/README.md b/cpp/.vcpkg-overlay/README.md index 17672e51f8b..e8e98953a3b 100644 --- a/cpp/.vcpkg-overlay/README.md +++ b/cpp/.vcpkg-overlay/README.md @@ -1,25 +1,14 @@ -# vcpkg overlay ports - -## Intel MPI +# vcpkg overlay port for Intel MPI This vcpkg overlay port contains scripts for installing Intel MPI on Windows (only). MSMPI, which is used by default with vcpkg, does not support the MPI3 standard. Using this port requires that Intel OneAPI binaries are already installed. On Unix systems the built-in OpenMPI or MPICH ports can be used. -This overlay port was adapted from the original at: - -https://github.com/arcaneframework/framework-ci - -## libaec - -The upstream vcpkg port for `libaec` downloads from `https://gitlab.dkrz.de`, -which can intermittently fail with HTTP 429 (rate limiting) errors in CI. -This overlay port redirects the download to the GitHub mirror at -https://github.com/MathisRosenhauer/libaec. - -## Usage - From the root of this repository it can be activated by e.g.: cmake -DCMAKE_TOOLCHAIN_FILE=%VCPKG_ROOT%/scripts/buildsystems/vcpkg.cmake -DVCPKG_OVERLAY_PORTS="cpp/.vcpkg-overlay" -B build-dir -S cpp/ + +This overlay port was adapted from the original at: + +https://github.com/arcaneframework/framework-ci diff --git a/cpp/.vcpkg-overlay/libaec/portfile.cmake b/cpp/.vcpkg-overlay/libaec/portfile.cmake deleted file mode 100644 index 65c34cc787a..00000000000 --- a/cpp/.vcpkg-overlay/libaec/portfile.cmake +++ /dev/null @@ -1,44 +0,0 @@ -vcpkg_from_github( - OUT_SOURCE_PATH SOURCE_PATH - REPO MathisRosenhauer/libaec - REF "v${VERSION}" - SHA512 55bd605590015e0f903a231268265051336c172c18935fdfecead5630454a99811f3f196117556ad1b62f3052d54894fd8e257a9ea9f45cc03f59c93cde00cda - HEAD_REF master -) - -string(COMPARE EQUAL "${VCPKG_LIBRARY_LINKAGE}" "static" BUILD_STATIC) - -vcpkg_cmake_configure( - SOURCE_PATH "${SOURCE_PATH}" - OPTIONS - -DBUILD_STATIC_LIBS=${BUILD_STATIC} - -Dlibaec_INSTALL_CMAKEDIR=share/${PORT} -) -vcpkg_cmake_install() -vcpkg_copy_pdbs() -vcpkg_cmake_config_fixup() -vcpkg_replace_string("${CURRENT_PACKAGES_DIR}/share/libaec/libaec-config.cmake" - "if(libaec_USE_STATIC_LIBS)" - "if(\"${BUILD_STATIC}\") # forced by vcpkg" -) - -# Compatibility with user's CMake < 3.18 (vcpkg claims support for >= 3.16): -# Make imported targets global so that libaec-config.cmake can create ALIAS targets. -set(_target_file "libaec_shared-targets") -if(BUILD_STATIC) - set(_target_file "libaec_static-targets") -endif() -file(READ "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" libaec_targets) -string(REGEX REPLACE " (SHARED|STATIC) IMPORTED" " \\1 IMPORTED \${libaec_maybe_global}" libaec_targets "${libaec_targets}") -file(WRITE "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" "set(libaec_maybe_global \"\") -if(CMAKE_VERSION VERSION_LESS 3.18) - set(libaec_maybe_global \"GLOBAL\") -endif() -${libaec_targets} -" -) - -file(REMOVE_RECURSE "${CURRENT_PACKAGES_DIR}/debug/include") - -file(INSTALL "${CURRENT_PORT_DIR}/usage" DESTINATION "${CURRENT_PACKAGES_DIR}/share/${PORT}") -vcpkg_install_copyright(FILE_LIST "${SOURCE_PATH}/LICENSE.txt") diff --git a/cpp/.vcpkg-overlay/libaec/usage b/cpp/.vcpkg-overlay/libaec/usage deleted file mode 100644 index a6aeb5f6a63..00000000000 --- a/cpp/.vcpkg-overlay/libaec/usage +++ /dev/null @@ -1,7 +0,0 @@ -libaec provides CMake targets: - - find_package(libaec CONFIG REQUIRED) - # libaec API - target_link_libraries(main PRIVATE libaec::aec) - # szip compatible API - target_link_libraries(main PRIVATE libaec::sz) diff --git a/cpp/.vcpkg-overlay/libaec/vcpkg.json b/cpp/.vcpkg-overlay/libaec/vcpkg.json deleted file mode 100644 index 69c73b3559d..00000000000 --- a/cpp/.vcpkg-overlay/libaec/vcpkg.json +++ /dev/null @@ -1,17 +0,0 @@ -{ - "name": "libaec", - "version": "1.1.7", - "description": "Adaptive Entropy Coding library", - "homepage": "https://github.com/MathisRosenhauer/libaec", - "license": "BSD-2-Clause", - "dependencies": [ - { - "name": "vcpkg-cmake", - "host": true - }, - { - "name": "vcpkg-cmake-config", - "host": true - } - ] -} From 77a7ce491db9095ff3e849935651ff34121b3a51 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 21 Jul 2026 09:33:26 +0200 Subject: [PATCH 20/20] Simplify further --- python/dolfinx/fem/bcs.py | 23 ++------- python/dolfinx/fem/element.py | 48 ++++-------------- python/dolfinx/fem/forms.py | 29 +++-------- python/dolfinx/fem/function.py | 92 ++++++---------------------------- 4 files changed, 36 insertions(+), 156 deletions(-) diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index d3534e3f24d..81884699bf9 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -98,8 +98,9 @@ class DirichletBC(Generic[Scalar]): The conditions are imposed on a linear system. """ - # Matched-precision built-ins (geometry == real scalar part). - _cpp_registry: ClassVar[dict] = { + # Matched-precision built-ins (geometry == real scalar part). Public: + # extend with additional (scalar, geometry) dtype pairs as needed. + cpp_types: ClassVar[dict] = { (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.DirichletBC_float32, (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.DirichletBC_float64, (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.DirichletBC_complex64, @@ -112,22 +113,6 @@ class DirichletBC(Generic[Scalar]): | _cpp.fem.DirichletBC_float64 ) - @classmethod - def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): - """Register a compiled C++ type for a (scalar, geometry) pair.""" - cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, scalar_dtype, geometry_dtype): - """Compiled C++ type registered for a (scalar, geometry) pair.""" - try: - return cls._cpp_registry[scalar_dtype, geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " - f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__(self, bc): """Initialise a Dirichlet boundary condition. @@ -234,7 +219,7 @@ def dirichletbc( geometry_dtype = np.dtype(value.function_space.mesh.geometry.x.dtype) else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - bctype = DirichletBC.get_cpp_type(dtype, geometry_dtype) + bctype = DirichletBC.cpp_types[dtype, geometry_dtype] # Unwrap value object, if required if isinstance(value, np.ndarray): diff --git a/python/dolfinx/fem/element.py b/python/dolfinx/fem/element.py index 181c78aee42..c7313b00b93 100644 --- a/python/dolfinx/fem/element.py +++ b/python/dolfinx/fem/element.py @@ -20,29 +20,14 @@ class CoordinateElement(Generic[Real]): """Coordinate element describing the geometry map for mesh cells.""" - # Built-in geometry types. - _cpp_registry: ClassVar[dict] = { + # Built-in geometry types. Public: extend with additional geometry + # dtypes as needed. + cpp_types: ClassVar[dict] = { np.dtype(np.float32): _cpp.fem.CoordinateElement_float32, np.dtype(np.float64): _cpp.fem.CoordinateElement_float64, } _cpp_object: _cpp.fem.CoordinateElement_float32 | _cpp.fem.CoordinateElement_float64 - @classmethod - def register_cpp_type(cls, cpp_type, geometry_dtype): - """Register a compiled C++ type for a geometry dtype.""" - cls._cpp_registry[geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, geometry_dtype): - """Compiled C++ type registered for a geometry dtype.""" - try: - return cls._cpp_registry[geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " - f"Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, cmap: _cpp.fem.CoordinateElement_float32 | _cpp.fem.CoordinateElement_float64 ): @@ -165,7 +150,7 @@ def coordinate_element( Returns: A coordinate element. """ - cpp_type = CoordinateElement.get_cpp_type(np.dtype(dtype)) + cpp_type = CoordinateElement.cpp_types[np.dtype(dtype)] return CoordinateElement(cpp_type(celltype, degree, variant)) @@ -181,36 +166,21 @@ def _(e: basix.finite_element.FiniteElement) -> CoordinateElement: Returns: A coordinate element. """ - cpp_type = CoordinateElement.get_cpp_type(e.dtype) + cpp_type = CoordinateElement.cpp_types[e.dtype] return CoordinateElement(cpp_type(e._e)) class FiniteElement(Generic[Real]): """A finite element.""" - # Built-in geometry types. - _cpp_registry: ClassVar[dict] = { + # Built-in geometry types. Public: extend with additional geometry + # dtypes as needed. + cpp_types: ClassVar[dict] = { np.dtype(np.float32): _cpp.fem.FiniteElement_float32, np.dtype(np.float64): _cpp.fem.FiniteElement_float64, } _cpp_object: _cpp.fem.FiniteElement_float32 | _cpp.fem.FiniteElement_float64 - @classmethod - def register_cpp_type(cls, cpp_type, geometry_dtype): - """Register a compiled C++ type for a geometry dtype.""" - cls._cpp_registry[geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, geometry_dtype): - """Compiled C++ type registered for a geometry dtype.""" - try: - return cls._cpp_registry[geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " - f"Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, cpp_object: _cpp.fem.FiniteElement_float32 | _cpp.fem.FiniteElement_float64, @@ -382,7 +352,7 @@ def finiteelement( the selected element. FiniteElement_dtype: Geometry type of the element. """ - CppElement = FiniteElement.get_cpp_type(np.dtype(FiniteElement_dtype)) + CppElement = FiniteElement.cpp_types[np.dtype(FiniteElement_dtype)] if ufl_e.is_mixed: elements = [ diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index a601e044a60..a54a6987378 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -37,8 +37,9 @@ class Form(typing.Generic[Scalar]): """A finite element form.""" - # Matched-precision built-ins (geometry == real scalar part). - _cpp_registry: typing.ClassVar[dict] = { + # Matched-precision built-ins (geometry == real scalar part). Public: + # extend with additional (scalar, geometry) dtype pairs as needed. + cpp_types: typing.ClassVar[dict] = { (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.Form_float32, (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.Form_float64, (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.Form_complex64, @@ -52,22 +53,6 @@ class Form(typing.Generic[Scalar]): ) _code: str | list[str] | None - @classmethod - def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): - """Register a compiled C++ type for a (scalar, geometry) pair.""" - cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, scalar_dtype, geometry_dtype): - """Compiled C++ type registered for a (scalar, geometry) pair.""" - try: - return cls._cpp_registry[scalar_dtype, geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " - f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, form: _cpp.fem.Form_complex64 @@ -247,7 +232,7 @@ def form_cpp_class( """ # Geometry is the real part of the scalar type (matched precision). scalar_dtype = np.dtype(dtype) - return Form.get_cpp_type(scalar_dtype, scalar_dtype.type(0).real.dtype) + return Form.cpp_types[scalar_dtype, scalar_dtype.type(0).real.dtype] _ufl_to_dolfinx_domain = { @@ -311,7 +296,7 @@ def mixed_topology_form( comm = mesh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = Form.get_cpp_type(np.dtype(dtype), mesh.geometry.x.dtype) + ftype = Form.cpp_types[np.dtype(dtype), mesh.geometry.x.dtype] ufcx_forms = [] modules = [] @@ -399,7 +384,7 @@ def _form(form): comm = msh.comm if jit_comm is None else jit_comm # Geometry type is fixed by the mesh. - ftype = Form.get_cpp_type(np.dtype(dtype), msh.geometry.x.dtype) + ftype = Form.cpp_types[np.dtype(dtype), msh.geometry.x.dtype] ufcx_form, module, code = jit.ffcx_jit( comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options @@ -464,7 +449,7 @@ def _zero_form(form): msh = V[0].mesh # Geometry type is fixed by the mesh. - ftype = Form.get_cpp_type(np.dtype(dtype), msh.geometry.x.dtype) + ftype = Form.cpp_types[np.dtype(dtype), msh.geometry.x.dtype] f = ftype( spaces=V, diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 759528df915..83ae7fc6b93 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -35,8 +35,9 @@ class Constant(ufl.Constant, Generic[Scalar]): """A constant with respect to a domain.""" - # Built-in scalar types. - _cpp_registry: typing.ClassVar[dict] = { + # Built-in scalar types. Public: extend with additional scalar dtypes + # as needed. + cpp_types: typing.ClassVar[dict] = { np.dtype(np.float32): _cpp.fem.Constant_float32, np.dtype(np.float64): _cpp.fem.Constant_float64, np.dtype(np.complex64): _cpp.fem.Constant_complex64, @@ -49,22 +50,6 @@ class Constant(ufl.Constant, Generic[Scalar]): | _cpp.fem.Constant_float64 ) - @classmethod - def register_cpp_type(cls, cpp_type, scalar_dtype): - """Register a compiled C++ type for a scalar dtype.""" - cls._cpp_registry[scalar_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, scalar_dtype): - """Compiled C++ type registered for a scalar dtype.""" - try: - return cls._cpp_registry[scalar_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{scalar_dtype}'. " - f"Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, domain, @@ -79,7 +64,7 @@ def __init__( c = np.asarray(c) super().__init__(domain, c.shape) try: - cpp_type = Constant.get_cpp_type(c.dtype) + cpp_type = Constant.cpp_types[c.dtype] except AttributeError as err: raise AttributeError("Constant value must have a dtype attribute.") from err self._cpp_object = cpp_type(c) @@ -129,7 +114,8 @@ class Expression(Generic[Scalar]): # Matched-precision built-ins (geometry == real scalar part). An # ``Expression`` is built through a factory, not a class constructor. - _cpp_registry: typing.ClassVar[dict] = { + # Public: extend with additional (scalar, geometry) dtype pairs. + cpp_types: typing.ClassVar[dict] = { (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.create_expression_float32, (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.create_expression_float64, (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.create_expression_complex64, @@ -145,22 +131,6 @@ class Expression(Generic[Scalar]): ) _code: str - @classmethod - def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): - """Register a C++ factory for a (scalar, geometry) pair.""" - cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, scalar_dtype, geometry_dtype): - """C++ factory registered for a (scalar, geometry) pair.""" - try: - return cls._cpp_registry[scalar_dtype, geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " - f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, e: ufl.core.expr.Expr, @@ -252,7 +222,7 @@ def __init__( geometry_dtype = expr_domain.ufl_cargo().geometry.x.dtype else: geometry_dtype = np.dtype(dtype).type(0).real.dtype - create_expression = Expression.get_cpp_type(np.dtype(dtype), geometry_dtype) + create_expression = Expression.cpp_types[np.dtype(dtype), geometry_dtype] _entity_maps = ( [entity_map._cpp_object for entity_map in entity_maps] @@ -378,8 +348,9 @@ class Function(ufl.Coefficient, Generic[Scalar]): """ - # Matched-precision built-ins (geometry == real scalar part). - _cpp_registry: typing.ClassVar[dict] = { + # Matched-precision built-ins (geometry == real scalar part). Public: + # extend with additional (scalar, geometry) dtype pairs. + cpp_types: typing.ClassVar[dict] = { (np.dtype(np.float32), np.dtype(np.float32)): _cpp.fem.Function_float32, (np.dtype(np.float64), np.dtype(np.float64)): _cpp.fem.Function_float64, (np.dtype(np.complex64), np.dtype(np.float32)): _cpp.fem.Function_complex64, @@ -394,22 +365,6 @@ class Function(ufl.Coefficient, Generic[Scalar]): _x: la.Vector[Scalar] - @classmethod - def register_cpp_type(cls, cpp_type, scalar_dtype, geometry_dtype): - """Register a compiled C++ type for a (scalar, geometry) pair.""" - cls._cpp_registry[scalar_dtype, geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, scalar_dtype, geometry_dtype): - """Compiled C++ type registered for a (scalar, geometry) pair.""" - try: - return cls._cpp_registry[scalar_dtype, geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for scalar '{scalar_dtype}' on geometry " - f"'{geometry_dtype}'. Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, V: FunctionSpace, @@ -439,7 +394,7 @@ def __init__( # Scalar type (dtype) is independent of the geometry type, which # is fixed by the mesh. geometry_dtype = V.mesh.geometry.x.dtype - cpp_type = Function.get_cpp_type(np.dtype(dtype), geometry_dtype) + cpp_type = Function.cpp_types[np.dtype(dtype), geometry_dtype] if x is not None: self._cpp_object = cpp_type(V._cpp_object, x._cpp_object) # type: ignore else: @@ -747,7 +702,7 @@ def functionspace( ) # Initialize the cpp.FunctionSpace. Geometry type is fixed by the mesh. - cpp_type = FunctionSpace.get_cpp_type(dtype) + cpp_type = FunctionSpace.cpp_types[dtype] cppV = cpp_type(mesh._cpp_object, element._cpp_object, cpp_dofmap) # type: ignore return FunctionSpace(mesh, ufl_e, cppV) @@ -756,30 +711,15 @@ def functionspace( class FunctionSpace(ufl.FunctionSpace, Generic[Real]): """A space on which Functions (fields) can be defined.""" - # Built-in geometry types. - _cpp_registry: typing.ClassVar[dict] = { + # Built-in geometry types. Public: extend with additional geometry + # dtypes as needed. + cpp_types: typing.ClassVar[dict] = { np.dtype(np.float32): _cpp.fem.FunctionSpace_float32, np.dtype(np.float64): _cpp.fem.FunctionSpace_float64, } _cpp_object: _cpp.fem.FunctionSpace_float32 | _cpp.fem.FunctionSpace_float64 _mesh: Mesh[Real] - @classmethod - def register_cpp_type(cls, cpp_type, geometry_dtype): - """Register a compiled C++ type for a geometry dtype.""" - cls._cpp_registry[geometry_dtype] = cpp_type - - @classmethod - def get_cpp_type(cls, geometry_dtype): - """Compiled C++ type registered for a geometry dtype.""" - try: - return cls._cpp_registry[geometry_dtype] - except KeyError: - raise NotImplementedError( - f"No compiled {cls.__name__} for geometry '{geometry_dtype}'. " - f"Register one with {cls.__name__}.register_cpp_type()." - ) from None - def __init__( self, mesh: Mesh[Real], @@ -826,7 +766,7 @@ def clone(self) -> FunctionSpace[Real]: A new function space that shares data """ # noqa: D301 # Geometry type is fixed by the mesh. - cpp_type = FunctionSpace.get_cpp_type(self._mesh.geometry.x.dtype) + cpp_type = FunctionSpace.cpp_types[self._mesh.geometry.x.dtype] Vcpp = cpp_type(self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap) # type: ignore return FunctionSpace(self._mesh, self.ufl_element(), Vcpp)