Inconsistency in custom line integral for DG elements

I would like to integrate discontinuous functions over a line for a given 2D domain (In practice I am working in 3D so I want to compute the surface integral but convenience I present the 2D case). Here one has to mark the internal edges and then to integrate with the appropriate labels. In the following example I generate a DG function where u = 0 for x[0] < 0.5 and u = 1 for x[0] >=0.5 as one can see from the attached picture

Then integrating u(‘-’) over the line x=0.5 the values should be 1.
In parallel there is some issue when I am using the command
mpirun -np N python example.py

The value of the computed integral is incorrect when N=9 or N =11. Here a minimal example is presented

from dolfin import *
parameters["ghost_mode"] = "shared_facet"


class DiscontinuousOverlLine(UserExpression):
    def __init__(self, **kwargs):
        try:
            super().__init__(**kwargs)
        except:
            print('super().init failed')
        
    def eval(self, values, x):        
        values[0] = 1
        if x[0] <0.5:
            values[0] = 0
        return

    def value_shape(self):
        return ()

M = 2**3
mesh = UnitSquareMesh(M,M)

interior_surface = CompiledSubDomain('near(x[0], 0.5, DOLFIN_EPS)')
interior_surface_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
interior_surface_marker.set_all(0)
interior_surface.mark(interior_surface_marker,1)


V = FunctionSpace(mesh, 'DG', 0)
ut = Function(V)
ut.interpolate(DiscontinuousOverlLine() )

utF = File('ut.pvd')
utF << ut

L =ut('-')*dS(1, subdomain_data=interior_surface_marker)
print('L=', assemble(L))

L =Constant(1)*dS(1, subdomain_data=interior_surface_marker, domain=mesh)
print('L=', assemble(L))

See Integrating over an interior surface - #4 by MiroK

Thank you for the provided link. I have just realized one should define a subdomain to define ‘-’ and ‘+’ restrictions. Many thanks