"""
Flow between axially moving concentric cylinders
"""
from dolfin import *

r0 = 0.2
r1 = 1
U0 = 0.0
U1 = 1.0
mesh = IntervalMesh(500, r0, r1)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
x = Expression("x[0]")
dpdz = Constant(0.0)
mu = Constant(0.01)

F = inner(grad(v), grad(u))*x*dx + 1./mu*dpdz*x*v*dx

def inner_bnd(x, on_bnd):
    return on_bnd and near(x[0], r0)

def outer_bnd(x, on_bnd):
    return on_bnd and near(x[0], 1)

bci = DirichletBC(V, Constant(U0), inner_bnd)
bco = DirichletBC(V, Constant(U1), outer_bnd)

u_ = Function(V)
solve(lhs(F) == rhs(F), u_, bcs=[bci, bco])
plot(u_, title="Computed velocity")

u_exact = (U1-U0)*ln(x)/ln(r1/r0) + U0 - (U1-U0)*ln(r0)/ln(r1/r0)
u_exact = project(u_exact, V)

plot(u_exact, mesh=mesh, title="Exact velocity")
plot(u_-u_exact, mesh=mesh, title="Error velocity")
print "Computed errornorm = ", errornorm(u_, u_exact, degree_rise=0)