JIT Expression - FEniCS Q&A

The modified nonlinear Poisson code below contains a nonlinear term that is implemented in two ways: once using UFL, as in the demo (hence supporting automated differentiation and fast evaluation); and then also using an expression.

from dolfin import *

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u_k = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])")

# define F and J
method = 2
if method == 1: # as in the demo-code
    F = inner((1 + u_k**2)*grad(u_k), grad(v))*dx - f*v*dx
    J = derivative(F, u_k, du)

elif method == 2: # using custom expressions, doing the same thing
    class CustomNonlinearity(Expression):
        def eval(self, value, x):
            value[0] = 1 + u_k(x)**2
    q = CustomNonlinearity()

    class DCustomNonlinearity(Expression):
        def eval(self, value, x):
            value[0] = 2 * u_k(x)
    Dq = DCustomNonlinearity()

    F = inner(q*grad(u_k), grad(v))*dx - f*v*dx
    J = inner(q*grad(du), grad(v)) * dx + inner(Dq * du * grad(u_k), grad(v)) * dx

# Compute solution
problem = NonlinearVariationalProblem(F, u_k, bc, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

# Plot solution and solution gradient
plot(u_k, title="Solution")
plot(grad(u_k), title="Solution gradient")
interactive()

In my code the nonlinear term is slightly more complicated, so option 1 is not possible. Now I would like to compile this expression to speed things up. But how can I pass the Function u_k (or its underlying C++ object) as an argument to the compiler?

I would also need an external C function in this JIT-expression. How can I tell the compiler where to look for this? Here is some code that might clarify a little bit:

from dolfin import *
code = '''
**line below --> how to tell compiler where to look for this?**
#include "external_C_function.h" # defines a function `f`

namespace dolfin {
    class MyFun : public Expression
    {
        public:
            MyFun(): Expression() {};
        void eval(Array<double>& values, const Array<double>& x) const {
            **line below --> how to access the function u_k?**
            argument = u_k(x)    --
            values[0] = f(argument);
        }
    };
}'''

f=Expression(code)

Many thanks