Dear all,

I want to solve a Laplace’s equation on a domain which is subdivided into two subdomains like this:

The idea is to build FE function spaces on these subdomains and use `MixedFunctionSpace`

to solve the full domain problem. For enforcing the solution continuity *weakly* at the interface, the idea is to use penalty methods, e.g., interior penalty method. I tried to write the code for this setting, but I cannot figure out correct integration measures and penalty terms for the interface. Here is the MWE:

```
from dolfin import *
# Create mesh
mesh = UnitSquareMesh(8, 8)
# Define domain boundaries
class Left(SubDomain):
def inside(self,x,on_boundary):
return x[0] < 0.0 + DOLFIN_EPS
class Right(SubDomain):
def inside(self,x,on_boundary):
return x[0] > 1.0 - DOLFIN_EPS
class Interface(SubDomain):
def inside(self,x,on_boundary):
return near(x[0], 0.5)
class Top(SubDomain):
def inside(self,x,on_boundary):
return x[1] > 1.0 - DOLFIN_EPS
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return x[1] < 0.0 + DOLFIN_EPS
# Define subdomains
domains = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
domains[c] = c.midpoint().x() > 0.5
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)
# Create submeshes
mesh_L = MeshView.create(domains, 0)
mesh_R = MeshView.create(domains, 1)
mesh_I = MeshView.create(facets, 1)
# Mark the boundaries and define integration measures(this is probably wrong)
mf_1 = MeshFunction("size_t", mesh_L, 1)
mf_1.set_all(0)
Left().mark(mf_1, 1)
Interface().mark(mf_1, 2)
mf_2 = MeshFunction("size_t", mesh_R, 1)
mf_2.set_all(0)
Right().mark(mf_2, 1)
Interface().mark(mf_2, 2)
mf = MeshFunction("size_t", mesh, 1)
mf.set_all(0)
Left().mark(mf, 1)
Right().mark(mf, 2)
Top().mark(mf, 3)
Bottom().mark(mf, 4)
Interface().mark(mf, 5)
dx1 = Measure('dx', domain=mesh_L)
dx2 = Measure('dx', domain=mesh_R)
ds1 = Measure('ds', domain=mesh_L, subdomain_data=mf_1)
ds2 = Measure('ds', domain=mesh_R, subdomain_data=mf_2)
dS = Measure('dS', domain=mesh, subdomain_data=mf, subdomain_id= 5)
# Construct function spaces
V_L = FunctionSpace(mesh_L, 'CG', 1)
V_R = FunctionSpace(mesh_R, 'CG', 1)
V = MixedFunctionSpace(V_L,V_R)
# Define boundary conditions
uD_l = Constant((1.0))
uD_r = Constant((0.0))
bc = [ DirichletBC(V.sub_space(0), uD_l, mf_1, 1) \
,DirichletBC(V.sub_space(1), uD_r, mf_2, 1) ]
# Define trial and test functions
u1,u2 = TrialFunctions(V)
v1,v2 = TestFunctions(V)
# Define model parameters, normals, and penalty parameters
f = Constant(0.0)
n = FacetNormal(mesh)
n1 = FacetNormal(mesh_L)
n2 = FacetNormal(mesh_R)
alpha = Constant(100.0)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-')) / 2.0
# probably defining it the wrong way
jump_u = (u1('+')*n1('+') + u2('-')*n2('-'))
jump_v = (v1('+')*n1('+') + v2('-')*n2('-'))
avg_grad_u = 0.5*(grad(u1)('+') + grad(u2)('-'))
avg_grad_v = 0.5*(grad(v1)('+') + grad(v2)('-'))
# Variational form
a = inner(grad(u1), grad(v1))*dx1 + inner(grad(u2), grad(v2))*dx2
L = f*v1*dx1 + f*v2*dx2
# This structure for penalty terms is probably also wrong
a+= - inner(avg_grad_u, jump_v )*dS \
- inner(jump_u , avg_grad_v)*dS \
+ (alpha / h_avg)*inner(jump_u , jump_v)*dS
# Compute solution
u = Function(V)
solve(a == L, u, bc)
```

This of course does not work. I am sure I am implementing the penalty terms incorrectly and also the integration measure taken for the penalty terms is also wrong. I have read somewhere that for `MixedFunctionSpace`

the measure `dS`

is not supported. But then I am not sure how to code the problem for this setting correctly.

P.S. I am aware that using Lagrange multipliers at the interface could also enforce the solution continuity at the interface, but that approach might not be useful for the bigger problem I am working on.

Any help will be highly appreciated!