adds required methods

cpp/kspace.cpp

@@ -626,6 +626,77 @@ void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restri
     }
 }
 
+template <field_backend be,
+          kspace_dealias_type dt>
+template <typename rnumber>
+void kspace<be, dt>::project_divfree(
+        typename fftw_interface<rnumber>::complex *__restrict__ xa,
+        typename fftw_interface<rnumber>::complex *__restrict__ ya,
+        typename fftw_interface<rnumber>::complex *__restrict__ za,
+        const bool maintain_energy)
+{
+    TIMEZONE("kspace::project_divfree");
+    this->CLOOP(
+                [&](const ptrdiff_t cindex,
+                    const ptrdiff_t xindex,
+                    const ptrdiff_t yindex,
+                    const ptrdiff_t zindex,
+                    const double k2){
+                if (k2 > 0)
+        {
+                    const typename fftw_interface<rnumber>::complex tval = {
+                     rnumber((this->kx[xindex]*((*(xa + cindex))[0]) +
+                              this->ky[yindex]*((*(ya + cindex))[0]) +
+                              this->kz[zindex]*((*(za + cindex))[0]) ) / k2),
+                     rnumber((this->kx[xindex]*((*(xa + cindex))[1]) +
+                              this->ky[yindex]*((*(ya + cindex))[1]) +
+                              this->kz[zindex]*((*(za + cindex))[1]) ) / k2)};
+                    double initial_size = 1;
+                    double projected_size = 1;
+                    if (maintain_energy)
+                    {
+                        initial_size = std::sqrt(
+                                ((*(xa + cindex))[0])*((*(xa + cindex))[0]) +
+                                ((*(ya + cindex))[0])*((*(ya + cindex))[0]) +
+                                ((*(za + cindex))[0])*((*(za + cindex))[0]) +
+                                ((*(xa + cindex))[1])*((*(xa + cindex))[1]) +
+                                ((*(ya + cindex))[1])*((*(ya + cindex))[1]) +
+                                ((*(za + cindex))[1])*((*(za + cindex))[1]));
+                    }
+                    for (int imag_part=0; imag_part<2; imag_part++)
+                    {
+                        xa[cindex][imag_part] -= tval[imag_part]*this->kx[xindex];
+                        ya[cindex][imag_part] -= tval[imag_part]*this->ky[yindex];
+                        za[cindex][imag_part] -= tval[imag_part]*this->kz[zindex];
+                    }
+                    if (maintain_energy)
+                    {
+                        projected_size = std::sqrt(
+                                ((*(xa + cindex))[0])*((*(xa + cindex))[0]) +
+                                ((*(ya + cindex))[0])*((*(ya + cindex))[0]) +
+                                ((*(za + cindex))[0])*((*(za + cindex))[0]) +
+                                ((*(xa + cindex))[1])*((*(xa + cindex))[1]) +
+                                ((*(ya + cindex))[1])*((*(ya + cindex))[1]) +
+                                ((*(za + cindex))[1])*((*(za + cindex))[1]));
+                        if (projected_size > 0)
+                        #pragma omp simd
+                            for (int imag_part=0; imag_part<2; imag_part++)
+                            {
+                                (*(xa + cindex))[imag_part] *= initial_size / projected_size;
+                                (*(ya + cindex))[imag_part] *= initial_size / projected_size;
+                                (*(za + cindex))[imag_part] *= initial_size / projected_size;
+                            }
+                    }
+           }
+        }
+    );
+    if (this->layout->myrank == this->layout->rank[0][0]) {
+        std::fill_n((rnumber*)(xa), 2, 0.0);
+        std::fill_n((rnumber*)(ya), 2, 0.0);
+        std::fill_n((rnumber*)(za), 2, 0.0);
+    }
+}
+
 template <field_backend be,
           kspace_dealias_type dt>
 template <typename rnumber>
@@ -1552,6 +1623,17 @@ template void kspace<FFTW, SMOOTH>::project_divfree<double>(
        typename fftw_interface<double>::complex *__restrict__ a,
        const bool);
 
+template void kspace<FFTW, SMOOTH>::project_divfree<float>(
+       typename fftw_interface<float>::complex *__restrict__ xa,
+       typename fftw_interface<float>::complex *__restrict__ ya,
+       typename fftw_interface<float>::complex *__restrict__ za,
+       const bool);
+template void kspace<FFTW, SMOOTH>::project_divfree<double>(
+       typename fftw_interface<double>::complex *__restrict__ xa,
+       typename fftw_interface<double>::complex *__restrict__ ya,
+       typename fftw_interface<double>::complex *__restrict__ za,
+       const bool);
+
 template void kspace<FFTW, SMOOTH>::force_divfree<float>(
        typename fftw_interface<float>::complex *__restrict__ a);
 template void kspace<FFTW, SMOOTH>::force_divfree<double>(

cpp/kspace.hpp

@@ -289,6 +289,12 @@ class kspace
         void project_divfree(
                 typename fftw_interface<rnumber>::complex *__restrict__ a,
                 const bool maintain_energy = false);
+        template <typename rnumber>
+        void project_divfree(
+                typename fftw_interface<rnumber>::complex *__restrict__ xa,
+                typename fftw_interface<rnumber>::complex *__restrict__ ya,
+                typename fftw_interface<rnumber>::complex *__restrict__ za,
+                const bool maintain_energy = false);
         // TODO: can the following be done in a cleaner way?
         template <typename rnumber>
         void force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a){