Here is a model problem for what I want to do. Consider a unit square \Omega with a smaller square \Omega_i fully contained in \Omega. Let \Omega_e = \Omega \setminus \Omega_i be the region inside the unit square but outside the smaller square. We want to solve Poisson’s equation in \Omega_e and \Omega_i:

I want the normal derivative to be smooth across \partial \Omega_i:

I also want to specify the discontinuity in \Theta across the boundary:

for a given function V_m defined on \partial \Omega_i.

I believe the weak form of this problem (not including that last condition on the discontinuity across \partial \Omega_i) is:

So far, I think I can implement this in fenics with the rudimentary knowledge I’ve pieced together from tutorials (in this example code I impose Dirichlet BC on the left and bottom of the outer square, just to make the solution nontrivial):

```
from fenics import *
from mshr import *
# Define constants for marking subdomains/boundaries
OMEGA_E = 1
OMEGA_I = 2
INTERFACE = 3
LEFTEDGE = 4
BOTTOMEDGE = 5
# Create domain and mark subdomains
Omega = Rectangle(Point(0, 0), Point(1, 1))
Omega_i = Rectangle(Point(0.2, 0.2), Point(0.4, 0.4))
Omega.set_subdomain(OMEGA_E, Omega - Omega_i)
Omega.set_subdomain(OMEGA_I, Omega_i)
mesh = generate_mesh(Omega, 256)
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())
# Mark boundaries
boundaries = MeshFunction("size_t", mesh, 1)
for f in facets(mesh):
p0 = Vertex(mesh, f.entities(0)[0])
p1 = Vertex(mesh, f.entities(0)[1])
x0, y0 = p0.x(0), p0.x(1)
x1, y1 = p1.x(0), p1.x(1)
on_vert_edge = lambda x, y: (near(x, .2) or near(x, .4)) and y > .2 and y < .4
on_horiz_edge = lambda x, y: (near(y, .2) or near(y, .4)) and x > .2 and x < .4
on_edge = lambda x, y: on_vert_edge(x, y) or on_horiz_edge(x, y)
if on_edge(x0, y0) and on_edge(x1, y1):
boundaries[f] = INTERFACE
elif x0 < DOLFIN_EPS and x1 < DOLFIN_EPS:
boundaries[f] = LEFTEDGE
elif y0 < DOLFIN_EPS and y1 < DOLFIN_EPS:
boundaries[f] = BOTTOMEDGE
# TEST
# Define a function of ones on the mesh:
# V = FunctionSpace(mesh, "CG", 1)
# f = Function(V)
# f.vector()[:] = 1.
# dStest = dS(subdomain_data=boundaries)
# perim = assemble(f*dStest(INTERFACE))
# print("Perimeter from integrating 1 on the interface: {} (expected 0.8)".format(perim))
# END TEST
V = FunctionSpace(mesh, "Lagrange", 1)
dx = Measure("dx")(subdomain_data=subdomains)
dS = Measure("dS")(subdomain_data=boundaries)
ds = Measure("ds")(subdomain_data=boundaries)
Theta = Function(V)
v = TestFunction(V)
n = FacetNormal(mesh)
LHS = inner(grad(Theta), grad(v)) * dx(OMEGA_E)
LHS += inner(grad(Theta), grad(v)) * dx(OMEGA_I)
LHS += inner(n('-'), grad(Theta('-'))) * v('-') * dS(INTERFACE)
LHS -= inner(n('+'), grad(Theta('+'))) * v('+') * dS(INTERFACE)
solve(LHS == 0, Theta, [DirichletBC(V, 1.0, boundaries, LEFTEDGE),
DirichletBC(V, 2.0, boundaries, BOTTOMEDGE)])
```

Now I need to implement that condition \Theta_e - \Theta_i = V_m \text{ on } \partial \Omega_i, and I have no idea where to start. One hint is that the paper I got this problem from specifies the weak form of that condition as

But I’m not sure how to apply that condition in fenics. I’m grateful for any advice that might help point me in the right direction!