convergence order of FEM on L-shape domain

Hi, I want to test the convergence of linear FEM method for Poisson equation

                              - \Delta u = f, with Dirichlet BC

on a L-shape domain with the origin as the corner.

I set exact solution u = r^{2/3}, then the right hand side f = -4/9*r^{-4/3}.

I try to use the following Python code, but it does not work.

Maybe, one reason is that the function f(x) is singular near the
origin, but f(x) is in L^{1}, the righthandside L = (f,v) is still
reasonable, I am not sure whether the numerical qudrature fails for
L = (f,v) in this case.

Thanks for any help ~

# ---------begin of the program-------
from dolfin import *

# the L-shape mesh can be find at
# http://fenicsproject.org/pub/data/meshes/
mesh = Mesh("lshape.xml.gz")
mesh.coordinates()[:] = (mesh.coordinates()[:] -0.5)*2.0

ue = Expression("pow(x[0]*x[0]+x[1]*x[1],1.0/3.0)")
f = Expression("-4.0/9.0*pow(x[0]*x[0]+x[1]*x[1],-2.0/3.0)")

def D_boundary(x, on_boundary):
    return on_boundary

V = FunctionSpace(mesh,'CG',1)
uh = TrialFunction(V)
vh = TestFunction(V)

a = inner(grad(uh),grad(vh))*dx
L = f*vh*dx

A = assemble(a)
b = assemble(L)
bc = DirichletBC(V, ue, D_boundary)
bc.apply(A,b)

#---------------------------------------------------
solver = KrylovSolver("gmres", "amg")
solver.parameters["absolute_tolerance"] = 1E-13
solver.parameters["relative_tolerance"] = 1E-13
solver.parameters["maximum_iterations"] = 1000000
solver.parameters["nonzero_initial_guess"] = True
#--------------------------------------------
uh = Function(V)
#uh.vector()[:] = uh.vector()[:]*0.8
solver.solve(A, uh.vector(), b)

EL2 = errornorm(ue,uh,"L2", degree_rise = 1)
EH1 = errornorm(ue,uh,"H1", degree_rise = 1)

print EL2, EH1
# ---------end of the program-------