Integral over single facet in meshview gives wrong result

[Jørgen Dokken]

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(2, 2)
marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
for f in facets(mesh):
    marker[f] = np.isclose(f.midpoint().y(), 1)
    if marker[f] == 1:
        break
print(marker.array())
submesh = create_meshview(marker, 1)

expr = Expression("x[0] + x[1]", degree=1)
V1 = FunctionSpace(mesh, "Lagrange", 1)
dx_sub = Measure("dx", domain=submesh)
u = interpolate(expr, V1)

# Now compute integral over boundary in three different ways
# Method 1: integrate over submesh
print(f"Value of submesh integral is {assemble(u*dx_sub)}")
# Method 3: integrate u over boundary (ds)
print(
    f"Value of boundary integral is {assemble(u*ds(domain=mesh, subdomain_data=marker, subdomain_id=1))}")

V_sub = FunctionSpace(submesh, "Lagrange", 1)
u_sub = interpolate(expr, V_sub)
print(f"SubMesh interpolated function {assemble(u_sub*dx)}")

returning

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]
Value of submesh integral is 0.375
Value of boundary integral is 0.625
SubMesh interpolated function 0.625

‌