# How to use bessel functions in Expression

I am using legacy FEniCS 2019.1.0 and I need to find the eigenvalues and eigenfunctions of a Matern covariance operator, which is defined using modified bessel functions. However, the following mwe, defining an expression using bessel functions, does not compile:

``` from fenics import * import numpy as np mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20) V = FunctionSpace(mesh, ‘P’, 1) u = TrialFunction(V) v = TestFunction(V) k = Expression(‘pow(sqrt(pow(x-x_i,2)+pow(x-y_i,2)),nu)*cyl_bessel_k(nu, sqrt(pow(x-x_i,2)+pow(x-y_i,2)))’, degree=1, nu=1, x_i=0, y_i=0) c = np.zeros((mesh.num_vertices(), mesh.num_vertices())) coords = mesh.coordinates() dof2v = dof_to_vertex_map(V) for i in range(len(coords)): k.x_i = Vertex(mesh, dof2v[i]).midpoint() k.y_i = Vertex(mesh, dof2v[i]).midpoint() c[i, :] = assemble(k * v * dx) vals, vecs = np.linalg.eig(c) ```

The error message is the following:
/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp: In member function ‘virtual void dolfin::dolfin_expression_4a2c2db2129efa596eb787cd10ed666d::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:
/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp:64:69: error: ‘cyl_bessel_k’ was not declared in this scope

Apparently cyl_bessel_k is not declared even though it is listed in the documentation of cmath: Standard library header <cmath> - cppreference.com

I have also looked into this (old) solution but the code provided there does not compile either, giving an extremely long error message, mostly repeating error: ‘namespace’ definition is not allowed here https://fenicsproject.org/qa/3498/use-bessel-functions-in-jit-compiled-expressions/

Is there a different way of approaching this problem? Or am I just using the wrong alias for bessel functions?

I would use the ufl formulation of bessel functions directly, instead of trying to generate code for it.
i.e.

``````from fenics import *

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)

x = SpatialCoordinate(mesh)
nu = 1
x_i, y_i = 0,0
k = sqrt((x-x_i)**2+(x-y_i)**2)**nu*bessel_K(nu, sqrt((x-x_i)**2+(x-y_i)**2))
print(assemble(k*dx))
``````

This gives you a more exact approximation of your function as it is evaluated at the quadrature points, and not interpolated into a Lagrange space of order `degree` (1 in your case).

``` from fenics import * import numpy as np import time mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20) V = FunctionSpace(mesh,‘P’,1) u = TrialFunction(V) v = TestFunction(V) x = SpatialCoordinate(mesh) def materncov(x, nu, x_i, y_i): k = sqrt((x-x_i)**2+(x-y_i)**2)**nu*bessel_K(nu, sqrt((x-x_i)**2+(x-y_i)**2)) return k c = np.zeros((mesh.num_vertices(), mesh.num_vertices())) coords = mesh.coordinates() dof2v = dof_to_vertex_map(V) for i in range(len(coords)): starttime = time.time() x_i = Vertex(mesh, dof2v[i]).midpoint() y_i = Vertex(mesh, dof2v[i]).midpoint() c[i, :] = assemble(materncov(x, 1, x_i, y_i) * v * dx) print(time.time() - starttime) vals, vecs = np.linalg.eig(c) ```

While that does work for one iteration, in order to find the eigenfunctions I need to compare each pair of points in the mesh, so x_i,y_i have to loop over all vertices. If I have to use assemble in each iteration of this loop, this works quickly for the first few iterations (taking ~3ms/iteration), but after a handful of iterations this code starts calling the JIT compiler in every iteration, taking 2+ secs per iteration. Hence my original plan of using Expressions, since then I could just update the parameters x_i,y_i and not have to call the compiler every time. Is there a fix to this?

You should use `Constant` for `x_i` and `y_i` to avoid recompilation:

``````from fenics import *
import time
import numpy as np
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)
V = FunctionSpace(mesh, "Lagrange", 1)
x = SpatialCoordinate(mesh)
nu = 1
x_i, y_i = Constant(0),Constant(0)
k = sqrt((x-x_i)**2+(x-y_i)**2)**nu*bessel_K(nu, sqrt((x-x_i)**2+(x-y_i)**2))
v = TestFunction(V)
res = k*v*dx

c = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for i in range(len(coords)):
starttime = time.time()
mp =  Vertex(mesh, dof2v[i]).midpoint()
x_i.assign(mp)
y_i.assign(mp)
c[i, :] = assemble(res)
print(time.time() - starttime)
``````

Reduces the runtime of assembly to ~0.003 seconds per vertex

EDIT
If you tested this with the previous approach, I would strongly suggest to clean your cache with `dijitso clean` to avoid having tons of pre-compiled forms to search through.