Hi everyone,

I am writing a time dependent coupled temperature, flow, and pressure solver based on the Cahn-Hilliard FEniCS demo. I have boundary conditions and source terms that depend on time as well, and I have been using classes to define these. For example,

```
class uExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = -np.exp(-c*self._t(x))*np.sin(2*np.pi*(1-x[1]))
value[1] = np.exp(-c*self._t(x))*np.sin(2*np.pi*x[0])
def value_shape(self):
return(2,)
```

to define the exact solution I am comparing against. Then,

```
t = Constant(0.0)
uexact = uexactDef(t)
```

to evaluate it the first time. I use this class for Dirichlet boundary conditions successfully, for example

```
bcu_wall = DirichletBC(W.sub(0), uexact, wall)
```

which works well. The timestepping is then done in a while loop, as in the Cahn-Hilliard demo, i.e.

```
t2 = +dt
t.assign(t2)
```

where I use t2 as a helper, (because t.assign(t2)) makes sense to me). This timesteps the boundary condition as well (and other time dependent parts). However, I am also imposing the normal stress at one of the boundaries, for which I need `uexact`

and `pexact`

. I want to do

```
def epsilon(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
[v[1].dx(0), v[1].dx(1)]]))
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-p*Identity(2)
b1 = - dot(dot(sigma(u_exact, p_exact), v), n) * ds(1)
```

where ds(1) is the boundary where I want to impose the condition, and `b1`

is the boundary term added to the variational form. However, this does not work. I get the error

```
error("Cannot determine geometric dimension from expression.")
```

I am currently using the alternative of defining a separate expression

```
u_Neumann = Expression(('-exp(-c*t)*sin(2*pi*(1-x[1]))', 'exp(-c*t)*sin(2*pi*x[0])'), c=c, t=0.0, degree=2, domain=mesh)
p_Neumann = Expression('-exp(-c*t)*(cos(2*pi*x[0])*cos(2*pi*(1-x[1])))', c=c, t=0.0, degree=2, domain=mesh)
```

and then

```
u_Neumann.t = t2
p_Neumann.t = t2
```

which then works as expected, namely

```
b1 = - dot(dot(sigma(u_Neumann, p_Neumann), v), n) * ds(1)
```

gives no errors and the correct answer but with *significant* slow down. My questions are:

- Is there a way to use the class
`uExactDef`

to impose the normal stress condition I need without errors? - Is there another way of imposing it instead of
`Expression`

? This makes my computations incredibly slow. - If I want everything to be consistent, would you suggest leaving everything as Expressions or classes?

NB I realise that I didn’t supply the MWE, but it is because the issue I am having is only relevant for the time dependent coupled case (as the time derivative is in the temperature problem). Therefore my MWE is actually quite long, and I think that my question is more related to the correct form of `uexact`

so that it can be used in `sigma(uexact, pexact)`

. If you require the MWE, I will copy it as well (as simplified as I can).