 # Converting a set of functions to a class

Hello,

I tried to convert the following set of functions

``````############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA

def euler_c_s( t0, c_s, dt, t, u):
temp = -0

# Iterating till the point at which we
# need approximation
while t0 < t:
temp = c_s

c_s = c_s + assemble(dt * func_c_s(t0, c_s, u))
t0 = t0 + dt
return c_s

def func_c_s(t, c_s, u):
return (u-c_s)*ds

def value_shape(self):
return ()
``````

to a class, derived from an UserExpression, by the following code

``````############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
class My_c_s(UserExpression):

def __init__(self, t, u, **kwargs):
super().__init__(**kwargs)
self.t = t
self.u = u

def euler_c_s(self, t0, c_s, dt, t, u):
temp = -0

# Iterating till the point at which we
# need approximation
while t0 < self.t:
temp = c_s

c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u))
t0 = t0 + dt
return c_s

def func_c_s(self, t, c_s, u):
return (self.u-c_s)*ds

def value_shape(self):
return ()
``````

Replacing the set of functions from above with the class `My_c_s(UserExpression)` in the running code, which also serves as a MWCE,

``````from fenics import *
import numpy as np
import dolfin as df

############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
class My_c_s(UserExpression):

def __init__(self, t, u, **kwargs):
super().__init__(**kwargs)
self.t = t
self.u = u

def euler_c_s(self, t0, c_s, dt, t, u):
temp = -0

# Iterating till the point at which we
# need approximation
while t0 < self.t:
temp = c_s

c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u))
t0 = t0 + dt
return c_s

def func_c_s(self, t, c_s, u):
return (self.u-c_s)*ds

def value_shape(self):
return ()

######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
def eval(self, value, x):

if x >= 0 and x <= 0.4:
value = u_avg_domain_t0

elif x >= 0.6 and x <= 1.0:
value = u_avg_domain_t0

else:
value = 0

def value_shape(self):
return ()

######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)

################################ geometry #####################################

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x, 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x, 0.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1.0)

class Obstacle_left(SubDomain):
def inside(self, x, on_boundary):
return (between(x, (0.0, 0.4)))

class Obstacle_right(SubDomain):
def inside(self, x, on_boundary):
return (between(x, (0.6, 1.0)))

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle_left = Obstacle_left()
obstacle_right = Obstacle_right()

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
domains.set_all(0)
obstacle_left.mark(domains, 1)
obstacle_right.mark(domains, 2)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, space_dim-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

##################### define constants #########################################
t = 0
t0 = 0
T = 2.0
num_steps = 10
dt = T / num_steps

u_avg_domain_t0 = 1
u_t0 = My_u_t0()
u = interpolate(u_t0, V)

c_s_0 = 0
c_sat = 0.9
#c_s = euler_c_s(t0, c_s_0, dt, t, u)
c_s = My_c_s(t, u, degree=1)

while t < T:
t = t+dt
#c_s = euler_c_s(t0, c_s_0, dt, t, u)
c_s.t = t
c_s.u = u

if c_s < c_sat:
print('t = %.1f: c_s = %.9f' % (t, c_s))
else:
print('Solution saturated.')
``````

results in the following error message

``````WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
File "20210412_c_s_class.py", line 136, in <module>
if c_s < c_sat:
File "/usr/local/lib/python3.6/dist-packages/ufl/conditional.py", line 46, in __bool__
error("UFL conditions cannot be evaluated as bool in a Python context.")
File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.
``````

I tried to evaluate `c_s` by integrating it over a computational domain by replacing

``````c_s = My_c_s(t, u, degree=1)
``````

with

``````c_s_help = My_c_s(t, u, degree=1)
c_s = assemble(c_s_help*dx(domain=mesh))
``````

but this only resulted in an additional error message

``````WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Traceback (most recent call last):
File "20210412_c_s_class.py", line 130, in <module>
c_s = assemble(c_s_help*dx(domain=mesh))
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 213, in assemble
assembler.assemble(tensor, dolfin_form)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate expression.
*** Reason:  Missing eval() function (must be overloaded).
*** Where:   This error was encountered inside Expression.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
``````

Can somebody tell me how the class `My_c_s(UserExpression)` should be implemented the right way?

Thanks!

Do you need any further information?

We likely need less information. Most of us don’t have time nor willpower to read long posts to track down the problem. So if you can produce a minimal working example which reproduces your error, you’re more likely to get help. See here.

1 Like

You´re right, @nate!

Basically the question is, how to convert the calculated concentration of the solution `c_s` to a type format, that can be evaluated in a boolean way by `if c_s < c_sat`?

