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 ################
def load_gmsh(path_to_mesh):
space_dim=2
mesh = Mesh()
with XDMFFile(path_to_mesh + "_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(path_to_mesh + "_mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
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()
gmsh.model.add("square")
gmsh.model.occ.addPoint(0,0,0,p,1)
gmsh.model.occ.addPoint(1,0,0,p,2)
gmsh.model.occ.addPoint(1,1,0,p,3)
gmsh.model.occ.addPoint(0,1,0,p,4)
gmsh.model.occ.addLine(1,2,1)
gmsh.model.occ.addLine(2,3,2)
gmsh.model.occ.addLine(3,4,3)
gmsh.model.occ.addLine(4,1,4)
gmsh.model.occ.addCurveLoop([1,2,3,4],1)
gmsh.model.occ.addPlaneSurface([1],1)
gmsh.model.occ.addPoint(2,0,0,p,5)
gmsh.model.occ.addPoint(2,1,0,p,6)
# old line 2 here
gmsh.model.occ.addLine(3,6,5)
gmsh.model.occ.addLine(6,5,6)
gmsh.model.occ.addLine(5,2,7)
gmsh.model.occ.addCurveLoop([2,5,6,7],2)
gmsh.model.occ.addPlaneSurface([2],2)
gmsh.model.occ.synchronize()
gmsh.model.add_physical_group(1,[1,3,4,5,6,7],1)
gmsh.model.add_physical_group(1,[2],2)
gmsh.model.add_physical_group(2,[1],3)
gmsh.model.add_physical_group(2,[2],4)
gmsh.model.mesh.generate(2)
gmsh.write(expanduser("./mesh_files/square.msh"))
gmsh.fltk.run()
gmsh.finalize()
msh = meshio.read(expanduser("./mesh_files/square.msh"))
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)
mesh,mf=load_gmsh(expanduser("./mesh_files/square"))
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)
```