I have a code that used to work with previous (2017) versions of FEniCS. Now, after a little adjustment because of code obsolescence, it returns the error

Generating mesh with CGAL 3D mesh generator

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

Solving linear variational problem.

Traceback (most recent call last):

File “MWE.py”, line 70, in

solve(a == L, u, bcs = bc0)

File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve

_solve_varproblem(*args, **kwargs)

File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 247, in _solve_varproblem

solver.solve()

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

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,

*** include aminimalrunning 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

*** -------------------------------------------------------------------------

The code is

```
from dolfin import *
import numpy as np
import mshr
cylinder = mshr.Cylinder(Point(0, 0, 0), Point(0, 0, -9), 8.0, 8.0)
domain = cylinder
mesh = mshr.generate_mesh(domain, 20)
V = FunctionSpace(mesh, 'CG', 1)
# Define boundary conditions
tol = 1E-14
class Top(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], 0.0, tol)
# Define subdomains.
# top
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= -6 - tol
# Other region.
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[2] <= -6 + tol
# Define the mesh
region = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
# Use the subdomains to mark the cells belonging to each subdomain
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(region, 0)
subdomain_1.mark(region, 1)
# k values of the layers
k_0 = 75
k_1 = 7
# Define the values of kappa in each subdomain
class K(UserExpression):
def __init__(self, region, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.region = region
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.region[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
kappa = K(region, k_0, k_1)
# Apply boundary conditions
u_D0 = Constant(10.0)
bc0 = DirichletBC(V, u_D0, Top())
# Define initial value
u_0 = Constant(5.0)
u_n = interpolate(u_0, V)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
F = u*v*dx + kappa*dot(grad(u), grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)
u = Function(V)
solve(a == L, u, bcs = bc0)
```

I suspect the error is due to the way the boundary conditions are implemented, which used to work but not anymore. I have no idea how to fix the problem.