Compiled Expression, expected an overload 'eval' or 'eval_data' error

Question

Hello,
still struggling to speed up my very non linear source term (for advection-diffusion).

I had this implementation of an user defined expression:

class SprayGridExpression(Expression):

    def set_spray_grid(self, spray_grid):
        self.spray_grid = spray_grid

    def eval(self, value, x):
        X = x[0]
        Y = x[1]

        value[0] = self.spray_grid(X, Y)  # Volumetric mass flow flux

where spray_grid is quite a complex object (actually a "functor").

I decided to reimplement all the class hierarchy in C++ to create spray_grid in C++ (spray_grid is basically a sum of three gaussian distributions that are then placed all over the domain specifying the location of them with the vector<double> x0s, y0s).
The Expression seems to compile correctly, but then I get the error:

Calling DOLFIN just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "test_dolfin.py", line 148, in <module>
    spray_exp = Expression(ccp_expression)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dolfin/functions/expression.py", line 666, in __new__
    return object.__new__(create_compiled_expression_class(cpp_base))
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dolfin/functions/expression.py", line 187, in create_compiled_expression_class
    {"__init__": __init__})
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dolfin/functions/expression.py", line 357, in __new__
    raise TypeError("expected an overload 'eval' or 'eval_cell' method")
TypeError: expected an overload 'eval' or 'eval_cell' method

This is the section in the code when using the cpp Expression:

with open('classes.hpp', 'r') as fh:
    ccp_expression = fh.read()

print(ccp_expression)

spray_exp = Expression(ccp_expression)

---

classes.hpp

#include <iostream>
#include <cmath>
#include <vector>


using namespace std; 

/* Linspace like function as in MATLAB or Python */
template <typename T>
vector<T> linspace(T a, T b, size_t N) {
      T h = (b - a) / static_cast<T>(N-1);
        vector<T> xs(N);
          typename vector<T>::iterator x;
            T val;
              for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
                      *x = val;
                return xs;
}

// Simple 1D Gaussian Function
class RadialGaussian { 
    protected: 
        double _sigma2, _mu;

        double exponent(double x) const{
            return this->exponent_fun(_sigma2, _mu, x);

        }
        double denominator() const {
            return this->denominator_fun(_sigma2);
        }

    public: 
        // Constructor
        RadialGaussian(double sigma2, double mu) : 
            _sigma2(sigma2), 
            _mu(mu)
        { }


        static double exponent_fun(double sigma2, double mu, double x) {
            return -pow(x - mu, 2)/(2*sigma2);
        }

        static double denominator_fun(double sigma2) {
            return sqrt(2*M_PI*sigma2);
        }

        // Acting as a functor
        double operator() (double x) {
            return 1/this->denominator()*exp(this->exponent(x));
        }

        // Accessors
        double getSigma2 () { return _sigma2; }
        double getMu () { return _mu; }
};

// 3D Gaussian "Hat"
class CartesianGaussian3D : public RadialGaussian {
    public:
        CartesianGaussian3D(double sigma2, double mu) : 
            RadialGaussian(sigma2, mu)
        { }

        CartesianGaussian3D(RadialGaussian gaussian) : 
            RadialGaussian(gaussian.getSigma2(), gaussian.getMu())
        { }

        double operator() (double x, double y) const {
            double r = sqrt(pow(x, 2) + pow(y, 2));
            return 1/this->denominator()*exp(this->exponent(r));
        }
};

// The Space function that represents a single nozzle
class SprayModel {
    protected: 
        vector<double> _sigma2s, _mus;

    public:
        SprayModel(vector<double> sigma2s, vector<double> mus) :
            _sigma2s(sigma2s), 
            _mus(mus)
        { }

        double operator() (double x, double y) const {
            double r = sqrt(pow(x, 2) + pow(y, 2));
            double phi;

            for(auto sigma2=_sigma2s.begin(), mu = _mus.begin(); sigma2 != _mus.end(); ++sigma2, ++mu) {
                // boost::tie(sigma2, mu) = tup;

                phi += 1.0/CartesianGaussian3D::denominator_fun(*sigma2)*exp(CartesianGaussian3D::exponent_fun(*sigma2, *mu, r));
            }

            return phi;
        }
};

// Just summing up all the nozzles at the different x0s,y0s locations
class SprayGrid {
    protected:
        SprayModel _spray;
        vector<double> _x0s, _y0s;

    public: 
        SprayGrid(SprayModel spray, vector<double> x0s, vector<double> y0s) :
            _spray(spray),
            _x0s(x0s),
            _y0s(y0s)
        { }

        double operator() (double x, double y) const {
            double phi;

            for(auto x0: _x0s){
                for(auto y0: _y0s) {
                    phi += _spray(x-x0, y-y0);
                }
            }

            return phi;
        }


};


// DOLFIN expression
class SprayGridExpression : public Expression
{
    vector<double> sigma2s = {0.0022592878436543907, 0.018770022717831417, 0.001473404233081996}, 
        mus = {0.003827877912948519, 0.042287499332578256, 0.09907663376990554};

    SprayGrid spray_grid = SprayGrid(SprayModel(sigma2s, mus), sigma2s, mus);

    public:

    SprayGridExpression() : Expression()
    {
    }

    void eval(Array<double>& values, const Array<double>& x) const{
        double xx = x[0], yy=x[1];

        values[0] = spray_grid(xx, yy);

    }
};