Use a `Function` as input for C++ `Expression`

6 months ago by
I'd like have a constitutive relation that depends on the independent variable, for example mass density as a function of temperature. In a simple case that I can easily express with UFL, that works fine, obviously:
from fenics import *

mesh = UnitIntervalMesh(400)

P2 = FiniteElement('CG', interval, 2)
V = FunctionSpace(mesh, P2)

w = TestFunction(V)
T = Function(V)

def rho_w(temperature):
    return 0.1*temperature

# just an example; should be allowed to vary over the domain
L = w*rho_w(T)*dx
m_vec = assemble(L)
Now, my expression for the density is a bit more involved, namely
\[ ρ_w(θ) = \sum_{i=0}^5 a_i (θ - 273.15)^i (-1 \times 10^7) + \sum_{i=0}^5 b_i (θ - 273.15)^i \]
with values for \(a_i\) and \(b_i\).

While it would be possible (just), to code this in UFL, I'd like to use a C++ expression here, like this
def rho_w(temperature):
    code = """
    class Rho : public Expression
        void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
            values[0] = 0.0;
            for (int i = 0; i < as.size(); ++i)
                values[0] += 1e-7 * as[i] * std::pow(theta - 273.15, i);
                values[0] += bs[i] * std::pow(theta - 273.15, i);

        std::vector<double> as = {4.8863e-7, -1.6528e-9, 1.8621e-12, 2.4266e-13, -1.5996e-15, 3.3703e-18};
        std::vector<double> bs = {1.0213e3, -7.7377e-1, 8.7696e-3, -9.2118e-5, 3.3534e-7, -4.4034e-10};
        double theta;
    return Expression(code, theta = temperature, degree = 0)
Unfortunately, this fails with
Traceback (most recent call last):
  File "", line 34, in <module>
    L = w*rho_w(T)*dx
  File "", line 31, in rho_w
    return Expression(code, theta = temperature, degree = 0)
  File "/usr/lib/python3/dist-packages/dolfin/functions/", line 177, in __init__
    setattr(self, member, value)
TypeError: in method 'Rho_theta_set', argument 2 of type 'double'
Aborted (core dumped)​
running Python3 on the stable docker image of FEniCS.

Any suggestions?

Edit The temperature is supposed to be variable over the domain. Using "Constant" was just for the sake of an example. The temperature could be from a previous solution step for example.
Community: FEniCS Project

1 Answer

6 months ago by
Try not to call rho_w with argument being of class Function.  As in your cpp code, theta is double, so call rho_w(273.15)
Well this was supposed to be just an example. In the "real" code, T will vary over the domain, as a solution of a previous `solve(a == L, T, bc)` for example.
written 6 months ago by Christoph Pohl  
I see, then I would go for UFL. Otherwise you must pass Function to the Expression and evaluate the function fo const Array<double>& x and I think it will be less efficient than UFL. In addition, with UFL you can later simply solve the implicit problem (not taking temperature from previous step, but as an unknown).

If you really want the cpp version, then just change theta to be Function and call theta(x) to evaluate.
written 6 months ago by Michal Habera  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »