# Seperated Domains within problem

Hi all,

A couple of days ago I posted implementing internal boundaries in fenics and got some really helpful feedback. However it has made me realise I now have another issue. For now, I am considering a 2D convection-diffusion problem. However what I would like is to solve the problem in two decoupled domain simaltaneously, with a method that would allow me to add coupling to the BC’s at a later date. Right now I use no flux boundary condtions at the extenral boundaries and this gives me the expected results for the steady state solution of the problem (expontetial function which builds up at the end we convect towards). This is image 1.

I then added an internal boundary which also has the no flux condtion, and what I expected to see was the domain split in half, with the original steady state solution essentially copied onto each half. However, what I see is some weird solution which ssems somewhat unphysical (image 2).

Having tried to look around on these forums, I understand that this is because I’m using continuous elements. I get the impression to get the behaviour I want I need to either use DG elements, or use the mixed dimesnional brach of fenics. I wanted to confirm that this is right and get the opinions of others on what to do next.

To add some context, the real problem I want to solve involves an elastic medium, which is deformed by a chemical undergoing reaction-diffusion on top, and this is in turn advected by the medium. I have a version of this code working already, but want to add internal boundaries which the chemical cannot move across, (but the elastic medium can interact accros the boundaries, which is why I want to do this all as one simulation). So I’m using this problem to try and understand how to do advection-diffusion with internal boundaries. Any advice on how you would proceed would be appreciated.

Here’s my code for the simplified problem:

``````from fenics import *
from ufl import Index

import numpy as np
from os.path import expanduser,join

# import gmsh
# import meshio

############### mesh conversion funcs ################
space_dim=2
mesh = Mesh()
with XDMFFile(path_to_mesh + "_mesh.xdmf") as infile:
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(path_to_mesh + "_mf.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
return mesh, mf

def create_meshio_mesh(mesh, cell_type, prune_z=True):
import meshio
cells = mesh.get_cells_type(cell_type)

cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh

############### Make mesh ################
p=0.01
gmsh.initialize()

# old line 2 here

gmsh.model.occ.synchronize()

gmsh.model.mesh.generate(2)
gmsh.write(expanduser("./mesh_files/square.msh"))
gmsh.fltk.run()
gmsh.finalize()

line_mesh = create_meshio_mesh(msh, "line", prune_z=True)
triangle_mesh = create_meshio_mesh(msh, "triangle", prune_z=True)

meshio.write(join(expanduser("./mesh_files/square_mf.xdmf")), line_mesh)
meshio.write(join(expanduser("./mesh_files/square_mesh.xdmf")), triangle_mesh)

############### Run simulation ################

set_log_level(30)
parameters["ghost_mode"] = "shared_facet"
T=10
dt=10**-2
num_steps=int(T/dt)

P1=FiniteElement('P',mesh.ufl_cell(),1)
func_space = FunctionSpace(mesh,P1)

rho=Function(func_space)
rho_prev=Function(func_space)
test_rho=TestFunction(func_space)

rho_0=interpolate(Constant(1),func_space)
assign(rho_prev,rho_0)
assign(rho,rho_0)

out_file1 = File(expanduser("./con_diff/output_rho.pvd"))

n=FacetNormal(mesh)
dx=Measure('dx',domain=mesh,subdomain_data=mf)
ds_ext=Measure('ds',domain=mesh,subdomain_data=mf,subdomain_id=1)
ds_int=Measure('dS',domain=mesh,subdomain_data=mf,subdomain_id=2)
j=Index()
v=Expression(('1','0'),degree=2,domain=mesh)
D=Constant(1)
dt=Constant(dt)

F =(
(rho-rho_prev)*test_rho*dx + dt*Dx(rho*v[j],j)*test_rho*dx
+ dt*D*Dx(rho,j)*Dx(test_rho,j)*(dx)
- dt*rho*v[j]*n[j]*test_rho*ds_ext - dt('+')*rho('+')*v[j]('+')*n[j]('+')*test_rho('+')*ds_int
)

for step in range(num_steps):
t=step*dt
solve(F==0,rho)
rho_prev.assign(rho)
out_file1 << (rho_prev,t)
print(t)
``````

Hi,

In order to get a function that is discontinuous across the interface you can combine two functions that are P1 on each submesh

From your code I see that you are assembling and solving the problem in the form F==0. Is this because you want to solve a non-linear problem eventually? The reason I ask is that there is a bug causing non-linear problems to fail with fenics-mixed-dim.

We can fix that by switching to fenics_ii, with the caveat that the syntax is a bit more involved.

sketch of fenics_ii code
``````from fenics import *
from ufl import Index
set_log_level(30)
parameters["ghost_mode"] = "shared_facet"
from xii import *

# Make 2d mesh
mesh = UnitSquareMesh(30,30)

marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

class Boundary(SubDomain):
def inside(self, x, on_boundary):
return x[0]<0.5

bound = Boundary()
bound.mark(marker, 1)
meshleft = SubMesh(mesh, marker, 1)
meshright = SubMesh(mesh, marker, 0)

# Run simulation
T=1
dt=0.1
num_steps=int(T/dt)

W = [FunctionSpace(mesh, 'CG', 1) for mesh in (meshleft, meshright)]

rho = ii_Function(W)
rhol, rhor = rho

rho_prev = ii_Function(W)
rhol_prev, rhor_prev = rho_prev

test_rho = list(map(TestFunction, W))

assign(rhol_prev,interpolate(Constant(1),W[0]))
assign(rhor_prev,interpolate(Constant(1),W[1]))

j=Index()
v=Expression(('1','0'),degree=2)
D=Constant(1)
dt=Constant(dt)

F = [0,0]
for i in range(len(F)):

rhoi = rho[i]
rhoi_prev = rho_prev[i]
test_rhoi = test_rho[i]
vi = interpolate(v, VectorFunctionSpace(W[i].mesh(), 'CG', 1))

dxi = Measure('dx', domain=W[i].mesh())
dsi = Measure('ds', domain=W[i].mesh())
n = FacetNormal(W[i].mesh())

# note: I changed some of the syntax here to fit fenics_ii, please double check
F[i] +=(
+ dt*D*rhoi.dx(j)*test_rhoi.dx(j)*dxi
- dt*rhoi*dot(vi,n)*test_rhoi*dsi
)

from xii.nonlin.jacobian import block_jacobian
# Solution
dF = block_jacobian(F, rho)

# Newton
eps = 1.0
tol = 1.0E-10
niter = 0
maxiter = 125

drho = ii_Function(W)

for step in range(num_steps):
t=step*dt
print('timestep: ',step)

# newton solver at this time step
while eps > tol and niter < maxiter:
niter += 1

A, b = list(map(ii_assemble, (dF, F)))
A, b = list(map(ii_convert, (A, b)))

solve(A, drho.vector(), b)

eps = sqrt(sum(x.norm('l2')**2 for x in drho.vectors()))

for i in range(len(W)):
rho[i].vector().axpy(-1, drho[i].vector())

rhol.assign(rho[0])
rhor.assign(rho[1])

File('shared/mixed-dim-comparisons/rhol.pvd') << rhol
File('shared/mixed-dim-comparisons/rhor.pvd') << rhor
``````

Hope that helps, feel free to ask more questions and I will try my best to help

1 Like