from dolfin import *
from numpy import cosh, cos

dpdx = Constant(-0.01)
mu = Constant(0.01)
a = 2
b = 1

mesh = RectangleMesh(-a, -b, a, b, 20, 20)

V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

F = inner(grad(u), grad(v))*dx + 1/mu*dpdx*v*dx

bc = DirichletBC(V, Constant(0), DomainBoundary())

u_ = Function(V)
solve(lhs(F) == rhs(F), u_, bcs=bc)

plot(u_, title="Velocity normal to plane")

factor = 16.*a*a/mu(0)/pi**3*(-dpdx(0))
class u_exact(Expression):
        
    def eval(self, value, x):
        u = 0
        for i in range(1, 100, 2):
            u += (-1)**((i-1)/2)*(1-cosh(i*pi*x[1]/2./a)/ \\
                 cosh(i*pi*b/2/a))*cos(i*pi*x[0]/2./a)/i**3
        value[0] = u * factor        

u_e = project(u_exact(), V, bcs=bc)
plot(u_e, title="Exact velocity normal to plane")

u_error = errornorm(u_, u_e, degree_rise=0)

print "Computed error = ", u_error

interactive()