Boundary integral values do not match

General
1

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?