Enforcing weak continuity at interface for a mixed dimensional problem

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!

Wrt. to MixedFunctionSpace, You should consider:

which comes from the paper:
https://dl.acm.org/doi/pdf/10.1145/3471138 by @cdaversin

Hi @dokken,

Thanks for your reply. The link you referred me to uses Lagrange multiplier living at the domain interface for enforcing the continuity, however I don’t want to use this approach. I am rather interested in using penalty terms for enforcing weak continuity at the interface. Is it doable using MixedFunctionSpace?

I think it should be possible, I guess you just need to take some care when adding the weak enforcement on that boundary (I do not use the MixedFunctionSpace myself), maybe @IngeborgGjerde knows?

Yes, this is something I cannot figure out. I think defining the correct integration measure at the interface and formulating the correct penalty terms for the jump between u1 and u2 at the interface should solve the problem. But so far I am unable to find anything relevant to this

I’m not sure this is possible with fenics-mixed-dim due to the use of the dS function for interface. As far as I know this is not supported.

Here is an attempt using fenics-mixed-dim, but I could not get it to work.
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)


# 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, n2, n = FacetNormal(mesh_L), FacetNormal(mesh_R), FacetNormal(mesh)
alpha = Constant(100.0)
h1, h2 = CellDiameter(mesh_L), CellDiameter(mesh_R)
h_avg = 0.5*(h1+h2)

# 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

a += inner(grad(u1), v1*n1)*ds1(2) # this term is fine
a += inner(grad(u1)('+'), v2('+')*n2('+'))*dS(5) # this term is not allowed
#a += inner(grad(u2), v1*n1)*dS(5) # this term is not allowed
a += inner(grad(u2), v2*n2)*ds2(2) # this term is fine
#a += (alpha / h_avg)*inner(jump_u , jump_v)*dS(5)

# Compute solution
u = Function(V)
solve(a == L, u, bc)

The problem can be implemented quite nicely using multiphenics:

from dolfin import *
from multiphenics import *

# Create mesh
mesh = UnitSquareMesh(16,16)

# Interface along x=0.5
class Interface(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[0], 0.5) 
    
class LeftBoundary(SubDomain):
    def inside(self,x,on_boundary):
        return x[0]<0.5 and on_boundary
    
class RightBoundary(SubDomain):
    def inside(self,x,on_boundary):
        return x[0]>=0.5 and on_boundary

# This interface gives rise to left and right subdomains
domains = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
    domains[c] = c.midpoint().x() > 0.5
    
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<0.5
    
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]>=0.5

mf = MeshFunction("size_t", mesh, 1)
mf.set_all(0)
Interface().mark(mf, 5)
LeftBoundary().mark(mf,2)
RightBoundary().mark(mf,3)


dx = Measure('dx')(subdomain_data=domains)
ds = Measure('ds')(subdomain_data=mf)
dS = Measure('dS')(subdomain_data=mf)

left = MeshRestriction(mesh, Left())
right = MeshRestriction(mesh, Right())

# Construct function spaces
V = FunctionSpace(mesh, 'CG', 1)
W = BlockFunctionSpace([V, V], restrict=[left, right])


# Define trial and test functions
u1u2 = BlockTrialFunction(W)
(u1, u2) = block_split(u1u2)
v1v2 = BlockTestFunction(W)
(v1, v2) = block_split(v1v2)


# model parameters
alpha = Constant(0.2/mesh.hmin())
f = Expression('sin(2*pi*x[0])*sin(2*pi*x[1])', degree=2)
u_bc = Constant(0)

n = FacetNormal(mesh)
a = [[0,0], [0,0]]

# Variational formulation
# - standard terms
a[0][0] = inner(grad(u1), grad(v1))*dx(0)
a[1][1] = inner(grad(u2), grad(v2))*dx(1)

# - dg terms on interface: 
# [[u]]{{v}}-[[v]]{{u}}+alpha*[[u]][[v]] = -0.5*(u1-u2)*(v1+v2)-0.5*(u1+u2)*(v1-v2)+alpha*(u1*v1+u2*v2)
# = (alpha-1)*u1*v1 +(alpha+1)*u2*v2 - alpha*u2*v1 - alpha*u1*v2 

# the restriction ('+') and ('-') is arbitrary, we do it so fenics does not complain
a[0][0] += +inner(grad(u1)('+'), v1('+')*n('+'))*dS(5) + (alpha-1)*u1('+')*v1('+')*dS(5)
a[1][1] += +inner(grad(u2)('+'), v2('+')*n('+'))*dS(5) + (alpha+1)*u2('+')*v2('+')*dS(5)

a[0][1] += inner(grad(u2)('+'), v1('+')*n('+'))*dS(5) -alpha*u2('+')*v1('+')*dS(5)
a[1][0] += inner(grad(u1)('+'), v2('+')*n('+'))*dS(5) -alpha*u1('+')*v2('+')*dS(5)

# - right hand side
L =  [f*v1*dx(0)                       , f*v2*dx(1)                   ]

bc1 = DirichletBC(W.sub(0), u_bc, mf, 2)
bc2 = DirichletBC(W.sub(1), u_bc, mf, 3)
bcs = BlockDirichletBC([bc1, bc2])

# solve
A = block_assemble(a)
F = block_assemble(L)
bcs.apply(A)
bcs.apply(F)

U = BlockFunction(W)
block_solve(A, U.block_vector(), F)

# add U1 and U2 to get full solution
U_sol = Function(V)
U_sol.vector()[:] = U[0].vector().get_local()[:] + U[1].vector().get_local()[:]
File("U.pvd") << U_sol

You can increase continuity at the interface by increasing alpha.

Hope that helps, let me know if you have any questions!

Hi @IngeborgGjerde ,

Thanks for looking into the problem. This is exactly what I have been trying to do and it seems like multiphenics offers a nice structured way to implement this.

Now one trivial question, which is not so trivial for me since I am not good at installing new libraries in Linux and I cannot see anything related to installation of multiphenics on its repository:

How can I neatly install multiphenics package in Ubuntu without messing up any of the dolfin package files? The last thing I would want is to see an error when calling from dolfin import * only because I have installed multiphenics incorrectly.

Thanks for your help! :slight_smile:

Great :smiley:

multiphenics can be smoothly installed via pip. In the terminal, you write

git clone https://github.com/multiphenics/multiphenics
cd multiphenics
python3 -m pip install .
2 Likes

Note that multiphenics marks itself as closed (archived), pointing to the new version that works with the next-generation version of FEniCS (dolfinx), https://github.com/multiphenics/multiphenicsx

2 Likes

Hello again @IngeborgGjerde,

Could you please rerun the code for the following changes:
f = Constant(0)
bc1 = DirichletBC(W.sub(0), Constant(1), mf, 2)
bc2 = DirichletBC(W.sub(1), Constant(0), mf, 3)
bcs = BlockDirichletBC([bc1, bc2])

The rest remains the same.

For these changes, the response should be a diffusion from left boundary to the right one, but I don’t see the correct response. Also, I noticed that there are probably some mistakes in your penalty terms.

# dg terms on interface:
# [[u]]{grad{v}} - [[v]]{grad{u}} + alpha*[[u]][[v]]

In total this should result in 12 terms, I think. But in your code you have 10.

I tried to put them my way, but I still don’t get the desired solution.
Could it be the solver issue?