Value appearing in solutions at points not defined on the mesh

Question

Hi,

I'm trying to solve time dependent Navier-Stokes equations using strictly Dirichlet boundary conditions in the unit square. When I plot the solution at each time step I appear to be getting values at points in the top right which are not even defined on my mesh. The code included below will generate the plot with the "ghost values".

from dolfin import *
from mshr import *

#Time date and Reynolds number
Re = 400.0
t_init = 0.0
t_final = 10.0
t_num = 10000
nu = 1.0/Re
dt = ( t_final - t_init ) / t_num
t = t_init

#Define mesh and subdomains
class Noslip(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class topBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] > 1.0 - DOLFIN_EPS and on_boundary

h = 10
mesh = UnitSquareMesh(h, h)

# Create mesh functions over the cell facets
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

sub_domains.set_all(2)

#Mark no slip bounday conditions
noslip = Noslip()
noslip.mark(sub_domains,0)

#Mark boundary condition at the top of the domain
top_bound = topBoundary()
top_bound.mark(sub_domains,1)

#plot(sub_domains,interactive = True)

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

#Define boundary conditions
g0 = Constant((0.0,0.0))
g1 = Constant((1.0,0.0))

bc0 = DirichletBC(W.sub(0),g0,sub_domains,0)
bc1 = DirichletBC(W.sub(0),g1,sub_domains,1)
bcs = [bc0,bc1]

#Define functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
w = Function(W)
u0 = Function(V)
uplot = Function(V)  # used for plotting in loop so get only 1 frame

#Define the weak formulations of the Navier Stokes equations using forward euler
LNSE = inner(u0,v)*dx

NSE = (inner(u,v) + dt*(.5*inner(grad(u)*u0,v) - .5*inner(grad(v)*u0,u)\
                       + nu*inner(grad(u),grad(v)) - div(v)*p) + q*div(u) )*dx

u_init = Expression(("0.0","0.0"))
u = interpolate(u_init,V)

# timestepping
while t < t_final * (1. + 1.e-10):
    t += dt
    u0.assign(u)
    solve(NSE == LNSE, w, bcs, solver_parameters=dict(linear_solver="lu"))
    print "The time is:"
    print(t)
    # split vel, pres out of w (deepcopy)
    (u, p) = w.split(True)
    uplot.assign(u)
    plot(u,interactive = True)
    #plot(uplot,interactive = True)