Evaluate UserExpression with eval_cell

1

Hi all,

I’m trying to evaluate a UserExpression subclass which uses a MeshFunction. It works great when assembling, but I can’t evaluate it at a specific point. Here’s the MWE and error message:


from fenics import *


class CoeffClass(UserExpression):
    def __init__(self, mesh, vm, **kwargs):
        super().__init__(kwargs)
        self.mesh = mesh
        self.vm = vm

    def eval_cell(self, value, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        if self.vm[cell] == 1:
            value[0] = 2
        else:
            value[0] = 3

    def value_shape(self):
        return ()


nx = ny = 2
mesh = UnitSquareMesh(nx, ny)
n = FacetNormal(mesh)
V = FunctionSpace(mesh, 'P', 1)
D0 = FunctionSpace(mesh, 'DG', 0)
volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 1)

D = CoeffClass(mesh=mesh, vm=volume_markers, degree=1)
# Doesn't work
# print(D(0.5, 0.5))
u = Expression("2*x[0]*x[0] + t", degree=2, t=1)
u = interpolate(u, V)
print(assemble(D*dot(grad(u), n)*ds))

ERROR:
Traceback (most recent call last):
File “test_function.py”, line 71, in
print(D(0.5, 0.5))
File “/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py”, line 191, in call
self._cpp_object.eval(values, x)
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
*** -------------------------------------------------------------------------

Thanks for the help !

2

Hi Rem,

as a workaround you can interpolate it in a function which you can evaluate:

D = CoeffClass(...)
d = interpolate(D, V)
print(d(0.5, 0.5))

BR

2 Likes
3

Hi again,
I found that workaround as well cheers.