Integrating over an interior surface

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
print(assemble(L0))