Plotting material properties in subdomains given by python class

Thank you in advance. How can I plot the material distribution given by a python class that is used in a variational statement?

For example, in the code below I have the following lines to define k and then I use k in the variational statement.

# Define thermal conductivity
class Thermal_Conductivity(Expression):
    def __init__(self, mesh, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if markers[cell.index] == 0:
            values[0] = k_1
        else:
            values[0] = k_2

k = Thermal_Conductivity(mesh, degree=1)

The code is given below:

from fenics import *
from mshr import * #library for meshing
import math as m   #Import like this to prevent overiding of fenics' UFL math functions

# Use compiler optimizations
parameters["form_compiler"]["cpp_optimize"] = True
# Allow approximating values for points that may be generated outside
# of domain (because of numerical inaccuracies)
parameters["allow_extrapolation"]  = True
# Choose refinement algorithm (needed for solver to choose right refinement algorithm)
parameters["refinement_algorithm"] = "plaza_with_parent_facets"

R_inner  =  0.009525
R_2      =  0.025
R_outer  =  0.06
h_p   =  0.50
h_2   =  0.30
T_g   = -70.0 #degrees C
T_inf =  20.0 #degrees C
#Conduction values
k_1  =  16.3
k_2  =  0.25
#Convection values
h_a    =  10.0  #Free convection
h_g    =  500.0 #Forced convection

#Generate mesh with sub-domains
domain = Rectangle(Point(R_inner,0.0), Point(R_outer,h_p/2.0))
h = h_2/2.0
pipe_vertices = [Point(R_inner,0.0), Point(R_2,0.0), Point(R_2,h), Point(R_inner,h)]
pipe   = Polygon(pipe_vertices)
domain.set_subdomain(1, pipe)
meshResolution = 30
mesh = generate_mesh(domain, meshResolution)

# Define boundaries
tol = 1.0e-2
class insideEdge(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_inner)

class outsideEdge(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_outer)

insideEdge  = insideEdge()
outsideEdge = outsideEdge()

boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
insideEdge.mark(boundaries, 1)
outsideEdge.mark(boundaries, 2)
ds = Measure('ds')[boundaries]

# Define function space
V = FunctionSpace(mesh, 'P', 1)

# Define subdomain markers
markers = MeshFunction('size_t', mesh, 2, mesh.domains())

# Define thermal conductivity
class Thermal_Conductivity(Expression):
    def __init__(self, mesh, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if markers[cell.index] == 0:
            values[0] = k_1
        else:
            values[0] = k_2

k = Thermal_Conductivity(mesh, degree=1)

# Define variational problem
h_g = Constant(h_g)
T_g = Constant(T_g)
h_a = Constant(h_a)
T_inf = Constant(T_inf)
T = TrialFunction(V)
q = TestFunction(V)
a = k*inner(grad(T), grad(q))*dx() + h_g*T*q*ds(1) + h_a*T*q*ds(2)
L = h_g*T_g*q*ds(1) + h_a*T_inf*q*ds(2)


# Compute solution (adaptive)
T   = Function(V)
M   = T*dx()
tol = 1.0e-5
problem = LinearVariationalProblem(a,L,T)
solver  = AdaptiveLinearVariationalSolver(problem,M)
solver.solve(tol)
solver.summary()
plot(mesh.root_node(), title='Original coarse mesh')
plot(mesh.leaf_node(), title='Fine mesh')
plot(T.root_node(), title='Temperature distribution (coarse)')
plot(T.leaf_node(), title='Temperature distribution (fine)')

# Plot heat flow
plot(-k*grad(T.leaf_node()), scale=1.5, title='Heat flow vector field')
interactive()