Dear Community,

I am trying to set nonlinear material properties, which depend on the solution in DOLFINx (version: 0.1.1.0). The task is pretty simple if polynomials are used. It is demonstrated in the following code (2D rectangle is cooling down during some time through heat exchange with the ambient):

```
import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp
mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
# Define initial conditions
u_n = dolfinx.Function(V)
with u_n.vector.localForm() as loc:
loc.set(1000)
#Define solution variable
uh = dolfinx.Function(V)
with uh.vector.localForm() as loc:
loc.set(1000)
v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)
## Material properties
# Definition using polynom
def heat_cond_polynom(u):
return 0.00004*u + 0.02
t = 0
dt=0.1
#Time loop
for i in range(1,6):
t += dt
print("time: " + str(np.round(t, 2)))
q_ = heat_cond_polynom(uh)
F = 0.54 * 0.0078 * uh * v * ufl.dx + q_*dt*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n + dt*f) * v * ufl.dx
# Boundary conditions (Robin on all sides)
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
marker = 1
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
F += dt * 0.03 * ufl.inner(uh-20, v) * ds(marker)
#Problem definition
problem = dolfinx.fem.NonlinearProblem(F, uh)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")
#Print max temperature
print(np.max(np.round(np.array(uh.vector.array),2)))
# Update solution at previous time step (u_n)
with uh.vector.localForm() as loc, u_n.vector.localForm() as loc_n:
loc.copy(result=loc_n)
```

output:

```
time: 0.1
Number of interations: 3
temperature: 556.22
time: 0.2
Number of interations: 3
temperature: 298.18
time: 0.3
Number of interations: 3
temperature: 164.14
time: 0.4
Number of interations: 2
temperature: 94.73
time: 0.5
Number of interations: 2
temperature: 58.7
```

However, i need to set material properties using data points rather than polynomials. Just defining e.g. the interp1d function instead of the polynom results in “ValueError: object arrays are not supported”. There is a number of discussions regarding this issue in dolfin such as RemDelaporteMathurin and freya. Combining approaches suggested in these threads the following code is obtained:

```
import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp
mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
# Define initial conditions
u_n = dolfinx.Function(V)
with u_n.vector.localForm() as loc:
loc.set(1000)
#Define solution variable
uh = dolfinx.Function(V)
with uh.vector.localForm() as loc:
loc.set(1000)
v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)
## Material properties
# Definition using interpolation
def heat_cond_interpolation(u):
heat_conduction_data = sp.interp1d([0,1000],[0.02,0.06])
return heat_conduction_data(u)
t = 0
dt=0.1
#Time loop
for i in range(1,6):
t += dt
print("time: " + str(np.round(t, 2)))
q_ = dolfinx.Function(V)
q_.vector.array[:] = heat_cond_interpolation(np.array(uh.vector.array))
F = 0.54 * 0.0078 * uh * v * ufl.dx + q_*dt*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n + dt*f) * v * ufl.dx
# Boundary conditions (Robin on all sides)
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
marker = 1
facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
F += dt * 0.03 * ufl.inner(uh-20, v) * ds(marker)
#Problem definition
problem = dolfinx.fem.NonlinearProblem(F, uh)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")
#Print max temperature
print(np.max(np.round(np.array(uh.vector.array),2)))
# Update solution at previous time step (u_n)
with uh.vector.localForm() as loc, u_n.vector.localForm() as loc_n:
loc.copy(result=loc_n)
```

which gives the following output and it is obvious that the temperatures as well as the number of iterations are not the same as they were while using the polynom:

```
time: 0.1
Number of interations: 1
temperature: 516.69
time: 0.2
Number of interations: 1
temperature: 269.04
time: 0.3
Number of interations: 1
temperature: 147.52
time: 0.4
Number of interations: 1
temperature: 85.95
time: 0.5
Number of interations: 1
temperature: 54.17
```

Basically, what i need is to reproduce the output of the first code using the material definition (Definition using interpolation) from the second code. There are other related topics such as RemDelaporteMathurin. However it does not work anymore since the “UserExpression” class has been deprecated in DOLFINx. I tried to come up with another solution using recommendations from Dokken, but have not succeeded. It has been more than two years since RemDelaporteMathurin posted his question. I have been hoping that during this time a simple solution to that issue appeared and is now available in DOLFINx.

Thank you

Sincerely

Anton