From a88e4c1942f3eb73da49fcbf444f27979520fff6 Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <clalesc1@jhu.edu>
Date: Thu, 24 Dec 2015 20:39:18 +0100
Subject: [PATCH] move get_grid_coords to interpolator

---
 bfps/cpp/rFFTW_interpolator.cpp | 31 +++++++++++++++++++++++++++----
 bfps/cpp/rFFTW_interpolator.hpp | 10 ++++++----
 bfps/cpp/rFFTW_particles.cpp    |  2 +-
 done.txt                        |  1 +
 todo.txt                        |  1 -
 5 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index 5e95d742..46673eda 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -121,12 +121,35 @@ int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
 }
 
 template <class rnumber, int interp_neighbours>
-void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
-        double t,
+void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
+        const int nparticles,
+        const int pdimension,
+        const double *x,
         int *xg,
-        double *xx,
+        double *xx)
+{
+    static double grid_size[] = {this->dx, this->dy, this->dz};
+    double tval;
+    std::fill_n(xg, nparticles*3, 0);
+    std::fill_n(xx, nparticles*3, 0.0);
+    for (int p=0; p<nparticles; p++)
+    {
+        for (int c=0; c<3; c++)
+        {
+            tval = floor(x[p*pdimension+c]/grid_size[c]);
+            xg[p*3+c] = MOD(int(tval), this->descriptor->sizes[2-c]);
+            xx[p*3+c] = (x[p*pdimension+c] - tval*grid_size[c]) / grid_size[c];
+        }
+    }
+}
+
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
+        const double t,
+        const int *xg,
+        const double *xx,
         double *dest,
-        int *deriv)
+        const int *deriv)
 {
     double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
     if (deriv == NULL)
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 0e9be896..c7efb9d8 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -71,22 +71,24 @@ class rFFTW_interpolator
 
         /* map real locations to grid coordinates */
         void get_grid_coordinates(
+                const int nparticles,
+                const int pdimension,
                 const double *__restrict__ x,
-                const int *__restrict__ xg,
-                const double *__restrict__ xx);
+                int *__restrict__ xg,
+                double *__restrict__ xx);
         /* interpolate field at an array of locations */
         void sample(
                 const int nparticles,
                 const double t,
                 const double *__restrict__ x,
-                const double *__restrict__ y,
+                double *__restrict__ y,
                 const int *deriv = NULL);
         /* interpolate 1 point */
         void operator()(
                 const double t,
                 const int *__restrict__ xg,
                 const double *__restrict__ xx,
-                const double *__restrict__ dest,
+                double *__restrict__ dest,
                 const int *deriv = NULL);
         int read_rFFTW(void *src);
 };
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 42a308f2..4d0ab85b 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -132,7 +132,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_fiel
     double *xx = new double[3*this->nparticles];
     double *yy =  new double[3*this->nparticles];
     std::fill_n(yy, 3*this->nparticles, 0.0);
-    this->get_grid_coordinates(x, xg, xx);
+    vec->get_grid_coordinates(this->nparticles, this->ncomponents, x, xg, xx);
     /* perform interpolation */
     for (int p=0; p<this->nparticles; p++)
         (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
diff --git a/done.txt b/done.txt
index 78739201..8dd974c1 100644
--- a/done.txt
+++ b/done.txt
@@ -1,2 +1,3 @@
 x 2015-12-04 make code py3 compatible                                @python3
 x 2015-12-23 decide on versioning system                                           +merge0
+x 2015-12-24 move get_grid coords to interpolator                    @optimization +v1.0
diff --git a/todo.txt b/todo.txt
index 71c2677e..dcc58366 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,7 +1,6 @@
 (B) FFTW interpolator doesn't need its own field            @optimization +v1.0
 (B) compute z polynomials only when needed                  @optimization +v1.0
 (B) get rid of temporal interpolation                       @optimization +v1.0
-(B) move get_grid coords to interpolator                    @optimization +v1.0
 (B) use less memory                                         @optimization
 (C) clean up machine_settings mess                          @design @documentation +v2.0
 (C) code overview                                           @documentation
-- 
GitLab