Including Bessel functions inside a cpp string to use with Expression

1

Hi,

I’m struggeling to create an expression containing Bessel functions. I need something like this: f = Expression("std::cyl_bessel_i(1,x[0])", degree=2) but I get an compilation error that the Bessel function is not found.

I also tried to use the Boost lib version of the Bessel function with f = Expression("boost::math::cyl_bessel_i(1,x[0])", degree=2) but this function is also not found.

MWE:

from dolfin import *
f = Expression("std::cyl_bessel_i(1,x[0])", degree=2)
print(f(1.23))

Does someone has a tip for me? Thanks!

2

You can try to import it from ufl package:

from ufl import bessel_I

With the following signature:

Signature: bessel_I(nu, f)
Docstring: UFL operator: regular modified cylindrical Bessel function.
3

Yes that works. Can I use this Bessel function from the ufl package inside an Expression?

4

I don’t think so …
But you can mix Expression and ufl operator into one ufl variable and use it as an Expression:

import matplotlib.pyplot as plt
import fenics as fe
from ufl import bessel_I

mesh = fe.IntervalMesh(100, 0, 10)
V = fe.FunctionSpace(mesh,'P',1)

x = fe.SpatialCoordinate(mesh)
bs = bessel_I(0,x[0])
expr = fe.Expression('cos(x[0]*b)', b=10, element=V.ufl_element())

expr_bessel = expr*bs
proj_expr = fe.project(expr_bessel, V)

fe.plot(proj_expr)
plt.show()
5

Thank you, yes that works.
However, I can only use the project method for the expr_bessel and can not use the interpolate function. But the project returns strange results if an radial angle is uses, see: Strange behavior of project vs interpolate for C++/CPP vs UFL Expression (azimuth, radial or polar angle)

Is there a way to use UFL expression such as Constant(1.2345)*bessel_I(0,x[0]) inside the interpolate function?

6

Apparently from the link you post, project do not capture discontinuities if your space is continuous.
If the expression you project is continuous, there is normally no problem.

To my knowledge, if you really want to use interpolation you have to use an Expression with cpp code like in your first post.

7

Yes, I ended up using a CPP Expression to allow for an interpolation.