Fix exterior facet integrals from meshtags (#4041)

* Add test that fails on main branch

* Add fix

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)