Update Expression for DirichletBC

1

Hi
I’m trying to implement a DirichletBC based on a UserExpression class which is time dependent.

In the following MWE, 2 different DirichletBC are set, one that only rely on the UserExpression and another one which is pondered by 2 (which is similar to what i’d like to achieve). In the latter case, the boundary condition is not updated for some reason. Does anyone have an explanation ?

Cheers!
Rem

from fenics import *

mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)
T = Expression("t", t=1, degree=1)


class Coeff(UserExpression):
    def __init__(self, T, **kwargs):
        super().__init__(kwargs)
        self._T = T

    def eval_cell(self, value, x, ufc_cell):
        value[0] = 2*self._T(x)*x[0]

    def value_shape(self):
        return ()


coeff = Coeff(T)

sm = MeshFunction("size_t", mesh, 1, 0)
right = CompiledSubDomain('x[0] > 0.75')
right.mark(sm, 2)
bc1 = DirichletBC(V, coeff, sm, 2)
bc2 = DirichletBC(V, 2*coeff, sm, 2)

F = dot(grad(u), grad(v))*dx

set_log_level(30)
for i in range(1, 10):
    T.t = i
    coeff._T = T

    solve(F == 0, u, bc1)
    a = u(0.5, 0.5)
    solve(F == 0, u, bc2)
    b = u(0.5, 0.5)

    print(a, b)
2

This is because the class of this two objects differs, and is handled differently by the DirichletBC.
You can see this with the following code:

import ufl
print(isinstance(Coeff(T), ufl.Coefficient))
print(isinstance(2*Coeff(T), ufl.Coefficient))

which returns

True
False

This means that the latter input is treated as something that has to be projected internally in the DirichletBC, and is thus not prone to updates. What you should do here is to add a scaling parameter to your UserExpression:

from fenics import *

mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)
T = Constant(1)


class Coeff(UserExpression):
    def __init__(self, T, **kwargs):
        super().__init__(kwargs)
        self._T = T
        self.alpha = 1

    def eval_cell(self, value, x, ufc_cell):
        value[0] = self.alpha*2*self._T(x)*x[0]

    def value_shape(self):
        return ()

coeff = Coeff(T)
sm = MeshFunction("size_t", mesh, 1, 0)
right = CompiledSubDomain('x[0] > 0.75')
right.mark(sm, 2)
bc1 = DirichletBC(V, coeff, sm, 2)

bc2 = DirichletBC(V, coeff, sm, 2)

F = dot(grad(u), grad(v))*dx

set_log_level(30)
for i in range(1, 10):
    coeff.alpha = 1
    coeff._T.assign(i)
    solve(F == 0, u, bc1)
    a = u(0.5, 0.5)
    coeff.alpha = 2
    solve(F == 0, u, bc2)
    b = u(0.5, 0.5)

    print(a, b)
1 Like
3

Cheers ! works perfectly !

4

Dear all,

thanks for this (solved) problem.
A naive question: how should we change the Coeff class for prescribing conditions to a vector?

Thanks

Alessandro

5

value is then multi-dimensional, as for instance shown in:

# Define velocity field
class uExactDef(UserExpression):
    def __init__(self, t, **kwargs):
        super().__init__(kwargs)
        self._t = t

    def eval_cell(self, value, x, ufc_cell):
        value[0] = -self._t(x) * np.sin(2 * np.pi * (1 - x[1]))
        value[1] = Constant(0.0)

    def value_shape(self):
        return (2,)

as shown in: Same time dependent temperature problem in two different ways: one works, one doesn't. Why?

6

Thanks!

It worked perfectly

Alessandro