Use Bessel functions in JIT-compiled expressions

Question

At least, there are two ways to use Bessel functions in an dolfin expression with Python: (a) by subclassing and then using scipy.special Bessel functions, and (b) by using UFL functions (for instance bessel_j). Both methods lead to inefficient procedures if the functions has to be evaluated or projected in a fine mesh.

As alternative, I'm considering to use a JIT-compiled expressions involving the c++ functions cyl_bessel_*. However, I couldn’t compile it successfully (in Fenics 1.3.0). The minimal example to reproduce the compiling error is:

from dolfin import *
code = '''
class MyFun : public Expression
{
public:
    MyFun(): Expression() {};
void eval(Array<double>& values, const Array<double>& x) const {
    values[0] = cyl_bessel_j(0,x[0]);
}};'''

f=Expression(code)
f(1.0)

Is it possible to call Bessel functions in a JIT-compiled code? In that case, is it correct the use of cyl_bessel_* functions in the code?

Answer 1

You could do this by including a header file and a using-statement:

from dolfin import *
code = '''
#include <boost/math/special_functions/bessel.hpp>
using boost::math::cyl_bessel_i;
using boost::math::cyl_bessel_j;
using boost::math::cyl_bessel_k;
using boost::math::cyl_neumann;

namespace dolfin {
    class MyFun : public Expression
    {
        public:
            MyFun(): Expression() {};
        void eval(Array<double>& values, const Array<double>& x) const {
            values[0] = cyl_bessel_j(0,x[0]);
        }
    };
}'''

f=Expression(code)
print f(1.0)

Comments

That's solved my question. Thank you very much!

---

If I create a UnitSquareMesh and then have some function say f=Expression("ax[0] + bx[1]").
Then how do I use the above code to calculate cyl_bessel_j(f) ?
I am using fenics and python for 1 month only.

---

From the above code, try to modify it like this:

from dolfin import *
code = '''
#include <boost/math/special_functions/bessel.hpp>
using boost::math::cyl_bessel_i;
using boost::math::cyl_bessel_j;
using boost::math::cyl_bessel_k;
using boost::math::cyl_neumann;

namespace dolfin {
    class MyFun : public Expression
    {
        double a,b;
        public:
            MyFun(): Expression() {};
        void eval(Array<double>& values, const Array<double>& x) const {
            double f = a*x[0]+b*x[1];
            values[0] = cyl_bessel_j(0,f);
        }

        void update(double _a, double _b) {
            a = _a;
            b = _b;
        }
    };
}'''

f=Expression(code)
f.update(2.0, 3.0)
print f((0.1,0.5))

---

Thank you Øyvind Evju for the answer! I could finally finish the code.
I want to highlight a small issue I faced. I used the following function as "f" and I could not
use the normal log() in that. log10(), logl(), logf() all would work but not the simple log function. I also included math.h but that did not work. I must have been wrong somewhere.

      double f = sqrt( logf(x[0]/x0)*logf(x[0]/x0) + 
                       logf(x[1]/y0)*logf(x[1]/y0) )/c0 ;

The whole code for reference:

  from dolfin import *
  import numpy
  from numpy import sqrt
  from math import log

  # System size definition 
  # system = (x1, x1+Lx) x (y1, y1+Ly)
  Lx = 2.0; Ly = 2.0
  x1 = 10.0; y1 = 10.0

  c0 = 1.0
  D0 = 1.0
  x0 = 1.5
  y0 = 1.5

  # Create mesh and define function space
  dx0 = dx1 = 0.1
  nx = int(Lx/dx0) 
  ny = int(Ly/dx1)
  mesh = RectangleMesh(Point(x1, y1), Point(x1+Lx, y1+Ly), nx, ny,"right/left")
  V = FunctionSpace(mesh, 'Lagrange', 1)

  code = '''
  #include <math.h>
  #include <boost/math/special_functions/bessel.hpp>
  using boost::math::cyl_bessel_i;
  using boost::math::cyl_bessel_j;
  using boost::math::cyl_bessel_k;
  using boost::math::cyl_neumann;

  namespace dolfin {
       class MyFun : public Expression
      {
              double c0,D0,x0,y0;
              public:
                        MyFun(): Expression() {};
              void eval(Array<double>& values, const Array<double>& x) const {
                      double f = D0*sqrt( logf(x[0]/x0)*logf(x[0]/x0) + 
                                       logf(x[1]/y0)*logf(x[1]/y0) )/c0 ;
                      values[0] = cyl_bessel_k(0,f);
              }

              void update(double _c0, double _D0, double _x0, double _y0) {
                    c0 = _c0;
                    D0 = _D0;
                    x0 = _x0;
                    y0 = _y0;
            }
      };
}'''

  f=Expression(code)
  f.update(c0,D0,x0,y0)

  ss = interpolate(f,V)
  plot(ss, interactive=True)