Hi, consider using restriction using `+`

and `-`

. Here `foo('+')`

restricts `foo`

to the positive cell and `FacetNormal('+')`

is the outer normal of a positive cell. Following the discussion here the `+`

, `-`

side can be changed by providing a form which has a `dx`

measure term which has specified `subdomain_data`

by some cell function `g`

. Then a positive cell for a facet is the one for which value of g is higher. This allows for the following

```
from dolfin import *
from mshr import *
omega1 = Rectangle(Point(0, 0), Point(0.5, 2))
omega2 = Rectangle(Point(0.5, 0), Point(1, 2))
domain = omega1 + omega2
domain.set_subdomain(1, omega2)
mesh = generate_mesh(domain, 32)
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, 'CG', 1)
# Make up something with discontinuous gradient
f = interpolate(Expression('std::max(0., x[0]-0.5)', degree=1), V)
dS = Measure('dS', domain=mesh, subdomain_data=interior_surface_marker, subdomain_id=1)
n = FacetNormal(mesh)
# +, - sides are by default assigned arbitrarily so assembling below
# will not yield meaningful result
L0 = dot(grad(f)('+'), n('+'))*dS
print(assemble(L0))
L1 = dot(grad(f)('-'), n('-'))*dS
print(assemble(L1))
mf = MeshFunction('size_t', mesh, 2, mesh.domains())
# Add dummy terms that will allow for controlling restrictions
# This is (0, 0) * (1, 0) * 2
L2 = dot(grad(f)('-'), n('-'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L2))
# This is (1, 0) . (-1, 0) * 2
L3 = dot(grad(f)('+'), n('+'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L3))
```