Hiya!

I’ve been toying around with bulk-surface PDEs and on my journey through the FEniCS jungle I ran into an inconsistency that I cannot explain. Here’s a MWE to illustrate my confusion:

```
from dolfinx.fem import functionspace, Function, form, locate_dofs_topological, assemble_scalar
from dolfinx.mesh import create_submesh, create_unit_square, locate_entities_boundary, meshtags
from mpi4py import MPI
from ufl import dx, Measure
import numpy as np
# Marker for the surface where "w" lives
def on_gamma(x):
return np.isclose(x[1], 0.0)
bulk_mesh = create_unit_square(MPI.COMM_WORLD, 20, 10)
dim = bulk_mesh.topology.dim
surface_entities = locate_entities_boundary(bulk_mesh, dim - 1, on_gamma)
surface_mesh, entity_map, vertex_map, geometry_map = create_submesh(bulk_mesh, dim - 1, surface_entities)
gamma_tags = meshtags(bulk_mesh, dim-1, surface_entities, np.full_like(surface_entities, 1))
element = ("Lagrange", 1)
Vb = functionspace(bulk_mesh, element)
Vs = functionspace(surface_mesh, element)
def u_func(x):
return np.square(x[0]) + np.square(x[1])
def w_func(x):
return x[0]
u = Function(Vb)
u.interpolate(u_func)
w = Function(Vs)
w.interpolate(w_func)
dS = Measure('ds', domain=bulk_mesh, subdomain_data=gamma_tags)
w_bulk = Function(Vb)
u_surface = Function(Vs)
parent_dofs = locate_dofs_topological(Vb, dim-1, surface_entities)
u_surface.x.array[:] = u.x.array[parent_dofs]
w_bulk.x.array[parent_dofs] = w.x.array[:]
print("∫u dx = ", assemble_scalar(form(u * dx)))
print("∫u dS = ", assemble_scalar(form(u * dS)))
print("∫u_surface dx = ", assemble_scalar(form(u_surface * dx)))
print("∫w dx = ", assemble_scalar(form(w * dx)))
print("∫w_bulk dS = ", assemble_scalar(form(w_bulk * dS)))
```

Let \Omega = [0,1]^2 and \Gamma = [0,1] \times \{0\} \subset \partial\Omega. I have a function u: \Omega \to \mathbb{R} and a function w: \Gamma \to \mathbb{R}. I define u_\text{surface} = \left . u \right |_\Gamma, and w_\text{bulk} = \begin{cases} u & \text{on } \Gamma\text{,}\\ 0 & \text{elsewhere.} \end{cases}. I would expect that \int u \, \mathrm{d}x > \int u \, \mathrm{d}S = \int u_\text{surface} \, \mathrm{d}x and \int w \, \mathrm{d}x = \int w_\text{bulk} \, \mathrm{d}S, but neither is the case. From what I can tell, only the integrals involving the custom measure `dS`

yield incorrect values. What causes this?