Compute PDE parameters via interpolation when using FENICS as solver

I am trying to solve the heat transfer equation considering the temperature-dependent thermal conductivity. I have experimental data for the thermal conductivity. I have been trying to interpolate this parameter using the temperature as input and the interpld in Python. However, I am getting the following error when doing this in FENICS.



I know this error is related to the format of T. How can I convert function type data in FENICS into numpy data, for example, so I can use interpld to interpolate the thermal conductivity? Or is there another way to do this? The Python script is below.


# FEATool Multiphysics -> FEniCS project conversion
# File featool-fenics.py created 21-Mar-2025 13:22:40
try:
    from fenics import *
except:
    raise ImportError("Could not find/import fenics.")
from timeit import default_timer as timer
import atexit
comm = MPI.comm_world
size = comm.Get_size()
rank = comm.Get_rank()

# import libraries
import pandas as pd
from scipy.interpolate import interp1d
   
# read data from txt file
tc_matrix = pd.read_csv("/home/christ96/Documents/MATLAB/FEA_TOOLBOX/Data/ThermalCond2.txt")
Tv  = tc_matrix.iloc[:,0].to_numpy()
TCv =  tc_matrix.iloc[:,1].to_numpy()


# Mesh and subdomains.
with HDF5File(comm, "/tmp/featool-fenics.h5", "r") as f:
    mesh = Mesh()
    f.read(mesh, "/mesh", False)
    subd = MeshFunction("size_t", mesh, mesh.topology().dim())
    f.read(subd, "/mesh")
    bdr = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    f.read(bdr, "/boundary")
dx = Measure("dx", domain=mesh, subdomain_data=subd)
ds = Measure("ds", domain=mesh, subdomain_data=bdr)
n = FacetNormal(mesh)
nz = n[0]
nr = n[1]


# Finite element spaces.
E1 = FiniteElement("P", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, E1)


#  Function spaces.
T = Function(V)
T_t = TestFunction(V)
T_n = Function(V)


# Model constants and expressions.
dt = 0.001
tmax = 0.35
t = 0.0
t_i = Expression("t", degree=1, t=0.0)
z = Expression("x[0]", degree=1)
r = Expression("x[1]", degree=1)
h_grid = CellDiameter(mesh) if "CellDiameter" in globals() else CellSize(mesh)


# thermal conductivity
def tc_inter(T, Tv, TCv):
    # interpolate
    f = interp1d(Tv, TCv,kind='cubic')
   
    tc_val = f(T)
   
    return tc_val

rho_ht_0 = Constant(0.8)
cp_ht_0 = Constant(7000)
u_ht_0 = Constant(0)
v_ht_0 = Constant(0)
q_ht_0 = Constant(5000000)
k_ht_0 = tc_inter(T, Tv, TCv)


# Mass forms.
m = ((r*rho_ht_0*cp_ht_0)*T*T_t)*dx


# Bilinear forms.
a = dt*(((r*rho_ht_0*cp_ht_0*u_ht_0)*T.dx(0)*T_t + (r*k_ht_0)*T.dx(0)*T_t.dx(0) + (r*rho_ht_0*cp_ht_0*v_ht_0)*T.dx(1)*T_t + (r*k_ht_0)*T.dx(1)*T_t.dx(1))*dx)


# Linear forms.
f = ((dt*(r*q_ht_0) + (r*rho_ht_0*cp_ht_0)*T_n)*T_t)*dx


# Boundary conditions.
dbc0 = DirichletBC(V, Constant(300), bdr, 1)
dbc1 = DirichletBC(V, Constant(300), bdr, 2)
dbc = [dbc0, dbc1]

g = dt*( ((rho_ht_0*cp_ht_0*(nz*u_ht_0+nr*v_ht_0)*T))*T_t*ds(0) + \
    ((rho_ht_0*cp_ht_0*(nz*u_ht_0+nr*v_ht_0)*T))*T_t*ds(3) )


# Initial conditions.
dvar = ["T"]
T_n = project(Constant(5000), V)
assign(T, T_n)


# Solver parameters.
param = {
    "newton_solver": {
        "linear_solver": "mumps" if "mumps" in linear_solver_methods() else "default",
        "preconditioner": "default" if "default" in krylov_solver_preconditioners() else "default",
        "maximum_iterations": 20,
        "relaxation_parameter": 1.0,
        "relative_tolerance": 1e-06,
        "absolute_tolerance": 1e-06,
        }
    }


# Solve.
f_data = HDF5File(comm, "/tmp/featool-fenics.h5", "a")
atexit.register(lambda: f_data.close())
t_start = timer()
f_data.write(T, "/" + dvar[0], 0.0)
n_ts = int(-(tmax // -dt))
for i_step in range(n_ts):
    t += dt
    print("Time step %i, time = %f" % (i_step+1, t)) if rank == 0 else None
    t_i.t = t
    f = ((dt*(r*q_ht_0) + (r*rho_ht_0*cp_ht_0)*T_n)*T_t)*dx
    solve(m + a - f - g == 0, T, dbc, solver_parameters=param)
    T_n.assign(T)
    f_data.write(T, "/" + dvar[0], t)
t_solve = comm.gather(timer() - t_start)
print("t_solve = %g (tot: %g)" % (max(t_solve), sum(t_solve))) if rank == 0 else None