Here’s my code:

```
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
L = 0.01
W = 0.01
nx = 200
ny = 200
mesh = RectangleMesh(Point(0, 0), Point(L, W), nx, ny)
V = FunctionSpace(mesh, 'P', 1)
def k(T):
return 174.9274 - 0.1067 * T + 5.0067e-5 * T**2 - 7.8349e-9 * T**3
mu = Constant(19300)
c = Constant(137)
q = Expression('1e7', degree=0)
#h = Expression('70000', degree=0)
T_ini = Expression('298.15', degree=0)
#T_env = Expression('323.15', degree=0)
f = Constant(0.0)
t_total = 100
num_steps = 5000
dt = t_total / num_steps
class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]-0.00)<DOLFIN_EPS
class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]-0.01)<DOLFIN_EPS
class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1]-0.00)<DOLFIN_EPS
class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1]-0.01)<DOLFIN_EPS
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bc_bottom = DirichletBC(V, T_ini, by0)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
#boundary_markers.set_all(9999)
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)
boundary_conditions = {0: {'Neumann': Constant(0)}, 1: {'Neumann': Constant(0)}, 2: {'Dirichlet': T_ini}, 3: {'Neumann': q}}
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
bcs.append(bc)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
T_n = interpolate(T_ini, V)
T = Function(V)
v = TestFunction(V)
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
#if boundary_conditions[i]['Neumann'] != 0:
q = boundary_conditions[i]['Neumann']
integrals_N.append(q*v*ds(i))
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
h, T_envi = boundary_conditions[i]['Robin']
integrals_R.append(h*(T - T_envi)*v*ds(i))
F = mu*c*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu*c*T_n*v*dx - dt*f*v*dx - dt*sum(integrals_N)
xdmffile_T = XDMFFile('10_10_100s_slab/temperature.xdmf')
timeseries_T = TimeSeries('10_10_100s_slab/temperature_series')
File('10_10_100s_slab/mesh.xml.gz') << mesh
xdmffile_T.parameters["flush_output"] = True
xdmffile_T.parameters["functions_share_mesh"] = True
progress = Progress('Time-stepping', num_steps)
t = 0
times = [0]
temperatures=[298.15]
plt.figure(figsize=(12,8))
for n in range(num_steps):
t += dt
solve(F == 0, T, bcs)
plot(T)
point = Point(0.005, 0.01)
temperature = T(point)
temperatures.append(temperature)
times.append(t)
xdmffile_T.write(T, t)
timeseries_T.store(T.vector(), t)
T_n.assign(T)
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
```