From 7e26002f28096d395fe1ed3645c0750bf36fd58c Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 3 Jun 2026 14:15:08 +0000 Subject: [PATCH 1/6] Fix expressions on facets --- src/scifem/bcs.py | 45 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 39 insertions(+), 6 deletions(-) diff --git a/src/scifem/bcs.py b/src/scifem/bcs.py index b657c53..5d6ea0a 100644 --- a/src/scifem/bcs.py +++ b/src/scifem/bcs.py @@ -3,10 +3,37 @@ import ufl import numpy.typing as npt import numpy as np +import itertools +from ffcx.ir.elementtables import ( + permute_quadrature_interval, + permute_quadrature_triangle, + permute_quadrature_quadrilateral, +) __all__ = ["interpolate_function_onto_facet_dofs"] +def build_quadrature_permutations(facet_type, points): + if facet_type == basix.CellType.interval: + return [permute_quadrature_interval(points, ref) for ref in range(2)] + elif facet_type == basix.CellType.triangle: + return itertools.chain.from_iterable( + [ + [permute_quadrature_triangle(points, ref, rot) for ref in range(3)] + for rot in range(2) + ] + ) + elif facet_type == basix.CellType.quadrilateral: + return itertools.chain.from_iterable( + [ + [permute_quadrature_quadrilateral(points, ref, rot) for ref in range(2)] + for rot in range(4) + ] + ) + else: + raise ValueError(f"Unsupported {facet_type=}") + + def interpolate_function_onto_facet_dofs( Q: dolfinx.fem.FunctionSpace, expr: ufl.core.expr.Expr, @@ -49,9 +76,10 @@ def interpolate_function_onto_facet_dofs( ) ref_top = ref_cmap.reference_topology ref_geom = ref_cmap.reference_geometry + facet_type = facet_types.pop() facet_cmap = basix.ufl.element( "Lagrange", - facet_types.pop(), + facet_type, 1, shape=(domain.topology.dim,), dtype=np.float64, @@ -71,9 +99,10 @@ def interpolate_function_onto_facet_dofs( else: assert np.allclose(reference_facet_points, ref_points) assert reference_facet_points is not None - # Create expression for BC - normal_expr = dolfinx.fem.Expression(expr, reference_facet_points) - + facet_points = build_quadrature_permutations(facet_type, reference_facet_points) + expressions = [] + for i, points in enumerate(facet_points): + expressions.append(dolfinx.fem.Expression(expr, points)) points_per_entity = [sum(ip.shape[0] for ip in ips) for ips in interpolation_points] offsets = np.zeros(domain.topology.dim + 2, dtype=np.int32) offsets[1:] = np.cumsum(points_per_entity[: domain.topology.dim + 1]) @@ -93,6 +122,9 @@ def interpolate_function_onto_facet_dofs( ) is_marked = np.zeros(num_facets_on_process, dtype=np.int8) is_marked[facets] = 1 + domain.topology.create_entity_permutations() + num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(domain.topology.cell_type, fdim) + facet_permutations = domain.topology.get_facet_permutations().reshape(-1, num_facets_per_cell) for i, cell in enumerate(all_connected_cells): values_per_entity[:] = 0.0 local_facets = c_to_f.links(cell) @@ -102,10 +134,11 @@ def interpolate_function_onto_facet_dofs( insert_pos = offsets[fdim] + reference_facet_points.shape[0] * j # Backwards compatibility entity = np.array([[cell, j]], dtype=np.int32) + perm = facet_permutations[cell, j] try: - normal_on_facet = normal_expr.eval(domain, entity) + normal_on_facet = expressions[perm].eval(domain, entity) except (AttributeError, AssertionError): - normal_on_facet = normal_expr.eval(domain, entity.flatten()) + normal_on_facet = expressions[perm].eval(domain, entity.flatten()) # NOTE: evaluate within loop to avoid large memory requirements values_per_entity[insert_pos : insert_pos + reference_facet_points.shape[0]] = ( normal_on_facet.reshape(-1, domain.geometry.dim) From ff548588b1fd06cbfc6e049f42854c5dbb6a542d Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 3 Jun 2026 15:13:38 +0000 Subject: [PATCH 2/6] Fix permutation logic by reverting it so we get the local expression --- src/scifem/bcs.py | 55 +++++++++++++++++++++++++++-------------------- tests/test_bcs.py | 8 +++---- 2 files changed, 36 insertions(+), 27 deletions(-) diff --git a/src/scifem/bcs.py b/src/scifem/bcs.py index 5d6ea0a..8820848 100644 --- a/src/scifem/bcs.py +++ b/src/scifem/bcs.py @@ -3,7 +3,7 @@ import ufl import numpy.typing as npt import numpy as np -import itertools + from ffcx.ir.elementtables import ( permute_quadrature_interval, permute_quadrature_triangle, @@ -17,19 +17,25 @@ def build_quadrature_permutations(facet_type, points): if facet_type == basix.CellType.interval: return [permute_quadrature_interval(points, ref) for ref in range(2)] elif facet_type == basix.CellType.triangle: - return itertools.chain.from_iterable( - [ - [permute_quadrature_triangle(points, ref, rot) for ref in range(3)] - for rot in range(2) - ] - ) + perms = [] + # FFCx order: rot is outer loop, ref is inner loop + for rot in range(3): + for ref in range(2): + # Counteract the mapping with the inverse permutation + rot_inv = (3 - rot) % 3 if ref == 0 else rot + perms.append(permute_quadrature_triangle(points, ref, rot_inv)) + return perms + elif facet_type == basix.CellType.quadrilateral: - return itertools.chain.from_iterable( - [ - [permute_quadrature_quadrilateral(points, ref, rot) for ref in range(2)] - for rot in range(4) - ] - ) + perms = [] + # FFCx order: rot is outer loop, ref is inner loop + for rot in range(4): + for ref in range(2): + # Counteract the mapping with the inverse permutation + rot_inv = (4 - rot) % 4 if ref == 0 else rot + perms.append(permute_quadrature_quadrilateral(points, ref, rot_inv)) + return perms + else: raise ValueError(f"Unsupported {facet_type=}") @@ -87,7 +93,6 @@ def interpolate_function_onto_facet_dofs( facet_cel = dolfinx.fem.CoordinateElement( dolfinx.cpp.fem.CoordinateElement_float64(facet_cmap.basix_element._e) ) - reference_facet_points = None for i, points in enumerate(interpolation_points[fdim]): geom = ref_geom[ref_top[fdim][i]] @@ -101,8 +106,8 @@ def interpolate_function_onto_facet_dofs( assert reference_facet_points is not None facet_points = build_quadrature_permutations(facet_type, reference_facet_points) expressions = [] - for i, points in enumerate(facet_points): - expressions.append(dolfinx.fem.Expression(expr, points)) + for i, perm in enumerate(facet_points): + expressions.append(dolfinx.fem.Expression(expr, perm)) points_per_entity = [sum(ip.shape[0] for ip in ips) for ips in interpolation_points] offsets = np.zeros(domain.topology.dim + 2, dtype=np.int32) offsets[1:] = np.cumsum(points_per_entity[: domain.topology.dim + 1]) @@ -114,7 +119,12 @@ def interpolate_function_onto_facet_dofs( all_connected_cells = dolfinx.mesh.compute_incident_entities( domain.topology, facets, domain.topology.dim - 1, domain.topology.dim ) - values = np.zeros(len(all_connected_cells) * offsets[-1] * domain.geometry.dim) + expr_value_size = expressions[0].value_size + + # Update array allocations to use the exact expression value size + values_per_entity = np.zeros((offsets[-1], expr_value_size), dtype=dolfinx.default_scalar_type) + values = np.zeros(len(all_connected_cells) * offsets[-1] * expr_value_size) + domain.topology.create_connectivity(domain.topology.dim, fdim) c_to_f = domain.topology.connectivity(domain.topology.dim, fdim) num_facets_on_process = ( @@ -141,18 +151,17 @@ def interpolate_function_onto_facet_dofs( normal_on_facet = expressions[perm].eval(domain, entity.flatten()) # NOTE: evaluate within loop to avoid large memory requirements values_per_entity[insert_pos : insert_pos + reference_facet_points.shape[0]] = ( - normal_on_facet.reshape(-1, domain.geometry.dim) + normal_on_facet.reshape(-1, expr_value_size) ) - values[ - i * offsets[-1] * domain.geometry.dim : (i + 1) * offsets[-1] * domain.geometry.dim - ] = values_per_entity.reshape(-1) + values[i * offsets[-1] * expr_value_size : (i + 1) * offsets[-1] * expr_value_size] = ( + values_per_entity.reshape(-1) + ) qh = dolfinx.fem.Function(Q) if hasattr(qh._cpp_object, "interpolate_f"): interpolate_func = qh._cpp_object.interpolate_f else: interpolate_func = qh._cpp_object.interpolate - interpolate_func(values.reshape(-1, domain.geometry.dim).T.copy(), all_connected_cells) + interpolate_func(values.reshape(-1, expr_value_size).T.copy(), all_connected_cells) qh.x.scatter_forward() - return qh diff --git a/tests/test_bcs.py b/tests/test_bcs.py index 4d0fa34..b6edab2 100644 --- a/tests/test_bcs.py +++ b/tests/test_bcs.py @@ -22,9 +22,9 @@ def right_facets(x): def test_normal_enforcement(cell_type: dolfinx.mesh.CellType): tdim = dolfinx.mesh.cell_dim(cell_type) if tdim == 2: - mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, cell_type=cell_type) elif tdim == 3: - mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3) + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3, cell_type=cell_type) V = dolfinx.fem.functionspace(mesh, ("BDM", 1)) @@ -51,9 +51,9 @@ def test_normal_enforcement(cell_type: dolfinx.mesh.CellType): def test_tangent_enforcement(cell_type: dolfinx.mesh.CellType): tdim = dolfinx.mesh.cell_dim(cell_type) if tdim == 2: - mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, cell_type=cell_type) elif tdim == 3: - mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3) + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3, cell_type=cell_type) else: raise ValueError(f"Unsupported {cell_type=}") From 87be35b68937a60c46515b2f4a96d023b13324af Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 3 Jun 2026 15:50:34 +0000 Subject: [PATCH 3/6] Extend interpolation to DG --- src/scifem/interpolation.py | 28 +++++++++++++++++----------- tests/test_interpolation.py | 5 ++++- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index bc6c720..223efe4 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -337,7 +337,10 @@ def interpolate_to_surface_submesh( """ Interpolate a function `u_volume` into the function `u_surface`. Note: - Does not work for DG as no dofs are associated with the facets + Does not work for DG as no dofs are associated with the facets in versions of DOLFINx + prior to https://github.com/FEniCS/dolfinx/pull/4140, which is included in version + 0.11.0 and later. + Args: u_volume: Function to interpolate data from u_surface: Function to interpolate data to @@ -362,18 +365,21 @@ def interpolate_to_surface_submesh( submesh.topology.create_entity_permutations() mesh.topology.create_entity_permutations() - cell_info = mesh.topology.get_cell_permutation_info() ft = V_surf.element.basix_element.cell_type - for i in range(integration_entities.shape[0]): - perm = np.arange(data.shape[1], dtype=np.int32) - V_vol.element.basix_element.permute_subentity_closure_inv( - perm, - cell_info[integration_entities[i, 0]], - ft, - int(integration_entities[i, 1]), - ) - data[i] = data[i][perm] + if Version(dolfinx.__version__) < Version("0.11.0.dev0"): + # Before the introduction of https://github.com/FEniCS/dolfinx/pull/4140 + # one needed to permute the data according to the facet permutations. + cell_info = mesh.topology.get_cell_permutation_info() + for i in range(integration_entities.shape[0]): + perm = np.arange(data.shape[1], dtype=np.int32) + V_vol.element.basix_element.permute_subentity_closure_inv( + perm, + cell_info[integration_entities[i, 0]], + ft, + int(integration_entities[i, 1]), + ) + data[i] = data[i][perm] if len(data.shape) == 3: # Data is now (num_cells, value_size,num_points) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 3d48732..2ffa93d 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -231,13 +231,16 @@ def test_discrete_curl(degree, use_petsc, cell_type): @pytest.mark.parametrize("degree", [1, 2, 3]) -@pytest.mark.parametrize("family", ["Lagrange"]) +@pytest.mark.parametrize("family", ["Lagrange", "DG"]) @pytest.mark.skipif( Version(dolfinx.__version__) < Version("0.10.0"), reason="Requires DOLFINx version >0.10.0" ) def test_interpolate_to_interface_submesh(family, degree): # Create a unit square comm = MPI.COMM_WORLD + + if Version(dolfinx.__version__) < Version("0.11.0.dev0") and family == "DG": + pytest.skip("Interpolation to surface submesh does not work for DG in DOLFINx < 0.11.0") domain = dolfinx.mesh.create_unit_square( comm, 48, 48, ghost_mode=dolfinx.mesh.GhostMode.shared_facet ) From b1a98ca2f283b909c7a32190f63300eadb686724 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 3 Jun 2026 15:50:46 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/scifem/interpolation.py | 2 +- tests/test_interpolation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index 223efe4..07336c3 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -369,7 +369,7 @@ def interpolate_to_surface_submesh( if Version(dolfinx.__version__) < Version("0.11.0.dev0"): # Before the introduction of https://github.com/FEniCS/dolfinx/pull/4140 - # one needed to permute the data according to the facet permutations. + # one needed to permute the data according to the facet permutations. cell_info = mesh.topology.get_cell_permutation_info() for i in range(integration_entities.shape[0]): perm = np.arange(data.shape[1], dtype=np.int32) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 2ffa93d..8fa9647 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -240,7 +240,7 @@ def test_interpolate_to_interface_submesh(family, degree): comm = MPI.COMM_WORLD if Version(dolfinx.__version__) < Version("0.11.0.dev0") and family == "DG": - pytest.skip("Interpolation to surface submesh does not work for DG in DOLFINx < 0.11.0") + pytest.skip("Interpolation to surface submesh does not work for DG in DOLFINx < 0.11.0") domain = dolfinx.mesh.create_unit_square( comm, 48, 48, ghost_mode=dolfinx.mesh.GhostMode.shared_facet ) From 69a7af638d6ef9f0450c7ac33872a4e792031358 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 3 Jun 2026 19:47:09 +0200 Subject: [PATCH 5/6] Add backward compat --- src/scifem/bcs.py | 61 ++++++++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 24 deletions(-) diff --git a/src/scifem/bcs.py b/src/scifem/bcs.py index 8820848..eef5f41 100644 --- a/src/scifem/bcs.py +++ b/src/scifem/bcs.py @@ -3,7 +3,7 @@ import ufl import numpy.typing as npt import numpy as np - +from packaging import Version from ffcx.ir.elementtables import ( permute_quadrature_interval, permute_quadrature_triangle, @@ -14,30 +14,43 @@ def build_quadrature_permutations(facet_type, points): - if facet_type == basix.CellType.interval: - return [permute_quadrature_interval(points, ref) for ref in range(2)] - elif facet_type == basix.CellType.triangle: - perms = [] - # FFCx order: rot is outer loop, ref is inner loop - for rot in range(3): - for ref in range(2): - # Counteract the mapping with the inverse permutation - rot_inv = (3 - rot) % 3 if ref == 0 else rot - perms.append(permute_quadrature_triangle(points, ref, rot_inv)) - return perms - - elif facet_type == basix.CellType.quadrilateral: - perms = [] - # FFCx order: rot is outer loop, ref is inner loop - for rot in range(4): - for ref in range(2): - # Counteract the mapping with the inverse permutation - rot_inv = (4 - rot) % 4 if ref == 0 else rot - perms.append(permute_quadrature_quadrilateral(points, ref, rot_inv)) - return perms - + if Version(dolfinx.__version__) < Version("0.11.0.dev0"): + # In older versions of dolfinx, the permutation is handled internally + # in the C++ code, so we can just return the original points + if facet_type == basix.CellType.interval: + num_permutations = 2 + elif facet_type == basix.CellType.triangle: + num_permutations = 6 + elif facet_type == basix.CellType.quadrilateral: + num_permutations = 8 + else: + raise ValueError(f"Unsupported {facet_type=}") + return [points for _ in range(num_permutations)] else: - raise ValueError(f"Unsupported {facet_type=}") + if facet_type == basix.CellType.interval: + return [permute_quadrature_interval(points, ref) for ref in range(2)] + elif facet_type == basix.CellType.triangle: + perms = [] + # FFCx order: rot is outer loop, ref is inner loop + for rot in range(3): + for ref in range(2): + # Counteract the mapping with the inverse permutation + rot_inv = (3 - rot) % 3 if ref == 0 else rot + perms.append(permute_quadrature_triangle(points, ref, rot_inv)) + return perms + + elif facet_type == basix.CellType.quadrilateral: + perms = [] + # FFCx order: rot is outer loop, ref is inner loop + for rot in range(4): + for ref in range(2): + # Counteract the mapping with the inverse permutation + rot_inv = (4 - rot) % 4 if ref == 0 else rot + perms.append(permute_quadrature_quadrilateral(points, ref, rot_inv)) + return perms + + else: + raise ValueError(f"Unsupported {facet_type=}") def interpolate_function_onto_facet_dofs( From 11c0204df42e8d674533bca38a4c5f3a037ae759 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 3 Jun 2026 21:03:59 +0200 Subject: [PATCH 6/6] Fix import --- src/scifem/bcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scifem/bcs.py b/src/scifem/bcs.py index eef5f41..63df756 100644 --- a/src/scifem/bcs.py +++ b/src/scifem/bcs.py @@ -3,7 +3,7 @@ import ufl import numpy.typing as npt import numpy as np -from packaging import Version +from packaging.version import Version from ffcx.ir.elementtables import ( permute_quadrature_interval, permute_quadrature_triangle,