Problem with boundary conditions: Missing eval() function (must be overloaded)

1

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

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.

2

Note that the eval_cell definition is not property indented, so it is just interpreted as a free function, and not a method of K. If that is indented, the error goes away. The warning about value_shape can be eliminated by adding the method

    def value_shape(self):
        return ()
1 Like