From 7ade28f6d42b2a4dd1f042157f7b5784b1381102 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 09:20:46 +0000 Subject: [PATCH 01/19] Make submesh expressions work (codim 0 for now, missing test for codim 1). --- cpp/dolfinx/fem/Expression.h | 60 ++++---- cpp/dolfinx/fem/assembler.h | 6 +- cpp/dolfinx/fem/pack.h | 121 +++++++++++++-- cpp/dolfinx/fem/utils.h | 9 +- python/dolfinx/fem/function.py | 12 +- .../wrappers/dolfinx_wrappers/assemble.h | 11 +- .../dolfinx/wrappers/dolfinx_wrappers/fem.h | 14 +- python/test/unit/mesh/cform.py | 144 ++++++++++++++++++ 8 files changed, 320 insertions(+), 57 deletions(-) create mode 100644 python/test/unit/mesh/cform.py diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h index eb524f0a234..4097d4547ae 100644 --- a/cpp/dolfinx/fem/Expression.h +++ b/cpp/dolfinx/fem/Expression.h @@ -1,4 +1,5 @@ -// Copyright (C) 2020-2025 Jack S. Hale, Michal Habera and Garth N. Wells +// Copyright (C) 2020-2026 Jack S. Hale, Michal Habera, Garth N. Wells and +// Jørgen S. Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -12,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -62,40 +64,32 @@ class Expression /// @param[in] fn Function for tabulating the Expression. /// @param[in] value_shape Shape of Expression evaluated at single /// point. + /// @param[in] entity_maps Maps between entities of different meshes. /// @param[in] argument_space Function space for Argument. Used to /// computed a 1-form expression, e.g. can be used to create a matrix /// that when applied to a degree-of-freedom vector gives the /// expression values at the evaluation points. - Expression(const std::vector>>& coefficients, - const std::vector>>& - constants, - std::span X, - std::array Xshape, - std::function - fn, - const std::vector& value_shape, - std::shared_ptr> argument_space - = nullptr) + Expression( + const std::vector>>& coefficients, + const std::vector>>& + constants, + std::span X, std::array Xshape, + std::function + fn, + const std::vector& value_shape, + const std::vector>& + entity_maps, + std::shared_ptr> argument_space + = nullptr) : _argument_space(argument_space), _coefficients(coefficients), _constants(constants), _fn(fn), _value_shape(value_shape), - _x_ref(std::vector(X.begin(), X.end()), Xshape) + _x_ref(std::vector(X.begin(), X.end()), Xshape), + _entity_maps(entity_maps) {}; - { - for (auto& c : _coefficients) - { - assert(c); - if (c->function_space()->mesh() - != _coefficients.front()->function_space()->mesh()) - { - throw std::runtime_error("Coefficients not all defined on same mesh."); - } - } - } - - /// Move constructor + /// Move constructorcreate_expression Expression(Expression&& e) = default; /// Destructor @@ -168,6 +162,13 @@ class Expression return _x_ref; } + /// @brief Maps between entities of different meshes. + const std::vector>& + entity_maps() const + { + return _entity_maps; + } + private: // Function space for Argument std::shared_ptr> _argument_space; @@ -189,5 +190,8 @@ class Expression // Evaluation points on reference cell std::pair, std::array> _x_ref; + + std::vector> + _entity_maps; }; } // namespace dolfinx::fem diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 352196d502f..9e2460b25c6 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -1,4 +1,4 @@ -// Copyright (C) 2018-2025 Garth N. Wells +// Copyright (C) 2018-2026 Garth N. Wells and Jørgen S. Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -115,7 +116,8 @@ void tabulate_expression(std::span values, const fem::Expression& e, std::vector>> c; std::ranges::transform(coefficients, std::back_inserter(c), [](auto c) -> const Function& { return *c; }); - fem::pack_coefficients(c, coffsets, entities, std::span(coeffs)); + fem::pack_coefficients(c, mesh, entities, e.entity_maps(), coffsets, + std::span(coeffs)); } std::vector constants = fem::pack_constants(e); diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 65777820ee8..287a150ebea 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -1,4 +1,4 @@ -// Copyright (C) 2013-2025 Garth N. Wells +// Copyright (C) 2013-2026 Garth N. Wells and Jørgen S. Dokken // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -332,17 +333,23 @@ void pack_coefficients(const Form& form, /// @brief Pack coefficient data over a list of cells or facets. /// /// Typically used to prepare coefficient data for an ::Expression. -/// @tparam T -/// @tparam U +/// @tparam T data type of coefficients +/// @tparam U floating point type of mesh geometry /// @param coeffs Coefficients to pack -/// @param offsets Offsets +/// @param mesh Mesh which the entities belong to /// @param entities Entities to pack over -/// @param[out] c Packed coefficients. +/// @param entity_maps Entity maps for the coefficient meshes +/// @param offsets Offsets +/// @param[in,out] c Packed coefficients. template void pack_coefficients( std::vector>> coeffs, - std::span offsets, fem::MDSpan2 auto entities, std::span c) + const mesh::Mesh& mesh, fem::MDSpan2 auto entities, + const std::vector>& + entity_maps, + std::span offsets, std::span c) { + assert(!offsets.empty()); const int cstride = offsets.back(); @@ -354,17 +361,105 @@ void pack_coefficients( { std::span cell_info = impl::get_cell_orientation_info(coeffs[coeff].get()); - - if constexpr (entities.rank() == 1) + // Get mesh of coefficient and check if entity map is required + auto mesh_c = coeffs[coeff].get().function_space()->mesh(); + if (mesh_c->topology() == mesh.topology()) { - impl::pack_coefficient_entity(std::span(c), cstride, coeffs[coeff].get(), - cell_info, entities, offsets[coeff]); + // If same mesh no mapping is needed + if constexpr (entities.rank() == 1) + { + impl::pack_coefficient_entity(std::span(c), cstride, + coeffs[coeff].get(), cell_info, entities, + offsets[coeff]); + } + else + { + auto cells = md::submdspan(entities, md::full_extent, 0); + impl::pack_coefficient_entity(std::span(c), cstride, + coeffs[coeff].get(), cell_info, cells, + offsets[coeff]); + } } else { - auto cells = md::submdspan(entities, md::full_extent, 0); - impl::pack_coefficient_entity(std::span(c), cstride, coeffs[coeff].get(), - cell_info, cells, offsets[coeff]); + // Helper function to get correct entity map + auto get_entity_map + = [mesh, &entity_maps](auto& mesh0) -> const mesh::EntityMap& + { + auto it = std::ranges::find_if( + entity_maps, + [mesh, mesh0](const mesh::EntityMap& em) + { + return ((em.topology() == mesh0->topology() + and em.sub_topology() == mesh.topology())) + or ((em.sub_topology() == mesh0->topology() + and em.topology() == mesh.topology())); + }); + + if (it == entity_maps.end()) + { + throw std::runtime_error( + "Incompatible mesh. argument entity_maps must be provided."); + } + return *it; + }; + // Find correct entity map and determine direction of the map + const mesh::EntityMap& emap = get_entity_map(mesh_c); + bool inverse = emap.sub_topology() == mesh_c->topology(); + + std::vector e_b; + if constexpr (entities.rank() == 1) + { + e_b = emap.sub_topology_to_topology( + std::span(entities.data_handle(), entities.size()), inverse); + md::mdspan e(e_b.data(), e_b.size()); + impl::pack_coefficient_entity(std::span(c), cstride, + coeffs[coeff].get(), cell_info, e, + offsets[coeff]); + } + else if (entities.rank() == 2) + { + const mesh::Topology topology = *mesh.topology(); + int tdim = topology.dim(); + // If codim is zero we extract the cells + int codim = tdim - mesh.topology()->dim(); + if (codim == 0) + { + auto cells = md::submdspan(entities, md::full_extent, 0); + assert(cells.stride(0) == 1 + && "Slice is not contiguous; cannot convert to std::span " + "safely!"); + e_b = emap.sub_topology_to_topology( + std::span(cells.data_handle(), cells.size()), inverse); + md::mdspan e(e_b.data(), e_b.size()); + impl::pack_coefficient_entity(std::span(c), cstride, + coeffs[coeff].get(), cell_info, e, + offsets[coeff]); + } + else if (codim == 1) + { + // Codim 1 mesh, need to map (cell, local index) to facets and then + // to cells of the submesh + auto c_to_f = topology.connectivity(tdim, tdim - 1); + assert(c_to_f); + std::vector parent_entities; + parent_entities.reserve(entities.extent(0)); + for (std::size_t e = 0; e < entities.extent(0); ++e) + { + auto pair = md::submdspan(entities, e, md::full_extent); + parent_entities.push_back(c_to_f->links(pair[0])[pair[1]]); + } + e_b = emap.sub_topology_to_topology(parent_entities, inverse); + md::mdspan e(e_b.data(), e_b.size()); + impl::pack_coefficient_entity(std::span(c), cstride, + coeffs[coeff].get(), cell_info, e, + offsets[coeff]); + } + } + else + { + throw std::runtime_error("Entities span has unsupported rank."); + } } } } diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 81fbc2f953d..5d4e87ecf8d 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -953,6 +953,8 @@ Expression create_expression( const ufcx_expression& e, const std::vector>>& coefficients, const std::vector>>& constants, + const std::vector>& + entity_maps, std::shared_ptr> argument_space = nullptr) { if (!coefficients.empty()) @@ -1008,7 +1010,7 @@ Expression create_expression( assert(tabulate_tensor); return Expression(coefficients, constants, std::span(X), Xshape, - tabulate_tensor, value_shape, argument_space); + tabulate_tensor, value_shape, entity_maps, argument_space); } /// @brief Create Expression from UFC input (with named coefficients and @@ -1019,6 +1021,8 @@ Expression create_expression( const std::map>>& coefficients, const std::map>>& constants, + const std::vector>& + entity_maps, std::shared_ptr> argument_space = nullptr) { // Place coefficients in appropriate order @@ -1057,7 +1061,8 @@ Expression create_expression( } } - return create_expression(e, coeff_map, const_map, argument_space); + return create_expression(e, coeff_map, const_map, entity_maps, + argument_space); } } // namespace dolfinx::fem diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 69c9cce7f2e..a8d062b29bb 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -117,6 +117,7 @@ def __init__( form_compiler_options: dict | None = None, jit_options: dict | None = None, dtype: npt.DTypeLike | None = None, + entity_maps: list[dolfinx.mesh.EntityMap] | None = None, ): """Create an Expression. @@ -131,6 +132,7 @@ def __init__( jit_options: Options controlling JIT compilation of C code. dtype: Type of the Expression values. If ``None``, the dtype is deduced from the UFL expression. + entity_maps: Maps between different meshes. Note: This wrapper is responsible for the FFCx compilation of the @@ -144,7 +146,10 @@ def __init__( # Get MPI communicator if comm is None: try: - domain = ufl.domain.extract_unique_domain(e) + domains = ufl.domain.extract_domains(e) + if len(domains) == 0: + raise RuntimeError("Could not extract MPI communicator for Expression.") + domain = domains[0] assert isinstance(domain, ufl.Mesh) mesh = domain.ufl_cargo() comm = mesh.comm @@ -199,11 +204,13 @@ def _create_expression(dtype): else: raise NotImplementedError(f"Type {dtype} not supported.") + _entity_maps = [entity_map._cpp_object for entity_map in entity_maps] if entity_maps is not None else [] ffi = module.ffi self._cpp_object = _create_expression(dtype)( ffi.cast("uintptr_t", ffi.addressof(self._ufcx_expression)), coeffs, constants, + _entity_maps, self.argument_space, ) @@ -262,7 +269,8 @@ def eval( raise TypeError("Passed values array does not have correct dtype.") constants = _cpp.fem.pack_constants(self._cpp_object) - coeffs = _cpp.fem.pack_coefficients(self._cpp_object, _entities) + breakpoint() + coeffs = _cpp.fem.pack_coefficients(self._cpp_object, mesh,_entities) _cpp.fem.tabulate_expression( values, self._cpp_object, constants, coeffs, mesh._cpp_object, _entities ) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h b/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h index dce4671a76e..082eef4c068 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/assemble.h @@ -212,6 +212,7 @@ void declare_assembly_functions(nanobind::module_& m) m.def( "pack_coefficients", [](const dolfinx::fem::Expression& e, + const dolfinx::mesh::Mesh& mesh, nb::ndarray entities) { std::vector coffsets = e.coefficient_offsets(); @@ -228,22 +229,22 @@ void declare_assembly_functions(nanobind::module_& m) if (entities.ndim() == 1) { dolfinx::fem::pack_coefficients( - c, coffsets, md::mdspan(entities.data(), entities.shape(0)), - std::span(coeffs)); + c, mesh, md::mdspan(entities.data(), entities.shape(0)), + e.entity_maps(), coffsets, std::span(coeffs)); } else { dolfinx::fem::pack_coefficients( - c, coffsets, + c, mesh, md::mdspan>( entities.data(), entities.shape(0), entities.shape(1)), - std::span(coeffs)); + e.entity_maps(), coffsets, std::span(coeffs)); } return dolfinx_wrappers::as_nbarray(std::move(coeffs), {entities.shape(0), cstride}); }, - nb::arg("expression"), nb::arg("entities"), + nb::arg("expression"), nb::arg("mesh"), nb::arg("entities"), "Pack coefficients for a Expression."); m.def( "pack_constants", diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index d208ff30b1b..018cffd35f0 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -683,6 +683,7 @@ void declare_objects(nb::module_& m, std::string type) nb::ndarray, nb::c_contig> X, std::uintptr_t fn_addr, const std::vector& value_shape, + const std::vector& entity_maps, std::shared_ptr> argument_space) { @@ -693,10 +694,11 @@ void declare_objects(nb::module_& m, std::string type) new (ex) dolfinx::fem::Expression( coefficients, constants, std::span(X.data(), X.size()), {X.shape(0), X.shape(1)}, tabulate_expression_ptr, value_shape, - argument_space); + ptr_to_ref_wrapper_vec(entity_maps), argument_space); }, nb::arg("coefficients"), nb::arg("constants"), nb::arg("X"), - nb::arg("fn"), nb::arg("value_shape"), nb::arg("argument_space")) + nb::arg("fn"), nb::arg("value_shape"), nb::arg("entity_maps"), + nb::arg("argument_space")) .def("X", [](const dolfinx::fem::Expression& self) { @@ -725,15 +727,17 @@ void declare_objects(nb::module_& m, std::string type) coefficients, const std::vector>>& constants, + 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, argument_space); + return dolfinx::fem::create_expression( + *p, coefficients, constants, ptr_to_ref_wrapper_vec(entity_maps), + argument_space); }, nb::arg("expression"), nb::arg("coefficients"), nb::arg("constants"), - nb::arg("argument_space").none(), + nb::arg("entity_maps"), nb::arg("argument_space").none(), "Create Expression from a pointer to ufc_form."); } diff --git a/python/test/unit/mesh/cform.py b/python/test/unit/mesh/cform.py new file mode 100644 index 00000000000..033feb25f0c --- /dev/null +++ b/python/test/unit/mesh/cform.py @@ -0,0 +1,144 @@ +import typing +from dataclasses import dataclass + +from mpi4py import MPI + +import numpy as np +import numpy.typing as npt + +import basix.ufl +import dolfinx +import ufl + + +@dataclass +class GeneralForm: + ufl_form: ufl.Form + ufcx_form: typing.Any + module: typing.Any + code: str + + +def compile_form( + comm: MPI.Intracomm, + form: ufl.Form, +) -> GeneralForm: + """Compile UFL form withtout associated DOLFINx data.""" + compiled_form = dolfinx.jit.ffcx_jit(comm, form) + return GeneralForm(form, *compiled_form) + + +def get_integration_domains(integral_type, subdomain): + """Get integration domains from subdomain data.""" + if subdomain is None: + return [] + else: + try: + if integral_type in ( + dolfinx.fem.IntegralType.exterior_facet, + dolfinx.fem.IntegralType.interior_facet, + ): + tdim = subdomain.topology.dim + subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) + subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) + domains = dolfinx.cpp.fem.compute_integration_domains( + integral_type, subdomain._cpp_object + ) + return [(s[0], np.array(s[1])) for s in domains] + except AttributeError: + return [(s[0], np.array(s[1])) for s in subdomain] + + +def form_cpp_creator( + dtype: npt.DTypeLike, +) -> ( + dolfinx.cpp.fem.Form_float32 + | dolfinx.cpp.fem.Form_float64 + | dolfinx.cpp.fem.Form_complex64 + | dolfinx.cpp.fem.Form_complex128 +): + """Return the wrapped C++ class of a variational form of a specific scalar type. + + Args: + dtype: Scalar type of the required form class. + + Returns: + Wrapped C++ form class of the requested type. + + Note: + This function is for advanced usage, typically when writing + custom kernels using Numba or C. + """ + if np.issubdtype(dtype, np.float32): + return dolfinx.cpp.fem.create_form_float32 + elif np.issubdtype(dtype, np.float64): + return dolfinx.cpp.fem.create_form_float64 + elif np.issubdtype(dtype, np.complex64): + return dolfinx.cpp.fem.create_form_complex64 + elif np.issubdtype(dtype, np.complex128): + return dolfinx.cpp.fem.create_form_complex128 + else: + raise NotImplementedError(f"Type {dtype} not supported.") + + +def create_form( + form: GeneralForm, + mesh: dolfinx.mesh.Mesh, + coefficient_map: dict[str, dolfinx.fem.Function], + constant_map: dict[str, dolfinx.fem.Constant], +): + sd = form.ufl_form.subdomain_data() + (domain,) = list(sd.keys()) # Assuming single domain + + # Subdomain markers (possibly empty list for some integral types) + subdomains = { + dolfinx.fem.forms._ufl_to_dolfinx_domain[key]: get_integration_domains( + dolfinx.fem.forms._ufl_to_dolfinx_domain[key], subdomain_data[0] + ) + for (key, subdomain_data) in sd.get(domain).items() + } + coefficients = {f"w{u.count()}": uh._cpp_object for (u, uh) in coefficient_map.items()} + constants = {f"c{c.count()}": ch._cpp_object for (c, ch) in constant_map.items()} + ftype = form_cpp_creator(dolfinx.default_scalar_type) + f = ftype( + form.module.ffi.cast("uintptr_t", form.module.ffi.addressof(form.ufcx_form)), + [], + coefficients, + constants, + subdomains, + mesh._cpp_object, + ) + return dolfinx.fem.Form(f, form.ufcx_form, form.code) + + +def test_compiled_form(): + # Create ufl form exclusively with basix.ufl and UFL + c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)) + domain = ufl.Mesh(c_el) + el = basix.ufl.element("Lagrange", "triangle", 2) + V = ufl.FunctionSpace(domain, el) + u = ufl.Coefficient(V) + w = ufl.Coefficient(V) + c = ufl.Constant(domain) + e = ufl.Constant(domain) + J = c * e * u * w * ufl.dx(domain=domain) + + # Compile form using dolfinx.jit.ffcx_jit + compiled_form = compile_form(MPI.COMM_WORLD, J) + + def create_and_integrate(N, compiled_form): + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) + el_2 = basix.ufl.element("Lagrange", "triangle", 2) + Vh = dolfinx.fem.functionspace(mesh, el_2) + uh = dolfinx.fem.Function(Vh) + uh.interpolate(lambda x: x[0]) + wh = dolfinx.fem.Function(Vh) + wh.interpolate(lambda x: x[1]) + eh = dolfinx.fem.Constant(mesh, 3.0) + ch = dolfinx.fem.Constant(mesh, 2.0) + form = create_form(compiled_form, mesh, {u: uh, w: wh}, {c: ch, e: eh}) + assert np.isclose(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(form), op=MPI.SUM), 1.5) + + # Create various meshes, that all uses this compiled form with a map from ufl to dolfinx functions and constants + for i in range(1, 4): + create_and_integrate(i, compiled_form) From 57e9fddeb2fff93c8bfa7e877f16d690d4c9b4ae Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 09:21:01 +0000 Subject: [PATCH 02/19] Minor cleanup --- python/dolfinx/fem/function.py | 3 +- python/test/unit/fem/test_expression.py | 49 ++++++++++++++++++++++++- 2 files changed, 49 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index a8d062b29bb..d7699937ca6 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -269,8 +269,7 @@ def eval( raise TypeError("Passed values array does not have correct dtype.") constants = _cpp.fem.pack_constants(self._cpp_object) - breakpoint() - coeffs = _cpp.fem.pack_coefficients(self._cpp_object, mesh,_entities) + coeffs = _cpp.fem.pack_coefficients(self._cpp_object, mesh._cpp_object, _entities) _cpp.fem.tabulate_expression( values, self._cpp_object, constants, coeffs, mesh._cpp_object, _entities ) diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py index c781eb4f261..b298a57a2a1 100644 --- a/python/test/unit/fem/test_expression.py +++ b/python/test/unit/fem/test_expression.py @@ -1,4 +1,4 @@ -# Copyright (C) 2019-2024 Michal Habera and Jørgen S. Dokken +# Copyright (C) 2019-2026 Michal Habera and Jørgen S. Dokken # # This file is part of DOLFINx (https://www.fenicsproject.org) # @@ -526,3 +526,50 @@ def test_rank1_blocked(): mask = np.ones(point_values.shape[2], dtype=bool) mask[offset::vs] = False np.testing.assert_allclose(point_values[i, j, mask], 0) + + +@pytest.mark.parametrize("qdegree", [1, 3, 5]) +@pytest.mark.parametrize( + "dtype", + [ + np.float32, + np.float64, + pytest.param(np.complex64, marks=pytest.mark.xfail_win32_complex), + pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex), + ], +) +def test_submesh_expression(dtype, qdegree): + xtype = dtype(0).real.dtype + mesh = create_unit_square(MPI.COMM_WORLD, 4, 3, dtype=xtype) + + V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) + u = dolfinx.fem.Function(V, dtype=dtype) + u.interpolate(lambda x: x[0] + 2.0 * x[1]) + + def mark_left_cells(x): + return x[0] <= 0.5 + 100 * np.finfo(xtype).resolution + + + left_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, mark_left_cells) + submesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, left_cells) + u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, ("Lagrange", 2)), dtype=dtype) + u_sub.interpolate(lambda x: x[1]**2 + x[0]**2) + + quadrature_points, _ = basix.make_quadrature(basix.CellType.triangle, qdegree) + quadrature_points = quadrature_points.astype(xtype) + + expr = dolfinx.fem.Expression(u*u_sub, quadrature_points, dtype=dtype, entity_maps=[entity_map]) + values = expr.eval(mesh, left_cells) + + num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local + values_sub = expr.eval(submesh, np.arange(num_sub_cells, dtype=np.int32)) + + tol = 50 * np.finfo(dtype).eps + np.testing.assert_allclose(values, values_sub, atol=tol) + + x = ufl.SpatialCoordinate(mesh) + u_exact = (x[0] + 2.0 * x[1]) * (x[1]**2 + x[0]**2) + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + expr_exact = dolfinx.fem.Expression(u_exact, quadrature_points, dtype=dtype) + values_exact = expr_exact.eval(mesh, np.arange(num_cells, dtype=np.int32)) + np.testing.assert_allclose(values, values_exact[left_cells], atol=tol) \ No newline at end of file From b9ddd8baee71604f51983f06a31d24656a150385 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:44:09 +0000 Subject: [PATCH 03/19] Add codim 1 test and fix bug --- cpp/dolfinx/fem/pack.h | 4 +- python/test/unit/fem/test_expression.py | 68 ++++++++++++++++++++++--- 2 files changed, 64 insertions(+), 8 deletions(-) diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 287a150ebea..41de76d24ab 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -421,10 +421,10 @@ void pack_coefficients( { const mesh::Topology topology = *mesh.topology(); int tdim = topology.dim(); - // If codim is zero we extract the cells - int codim = tdim - mesh.topology()->dim(); + int codim = tdim - mesh_c->topology()->dim(); if (codim == 0) { + // If codim is zero we extract the cells and map them auto cells = md::submdspan(entities, md::full_extent, 0); assert(cells.stride(0) == 1 && "Slice is not contiguous; cannot convert to std::span " diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py index b298a57a2a1..df899a8a2ec 100644 --- a/python/test/unit/fem/test_expression.py +++ b/python/test/unit/fem/test_expression.py @@ -15,7 +15,7 @@ from basix.ufl import quadrature_element from dolfinx import fem, la from dolfinx.fem import Constant, Expression, Function, form, functionspace -from dolfinx.mesh import create_unit_square +from dolfinx.mesh import create_rectangle, create_unit_square @pytest.mark.parametrize( @@ -549,16 +549,17 @@ def test_submesh_expression(dtype, qdegree): def mark_left_cells(x): return x[0] <= 0.5 + 100 * np.finfo(xtype).resolution - left_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, mark_left_cells) submesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, left_cells) u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, ("Lagrange", 2)), dtype=dtype) - u_sub.interpolate(lambda x: x[1]**2 + x[0]**2) + u_sub.interpolate(lambda x: x[1] ** 2 + x[0] ** 2) quadrature_points, _ = basix.make_quadrature(basix.CellType.triangle, qdegree) quadrature_points = quadrature_points.astype(xtype) - expr = dolfinx.fem.Expression(u*u_sub, quadrature_points, dtype=dtype, entity_maps=[entity_map]) + expr = dolfinx.fem.Expression( + u * u_sub, quadrature_points, dtype=dtype, entity_maps=[entity_map] + ) values = expr.eval(mesh, left_cells) num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local @@ -568,8 +569,63 @@ def mark_left_cells(x): np.testing.assert_allclose(values, values_sub, atol=tol) x = ufl.SpatialCoordinate(mesh) - u_exact = (x[0] + 2.0 * x[1]) * (x[1]**2 + x[0]**2) + u_exact = (x[0] + 2.0 * x[1]) * (x[1] ** 2 + x[0] ** 2) num_cells = mesh.topology.index_map(mesh.topology.dim).size_local expr_exact = dolfinx.fem.Expression(u_exact, quadrature_points, dtype=dtype) values_exact = expr_exact.eval(mesh, np.arange(num_cells, dtype=np.int32)) - np.testing.assert_allclose(values, values_exact[left_cells], atol=tol) \ No newline at end of file + np.testing.assert_allclose(values, values_exact[left_cells], atol=tol) + + +@pytest.mark.parametrize("qdegree", [1, 3, 5]) +@pytest.mark.parametrize( + "dtype", + [ + np.float32, + np.float64, + pytest.param(np.complex64, marks=pytest.mark.xfail_win32_complex), + pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex), + ], +) +def test_submesh_codim_one(dtype, qdegree): + xtype = dtype(0).real.dtype + mesh = create_rectangle( + MPI.COMM_WORLD, np.array([[0, 0], [2.1, 2.0]], dtype=xtype), [5, 3], dtype=xtype + ) + el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(), dtype=xtype) + V = dolfinx.fem.functionspace(mesh, el) + u = dolfinx.fem.Function(V, dtype=dtype) + u.interpolate(lambda x: x[0] + 2.0 * x[1]) + + tol = 50 * np.finfo(xtype).resolution + + def mark_left_facets(x): + return np.isclose(x[0], 1.0, atol=tol) + + left_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim - 1, mark_left_facets) + submesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, left_facets) + sub_el = basix.ufl.element( + "Lagrange", submesh.basix_cell(), 2, shape=(submesh.geometry.dim,), dtype=xtype + ) + u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, sub_el), dtype=dtype) + u_sub.interpolate(lambda x: (x[1] ** 2, -(x[0] ** 2))) + + quadrature_points, _ = basix.make_quadrature(basix.CellType.interval, qdegree) + quadrature_points = quadrature_points.astype(xtype) + + n_h = ufl.FacetNormal(mesh) + expr = u * (ufl.dot(u_sub, n_h) + n_h[0] * u_sub[1]) + + mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) + expr = dolfinx.fem.Expression(expr, quadrature_points, dtype=dtype, entity_maps=[entity_map]) + + entities = dolfinx.fem.compute_integration_domains( + dolfinx.fem.IntegralType.exterior_facet, mesh.topology, left_facets + ) + + values = expr.eval(mesh, entities.reshape(-1, 2)) + + x = ufl.SpatialCoordinate(mesh) + expr_exact = (x[0] + 2.0 * x[1]) * (x[1] ** 2 - x[0] ** 2) + expr_exact = dolfinx.fem.Expression(expr_exact, quadrature_points, dtype=dtype) + values_exact = expr_exact.eval(mesh, entities.reshape(-1, 2)) + np.testing.assert_allclose(values, values_exact, atol=tol) From 5c10ffa344f21fe4c5c4af43dea5ce553c0385b3 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:45:43 +0000 Subject: [PATCH 04/19] Change fffcx ref --- .github/workflows/fenicsx-refs.env | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fenicsx-refs.env b/.github/workflows/fenicsx-refs.env index 22c7e44a7e7..9bc01fe49cd 100644 --- a/.github/workflows/fenicsx-refs.env +++ b/.github/workflows/fenicsx-refs.env @@ -3,4 +3,4 @@ basix_ref=main ufl_repository=FEniCS/ufl ufl_ref=main ffcx_repository=FEniCS/ffcx -ffcx_ref=main +ffcx_ref=dokken/submesh-expression From 322b427dc9c4e75391720f82f64dfc881f29531d Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:47:45 +0000 Subject: [PATCH 05/19] Ruff formatting --- python/dolfinx/fem/function.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index d7699937ca6..2e0281d86be 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -204,7 +204,11 @@ def _create_expression(dtype): else: raise NotImplementedError(f"Type {dtype} not supported.") - _entity_maps = [entity_map._cpp_object for entity_map in entity_maps] if entity_maps is not None else [] + _entity_maps = ( + [entity_map._cpp_object for entity_map in entity_maps] + if entity_maps is not None + else [] + ) ffi = module.ffi self._cpp_object = _create_expression(dtype)( ffi.cast("uintptr_t", ffi.addressof(self._ufcx_expression)), From 0cb61fd9faa286d8ed170b78876c52ba6711e9bb Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:48:47 +0000 Subject: [PATCH 06/19] Mypy --- python/dolfinx/fem/function.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 2e0281d86be..f45b1a63438 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -26,7 +26,7 @@ if typing.TYPE_CHECKING: from mpi4py import MPI as _MPI - from dolfinx.mesh import Mesh + from dolfinx.mesh import Mesh, EntityMap class Constant(ufl.Constant): @@ -117,7 +117,7 @@ def __init__( form_compiler_options: dict | None = None, jit_options: dict | None = None, dtype: npt.DTypeLike | None = None, - entity_maps: list[dolfinx.mesh.EntityMap] | None = None, + entity_maps: list[EntityMap] | None = None, ): """Create an Expression. From bfc54828c3fdb5d2f1fcd7844fd4ad7bd90a10e3 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:50:49 +0000 Subject: [PATCH 07/19] Sort imports --- python/dolfinx/fem/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index f45b1a63438..1f23bafb8c9 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -26,7 +26,7 @@ if typing.TYPE_CHECKING: from mpi4py import MPI as _MPI - from dolfinx.mesh import Mesh, EntityMap + from dolfinx.mesh import EntityMap, Mesh class Constant(ufl.Constant): From 97155c21635c09a4b2531fc6be20ea725d6ae14e Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 11:53:31 +0000 Subject: [PATCH 08/19] Remove file --- python/test/unit/mesh/cform.py | 144 --------------------------------- 1 file changed, 144 deletions(-) delete mode 100644 python/test/unit/mesh/cform.py diff --git a/python/test/unit/mesh/cform.py b/python/test/unit/mesh/cform.py deleted file mode 100644 index 033feb25f0c..00000000000 --- a/python/test/unit/mesh/cform.py +++ /dev/null @@ -1,144 +0,0 @@ -import typing -from dataclasses import dataclass - -from mpi4py import MPI - -import numpy as np -import numpy.typing as npt - -import basix.ufl -import dolfinx -import ufl - - -@dataclass -class GeneralForm: - ufl_form: ufl.Form - ufcx_form: typing.Any - module: typing.Any - code: str - - -def compile_form( - comm: MPI.Intracomm, - form: ufl.Form, -) -> GeneralForm: - """Compile UFL form withtout associated DOLFINx data.""" - compiled_form = dolfinx.jit.ffcx_jit(comm, form) - return GeneralForm(form, *compiled_form) - - -def get_integration_domains(integral_type, subdomain): - """Get integration domains from subdomain data.""" - if subdomain is None: - return [] - else: - try: - if integral_type in ( - dolfinx.fem.IntegralType.exterior_facet, - dolfinx.fem.IntegralType.interior_facet, - ): - tdim = subdomain.topology.dim - subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) - subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) - domains = dolfinx.cpp.fem.compute_integration_domains( - integral_type, subdomain._cpp_object - ) - return [(s[0], np.array(s[1])) for s in domains] - except AttributeError: - return [(s[0], np.array(s[1])) for s in subdomain] - - -def form_cpp_creator( - dtype: npt.DTypeLike, -) -> ( - dolfinx.cpp.fem.Form_float32 - | dolfinx.cpp.fem.Form_float64 - | dolfinx.cpp.fem.Form_complex64 - | dolfinx.cpp.fem.Form_complex128 -): - """Return the wrapped C++ class of a variational form of a specific scalar type. - - Args: - dtype: Scalar type of the required form class. - - Returns: - Wrapped C++ form class of the requested type. - - Note: - This function is for advanced usage, typically when writing - custom kernels using Numba or C. - """ - if np.issubdtype(dtype, np.float32): - return dolfinx.cpp.fem.create_form_float32 - elif np.issubdtype(dtype, np.float64): - return dolfinx.cpp.fem.create_form_float64 - elif np.issubdtype(dtype, np.complex64): - return dolfinx.cpp.fem.create_form_complex64 - elif np.issubdtype(dtype, np.complex128): - return dolfinx.cpp.fem.create_form_complex128 - else: - raise NotImplementedError(f"Type {dtype} not supported.") - - -def create_form( - form: GeneralForm, - mesh: dolfinx.mesh.Mesh, - coefficient_map: dict[str, dolfinx.fem.Function], - constant_map: dict[str, dolfinx.fem.Constant], -): - sd = form.ufl_form.subdomain_data() - (domain,) = list(sd.keys()) # Assuming single domain - - # Subdomain markers (possibly empty list for some integral types) - subdomains = { - dolfinx.fem.forms._ufl_to_dolfinx_domain[key]: get_integration_domains( - dolfinx.fem.forms._ufl_to_dolfinx_domain[key], subdomain_data[0] - ) - for (key, subdomain_data) in sd.get(domain).items() - } - coefficients = {f"w{u.count()}": uh._cpp_object for (u, uh) in coefficient_map.items()} - constants = {f"c{c.count()}": ch._cpp_object for (c, ch) in constant_map.items()} - ftype = form_cpp_creator(dolfinx.default_scalar_type) - f = ftype( - form.module.ffi.cast("uintptr_t", form.module.ffi.addressof(form.ufcx_form)), - [], - coefficients, - constants, - subdomains, - mesh._cpp_object, - ) - return dolfinx.fem.Form(f, form.ufcx_form, form.code) - - -def test_compiled_form(): - # Create ufl form exclusively with basix.ufl and UFL - c_el = basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)) - domain = ufl.Mesh(c_el) - el = basix.ufl.element("Lagrange", "triangle", 2) - V = ufl.FunctionSpace(domain, el) - u = ufl.Coefficient(V) - w = ufl.Coefficient(V) - c = ufl.Constant(domain) - e = ufl.Constant(domain) - J = c * e * u * w * ufl.dx(domain=domain) - - # Compile form using dolfinx.jit.ffcx_jit - compiled_form = compile_form(MPI.COMM_WORLD, J) - - def create_and_integrate(N, compiled_form): - mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) - el_2 = basix.ufl.element("Lagrange", "triangle", 2) - Vh = dolfinx.fem.functionspace(mesh, el_2) - uh = dolfinx.fem.Function(Vh) - uh.interpolate(lambda x: x[0]) - wh = dolfinx.fem.Function(Vh) - wh.interpolate(lambda x: x[1]) - eh = dolfinx.fem.Constant(mesh, 3.0) - ch = dolfinx.fem.Constant(mesh, 2.0) - form = create_form(compiled_form, mesh, {u: uh, w: wh}, {c: ch, e: eh}) - assert np.isclose(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(form), op=MPI.SUM), 1.5) - - # Create various meshes, that all uses this compiled form with a map from ufl to dolfinx functions and constants - for i in range(1, 4): - create_and_integrate(i, compiled_form) From c32e457809685ae34f8d3349f1482fd374e0819d Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 12:37:39 +0000 Subject: [PATCH 09/19] Fix C++ demo and test and minor logic change for when to check for coordinate element hash --- cpp/demo/hyperelasticity/main.cpp | 2 +- cpp/dolfinx/fem/Expression.h | 16 ++++++++++++++-- cpp/dolfinx/fem/assembler.h | 13 +++++++++++++ cpp/dolfinx/fem/utils.h | 16 ++-------------- cpp/test/fem/form.cpp | 9 ++++++--- python/dolfinx/wrappers/dolfinx_wrappers/fem.h | 6 ++++-- 6 files changed, 40 insertions(+), 22 deletions(-) diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp index 35349f1cfbc..e0d40e051b0 100644 --- a/cpp/demo/hyperelasticity/main.cpp +++ b/cpp/demo/hyperelasticity/main.cpp @@ -334,7 +334,7 @@ int main(int argc, char* argv[]) // Compute Cauchy stress. Construct appropriate Basix element for // stress. fem::Expression sigma_expression = fem::create_expression( - *expression_hyperelasticity_sigma, {{"u", u}}, {}); + *expression_hyperelasticity_sigma, {{"u", u}}, {}, {}); constexpr auto family = basix::element::family::P; auto cell_type diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h index 4097d4547ae..f6087a92c42 100644 --- a/cpp/dolfinx/fem/Expression.h +++ b/cpp/dolfinx/fem/Expression.h @@ -82,14 +82,17 @@ class Expression const std::vector& value_shape, const std::vector>& entity_maps, + + std::uint64_t coordinate_element_hash, std::shared_ptr> argument_space = nullptr) : _argument_space(argument_space), _coefficients(coefficients), _constants(constants), _fn(fn), _value_shape(value_shape), _x_ref(std::vector(X.begin(), X.end()), Xshape), - _entity_maps(entity_maps) {}; + _entity_maps(entity_maps), + _coordinate_element_hash(coordinate_element_hash) {}; - /// Move constructorcreate_expression + /// Move constructor Expression(Expression&& e) = default; /// Destructor @@ -169,7 +172,16 @@ class Expression return _entity_maps; } + /// @brief Hash for coordinate element used to create the expression. + std::uint64_t coordinate_element_hash() const + { + return _coordinate_element_hash; + } + private: + // Hash for coordinate element used to create the expression + std::uint64_t _coordinate_element_hash; + // Function space for Argument std::shared_ptr> _argument_space; diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 9e2460b25c6..732fd988fe5 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -71,6 +71,12 @@ void tabulate_expression( std::pair>, std::size_t>> element) { + // Check that domain is the same as mesh of the expression + if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash()) + { + throw std::runtime_error( + "Expression was created on a different mesh. Cannot tabulate."); + } auto [X, Xshape] = e.X(); impl::tabulate_expression(values, e.kernel(), Xshape, e.value_size(), coeffs, constants, mesh, entities, element); @@ -96,6 +102,13 @@ template void tabulate_expression(std::span values, const fem::Expression& e, const mesh::Mesh& mesh, fem::MDSpan2 auto entities) { + // Check that domain is the same as mesh of the expression + if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash()) + { + throw std::runtime_error( + "Expression was created on a different mesh. Cannot tabulate."); + } + std::optional< std::pair>, std::size_t>> element = std::nullopt; diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 5d4e87ecf8d..386aac59ae6 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -957,19 +957,6 @@ Expression create_expression( entity_maps, std::shared_ptr> argument_space = nullptr) { - if (!coefficients.empty()) - { - assert(coefficients.front()); - assert(coefficients.front()->function_space()); - std::shared_ptr> mesh - = coefficients.front()->function_space()->mesh(); - if (mesh->geometry().cmap().hash() != e.coordinate_element_hash) - { - throw std::runtime_error( - "Expression and mesh geometric maps do not match."); - } - } - if (e.rank > 0 and !argument_space) { throw std::runtime_error("Expression has Argument but no Argument " @@ -1010,7 +997,8 @@ Expression create_expression( assert(tabulate_tensor); return Expression(coefficients, constants, std::span(X), Xshape, - tabulate_tensor, value_shape, entity_maps, argument_space); + tabulate_tensor, value_shape, entity_maps, + e.coordinate_element_hash, argument_space); } /// @brief Create Expression from UFC input (with named coefficients and diff --git a/cpp/test/fem/form.cpp b/cpp/test/fem/form.cpp index 348ae15c5eb..fe52fa1128e 100644 --- a/cpp/test/fem/form.cpp +++ b/cpp/test/fem/form.cpp @@ -35,7 +35,7 @@ void test_expression_cmap_compat(const auto& V) // Create Expression that expects P1 geometry dolfinx::fem::Expression expr1 = dolfinx::fem::create_expression(*expression_expr_Q6_P1, - {{"u1", u}}, {}); + {{"u1", u}}, {}, {}); auto [Xc, Xshape] = expr1.X(); std::vector grad_e(3 * Xshape[0] * cells.size()); fem::tabulate_expression(std::span(grad_e), expr1, *mesh, @@ -43,8 +43,11 @@ void test_expression_cmap_compat(const auto& V) // Create Expression that expects P2 geometry. Should throw because // mesh is P1. - CHECK_THROWS(dolfinx::fem::create_expression(*expression_expr_Q6_P2, - {{"u2", u}}, {})); + dolfinx::fem::Expression expr2 + = dolfinx::fem::create_expression(*expression_expr_Q6_P2, + {{"u2", u}}, {}, {}); + CHECK_THROWS(fem::tabulate_expression( + std::span(grad_e), expr2, *mesh, md::mdspan(cells.data(), cells.size()))); } } // namespace diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h index 018cffd35f0..84799555cbc 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/fem.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/fem.h @@ -684,6 +684,7 @@ void declare_objects(nb::module_& m, std::string type) std::uintptr_t fn_addr, const std::vector& value_shape, const std::vector& entity_maps, + std::uint64_t coordinate_element_hash, std::shared_ptr> argument_space) { @@ -694,11 +695,12 @@ void declare_objects(nb::module_& m, std::string type) new (ex) dolfinx::fem::Expression( coefficients, constants, std::span(X.data(), X.size()), {X.shape(0), X.shape(1)}, tabulate_expression_ptr, value_shape, - ptr_to_ref_wrapper_vec(entity_maps), argument_space); + ptr_to_ref_wrapper_vec(entity_maps), coordinate_element_hash, + argument_space); }, nb::arg("coefficients"), nb::arg("constants"), nb::arg("X"), nb::arg("fn"), nb::arg("value_shape"), nb::arg("entity_maps"), - nb::arg("argument_space")) + nb::arg("coordinate_element_hash"), nb::arg("argument_space")) .def("X", [](const dolfinx::fem::Expression& self) { From 2858d1235e219886179de436eeaf0d7a79f75d87 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 12:54:12 +0000 Subject: [PATCH 10/19] Try fixing create_expression --- cpp/dolfinx/fem/utils.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 386aac59ae6..763ea8cd8ea 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -996,9 +996,10 @@ Expression create_expression( throw std::runtime_error("Type not supported."); assert(tabulate_tensor); + std::uint64_t e_hash = e.coordinate_element_hash; return Expression(coefficients, constants, std::span(X), Xshape, - tabulate_tensor, value_shape, entity_maps, - e.coordinate_element_hash, argument_space); + tabulate_tensor, value_shape, entity_maps, e_hash, + argument_space); } /// @brief Create Expression from UFC input (with named coefficients and From 0720d7c02b13304920d34f1e7a2fa5229dfd7eec Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 12:55:50 +0000 Subject: [PATCH 11/19] Reorder attributes --- cpp/dolfinx/fem/Expression.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h index f6087a92c42..4a96e071547 100644 --- a/cpp/dolfinx/fem/Expression.h +++ b/cpp/dolfinx/fem/Expression.h @@ -179,9 +179,6 @@ class Expression } private: - // Hash for coordinate element used to create the expression - std::uint64_t _coordinate_element_hash; - // Function space for Argument std::shared_ptr> _argument_space; @@ -203,7 +200,11 @@ class Expression // Evaluation points on reference cell std::pair, std::array> _x_ref; + // Map between different mesh topologies std::vector> _entity_maps; + + // Hash for coordinate element used to create the expression + std::uint64_t _coordinate_element_hash; }; } // namespace dolfinx::fem From 6067f80cb9c921350f916c35551b66596cfe6e40 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 13:49:56 +0000 Subject: [PATCH 12/19] Add docs --- cpp/dolfinx/fem/Expression.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/dolfinx/fem/Expression.h b/cpp/dolfinx/fem/Expression.h index 4a96e071547..c637f9c0869 100644 --- a/cpp/dolfinx/fem/Expression.h +++ b/cpp/dolfinx/fem/Expression.h @@ -65,6 +65,8 @@ class Expression /// @param[in] value_shape Shape of Expression evaluated at single /// point. /// @param[in] entity_maps Maps between entities of different meshes. + /// @param[in] coordinate_element_hash Hash for coordinate element used + /// to create the expression. /// @param[in] argument_space Function space for Argument. Used to /// computed a 1-form expression, e.g. can be used to create a matrix /// that when applied to a degree-of-freedom vector gives the From 4487a75447f3c73f570ad4d6f370442fb1ede067 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 14:07:14 +0000 Subject: [PATCH 13/19] Fix import for docs --- python/dolfinx/fem/function.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 1f23bafb8c9..c64cf90d8c9 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -26,7 +26,8 @@ if typing.TYPE_CHECKING: from mpi4py import MPI as _MPI - from dolfinx.mesh import EntityMap, Mesh + from dolfinx.mesh import EntityMap as _EntityMap + from dolfinx.mesh import Mesh class Constant(ufl.Constant): @@ -117,7 +118,7 @@ def __init__( form_compiler_options: dict | None = None, jit_options: dict | None = None, dtype: npt.DTypeLike | None = None, - entity_maps: list[EntityMap] | None = None, + entity_maps: list[_EntityMap] | None = None, ): """Create an Expression. From de056d586d2c6e0fb3d96f2d480f9469b085c8d8 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 18:06:57 +0000 Subject: [PATCH 14/19] Copy data as it is strided when extracted from python --- cpp/dolfinx/fem/pack.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 41de76d24ab..3706f3fb2c8 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -426,11 +426,11 @@ void pack_coefficients( { // If codim is zero we extract the cells and map them auto cells = md::submdspan(entities, md::full_extent, 0); - assert(cells.stride(0) == 1 - && "Slice is not contiguous; cannot convert to std::span " - "safely!"); - e_b = emap.sub_topology_to_topology( - std::span(cells.data_handle(), cells.size()), inverse); + std::vector contiguous_cells(cells.extent(0)); + for (std::size_t i = 0; i < cells.extent(0); ++i) + contiguous_cells[i] = cells(i); + e_b = emap.sub_topology_to_topology(std::span(contiguous_cells), + inverse); md::mdspan e(e_b.data(), e_b.size()); impl::pack_coefficient_entity(std::span(c), cstride, coeffs[coeff].get(), cell_info, e, From 9322220a5f1601241186201491e7dcc98401d39c Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 18:48:02 +0000 Subject: [PATCH 15/19] Fix for parallel test --- python/test/unit/fem/test_expression.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py index df899a8a2ec..0657bb5c922 100644 --- a/python/test/unit/fem/test_expression.py +++ b/python/test/unit/fem/test_expression.py @@ -538,7 +538,7 @@ def test_rank1_blocked(): pytest.param(np.complex128, marks=pytest.mark.xfail_win32_complex), ], ) -def test_submesh_expression(dtype, qdegree): +def test_submesh_codim_zero(dtype, qdegree): xtype = dtype(0).real.dtype mesh = create_unit_square(MPI.COMM_WORLD, 4, 3, dtype=xtype) @@ -562,7 +562,8 @@ def mark_left_cells(x): ) values = expr.eval(mesh, left_cells) - num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local + sub_cellmap = submesh.topology.index_map(submesh.topology.dim) + num_sub_cells = sub_cellmap.size_local + sub_cellmap.num_ghosts values_sub = expr.eval(submesh, np.arange(num_sub_cells, dtype=np.int32)) tol = 50 * np.finfo(dtype).eps @@ -570,7 +571,8 @@ def mark_left_cells(x): x = ufl.SpatialCoordinate(mesh) u_exact = (x[0] + 2.0 * x[1]) * (x[1] ** 2 + x[0] ** 2) - num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + cell_map = mesh.topology.index_map(mesh.topology.dim) + num_cells = cell_map.size_local + cell_map.num_ghosts expr_exact = dolfinx.fem.Expression(u_exact, quadrature_points, dtype=dtype) values_exact = expr_exact.eval(mesh, np.arange(num_cells, dtype=np.int32)) np.testing.assert_allclose(values, values_exact[left_cells], atol=tol) From 05d71d4c2640fc864637442c8b14770766569535 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 1 Apr 2026 18:54:47 +0000 Subject: [PATCH 16/19] Further parallel fixes --- python/test/unit/fem/test_expression.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/python/test/unit/fem/test_expression.py b/python/test/unit/fem/test_expression.py index 0657bb5c922..c5804480150 100644 --- a/python/test/unit/fem/test_expression.py +++ b/python/test/unit/fem/test_expression.py @@ -557,13 +557,17 @@ def mark_left_cells(x): quadrature_points, _ = basix.make_quadrature(basix.CellType.triangle, qdegree) quadrature_points = quadrature_points.astype(xtype) + sub_cellmap = submesh.topology.index_map(submesh.topology.dim) + num_sub_cells = sub_cellmap.size_local + sub_cellmap.num_ghosts + parent_cells = entity_map.sub_topology_to_topology( + np.arange(num_sub_cells, dtype=np.int32), inverse=False + ) + expr = dolfinx.fem.Expression( u * u_sub, quadrature_points, dtype=dtype, entity_maps=[entity_map] ) - values = expr.eval(mesh, left_cells) + values = expr.eval(mesh, parent_cells) - sub_cellmap = submesh.topology.index_map(submesh.topology.dim) - num_sub_cells = sub_cellmap.size_local + sub_cellmap.num_ghosts values_sub = expr.eval(submesh, np.arange(num_sub_cells, dtype=np.int32)) tol = 50 * np.finfo(dtype).eps @@ -575,7 +579,7 @@ def mark_left_cells(x): num_cells = cell_map.size_local + cell_map.num_ghosts expr_exact = dolfinx.fem.Expression(u_exact, quadrature_points, dtype=dtype) values_exact = expr_exact.eval(mesh, np.arange(num_cells, dtype=np.int32)) - np.testing.assert_allclose(values, values_exact[left_cells], atol=tol) + np.testing.assert_allclose(values, values_exact[parent_cells], atol=tol) @pytest.mark.parametrize("qdegree", [1, 3, 5]) From 4a7c72e491c06c8f097d0a2aa8c6d8af780cbdd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Fri, 17 Apr 2026 15:47:06 +0200 Subject: [PATCH 17/19] Revert to main branch of FFCx MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- .github/workflows/fenicsx-refs.env | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fenicsx-refs.env b/.github/workflows/fenicsx-refs.env index 9bc01fe49cd..22c7e44a7e7 100644 --- a/.github/workflows/fenicsx-refs.env +++ b/.github/workflows/fenicsx-refs.env @@ -3,4 +3,4 @@ basix_ref=main ufl_repository=FEniCS/ufl ufl_ref=main ffcx_repository=FEniCS/ffcx -ffcx_ref=dokken/submesh-expression +ffcx_ref=main From e7161a865dd8161ba53fb7b63f2c7c82fad934cf Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Thu, 23 Apr 2026 07:13:32 +0000 Subject: [PATCH 18/19] Update API --- cpp/dolfinx/fem/assembler.h | 4 ++-- cpp/dolfinx/mesh/Geometry.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 9a6a41ea765..9f02223844e 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -72,7 +72,7 @@ void tabulate_expression( element) { // Check that domain is the same as mesh of the expression - if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash()) + if (e.coordinate_element_hash() != mesh.geometry().cmap(0).hash()) { throw std::runtime_error( "Expression was created on a different mesh. Cannot tabulate."); @@ -103,7 +103,7 @@ void tabulate_expression(std::span values, const fem::Expression& e, const mesh::Mesh& mesh, fem::MDSpan2 auto entities) { // Check that domain is the same as mesh of the expression - if (e.coordinate_element_hash() != mesh.geometry().cmaps()[0].hash()) + if (e.coordinate_element_hash() != mesh.geometry().cmap(0).hash()) { throw std::runtime_error( "Expression was created on a different mesh. Cannot tabulate."); diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index 102023258e9..650ece0b7d3 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -120,7 +120,7 @@ class Geometry /// map element in the geometry. /// @param[in] i Index of the requested degree-of-freedom map. The /// degree-of-freedom map corresponds to the geometry element - /// `cmaps()[i]`. + /// `cmap(i)`. /// @return A dofmap array, with shape `(num_cells, dofs_per_cell)`. md::mdspan> dofmap(std::size_t i) const @@ -152,7 +152,7 @@ class Geometry /// @brief The elements that describes the geometry map. /// - /// The coordinate element `cmaps()[i]` corresponds to the + /// The coordinate element `cmap(i)` corresponds to the /// degree-of-freedom map `dofmap(i)`. /// @param i The coordinate map to fetch /// @return The coordinate/geometry elements. From ecfa3a9175a1efb59079345b2469703ab734ab89 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Tue, 5 May 2026 12:19:07 +0200 Subject: [PATCH 19/19] Add some typing --- python/dolfinx/fem/function.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 3df7c7e5388..5a30ef2716b 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -112,6 +112,16 @@ class Expression(Generic[Scalar]): """ + _ufl_expression: ufl.core.expr.Expr + _argument_space: FunctionSpace | None + _cpp_object: ( + _cpp.fem.Expression_complex64 + | _cpp.fem.Expression_complex128 + | _cpp.fem.Expression_float32 + | _cpp.fem.Expression_float64 + ) + _code: str + def __init__( self, e: ufl.core.expr.Expr,