From f0431096d3c2abba58c0435bd78e2349c1bf8299 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 15 Jan 2026 08:31:39 +0000 Subject: [PATCH 1/3] Add test that fails on main branch --- python/test/unit/fem/test_assemble_domains.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/python/test/unit/fem/test_assemble_domains.py b/python/test/unit/fem/test_assemble_domains.py index 5b09171b5a..0c5bc1776c 100644 --- a/python/test/unit/fem/test_assemble_domains.py +++ b/python/test/unit/fem/test_assemble_domains.py @@ -349,3 +349,21 @@ def create_forms(dx, ds, dS): assert np.allclose(A.data, A_mt.data) assert np.allclose(b.array, b_mt.array) + + +def test_assemble_exterior_facet(): + """Check special handling of packing of integration entities for exterior facets, + which for any other co-dimensional entity is just a one-sided integral. + """ + domain = create_unit_square(MPI.COMM_WORLD, 2, 2) + fdim = domain.topology.dim - 1 + tag = 1 + + # Only tag interior facets + facets = locate_entities(domain, fdim, lambda x: np.isclose(x[0], 0.5)) + ft = meshtags(domain, fdim, facets, np.full_like(facets, tag)) + ds = ufl.Measure("ds", domain=domain, subdomain_data=ft) + + # Check that integral is 0 + value = domain.comm.allreduce(assemble_scalar(form(1.0 * ds(tag))), op=MPI.SUM) + assert np.isclose(value, 0.0) From 5dcb3f448c726d397bc244287af61bc4907437bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 15 Jan 2026 08:31:57 +0000 Subject: [PATCH 2/3] Add fix --- python/dolfinx/fem/forms.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index cdb4fdeecd..3c994db79c 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -176,14 +176,24 @@ def get_integration_domains( subdomain._cpp_object.topology.create_connectivity(tdim - 2, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 2) + # Special handling for exterior facets, compared to other + # one-sided entity integrals + if integral_type is IntegralType.exterior_facet: + exterior_facets = _cpp.mesh.exterior_facet_indices(subdomain.topology) + # Compute integration domains only for each subdomain id in # the integrals. If a process has no integral entities, # insert an empty array. for id in subdomain_ids: + entities = subdomain.find(id) + if integral_type is IntegralType.exterior_facet: + # Compute intersection of tag an exterior facets + entities = np.intersect1d(entities, exterior_facets) + integration_entities = _cpp.fem.compute_integration_domains( integral_type, subdomain._cpp_object.topology, - subdomain.find(id), + entities, ) domains.append((id, integration_entities)) return [(s[0], np.array(s[1])) for s in domains] From e8f826be361536f1a632417bfada0a088c4fded2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 4 Feb 2026 15:22:50 +0100 Subject: [PATCH 3/3] Ruff formatting --- python/test/unit/fem/test_element_integrals.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/python/test/unit/fem/test_element_integrals.py b/python/test/unit/fem/test_element_integrals.py index 7397833710..e6c543700e 100644 --- a/python/test/unit/fem/test_element_integrals.py +++ b/python/test/unit/fem/test_element_integrals.py @@ -223,9 +223,11 @@ def test_facet_integral(cell_type, dtype): elif cell_type == CellType.tetrahedron: s = 2**0.5 * 3 ** (1 / 3) # side length v.interpolate( - lambda x: (x[0] - s / 2) ** 2 - + (x[1] - s / 2 / np.sqrt(3)) ** 2 - + (x[2] - s * np.sqrt(2 / 3) / 4) ** 2 + lambda x: ( + (x[0] - s / 2) ** 2 + + (x[1] - s / 2 / np.sqrt(3)) ** 2 + + (x[2] - s * np.sqrt(2 / 3) / 4) ** 2 + ) ) elif cell_type == CellType.hexahedron: v.interpolate(lambda x: x[0] * (1 - x[0]) + x[1] * (1 - x[1]) + x[2] * (1 - x[2]))