# 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('-'))

# Variational form
L = f*v1*dx1 + f*v2*dx2

# This structure for penalty terms is probably also wrong
a+= - inner(avg_grad_u, jump_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)

# Variational form
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

# - 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)

# - 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!

2 Likes

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.

Great

`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 .
``````
3 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?

I was able to implement this. Thank you @IngeborgGjerde

However, there is still something strange `-` the continuity at the interface is scaling with the penalty parameter and yet I never get the exact nodal value of solution at the interface. There is always some numerical error no matter what the penalty parameter is (for very small and very big penalty parameters, the solution is of course not stable and that’s natural).

But let’s say if the solution is stable for some value of penalty parameter, should the solution not be exact at the interface or am I wrong?

Hi again,

I think the method works as expected. Comparing the fully coupled and this interface coupled solution method, we find that the difference between the two shrinks as the mesh size decreases and the penalty increases. The decoupled solution does not need to be exact at the interface as long as it converges.

``````alpha=1:    N     error
--------------
8    1.4e-03
16    8.0e-04
32    4.8e-04
64    2.8e-04
128    1.5e-04

alpha=10:    N     error
--------------
8    3.6e-04
16    1.5e-04
32    7.2e-05
64    3.5e-05
128    1.7e-05

alpha=100:    N     error
--------------
8    2.3e-04
16    5.7e-05
32    1.6e-05
64    5.0e-06
128    2.0e-06
``````
Code
``````from dolfin import *
from multiphenics import *

set_log_level(40)

def solve_decoupled(mesh, f, u_bc, alpha):

# 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)

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(alpha/mesh.hmin())

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

# Variational formulation
# - standard terms

# - 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)

# - 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
sol = Function(V)
sol.vector()[:] = U[0].vector().get_local()[:] + U[1].vector().get_local()[:]
File("U.pvd") << sol

return sol

def solve_coupled(mesh, f, u_bc):

V = FunctionSpace(mesh, 'CG', 1)

u = TrialFunction(V)
v = TestFunction(V)

L = f*v*dx

bc = DirichletBC(V, u_bc, 'on_boundary')
sol = Function(V)
solve(a==L, sol, bc)
return sol

# Experiment: compare decoupled and coupled solution for different alpha,
# make sure that error decreases in alpha and N
for alpha in [1, 10, 100]:

f = Expression('sin(2*pi*x[0])*sin(2*pi*x[1])', degree=2)
u_bc = Constant(0)

print(f'\nalpha={alpha}:    N     error')
print('           --------------')
for N in [8, 16, 32, 64, 128]:
mesh = UnitSquareMesh(N, N)
U_sol = solve_decoupled(mesh, f, u_bc, alpha)
u_sol = solve_coupled(mesh, f, u_bc)

# compare
error = errornorm(U_sol, u_sol, 'L2')
print(f'    {N:10.0f} {error:10.1e}')

``````

Hope that helps

Hi again ,

Thanks a lot for your time. With your help, I have almost managed to do what I wanted to do!

In your previous code, I wanted to mention a few things though:

1 -

This is probably incorrect because values of nodal DoFs at the interface for left and right subdomains should be the same (we make sure this happens through dg-type coupling). Moreover, the DoFs mapping for both functions `U[0]` and `U[1]` have the same position/index for the interface nodes, so adding them to get the full domain solution could be misleading, since it will double the solution at the interface nodes. Consider something like this:

`[0. 1.25 1.625 0. 1. 1.5 3. 3.25 0. ]`
`[2. 1.25 1.625 2.5 0. 0. 0. 3.25 4. ]`

These are the `U[0]` and `U[1]` solution vectors for the toy problem I solved, where you see `1.25,1.625,3.25` appear in both the solution vectors and these are indeed the interface DoFs. In my opinion, the full domain solution should be stitched together using these two subdomain solutions (I have not yet seen any function for this) rather than being a sum of these two subdomain solutions.

2 -

I am not sure which type of dg method that is and how you arrived to these terms.
The terms I get upon simplifying the interior penalty type dg coupling are:

``````a[0][0] += -0.5*dot(grad(v1)('-'), u1('-')*n('-'))*dS(3) -0.5*dot(grad(u1)('-'), v1('-')*n('-'))*dS(3) \
+ (alpha/h_avg)*dot(u1('-')*n('-'),v1('-')*n('-'))*dS(3)