Fix exterior facet integrals from meshtags (#4041)

* Add test that fails on main branch

* Add fix

* Ruff formatting

Author: jorgensd

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]

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)

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]))