Hello. I keep encountering a weird bug that I cannot explain or resolve. I am trying to model Stokes flow in a rectangular domain with a circular cutout. The boundary condition is Dirichlet, and inferred from experimental data that needs to be interpolated. I do this by constructing a UserExpression that contains two SciPy RectBivariateSpline objects that return the interpolated velocity components at any point in the domain. It is necessary to perform cubic 2D interpolation of the data because in my use case, the position of the inner circular cut out is not known a priori, and when post processing the data I need to calculate derivatives so the boundary data needs to be smooth.

The problem I encounter is that dolfin 2019.1.0 randomly crashes with the “Missing eval() function (must be overloaded)” error when I do this. It most often (but not always) does so after restarting the kernel, and often also on repeated execution. Sometimes the problem resolves itself when I just try to re-run it, but often it doesn’t and I can’t reliably use it.

Based on THIS topic I suspect it has something to do with the way interpolation is handled, but I’m not really sure how to work around it given the constraints mentioned above.

I’d be really grateful for an explanation of the error and suggestions how to resolve it! MWE below, with an example data file to execute it HERE: EDIT see comment below for updated MWE without the need to download a file.

EDIT: condensed the MWE slightly. Also to be clear, the error occurs in the line ‘solve(a == L, w, bcs)’:

EDIT2: condensed the MWE even further, the error also occurs if the markers of the inner and outer boundary are not distinct.

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

```
#remove boundary distinctions
from dolfin import *
import mshr
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp
import pandas as pd
'''Classes'''
class Velocity_from_PIV(UserExpression):
# Fenics object to evaluate/interpolate PIV velocity data anywhere
def __init__(self, data, **kwargs):
super().__init__(kwargs)
self.x=data['x'].drop_duplicates().values #vector of x coords (ascending)
self.y=data['y'].drop_duplicates().values #vector of y coords (ascending)
self.dx=data['dx'].values.reshape(self.x.size,self.y.size) #matrix of values
self.dy=data['dy'].values.reshape(self.x.size,self.y.size)
#Splines (default: cubic, no smoothing)
self.splx = interp.RectBivariateSpline(self.x,self.y,self.dx)
self.sply = interp.RectBivariateSpline(self.x,self.y,self.dy)
def eval_cell(self, value, x, ufc_cell):
value[0] = self.splx.ev(x[0],x[1])
value[1] = self.sply.ev(x[0],x[1])
def value_shape(self):
return (2,)
'''Methods'''
def main():
#load data
data = pd.read_csv('Data.txt',sep='\t',names=['i','j','x','y','dx','dy'])
#Calculate mesh
mesh, boundary_markers = meshing()
plot(mesh)
#Define function space and boundary conditions
W, bcs = defineBCs(mesh,boundary_markers,data)
#Perform FEM
u, p = calc_fem(W,bcs)
return u, p
def meshing():
"""Create the mesh and determine inner and outer boundaries"""
# For sake of MWE define rectangular dimensions manually:
xmin, xmax, ymin, ymax, PIVres = np.array([0.0, 0.7424643186232591, 0.0, 0.39774874211960304, 0.02651658280797354])
# x,y coords and radius of circular cutout:
centrex, centrey, radius = np.array([0.27,0.163,0.11])
#define interior and exterior boundary
outerb = mshr.Rectangle(Point(xmin,ymin),Point(xmax,ymax))
innerb = mshr.Circle(Point(centrex,centrey),radius)
mesh = mshr.generate_mesh(outerb-innerb,1/PIVres)
#label the boundary elements
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_markers.set_all(0)
return mesh, boundary_markers
def defineBCs(mesh,boundary_markers,data):
'''define FEM space, BCs'''
#Create fenics UserExpression to access data
data_fenics = Velocity_from_PIV(data)
# Define function space, choose Taylor-Hood elements
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 3)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
#Define Dirichlet BC on each part of the boundary
bcs = DirichletBC(W.sub(0), data_fenics, boundary_markers,0)
return W, bcs
def calc_fem(W,bcs):
"""obtain u and p from Stokes FEM"""
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
#Stokes equations in weak form
a = (inner(grad(u), grad(v)) - div(v)*p-div(u)*q -(1E-10)*p*q)*dx
L = -Constant(0)*q*dx
# Compute solution
w = Function(W)
solve(a == L, w, bcs)
#extract velocity and pressure fields
(u, p) = w.split()
return u, p
u, p = main()
```