Hi all,
I have encountered a problem when solving a pde due to marking of the boundaries which seems like a bug, you can judge by yourself: the pde is similar to http://fenicsproject.org/documentation/tutorial/materials.html
but with Dirichlet and Neumann boundary conditions.
The pde is of the type:
$$\nabla\cdot(\sigma(x,y)\nabla u(x,y))=0$$
with $\sigma(x,y) = 1$ in $\Omega_0$ and $\sigma(x,y)$ = 10 in $\Omega_1$, and the union of $\Omega_0$ and $\Omega_1$ is the unit square. $\Omega_1$ is defined as a ball of center (0.3,0.5) and radius 0.2 .
The boundary is split in 5 parts: top, right, bottom, top left and bottom left. I mark the boundaries in the following way:
top is marked as 3
right is marked as 4
bottom is marked as 5
top left is marked as 6
bottom left is marked as 7
Then I put Dirichlet conditions on the boundary "top left" and nonhomogeneous Neumann conditions on all other boundaries.
If I use the code corresponding to what I just described, then it works fine. If on the other hand I replace the mark 7 for bottom left by the mark 2 then I get a solution which is equal to zero everywhere! So it seems I cannot mark a boundary with the number 2.
If you want to try it yourself you can find the code below. The one I am posting is working well. Then replace 7 by 2 in
leftbot.mark( boundaries, 7)
and in
g*v*ds(7)
at the end of the code and you should get a solution equal to zero everywhere which is not correct (equal to zero).
from dolfin import *
import numpy as np
import pdb
#----------
# Some data
#----------
N = 31
sigma0 = 1
sigma1 = 10
#----------------------
# Define mesh for Omega
#----------------------
mesh = RectangleMesh(0, 0, 1, 1, N, N)
V = FunctionSpace(mesh, "Lagrange", 1)
phi = Function(V)
#------------------------
# Omega = Omega1 + Omega0
#------------------------
psi = Expression('(x[0]-0.3)*(x[0]-0.3) + (x[1]-0.5)*(x[1]-0.5) - 0.2*0.2')
phi = project(psi,V)
class Omega1(SubDomain):
def inside(self, x, on_boundary):
return True if x[0]>=0 and x[0]<=1 and x[1]>=0 and x[1]<=1 and phi(x) < 0 else False
class Top(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and abs(x[1]-1.0) < tol
class Right(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and abs(x[0]-1.0)< tol
class Bottom(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and abs(x[1]) < tol
class LeftTop(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and x[1]> 0.5 and abs(x[0]) < tol
class LeftBot(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and x[1]<= 0.5 and abs(x[0]) < tol
#--------------------------------
# Initialize sub-domain instances
#--------------------------------
lefttop = LeftTop()
leftbot = LeftBot()
top = Top()
right = Right()
bottom = Bottom()
omega1 = Omega1()
#----------------------------------------------
# Initialize mesh function for boundary domains
#----------------------------------------------
boundaries = FacetFunction("uint", mesh)
domains = CellFunction("uint", mesh)
domains.set_all(0)
omega1.mark(domains, 1)
top.mark( boundaries, 3)
right.mark( boundaries, 4)
bottom.mark( boundaries, 5)
lefttop.mark( boundaries, 6)
leftbot.mark( boundaries, 7)
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
#--------------------
# Boundary conditions
#--------------------
f0 = Expression("0")
f1 = Expression("0")
g = Expression("1")
bcs = [DirichletBC(V, 0.0, boundaries, 6) ]
#----------
# Solve pde
#----------
y = TrialFunction(V)
v = TestFunction(V)
a = sigma0*inner(grad(y), grad(v))*dx(0) + sigma1*inner(grad(y), grad(v))*dx(1)
L = f0*v*dx(0) + f1*v*dx(1) + g*v*ds(3) + g*v*ds(4) + g*v*ds(5) + g*v*ds(7)
y = Function(V)
solve(a == L, y, bcs)
#--------------
# Plot solution
#--------------
plot(y)
interactive()
