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

General
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