Problem in defining User-defined expressions by subclassing

variational formulation
1

Hello all,

I am trying to define my load as a time dependent expression as:
This load has constant magnitude = gt

gt = Expression('A*sin(omega*t)', A = 1.0,
                 omega = 1.0, t = 0.0, degree = 2)  
class IC(UserExpression):
    def eval(self, value, x):
        value[0] = 0.0
        value[1] = 0.0
        if x[2] == z1 and  0.051505 <= x[1] <= 0.098495 and  gt <= x[0] <= 1.1*gt:
            value[2] == -rho*gamma
        else :
            value[2] == 0.0
    def value_shape(self):
        return (3,)
tr = IC(element = V.ufl_element())

when I am calling tr which is traction on the boundary I am getting this error:

ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

So I can see the error is in this line inside the UserExpression :

if x[2] == z1 and  0.051505 <= x[1] <= 0.098495 and  gt <= x[0] <= 1.1*gt

any hint would be much appreciated. Thanks.

2

There are several issues with your code. The first one is that you have not defined an initializer.
You Need:

        def __init__(self, **kwargs):
            super().__init__(**kwargs)

Or else

Has no effect. Secondly,

This should be assigned (using one equal sign) not two. Two equal signs evaluate this as a boolean.
Thirdly,

Cannot be evaluated (Which is your error message). I would instead use the initializer to define a lambda function that defined gt.

3

Thanks Jorgen, because It’s my first time that I am dealing with this concepts, Could you please give a reference or some example that walk me through the process to understand the structure of user-defined expression.
P.S. I’ve read Dolfin tutorial but I could’nt understand how should I modify that for my time-dependent expression.

4

I was thinking something like this:

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(20,20,2)
V = VectorFunctionSpace(mesh, "CG", 1)
class IC(UserExpression):
    def __init__(self, t, A, omega, z1, rho, gamma, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.A = A
        self.omega = omega
    def eval(self, value, x):
        value[0] = 0.0
        value[1] = 0.0
        gt = self.A*np.sin(self.omega*self.t)
        if x[2] == z1 and  0.1 <= x[1] <= 0.3 and  gt <= x[0] <= 1.1*gt:
            value[2] = -rho*gamma
        else :
            value[2] = 0.0
    def value_shape(self):
        return (3,)
A = 1
omega = 1
t = 0
z1 = 0.5
rho = 0.1
gamma = 0.2
tr = IC(t, A, omega, z1, rho, gamma, element = V.ufl_element())
print(assemble(inner(tr, tr)*dx(domain=mesh)))
File("u.pvd") << project(tr, V)
tr.t = 2
print(assemble(inner(tr, tr)*dx(domain=mesh)))

Note that I have had to use random parameters as you havent supplied me with values, and to show that it is time dependent, consider the output:

5.99999999999999e-07
2.183333333333337e-06
1 Like
5

Hello @dokken

I am trying to make a function in parts. The first would be a gaussian distribution at a certain height H of the mesh and the other would be constant after that height H. However, I have not been successful.

from dolfin import *

mesh = RectangleMesh(Point(0.0, 0.0), Point(100.0, 50.0), 10, 10)

V = FunctionSpace(mesh, "Lagrange", 1)
mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L

class porosity(UserExpression):
    def __init__(self, phi_0, phi_air, H, **kwargs):
        super().__init__(**kwargs)
        self.phi_0 = phi_0
        self.phi_air = phi_air
        self.H = H
    def eval(self, value, x):
        if x[1] <= H + DOLFIN_EPS:
            value[0] = self.phi_0
        else:
            value[0] = self.phi_air
    def value_shape(self):
        return ()

mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L
phi_0 = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi)
phi_air = Expression('1.0', degree=1)
phi = porosity(phi_0, phi_air, H, degree=1)
File("u.pvd") << project(phi, V)

Any suggestions?

6

@dokken

I made another code and this one seems to work but I don’t know if it does what I require.

from dolfin import *

mesh = RectangleMesh(Point(0.0, 0.0), Point(100.0, 50.0), 10, 10)

V = FunctionSpace(mesh, "Lagrange", 1)
mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L

class porosity(UserExpression):
    def __init__(self, mu_phi, x0_phi, H, **kwargs):
        super().__init__(**kwargs)
        self.mu_phi = mu_phi
        self.x0_phi = x0_phi
        self.H = H
    def eval(self, value, x):
        phi0 = self.np.exp(-self.mu_phi * pow(x[0] - self.x0_phi, 2.0))
        if x[1] <= H + DOLFIN_EPS:
            value[0] = phi0
        else:
            value[0] = 1.0
    def value_shape(self):
        return ()

mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L
phi = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi, H=H)
File("u.pvd") << project(phi, V)
7

Mixing expression and UserExpression is not adviced, choose one or the other.

8

Okay,

I did it this way, but in the corner there seems to be problems.

image
image1920×1080 281 KB

This is the code

from dolfin import *
import numpy as np

mesh = RectangleMesh(Point(0.0, 0.0), Point(100.0, 50.0), 10, 10)

V = FunctionSpace(mesh, "Lagrange", 1)
mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L

class porosity(UserExpression):
    def __init__(self, mu_phi, x0_phi, H, **kwargs):
        super().__init__(**kwargs)
        self.mu_phi = mu_phi
        self.x0_phi = x0_phi
        self.H = H
    def eval(self, value, x):
        phi0 = np.exp(-self.mu_phi * pow(x[0] - self.x0_phi, 2.0))
        if x[1] <= self.H + DOLFIN_EPS:
            value[0] = phi0
        else:
            value[0] = 1.0
    def value_shape(self):
        return ()

mu_phi = 0.5
L = 100
H = 25
x0_phi = 0.5 * L
phi = porosity(mu_phi=mu_phi, x0_phi=x0_phi, H=H, element = V.ufl_element())
File("u.pvd") << project(phi, V)