Dear community,

In order to calculate a concentration at the borders of the domain, named c_s, by the Euler’s method, I wrote the following code

```
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):
return ((u_prev-c_s)*ds)
def euler_c_s( t0, c_s, dt, t, u_prev ):
temp = -0
# Iterating till the point at which we
# need approximation
while t0 < t:
temp = c_s
c_s = c_s + dt * func_c_s(t0, c_s, u_prev)
t0 = t0 + dt
return c_s
def value_shape(self):
return ()
```

, which is run by

```
while t < T:
t = t+dt
c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
```

,whereas `c_s_0`

is the intial value for t=0 of c_s and described by

```
class My_c_s_t0(UserExpression):
def eval(self, value, x):
# boundaries
if x[0] == 0:
value[0] = 0.2
elif x[0] == 1:
value[0] = 0.2
elif x[1] == 0:
value[0] = 0.2
elif x[1] == 1:
value[0] = 0.2
# domain without boundaries
else:
value[0] = 0
def value_shape(self):
return ()
```

```
c_s_t0 = My_c_s_t0()
c_s_0 = interpolate(c_s_t0, V)
```

and `u_prev`

, which is usually determined by `solve(a == L, u)`

, but here, in order to reduce the complexity of the MWCE given by

```
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 ()
```

```
u_t0 = My_u_t0()
u_prev = interpolate(u_t0, V)
```

When running the MWCE, I get the following error code

```
File "20201217_c_s.py", line 150, in <module>
c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
File "20201217_c_s.py", line 43, in euler_c_s
c_s = c_s + dt * func_c_s(t0, c_s, u_prev)
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'
```

**I´m unsure, if the error message is a result of a wrong implementation of the Euler´s method** (Before, I only used Euler´s method for problems with scalar initial values) **or** `func_c_s`

**can’t handle the matrix of** `c_s`

**in its current form?**

Thanks a lot for your time,

Cou

## Full MWCE

```
from fenics import *
import numpy as np
############### define concentration c_s (on boundaries) for t=0 ##############
class My_c_s_t0(UserExpression):
def eval(self, value, x):
# boundaries
if x[0] == 0:
value[0] = 0.2
elif x[0] == 1:
value[0] = 0.2
elif x[1] == 0:
value[0] = 0.2
elif x[1] == 1:
value[0] = 0.2
# domain without boundaries
else:
value[0] = 0
def value_shape(self):
return ()
############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):
return ((u_prev-c_s)*ds)
def euler_c_s( t0, c_s, dt, t, u_prev ):
temp = -0
# Iterating till the point at which we
# need approximation
while t0 < t:
temp = c_s
c_s = c_s + dt * func_c_s(t0, c_s, u_prev)
t0 = t0 + dt
return c_s
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)
################################ 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()
c_s_t0 = My_c_s_t0()
u_prev = interpolate(u_t0, V)
c_s_0 = interpolate(c_s_t0, V)
c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))
while t < T:
t = t+dt
c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))
```