# 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] >= 0 and x[0] <= 0.4:
value[0] = u_avg_domain_t0

elif x[0] >= 0.6 and x[0] <= 1.0:
value[0] = u_avg_domain_t0

else:
value[0] = 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.0)

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

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

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

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

class Obstacle_right(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (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] >= 0 and x[0] <= 0.4:
value[0] = u_avg_domain_t0

elif x[0] >= 0.6 and x[0] <= 1.0:
value[0] = u_avg_domain_t0

else:
value[0] = 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] >= 0 and x[0] <= 0.4:
value[0] = u_avg_domain_t0

elif x[0] >= 0.6 and x[0] <= 1.0:
value[0] = u_avg_domain_t0

else:
value[0] = 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!