C++ Code Snippet for user expression error "Missing eval()"

General
1

I am trying write a custom user expression with a c++ code snippet by following the simple tutorial in tutorial.

Here is a simple working example using the poisson demo from poisson.
I am using dolphin version 2019.2.0.dev0.

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

materials = MeshFunction('double', mesh,2)

cppcode = """
class K : public Expression
{
public:

  void eval(Array<double>& values,
            const Array<double>& x,
            const ufc::cell& cell) const
  {
    values[0] = 1;
  }

};
"""

kappa = UserExpression(cppcode, degree=0)


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",degree=2)
g = Expression("sin(5*x[0])",degree=2)
a = kappa*inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u)

I get the following error, despite my best efforts.

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.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------

I have simplified the c++ snippet as much as possible. I also had to change the demo from Expression to UserExpression. For some reason it cannot find the Eval method. I have also tried the examples with the pybind11 signatures, but I get a different type error about Dolfin::Expression.

If anyone can point me in the right direction about simple c++ snippets, that would be an immense help. Eventually I need to evaluate a BC term that will need to include facet information. The expression will be very similar to FacetArea() with some slight modifications.

2

Hi jhcollins,

I don’t have a lot of experience with implementing the UserExpression through a C++ code snippet like that, but I tried with this Python syntax (which I’m more used to):

class Kappa(UserExpression):
    kappa_value = 0.001
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = self.kappa_value
    def value_shape(self):
        return ()

and then defining kappa = Kappa(mesh=mesh). Using this in the variational form seems to be working, I tried manipulating kappa_value and the results change according to expectations.

I didn’t really get what you mean with “the examples with pybind11 signatures”, is this what you meant you had tried? Or are you referring to something like this? Anyways, I’m pretty sure using the expressions through Python code like I did should provide identical functionality. Also, as a sidenote, note that UserExpression assumes the expression is a scalar value when not implementing value_shape() (just in case you will implement e.g. a vector function in the future).

If you need to evaluate a BC term with this type of expression, you can evaluate e.g. values along the normal direction of a boundary using the expression I provided by implementing eval_cell() instead of eval(). When overloading eval_cell(), you can get specific cell information like the cell normal facet during evaluation of the function. See this demo on mixed formulation for Poisson for an example of how to do that.

Here is your full code with the Python syntax expression instead of the C++ code snippet:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

materials = MeshFunction('double', mesh,2)

class Kappa(UserExpression):
    kappa_value = 0.001
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = self.kappa_value
    def value_shape(self):
        return ()
kappa = Kappa(mesh=mesh)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",degree=2)
g = Expression("sin(5*x[0])",degree=2)
a = kappa*inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u)

Hope this helps.

Best,
Halvor

3

If you insist on using C++ expressions, consider
https://bitbucket.org/fenics-project/dolfin/src/1c52e837eb54cc34627f90bde254be4aa8a2ae17/python/demo/documented/tensor-weighted-poisson/demo_tensor-weighted-poisson.py?at=master#demo_tensor-weighted-poisson.py-115
with the mwe:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

materials = MeshFunction('double', mesh,2)

cppcode = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;

#include <dolfin/function/Expression.h>

class K : public dolfin::Expression
{
public:

  K() : dolfin::Expression() {}

  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x,
            const ufc::cell& cell) const
  {
    values[0] = 1;
  }

};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<K, std::shared_ptr<K>, dolfin::Expression>
    (m, "K")
    .def(py::init<>());
}
"""

kappa = CompiledExpression(compile_cpp_code(cppcode).K(), degree=0)


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",degree=2)
g = Expression("sin(5*x[0])",degree=2)
a = kappa*inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u)