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 = -np.exp(-c*self._t(x))*np.sin(2*np.pi*(1-x)) value = np.exp(-c*self._t(x))*np.sin(2*np.pi*x) 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
pexact . I want to do
def epsilon(v): return sym(as_tensor([[v.dx(0), v.dx(1)], [v.dx(0), v.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))', 'exp(-c*t)*sin(2*pi*x)'), c=c, t=0.0, degree=2, domain=mesh) p_Neumann = Expression('-exp(-c*t)*(cos(2*pi*x)*cos(2*pi*(1-x)))', c=c, t=0.0, degree=2, domain=mesh)
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
uExactDefto 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).