How to compute f(T(x),x) for some function T

variational formulation
1

I obtained a function T (here, the temperature) from a former finite element computation, and I would like to use f(T(x),x) as a coefficient in a further finite element computation (here, linear elasticity) on the same mesh, where x is the spatial variable and f is some nonlinear function (here, the thermal expansion, given by a material dependent, piecewise linear interpolation of values given in a textbook).

I think I solved this by employing a custom Python class derived from a UserExpression (didn’t validate the results). However, this is extremely slow. Is there a better, possibly canonical solution for this problem?

I am using fenics version 2018.1.0.

2

Are you asking how to do compute T(x)?
To evaluate a function at a point, say x=(0,0,0), you can just do T(0,0,0). Then you can pass this to f.

3

Let me state my question more precisely:

  1. My function x \mapsto T(x) is an element of FunctionSpace(mesh,‘P’,1).
  2. My function (T,x) \mapsto f(T,x) is currently implemented in the form of two scipy.interpolate.interp1d functions T \mapsto f_1(T), f_2(T), each for one of the materials involved, that is f(T,x)=f_1(T) if x lies in material 1, while f(T,x)=f_2(T) if x lies in material 2.

Now I need to implement the function x \mapsto f(T(x),x) in a way that fenics accepts as a coefficient in a bilinear form, for instance as an element of FunctionSpace(mesh,‘P’,1).

Even without the explicit x dependence I wouldn’t know how to do this. For instance, simply using f_1(T) in the bilinear form with f_1 as above doesn’t work:

ValueError: object arrays are not supported

4

If there are not too many data points, you might be able to get away with using nested conditionals for the piecewise linear functions:

from dolfin import *

def interpUFL(T,Tvalues,fvalues):
    if((len(Tvalues)<2) or (len(Tvalues) != len(fvalues))):
        print("ERROR: Invalid data.")
    t = (Tvalues[1]-T)/(Tvalues[1]-Tvalues[0])
    combo = t*fvalues[0] + (1.0-t)*fvalues[1]
    if(len(Tvalues)==2):
        return combo
    return conditional(gt(t,0),combo,interpUFL(T,Tvalues[1:],fvalues[1:]))

# Test:
from matplotlib import pyplot as plt
mesh = UnitIntervalMesh(100)
x = SpatialCoordinate(mesh)
Tvalues = [0.0,0.1,0.7,1.0]
fvalues = [0.0,1.0,0.0,2.0]
plot(interpUFL(x[0],Tvalues,fvalues))
plt.show()
5

Thank you for your suggestion!

However, I would expect this solution to be quite slow. In the meantime I found a different method which is both simple and efficient, namely to work directly with the underlying (numpy) arrays.

Suppose that f is implemented as a function working with numpy arrays and T \in V. Then

alpha = dolfin.Function(V)
alpha.vector()[:] = f(T.vector()[:], V.tabulate_dof_coordinates())

does the job.

3 Likes