Neumann boundary condition on complex domain

Question

Hi FEniCS community,

I am still quite new to FEniCS, but after my first steps I am now trying to solve the diffusion equation with Neumann boundary conditions (no-flux on the boundaries) on a domain with a hole in it. I created this domain using mesh2d in matlab and then imported it into FEniCS which seems to work out properly.

I already managed to apply Neumann BCs on the boundaries of the domain, but I cannot figure out how to specify that the circle is also a boundary of the domain.
I found some similar code for defining the boundary as this:

 class Boundary(SubDomain):
 def inside(x, on_boundary):        
     return on_boundary or on_circle

where i have on_circle return if x lies inside of the circular hole.
When I declare ds as follows I get an error

boundary = Boundary()
ds = Measure("ds")[boundary]

It seems like I am not using the function 'Measure' correctly.

So, how can I tell 'ds' that it should also take the boundary of the hole?

I also attached my code below, but I am not sure if I can also add the .xml file for the mesh.
Also note that for no-flux BCs I would pick g = 0. When picking g=0.5 it becomes clear that there is no flux around the hole.

from dolfin import *
import time

center = [1.5,2]
radius = 0.8

def on_circle(x, on_boundary):    
    oncircle = ((x[0]-center[0])*(x[0]-center[0]) + (x[1]-center[1])*(x[1]-center[1])) <= radius
    return (on_boundary or oncircle) and x[0]>DOLFIN_EPS

class Boundary(SubDomain):
    def inside(x, on_boundary):        
        return on_boundary or on_circle
        #return on_circle

mesh = Mesh('test.xml')        

boundary = Boundary()
#ds = Measure("ds")[boundary]

V = FunctionSpace(mesh, "Lagrange", 1)

# Time variables
dt = Constant(0.3); t = float(dt); T = 5
g_expr = "0.5"
g = Expression(g_expr)

# initial values
init_expr = "x[0]<DOLFIN_EPS"
init = Expression(init_expr)
u0 = interpolate(init, V);

u1 = Function(V)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = u*v*dx(0) + dt*inner(grad(u), grad(v))*dx(0)
L = u0*v*dx(0) + dt*f*v*dx(0) + dt*g*v*ds(0) #for neumann condition with g on boundary

u_plot = plot(u0, interactive=True,axes=True,rescale=False,scalarbar=True)
while t <= T:    
    solve(a == L, u1)
    u0.assign(u1)
    t += float(dt)    
    time.sleep(0.2)
    u_plot.update(u1)
plot(mesh,interactive=True,rescale=False,scalarbar=True)

And here is also a picture of the mesh:

I would really appreciate any help!
Thank you!

Answer 1

Hi Johannes,

class Boundary(SubDomain):
    def inside(x, on_boundary):        
        return on_boundary

is true for all points/edges on the outer boundary (i.e. the rectangle) and the inner boundary (i.e. the circle).

If you want to impose the same bc on both the outer and inner boundary then you don't have to explicitly define on_circle...

Best,
Umbe

Comments

Thank you very much for your quick reply, Umbe!
Does this mean that I don't even need to add

boundary = Boundary()
ds = Measure("ds")[boundary]

?

---

Just tried it and now it seems to work!
Thank you so much, Umbe!
This was really causing me headaches

---

Hi Johannes,

Yes, if you are integrating over all boundaries the symbol ds is already defined.
I also believe that the definition of ds is exactly the one you wrote above above.

Best,

Umbe