Hi, let me mention an alternative to previous answers which avoids the multiscale formulation. As such you can readily implement it in standard FEniCS. In addition, compared to the suggested mortaring approach you will get a better representation of flux/current. Moreover, standard preconditioners will work (for the multiscale preconditioners see here for mortaring and here mixed mortaring approach). Unless you have a use for the Lagrange multiplier further in your solver, as is done eg. here, the only downside is that a larger problem is to be solved.

```
from dolfin import *
# Let Gamma be the interface between Omega_1, Omega_2 and consider
# the problem
#
# -div(kappa_i*grad(u_i)) = f_i in Omega_i
#
# such that on Gamma
#
# kappa_0*dot(grad(u_0), n_0) + kappa_1*dot(grad(u_1), n_1) = 0, (1)
# u_0 - u_1 = J. (2)
#
# Here n_i is the outer normal of Omega_i. Condition (1) is the flux/current
# continuity while (2) forces jumps in the potentials/voltages.
#
# Let us introduce fluxes sigma_i = kappa_i*grad(u_i) as explicit unknowns.
# By (1) the fluxes are continuous accross the interfaces. With the new
# unknowns the system reads
#
# kappa^{-1}_i sigma_i - grad(u_i) = 0 in Omega_i (3)
# -div(sigma_i) = f_i in Omega_i
#
# while the interface conditions are
#
# dot(sigma_0, n_0) + dot(sigma_1, n_1) = 0 (4)
# u_0 - u_1 = J.
#
# We can enforce (3) by construction if sigma, such that sigma|_Omega_i = sigma_i
# is sought after in H-div conforming (vector valued) finite element space.
# Let then tau be a test function from such a space. Then integrating (3)
# on each subdomain leads to
#
# inner(kappa^{-1}_i*sigma_i, tau_i)*dx(i) + inner(div(tau_i), u_i)*dx(i) +
#
# an interface term -inner(u_i, dot(tau_i, n_i))*dS(1). However, (4) for
# tau by construction the equations can be added yielding the interface
# term
#
# -inner((u_0 - u1), dot(tau, n_0))*dS(1).
#
# Plugging u_0 - u_1 = J, the interface term contributes to the right
# hand side. Let's implement the idea:
mesh = UnitSquareMesh(32, 32)
left = CompiledSubDomain('x[0] < 0.5+DOLFIN_EPS')
thin = CompiledSubDomain('near(x[0], 0.5)')
# Conductivities for subdomains
kappa_1, kappa_2 = Constant(1), Constant(2)
# Forcing on subdomains
f_1, f_2 = Constant(1), Constant(2)
# Interface forcing
J = Constant(3)
tdim = mesh.topology().dim()
cell_f = MeshFunction('size_t', mesh, tdim, 2)
left.mark(cell_f, 1)
# Interface
facet_f = MeshFunction('size_t', mesh, tdim-1, 0)
thin.mark(facet_f, 1)
cell = mesh.ufl_cell()
# Flux and potential
elm = MixedElement([FiniteElement('Raviart-Thomas', cell, 1),
FiniteElement('Discontinuous Lagrange', cell, 0)])
W = FunctionSpace(mesh, elm)
sigma, u = TrialFunctions(W)
tau, v = TestFunctions(W)
dx = Measure('dx', domain=mesh, subdomain_data=cell_f)
dS = Measure('dS', domain=mesh, subdomain_data=facet_f)
n = FacetNormal(mesh)
a = inner((1./kappa_1)*sigma, tau)*dx(1) + inner((1./kappa_2)*sigma, tau)*dx(2)
a += inner(div(tau), u)*dx() + inner(div(sigma), v)*dx()
L = inner(-f_1, v)*dx(1) + inner(-f_2, v)*dx(2)
L += inner(avg(J), dot(avg(tau), n('+')))*dS(1)
wh = Function(W)
solve(a == L, wh)
sigma, u = wh.split(deepcopy=True)
File('u.pvd') << u
File('sigma.pvd') << sigma
```