reimplement sample as method of rFFTW_interpolator

bfps/cpp/rFFTW_interpolator.cpp

@@ -143,6 +143,36 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
     }
 }
 
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::sample(
+        const int nparticles,
+        const int pdimension,
+        const double t,
+        const double *__restrict__ x,
+        double *__restrict__ y,
+        const int *deriv)
+{
+    /* get grid coordinates */
+    int *xg = new int[3*nparticles];
+    double *xx = new double[3*nparticles];
+    double *yy =  new double[3*nparticles];
+    std::fill_n(yy, 3*nparticles, 0.0);
+    this->get_grid_coordinates(nparticles, pdimension, x, xg, xx);
+    /* perform interpolation */
+    for (int p=0; p<nparticles; p++)
+        this->operator()(t, xg + p*3, xx + p*3, yy + p*3, deriv);
+    MPI_Allreduce(
+            yy,
+            y,
+            3*nparticles,
+            MPI_DOUBLE,
+            MPI_SUM,
+            this->descriptor->comm);
+    delete[] yy;
+    delete[] xg;
+    delete[] xx;
+}
+
 template <class rnumber, int interp_neighbours>
 void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
         const double t,

bfps/cpp/rFFTW_interpolator.hpp

@@ -79,6 +79,7 @@ class rFFTW_interpolator
         /* interpolate field at an array of locations */
         void sample(
                 const int nparticles,
+                const int pdimension,
                 const double t,
                 const double *__restrict__ x,
                 double *__restrict__ y,

bfps/cpp/rFFTW_particles.cpp

@@ -154,7 +154,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double
     switch(particle_type)
     {
         case VELOCITY_TRACER:
-            this->sample_vec_field(this->vel, t, x, y);
+            this->vel->sample(this->nparticles, this->ncomponents, t, x, y);
             break;
     }
 }
@@ -255,29 +255,6 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
 }
 
 
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_grid_coordinates(double *x, int *xg, double *xx)
-{
-    static double grid_size[] = {this->dx, this->dy, this->dz};
-    double tval;
-    std::fill_n(xg, this->nparticles*3, 0);
-    std::fill_n(xx, this->nparticles*3, 0.0);
-    for (int p=0; p<this->nparticles; p++)
-    {
-        for (int c=0; c<3; c++)
-        {
-            tval = floor(x[p*this->ncomponents+c]/grid_size[c]);
-            xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
-            xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
-        }
-        /*xg[p*3+2] -= this->fs->rd->starts[0];
-        if (this->myrank == this->fs->rd->rank[0] &&
-            xg[p*3+2] > this->fs->rd->subsizes[0])
-            xg[p*3+2] -= this->fs->rd->sizes[0];*/
-    }
-}
-
 template <int particle_type, class rnumber, int interp_neighbours>
 void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data_file_id)
 {

bfps/cpp/rFFTW_particles.hpp

@@ -90,7 +90,6 @@ class rFFTW_particles
         void get_rhs(double *__restrict__ x, double *__restrict__ rhs);
         void get_rhs(double t, double *__restrict__ x, double *__restrict__ rhs);
 
-        void get_grid_coordinates(double *__restrict__ x, int *__restrict__ xg, double *__restrict__ xx);
         void sample_vec_field(
             rFFTW_interpolator<rnumber, interp_neighbours> *vec,
             double t,