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

205

views

0

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:

\[ ρ_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

Any suggestions?

```
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
T.assign(Constant(273.15))
L = w*rho_w(T)*dx
m_vec = assemble(L)
print(sum(m_vec.array()))
```

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
{
public:
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 "temperature_dependent_density.py", line 34, in <module>
L = w*rho_w(T)*dx
File "temperature_dependent_density.py", line 31, in rho_w
return Expression(code, theta = temperature, degree = 0)
File "/usr/lib/python3/dist-packages/dolfin/functions/expression.py", 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

0

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

If you really want the cpp version, then just change theta to be

`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.