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!
