Strange error somehow connected to boundary condition "in the interior"

Asked by Peter Franek on 2013-06-04

Can you help me please? I isolated the problem as much as possible. Solving a Stokes equation on the square with a prescribed force, a nonslip bc on the boundary and an additional boundary condition u=0 in some rectangle in the interior. Here is the code. It works if the rectangle is very "thin", such as [0.5, 0.6]\times [0.1, 0.8], but it returns an error when the rectangle is "thick", such as [0.2, 0.7] \times [0.1, 0.8]. In similar situation, the chorin projection demo in Navier Stokes causes no problem. Also the stokes demo online works ok, if I add this rectangle-bc. Any hint? Here is the complete code.

from dolfin import *
from boundaryparts import *

mesh=UnitSquare(10,10)

f = Expression(("sin(pi*x)*cos(pi*x+pi/2)", "5*x*x*sin((1-x))*(1-x)"))

# Function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
mixed = MixedFunctionSpace([V,Q,R]) # (vector, scalar, constant)

# Trial and test functions
A=TrialFunction(mixed)
B=TestFunction(mixed)
u,p,c = split(A)
v,q,d = split(B)

def u_boundary(x, on_boundary):
return on_boundary
# boundary condition for velocity
nonslip = DirichletBC(mixed.sub(0),('0', '0'), u_boundary)

############# CHANGING THIS HAS INFLUENCE ON SUCCESS ##########################
obstacle_region = BoundaryPart2D([0.1,0.1], [0.9,0.5])

obstacle = DirichletBC(mixed.sub(0), ('0', '0'), obstacle_region)
bcu = [nonslip, obstacle]

# Variational problem for stacionary Stokes flow
a = inner(grad(u),grad(v))*dx + p*div(v)*dx + q*div(u)*dx + (p*d+c*q)*dx
L = inner(f, v)*dx

Vel = Function(mixed) # solution for velocity and pressure
solve(a == L, Vel, bcu)
(vel, pres, a)=split(Vel)
plot(vel, title="Stokes flow", rescale=True)
interactive()

Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-06-04