No facet matching domain for velocity boundary in cylinder problem

Question

Hello every body
I have just started to learn FEniCS and I am an amateur. I'm trying to implement the Flow past a cylinder (2D) problem with the Stokes equations. I encountered this error for tentative velocity and computing velocity update:
warning: found no facet matching domain for boundary
and the velocity plot disappear after a while. could you please guide me how to solve it?

here is my code:

from dolfin import*

Set parameter values

dt = 0.01
T = 1
nu = 0.001
tol= 1e-10

mesh domain

mesh = Mesh("cylinder.xml")

## boundary condition

class Left(SubDomain):
def inside(self, x, on_boundary):
  return near(x[0], 0.0,tol)


class Right(SubDomain):
def inside(self, x, on_boundary):
    return near(x[0], 2.4,tol)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
    return near(x[1], 0.0,tol)


class Top(SubDomain):
def inside(self, x, on_boundary):
    return near(x[1], 0.4,tol)

class Cylinder(SubDomain):
def inside(self,x,on_boundary):
    r = sqrt((x[0] - 0.2)**2 + (x[1] -0.2)**2)
    return on_boundary and near(r, 0.05,tol)

left = Left()
top = Top()
right = Right()
bottom = Bottom()
cylindr = Cylinder()

Define function spaces (P2-P1)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)

Define trial and test functions

u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

Define time-dependent pressure and velocity boundary condition

v_in = Expression(('(6*(x[1]*(0.41-x[1]))*sin(pi*t/8.0))/0.0283','0'), t=0.0)

Define boundary conditions

inflowv= DirichletBC(V,v_in, left, method='pointwise')
noslip1=DirichletBC(V,(0,0), top, method='pointwise')
noslip2=DirichletBC(V,(0,0), bottom, method='pointwise')
noslip3=DirichletBC(V,(0,0), right, method='pointwise')
noslip4=DirichletBC(V,(0,0), cylindr, method='pointwise')
outflowp = DirichletBC(Q, 0, right, method='pointwise')

bcu = [inflowv,noslip1,noslip2,noslip3,noslip4]
bcp = [outflowp]

Create functions

u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

Define coefficients

k = Constant(dt)
f = Constant((0, 0))

Tentative velocity step

F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
 nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx

a1 = lhs(F1)
L1 = rhs(F1)

Pressure update

a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx

Velocity update

a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

Assemble matrices

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

Create files for storing solution

ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

Time-stepping

t = dt
while t < T + DOLFIN_EPS:

# Update pressure boundary condition

v_in.t = t

# Compute tentative velocity step
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")
end()

# Pressure correction
begin("Computing pressure correction")
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2, "gmres", "default")
end()

# Velocity correction
begin("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")
end()

# Plot solution
plot(p1, title="Pressure", rescale=True)
plot(u1, title="Velocity", rescale=True)

# Save to file
ufile << u1
pfile << p1

# Move to next time step
u0.assign(u1)
t += dt
print "t =", t

Hold plot

interactive()

mwelland · Answer 1 · Jul 17, 2014

Hi,
First off, next time you post a question, please use the "Code sample" button (see other questions for an example) so that the question is easier to read. Also, please post only a minimum working example again for the same reason. (Usually in writing the MWE, you solve your own problem). The way it is here ^ is pretty confusing, and a few lines, like when you define a mesh and them immediately load a different mesh, seem strange.

I get that error when no facets match a boundary I define. In your case I suspect it is with this block (since the others look simple enough):

class Cylinder(SubDomain):
def inside(self,x,on_boundary):
r = sqrt((x[0] - 0.2)2 + (x[1] -0.2)2)
return on_boundary and near(r,0.04, tol)

Without seeing your mesh file, there isn't much I can help you with, but to solve it I would suggest you try systematically eliminating parts of your code and seeing where the error starts exactly. You could also try increasing 'tol' in the above snipit. If the cylinder in your mesh is actually comprised of straight lines (a polygon, which I suspect it is), such a small tolerance is going to be nigh impossible to satisfy for any facet)

Thanks for replying. I am sorry for the mistakes and bad writting. I changed the code above as you suggested. In fact, I have commented the first lines but after copy, they became parts of codes. I should eliminate them first.
About the code, I have warning: found no facet matching domain for boundary. it can calculate pressure correctly but I really don't understand why it cant compute tentative velocity and velocity update. I have tried "pointwise" boundary condition but the velocity disappears again. could you please guide me how to debug it?
Best Regards
Saeedeh — saeedeh, Jul 18, 2014