solving Mixed Formulation for vector laplacian

I'm trying to implement the example of solving a vector laplacian (in a poisson's equation) on a L shaped domain. There are two parts to the example- one with linear elements for both subspaces and one with RT subspace for u. The equations and figures are from this paper:
www.ams.org/bull/2010-47-02/.../S0273-0979-10-01278-4.pdf
(page 296)
While the solution with both elements from linear space works nice, when i use RT elements, the program indicates 'shape mismatch' error. What may be wrong?

The code is:

"""
Solve
-u'' = f(x); f(x) = (pi/2)^2 cos(pi*x/2)
-1 < x < 1, u(-1) = u(1) = 0
14 subintervals
"""
def u0_boundary(x, on_boundary):
      return on_boundary

from dolfin import *

mesh = Mesh("/home/pratik/Install/share/dolfin/data/meshes/lshape.xml.gz")
LIN = FunctionSpace(mesh, 'CG', 1)
RT = VectorFunctionSpace(mesh, 'RT', 1)
W = LIN * RT

(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

f = Expression(("-1", "0"))
a = (sigma*tau - inner(u,grad(tau)) + inner(grad(sigma),v) + inner(curl(u),curl(v)))*dx
L = -inner(f,v)*dx

# Define function G such that G \cdot n = g
class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        values[0] = 0
        values[1] = 0
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh)

bc = DirichletBC(W.sub(1), G, u0_boundary)

w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split(True)
#plot(sigma)
plot(u)
#plot(mesh)
interactive()