add interpolation operator to interpolator

bfps/NavierStokes.py

@@ -234,14 +234,23 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         return None
     def add_particle_fields(
             self,
+            interp_type = 'spline',
             kcut = None,
             neighbours = 1,
+            smoothness = 1,
             name = 'particle_field'):
         self.particle_variables += ('interpolator<{0}, {1}> *vel_{2}, *acc_{2};\n' +
                                     '{0} *{2}_tmp;\n').format(self.C_dtype, neighbours, name)
-        self.particle_start += ('vel_{0} = new interpolator<{1}, {2}>(fs);\n' +
-                                'acc_{0} = new interpolator<{1}, {2}>(fs);\n' +
-                                '{0}_tmp = new {1}[acc_{0}->unbuffered_descriptor->local_size];\n').format(name, self.C_dtype, neighbours)
+        if interp_type == 'spline':
+            beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
+        elif interp_type == 'Lagrange':
+            beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
+        self.particle_start += ('vel_{0} = new interpolator<{1}, {2}>(fs, {3});\n' +
+                                'acc_{0} = new interpolator<{1}, {2}>(fs, {3});\n' +
+                                '{0}_tmp = new {1}[acc_{0}->unbuffered_descriptor->local_size];\n').format(name,
+                                                                                                           self.C_dtype,
+                                                                                                           neighbours,
+                                                                                                           beta_name)
         self.particle_end += ('delete vel_{0};\n' +
                               'delete acc_{0};\n' +
                               'delete[] {0}_tmp;\n').format(name)

bfps/cpp/interpolator.cpp

@@ -28,11 +28,13 @@
 
 template <class rnumber, int interp_neighbours>
 interpolator<rnumber, interp_neighbours>::interpolator(
-        fluid_solver_base<rnumber> *fs)
+        fluid_solver_base<rnumber> *fs,
+        base_polynomial_values BETA_POLYS)
 {
     int tdims[4];
     this->unbuffered_descriptor = fs->rd;
     this->buffer_size = (interp_neighbours+1)*this->unbuffered_descriptor->slice_size;
+    this->compute_beta = BETA_POLYS;
     tdims[0] = (interp_neighbours+1)*2*this->unbuffered_descriptor->nprocs + this->unbuffered_descriptor->sizes[0];
     tdims[1] = this->unbuffered_descriptor->sizes[1];
     tdims[2] = this->unbuffered_descriptor->sizes[2];
@@ -116,6 +118,39 @@ int interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
     return EXIT_SUCCESS;
 }
 
+template <class rnumber, int interp_neighbours>
+void interpolator<rnumber, interp_neighbours>::operator()(
+        int *xg,
+        double *xx,
+        double *dest,
+        int *deriv)
+{
+    double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
+    if (deriv == NULL)
+    {
+        this->compute_beta(0, xx[0], bx);
+        this->compute_beta(0, xx[1], by);
+        this->compute_beta(0, xx[2], bz);
+    }
+    else
+    {
+        this->compute_beta(deriv[0], xx[0], bx);
+        this->compute_beta(deriv[1], xx[1], by);
+        this->compute_beta(deriv[2], xx[2], bz);
+    }
+    std::fill_n(dest, 3, 0);
+    rnumber *field = this->f + this->buffer_size;
+    for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
+    for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
+    for (int ix = -interp_neighbours; ix <= interp_neighbours+1; ix++)
+        for (int c=0; c<3; c++)
+            dest[c] += field[((ptrdiff_t(    xg[2]+iz                             ) *this->descriptor->sizes[1] +
+                               ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1])))*this->descriptor->sizes[2] +
+                               ptrdiff_t(MOD(xg[0]+ix, this->descriptor->sizes[2])))*3+c]*(bz[iz+interp_neighbours]*
+                                                                                           by[iy+interp_neighbours]*
+                                                                                           bx[ix+interp_neighbours]);
+}
+
 template class interpolator<float, 1>;
 template class interpolator<float, 2>;
 template class interpolator<float, 3>;

bfps/cpp/interpolator.hpp

@@ -48,14 +48,17 @@ class interpolator
 {
     public:
         ptrdiff_t buffer_size;
+        base_polynomial_values compute_beta;
         field_descriptor<rnumber> *descriptor;
         field_descriptor<rnumber> *unbuffered_descriptor;
         rnumber *f;
 
         interpolator(
-                fluid_solver_base<rnumber> *FSOLVER);
+                fluid_solver_base<rnumber> *FSOLVER,
+                base_polynomial_values BETA_POLYS);
         ~interpolator();
 
+        void operator()(int *xg, double *xx, double *dest, int *deriv = NULL);
         /* destroys input */
         int read_rFFTW(void *src);
 };

tests/test_base.py

@@ -93,34 +93,28 @@ def launch(
     c.parameters['niter_part'] = 1
     c.parameters['famplitude'] = 0.2
     if c.parameters['nparticles'] > 0:
-        c.add_particle_fields(name = 'regular', neighbours = 6)
+        c.add_particle_fields(name = 'regular', neighbours = opt.neighbours, smoothness = opt.smoothness)
         c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
         c.add_particles(
                 kcut = 'fs->kM/2',
                 integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
                 fields_name = 'filtered')
-        for integr_steps in range(1, 7):
-            c.add_particles(
-                    integration_steps = 1,
-                    neighbours = opt.neighbours,
-                    smoothness = opt.smoothness,
-                    fields_name = 'regular')
-        for info in [(2, 1, 0, 'spline', 'AdamsBashforth'),
-                     (2, 1, 1, 'spline', 'AdamsBashforth'),
-                     (2, 6, 1, 'spline', 'AdamsBashforth'),
-                     (2, 6, 2, 'spline', 'AdamsBashforth'),
-                     (2, 6, 1, 'spline', 'Heun'),
-                     (4, 6, 1, 'spline', 'cRK4'),
-                     (2, 1, 1, 'Lagrange', 'AdamsBashforth'),
-                     (2, 6, 1, 'Lagrange', 'AdamsBashforth'),
-                     (2, 6, 1, 'Lagrange', 'Heun'),
-                     (4, 6, 1, 'Lagrange', 'cRK4')]:
+        #for integr_steps in range(1, 7):
+        #    c.add_particles(
+        #            integration_steps = integr_steps,
+        #            neighbours = opt.neighbours,
+        #            smoothness = opt.smoothness,
+        #            fields_name = 'regular')
+        for info in [(2, 1, 'spline', 'Heun'),
+                     (4, 1, 'spline', 'cRK4'),
+                     (2, 1, 'Lagrange', 'Heun'),
+                     (4, 1, 'Lagrange', 'cRK4')]:
             c.add_particles(
                     integration_steps = info[0],
-                    neighbours = info[1],
-                    smoothness = info[2],
-                    interp_type = info[3],
-                    integration_method = info[4],
+                    neighbours = opt.neighbours,
+                    smoothness = info[1],
+                    interp_type = info[2],
+                    integration_method = info[3],
                     fields_name = 'regular')
     c.fill_up_fluid_code()
     c.finalize_code()