I’m having some issues with how to use internal boundaries, which have been created in gmsh. I’m using legacy fenics as the code I’m working was previously written in fenics and I’m not keen to migrate it all to fenicsx, as I’m a little sure about the syntax modifications needed, however if there’s not way to do this in legacy fenics then I’m happy to move my code over.
Essentially I first create a mesh in gmsh, my MWE here is a simple rectangle with a line down the centre. I put this internal line into physical group 2 and all the external curves into physical group 1.
I then convert this to xdmf and import into fenics. In my real problem I want to solve a problem on this mesh with robin boundary condtions on the internal curves. However, in this toy problem I’ve just attempted to apply a more simple neumann one.
This is where the issue arises when I specfy ds=Measure(“ds”,domain=mesh,subdomain_data=mf,subdomain_id=1), and use this surface element in my weak form, the boundary condtion behaves as expeced. However when I use ds=Measure(“ds”,domain=mesh,subdomain_data=mf,subdomain_id=2), nothing appears to happen. I’ve also tried removing the subdomain_id keyword arguement and using ds(1) and ds(2) instead, with the same results.
I am able to specfiy the Dirichlet condtions for both the internal and external boundaries with no issues, so I suspect it is a problem with the way I’m implementing the neumann ones, rather than the way I’m making the mesh.
From some looking online I wonder if I maybe need to use DG elements, but when I try this the soltuion doesn’t seem to evolve at all. Any help is appreciated, I’m happy to try and read tutorials and work out hwo to do this myself, but I don’t really know where t start.
My MWE:
from fenics import *
import gmsh
import numpy as np
import meshio
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
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
p=0.05
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,2],3)
gmsh.model.mesh.generate(2)
gmsh.write("./mesh_files/square.msh")
# gmsh.fltk.run()
gmsh.finalize()
msh = meshio.read("./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(".mesh_files/square_mf.xdmf", line_mesh)
meshio.write(".mesh_files/square_mesh.xdmf", triangle_mesh)
set_log_level(30)
T=10
dt=10**-2
num_steps=int(T/dt)
mesh, mf = load_gmsh(".mesh_files/square")
V = FunctionSpace(mesh, "Lagrange", 1)
u0 = Expression('exp(-50*pow(x[0]-0.4, 2) - 50*pow(x[1]-0.5, 2))',
degree=2)
# bc=DirichletBC(V,0.0,mf,2) # works for both inernal and external boundaries depending on whehter last arguement is 1 or 2
u=Function(V)
v=TestFunction(V)
u_prev=Function(V)
u=interpolate(u0,V)
u_prev=interpolate(u0,V)
out_file = File("./diffusion/output.pvd")
n=FacetNormal(mesh)
ds_int=Measure("ds",domain=mesh,subdomain_data=mf,subdomain_id=2)
F = (u-u_prev)*v*dx + dt*dot(grad(u), grad(v))*dx - dt*(0.01)*v*ds
for step in range(num_steps):
t=dt*step
solve(F==0,u)
u_prev.assign(u)
out_file << (u,t)
print(t)