I am trying to compute 2D stokes flow past a confined cylinder using 1st order velocities.
The following code works on my machine (FEniCS 1.6).
from dolfin import *
from mshr import *
L1 = 4.; L2 = 7.; w = 4.
noslip = Constant((0.,0.))
inflow = Expression(("1.5*(1.-pow(x[1]*2./w,2.))","0."),w=w)
r1 = Rectangle(Point(-L1,-w/2.),Point(L2,w/2.))
c1 = Circle(Point(0.,0.),1.,8)
domain = r1-c1
mesh = generate_mesh(domain, 45)
def upper_lower(x, on_boundary):
return x[1] < -w/2.+1e3*DOLFIN_EPS or w/2.-1e3*DOLFIN_EPS < x[1]
def left(x, on_boundary):
return x[0] < -L1+1e3*DOLFIN_EPS
def right(x, on_boundary):
return L2 - 1e3*DOLFIN_EPS < x[0]
def cylbnd(x, on_boundary):
return x[0]*x[0]+x[1]*x[1] < 1.5 and on_boundary
class right2(SubDomain):
def inside(self, x, on_boundary):
return L2 - 1e3*DOLFIN_EPS < x[0]
Q = VectorFunctionSpace(mesh, "CG", 1)
V = FunctionSpace(mesh,'CG',1)
P = FunctionSpace(mesh,'DG',0)
W = MixedFunctionSpace([Q, V, P])
U = Function(W)
#boundaries
bc0 = DirichletBC(W.sub(0), noslip, cylbnd)
bc1 = DirichletBC(W.sub(0), noslip, upper_lower)
bc2 = DirichletBC(W.sub(1), Constant(0.), right)
bc3 = DirichletBC(W.sub(0), inflow, left)
bcs = [bc0, bc1, bc2, bc3]
(u, p, pDG) = TrialFunctions(W)
(u_test, p_test, pDG_test) = TestFunctions(W)
a = (inner(sym(grad(u)), sym(grad(u_test))) - div(u_test)*pDG \
+ p_test*div(u) + (p-pDG)*pDG_test)*dx
L = (inner(noslip, u_test))*dx
solve(a==L, U, bcs)
v, p, pDG = U.split()
plot(v[0], interactive=True)
plot(p, interactive=True)
but it fails, if I set the number of circle segments to 16. Why?
The problem disappears, when I use my own mesh adaptation library to generate the mesh. This has a true circle representation via analytical distance functions.
