Question on solution of Poisson equation


I want to solve the Poisson equation on a domain like this:

Where I want the normal component of the gradient of the solution to bei zero on red and apply different constant Dirichlet BCs on blue and white.

When I project the negative gradient of the solution I obtain following field:

What I noticed is that the velocity spike at the node between the inflow and free slip boundaries increases with increasing mesh density. My question is whether this is physical, or something is wrong somewhere.

an mwe:

from dolfin import *
from matplotlib import pyplot as plt

class Boundary(SubDomain):
    def inside(self,x, on_boundary):
        return on_boundary

class Inflow(SubDomain):
    def inside(self,x, on_boundary):
        return (.3 <= x[0] <= .7) and x[1] < .001 and on_boundary

class Top(SubDomain):
    def inside(self,x, on_boundary):
        return near(x[1],1) and on_boundary

mesh = UnitSquareMesh(25,25,diagonal='left')
CG1  = FunctionSpace(mesh,'CG',1)
DG0  = VectorFunctionSpace(mesh,'DG',0)

p = TrialFunction(CG1)
q = TestFunction(CG1)
n = FacetNormal(mesh)

sub_domains = MeshFunction('size_t', mesh, 1)
Boundary().mark(sub_domains, 3)
Top().mark(sub_domains, 2)
Inflow().mark(sub_domains, 1)

ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)

a   = dot(grad(p), grad(q))*dx
L   = Constant(0)*q*dx
bc1 = DirichletBC(CG1, Constant(1),sub_domains,1)
bc2 = DirichletBC(CG1, Constant(0),sub_domains,2)
bc  = [bc1, bc2]
p   = Function(CG1)

solve(a == L, p,bc,solver_parameters={'linear_solver': 'mumps'})

u = project(grad(-p), DG0)

c = plot(p)
plot(mesh, linewidth=.2)

c = plot(u)
plot(mesh, linewidth=.2)