# Refine the mesh uniformly
mesh.refine()

# Represent a function using higher order (accuracy)
V2 = FunctionSpace(mesh, "CG", 2)
exact = Expression("16*x[0]*(1.0 - x[0])*x[1]*(1.0 - x[1])", element=V2.ufl_element())

# Compute the L2 norm
e = exact - U
M = e*e*dx
e_L2 = sqrt(assemble(M, mesh=mesh))
print "error in L2 norm: ", e_L2