Here is my revised MWCE

``````from fenics import *

################ Euler's method to calculate the differential #################
################## of the concentration of the solution c_s ###################
# dc_s / dt = integral^{Gamma} (u - c_s) dA

class My_c_s(UserExpression):

def __init__(self, t, u, **kwargs):
super().__init__(**kwargs)
self.t = t
self.u = u

def euler_c_s(self, t0, c_s, dt, t, u):
temp = -0

# Iterating till the point at which we
# need approximation
while t0 < self.t:
temp = c_s

c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u))
t0 = t0 + dt
return c_s

def func_c_s(self, t, c_s, u):
return (self.u-c_s)*ds

def value_shape(self):
return ()

######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
def eval(self, value, x):

if x >= 0 and x <= 0.4:
value = u_avg_domain_t0

elif x >= 0.6 and x <= 1.0:
value = u_avg_domain_t0

else:
value = 0

def value_shape(self):
return ()

######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)

##################### define constants #########################################
t = 0   # time
t0 = 0  # t0 used by Euler's method
T = 2.0 # total time of the simulation
num_steps = 10
dt = T / num_steps

u_avg_domain_t0 = 1  # concentration of the corresponding domain used by My_u_t0
u_t0 = My_u_t0()
u = interpolate(u_t0, V)

c_s_0 = 0   # concentration of c_s for t=0
c_sat = 0.9 # saturation concentration of c_s
c_s = My_c_s(t, u, degree=1)

while t < T:
t = t+dt
c_s.t = t

if c_s < c_sat:
print('t = %.1f: c_s = %.9f' % (t, c_s))
else:
print('Solution saturated.')
``````

and the resulting error message

``````WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
File "20210426_c_s_class.py", line 80, in <module>
if c_s < c_sat:
File "/usr/local/lib/python3.6/dist-packages/ufl/conditional.py", line 46, in __bool__
error("UFL conditions cannot be evaluated as bool in a Python context.")
File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.
``````

I´m running FEniCS via Docker with the lastest stable version.

Thanks!

1 Like

`c_s` is a `My_c_s` which is an expression which is a field defined everywhere in your domain. `c_sat` is just a number. It doesn’t make sense to perform the boolean comparison `if c_s < c_sat`.

1 Like

Thanks for clarifying that, @nate!

The only reason I tried to convert the set of functions, which also calculate `c_s` using the Euler’s method, to a class, is because I know how to update `c_s` by `c_s.t = t` inside the for loop for every new time step.

Can somebody explain me how to update the set of functions for every time step after “calling” it initially by `c_s = euler_c_s(t0, c_s_0, dt, t, u)`?

That’s the revised MWCE, which only differs in using a set of functions to calculate `c_s` instead of the class `My_c_s`,

``````
################ Euler's method to calculate the differential #################
################## of the concentration of the solution c_s ###################
# dc_s / dt = integral^{Gamma} (u - c_s) dA

def euler_c_s( t0, c_s, dt, t, u):
temp = -0

# Iterating till the point at which we
# need approximation
while t0 < t:
temp = c_s

c_s = c_s + assemble(dt * func_c_s(t0, c_s, u))
t0 = t0 + dt
return c_s

def func_c_s(t, c_s, u):
return (u-c_s)*ds

def value_shape(self):
return ()

######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
def eval(self, value, x):

if x >= 0 and x <= 0.4:
value = u_avg_domain_t0

elif x >= 0.6 and x <= 1.0:
value = u_avg_domain_t0

else:
value = 0

def value_shape(self):
return ()

######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)

##################### define constants #########################################
t = 0   # time
t0 = 0  # t0 used by Euler's method
T = 2.0 # total time of the simulation
num_steps = 10
dt = T / num_steps

u_avg_domain_t0 = 1  # concentration of the corresponding domain used by My_u_t0
u_t0 = My_u_t0()
u = interpolate(u_t0, V)

c_s_0 = 0   # concentration of c_s for t=0
c_sat = 0.9 # saturation concentration of c_s
c_s = euler_c_s(t0, c_s_0, dt, t, u)
print('t = %.1f: c_s = %.9f' % (t, c_s))

while t < T:
t = t+dt
c_s.t = t
if c_s < c_sat:
print('t = %.1f: c_s = %.9f' % (t, c_s))
else:
print('Solution saturated.')
``````

which is creating to following error message

``````WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
t = 0.0: c_s = 0.000000000
Traceback (most recent call last):
File "20210426_c_s_func.py", line 72, in <module>
c_s.t = t
AttributeError: 'int' object has no attribute 't'
``````

Thanks!