Heat diffusion simulation in a 3D piston using FENICS

Question

I'm trying to simulate the heat diffusion in a 3D piston. I marked the boundaries on GMSH.

I have used a Dirichlet BC of 300 on the top face of piston. But the results look abnormal. There is a change in temperature along the edge on the top face, which should not be case.

Here is code snippet

from __future__ import print_function
from fenics import *
import numpy as np


mesh = Mesh("piston-3d.xml");

V = FunctionSpace(mesh, 'P', 1)

subdomains = MeshFunction("size_t", mesh, "piston-3d_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "piston-3d_facet_region.xml")


#Define neumann bcs
g_1= Constant(300.0)
g_2 = Constant(100.0)
g_3 = Constant(30.0)
g_4 = Constant(0.0) // On the cross-sectional faces, Neumann
g_5 = Constant(0.0) // On the cross-sectional faces, Neumann


bc1 = DirichletBC(V, g_1, boundaries, 1) // Top face of piston
bc2 = DirichletBC(V, g_2, boundaries, 2) // Liner face of the piston
bc3 = DirichletBC(V, g_3, boundaries, 3) // Inner side of piston and bottom

bcs = [bc1, bc2, bc3]

#define measures
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)


u_D = Constant(100.0)

# Define initial value
u_n = interpolate(u_D, V)


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = 0.0

a = dot(grad(u), grad(v))*dx 

L = f*v*dx - (g_4*v*ds(4) + g_5*v*ds(5))

# Time-stepping
u = Function(V)

# Compute solution
solve(a == L, u, bcs)

vtkfile = File("solution_steady.pvd")   
vtkfile << u

My question is why is the temperature of Top face not remaining constant (300)? Is it because of the boundary sharing?