To make the discussion a bit more concrete, consider the problem
$$ \begin{alignat}{2} - \nabla \cdot (a \nabla u) &= f &&\mbox{in }\Omega \\ u &= 0 &&\mbox{on }\partial\Omega \\ a \nabla u \cdot n &= k[\![u]\!] \quad &&\mbox{on }\Gamma, \end{alignat} $$
where $\Gamma$ is a curve dividing $\Omega$ into two subdomains $\Omega_1$ and $\Omega_2$, $n$ is the unit normal vector on $\Gamma$ pointing into $\Omega_2$, $[\![u]\!]$ is the jump in $u$ going from $\Omega_1$ to $\Omega_2$, and $k\geq 0$.
The variational formulation of the problem in the space $V = H^1_{0,\partial\Omega}(\Omega_1)+H^1_{0,\partial\Omega}(\Omega_2)$ reads
$$ \int_\Omega a \nabla u \cdot v\, dx + \int_\Gamma k [\![u]\!] [\![v]\!]\, ds = \int_\Omega f v \, dx$$
for all $v\in V$. The bilinear form is bounded, symmetric and positive definite for reasonable $a$ and $k$, so the problem is well-posed in $V$.
The ideal way to obtain a numerical solution to the problem would be to define Lagrange finite element spaces $V_{h,1}\subset H^1_{0,\partial\Omega_1}$ and $V_{h,1}\subset H^1_{0,\partial\Omega_2}$, and write down the variation problem on $V_h = V_{h,1} + V_{h,2} \subset V$ and solve. Unfortunately, this is not yet possible in FEniCS. But there are several ways to get around this limitation, at cost of having to solve a larger linear system.
One way is to adopt a discontinuous Galerkin formulation, as in the DG-Poisson demo, replacing the standard DG terms on the facets that are on $\Gamma$ with the the term $\int k[\![u]\!][\![v]\!]\, ds$.
A closely related method is to solve a mixed form of the Poisson problem, using discontinuous elements for the scalar field, as in the mixed-Poisson demo.
A third way to solve the problem is to let $V_{h,1}$ and $V_{h,2}$ be defined everywhere and implement $V_h$ as a mixed space $V_{h,1}\times V_{h,2}$. This leads to an excess number of dofs that have to be eliminated after assembling the bilinear form. See code example below.
from dolfin import *
from matplotlib import pyplot
parameters["plotting_backend"] = "matplotlib"
N = 32
mesh = UnitSquareMesh(N, N)
x = SpatialCoordinate(mesh)
# indicator functions for the two subdomains
I0 = conditional(x[0] < 0.5, 1, 0)
I1 = conditional(x[0] > 0.5, 1, 0)
# material parameters and source term
c0, c1 = Constant(1.0), Constant(1.0)
c = I0 * c0 + I1 * c1
k = Constant(8.0)
# define exact solution and source term
U0 = 4*x[0]**2
U1 = (1 + 4 * (c0/k)) * 4 * x[0] * (1 - x[0]) \
+ 4 * (c0/c1) * (2*x[0] - 1) * (1 - x[0])
U = (I0 * U0 + I1 * U1) * sin(pi * x[1])
f = - c * div(grad(U))
# define subdomains and cell measure
cell_domains = CellFunction("size_t", mesh)
right = CompiledSubDomain("x[0]>0.5-DOLFIN_EPS")
right.mark(cell_domains, 1)
dx = Measure("dx", domain = mesh, subdomain_data = cell_domains)
# define interface and facet measure
facet_domains = FacetFunction("size_t", mesh)
interface = CompiledSubDomain("fabs(x[0]-0.5)<DOLFIN_EPS")
interface.mark(facet_domains, 1)
dS = Measure("dS", domain = mesh, subdomain_data = facet_domains)
# define variational forms
W = VectorFunctionSpace(mesh, "CG", 1)
u, v = TrialFunction(W), TestFunction(W)
# dirichlet boundary on the exterior
bc = DirichletBC(W, (0., 0.), "on_boundary")
# cell integrals
a = c0 * inner(grad(u[0]), grad(v[0])) * dx(0) \
+ c1 * inner(grad(u[1]), grad(v[1])) * dx(1)
L = f * v[0] * dx(0) + f * v[1] * dx(1)
# interface integral
jmp_u = avg(u[0]) - avg(u[1])
jmp_v = avg(v[0]) - avg(v[1])
a += k * jmp_u * jmp_v * dS(1)
# assemble and eliminate excess dofs
A, b = assemble_system(a, L, bc)
A.ident_zeros()
# solve and compute L^2-error
w = Function(W)
solve(A, w.vector(), b)
L2_error = assemble((U-w[0])**2 * dx(0) + (U-w[1])**2 * dx(1))**.5
print "||u - uh; L^2|| = {0:1.3e}".format(L2_error)
# project to DG0 and plot
V = FunctionSpace(mesh, "DG", 0)
u, v = TrialFunction(V), TestFunction(V)
u_dg = Function(V)
solve(u*v * dx == v*w[0]*dx(0) + v*w[1]*dx(1), u_dg)
plot(u_dg); pyplot.show()
