Fenics solver and user-defined functions

I am trying to solve using the Fenics solver the heat transfer equation considering temperature-dependent parameters. However, I am getting the following error:



It seems that the user-defined function density  needs to be included somehow in the solver. Any thoughts on how this can be fixed? I am using FEATool Multiphysics as a Matlab toolbox on Linux.

Hard to say why the density was not defined without more information, but the workaround would be to manually edit the converted Python code (using the "Script" tab in the "FEniCS" solver settings dialog box) and add the variable manually before calling the solver.

All the user-defined functions are added in the Python script (see the image below), and still, the issue appears.

From the attached screenshot it doesn't look like you have defined "density" but are trying to used in defining "rho_ht_0", but it is hard to say as you haven't attached a complete script or model.

The complete script is shown below.

# FEATool Multiphysics -> FEniCS project conversion
# File featool-fenics.py created 18-Mar-2025 09:39:41
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()


# 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)

u_ht_0 = Constant(0)
v_ht_0 = Constant(0)
rho_ht_0 = density(T)
cp_ht_0 = cp_inter(T)
k_ht_0 = tc_inter(T)
q_ht_0 = input_power(T)


# 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(Expression("T0f(x[1])", degree=1), 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

As noted in the first and subsequent replies, the expressions for "density", as well as "cp_inter", "tc_inter", and "power" have not been defined in the conversion, so you will have to add them manually by editing the script before running it.

How can a matlab function be called from the fenics script? The fenics script has the matlabengine library, but there is a problem with the data type of the argument. I think the problem is related how fenics treats the variable T. I would appreciate it if some guidance could be provided. The fenics script is below.

# FEATool Multiphysics -> FEniCS project conversion
# File featool-fenics.py created 20-Mar-2025 15:57:12
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 matlab.engine
eng = matlab.engine.start_matlab()


# 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)

cp_ht_0 = Constant(7000)
k_ht_0 = Constant(4)
u_ht_0 = Constant(0)
v_ht_0 = Constant(0)
q_ht_0 = Constant(5000000)
rho_ht_0 = eng.density(T)


# 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

The error can be seen in the figure below.
 

Not quite sure this is a good idea, but the specific error here indicates that T in your engine call is a FEniCS object (dolfin.function..) which naturally Matlab doesn't know how to handle. So T would need to be converted to a datatype that Matlab can manage. You might want to also check with the Mathworks and/or FEniCS support/and forum for this as it is quite out of scope here.

That said calling Matlab from Python from a running Matlab toolbox does not seem like it would be a good idea. Moreover, variables involving dependent variables in FEniCS are called dynamically and may be evaluated many times in each solver iteration, making it very inefficient (as calling Matlab is surely quite expensive). The best approach here would probably be to convert your Matlab function to Python if its not too complicated.

Thank you for the recommendation.