Clean up kernel declarations

examples/src/example1.cpp

@@ -59,7 +59,7 @@ void nbody(vec&  src_coord, vec&  src_value,
 void fmm_test(size_t N, int mult_order, MPI_Comm comm){
 
   // Set kernel.
-  const pvfmm::Kernel<double>& kernel_fn=pvfmm::laplace_grad_d;
+  const pvfmm::Kernel<double>& kernel_fn=pvfmm::LaplaceKernel<double>::gradient();
 
   // Create target and source vectors.
   vec  trg_coord=point_distrib<double>(RandUnif,N,comm);

examples/src/example2.cpp

@@ -28,7 +28,7 @@ void fn_output(const double* coord, int n, double* out){
 void fmm_test(size_t N, int mult_order, int cheb_deg, double tol, MPI_Comm comm){
 
   // Set kernel.
-  const pvfmm::Kernel<double>& kernel_fn=pvfmm::laplace_potn_d;
+  const pvfmm::Kernel<double>& kernel_fn=pvfmm::LaplaceKernel<double>::gradient();
 
   // Construct tree.
   size_t max_pts=100;

examples/src/fmm_cheb.cpp

@@ -230,37 +230,35 @@ void fmm_test(int test_case, size_t N, size_t M, bool unif, int mult_order, int
       fn_input_=fn_input_t1<Real_t>;
       fn_poten_=fn_poten_t1<Real_t>;
       fn_grad_ =fn_grad_t1<Real_t>;
-      mykernel     =&pvfmm::LaplaceKernel<Real_t>::potn_ker();
-      //mykernel_grad=&pvfmm::LaplaceKernel<Real_t>::grad_ker();
+      mykernel     =&pvfmm::LaplaceKernel<Real_t>::potential();
+      //mykernel_grad=&pvfmm::LaplaceKernel<Real_t>::gradient();
       bndry=pvfmm::Periodic;
       break;
     case 2:
       fn_input_=fn_input_t2<Real_t>;
       fn_poten_=fn_poten_t2<Real_t>;
       fn_grad_ =fn_grad_t2<Real_t>;
-      mykernel     =&pvfmm::LaplaceKernel<Real_t>::potn_ker();
-      //mykernel_grad=&pvfmm::LaplaceKernel<Real_t>::grad_ker();
+      mykernel     =&pvfmm::LaplaceKernel<Real_t>::potential();
+      //mykernel_grad=&pvfmm::LaplaceKernel<Real_t>::gradient();
       bndry=pvfmm::FreeSpace;
       break;
     case 3:
       fn_input_=fn_input_t3<Real_t>;
       fn_poten_=fn_poten_t3<Real_t>;
-      mykernel     =&pvfmm::ker_stokes_vel;
-      //mykernel_grad=&pvfmm::ker_stokes_grad;
+      mykernel     =&pvfmm::StokesKernel<Real_t>::velocity();
+      //mykernel_grad=&pvfmm::StokesKernel<Real_t>::vel_grad();
       bndry=pvfmm::FreeSpace;
       break;
     case 4:
       fn_input_=fn_input_t4<Real_t>;
       fn_poten_=fn_poten_t4<Real_t>;
-      mykernel     =&pvfmm::ker_biot_savart;
-      //mykernel_grad=&pvfmm::ker_biot_savart_grad;
+      mykernel     =&pvfmm::BiotSavartKernel<Real_t>::potential();
       bndry=pvfmm::FreeSpace;
       break;
     case 5:
       fn_input_=fn_input_t5<Real_t>;
       fn_poten_=fn_poten_t5<Real_t>;
-      mykernel     =&pvfmm::ker_helmholtz;
-      //mykernel_grad=&pvfmm::ker_helmholtz_grad;
+      mykernel     =&pvfmm::HelmholtzKernel<Real_t>::potential();
       bndry=pvfmm::FreeSpace;
       break;
     default:
@@ -543,12 +541,13 @@ int main(int argc, char **argv){
   omp_set_num_threads( atoi(commandline_option(argc, argv,  "-omp",     "1", false, "-omp  <int> = (1)    : Number of OpenMP threads."          )));
   size_t   N=(size_t)strtod(commandline_option(argc, argv,    "-N",     "1",  true, "-N    <int>          : Number of point sources."           ),NULL);
   size_t   M=(size_t)strtod(commandline_option(argc, argv,    "-M",     "1", false, "-M    <int>          : Number of points per octant."       ),NULL);
-  bool  unif=              (commandline_option(argc, argv, "-unif",    NULL, false, "-unif                : Uniform point distribution."        )!=NULL);
+  bool  unif=   (1==strtoul(commandline_option(argc, argv, "-unif",     "0", false, "-unif <0/1> =  (0)   : Uniform point distribution."        ),NULL,10));
   int      m=       strtoul(commandline_option(argc, argv,    "-m",    "10", false, "-m    <int> = (10)   : Multipole order (+ve even integer)."),NULL,10);
   int      q=       strtoul(commandline_option(argc, argv,    "-q",    "14", false, "-q    <int> = (14)   : Chebyshev order (+ve integer)."     ),NULL,10);
   int      d=       strtoul(commandline_option(argc, argv,    "-d",    "15", false, "-d    <int> = (15)   : Maximum tree depth."                ),NULL,10);
   double tol=        strtod(commandline_option(argc, argv,  "-tol",  "1e-5", false, "-tol <real> = (1e-5) : Tolerance for adaptive refinement." ),NULL);
-  bool  adap=              (commandline_option(argc, argv, "-adap",    NULL, false, "-adap                : Adaptive tree refinement."          )!=NULL);
+  bool  adap=   (1==strtoul(commandline_option(argc, argv, "-adap",     "0", false, "-adap <0/1> =  (0)   : Adaptive tree refinement."          ),NULL,10));
+  bool    sp=   (1==strtoul(commandline_option(argc, argv,   "-sp",     "0", false, "-sp   <0/1> =  (0)   : Single Precision."                  ),NULL,10));
   int   test=       strtoul(commandline_option(argc, argv, "-test",     "1", false,
        "-test <int> = (1)    : 1) Laplace, Smooth Gaussian, Periodic Boundary\n\
                                2) Laplace, Discontinuous Sphere, FreeSpace Boundary\n\
@@ -560,7 +559,8 @@ int main(int argc, char **argv){
 
   // Run FMM with above options.
   pvfmm::Profile::Tic("FMM_Test",&comm,true);
-  fmm_test<double>(test, N,M,unif, m,q, d, adap,tol, comm);
+  if(sp) fmm_test<float >(test, N,M,unif, m,q, d, adap,tol, comm);
+  else   fmm_test<double>(test, N,M,unif, m,q, d, adap,tol, comm);
   pvfmm::Profile::Toc();
 
   //Output Profiling results.

include/kernel.hpp

@@ -156,132 +156,29 @@ Kernel<T> BuildKernel(const char* name, int dim, std::pair<int,int> k_dim,
 #endif
 namespace pvfmm{ // Predefined Kernel-functions
 
-////////////////////////////////////////////////////////////////////////////////
-////////                   LAPLACE KERNEL                               ////////
-////////////////////////////////////////////////////////////////////////////////
-
-/**
- * \brief Green's function for the Poisson's equation. Kernel tensor
- * dimension = 1x1.
- */
-template <class T, int newton_iter=0>
-void laplace_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-// Laplace double layer potential.
-template <class T, int newton_iter=0>
-void laplace_dbl_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-// Laplace grdient kernel.
-template <class T, int newton_iter=0>
-void laplace_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-
-//#ifdef PVFMM_QUAD_T
-//const Kernel<QuadReal_t> laplace_potn_q=BuildKernel<QuadReal_t, laplace_poten, laplace_dbl_poten>("laplace"     , 3, std::pair<int,int>(1,1));
-//const Kernel<QuadReal_t> laplace_grad_q=BuildKernel<QuadReal_t, laplace_grad                    >("laplace_grad", 3, std::pair<int,int>(1,3));
-//#endif
-
-const Kernel<double    > laplace_potn_d=BuildKernel<double    , laplace_poten<double,2>, laplace_dbl_poten<double,2> >("laplace"     , 3, std::pair<int,int>(1,1));
-const Kernel<double    > laplace_grad_d=BuildKernel<double    , laplace_grad <double,2>                              >("laplace_grad", 3, std::pair<int,int>(1,3),
-  &laplace_potn_d, &laplace_potn_d, NULL, &laplace_potn_d, &laplace_potn_d, NULL, &laplace_potn_d, NULL);
-
-const Kernel<float     > laplace_potn_f=BuildKernel<float     , laplace_poten<float,1>, laplace_dbl_poten<float,1> >("laplace"     , 3, std::pair<int,int>(1,1));
-const Kernel<float     > laplace_grad_f=BuildKernel<float     , laplace_grad <float,1>                             >("laplace_grad", 3, std::pair<int,int>(1,3),
-  &laplace_potn_f, &laplace_potn_f, NULL, &laplace_potn_f, &laplace_potn_f, NULL, &laplace_potn_f, NULL);
-
 template<class T>
 struct LaplaceKernel{
-  inline static const Kernel<T>& potn_ker();
-  inline static const Kernel<T>& grad_ker();
+  inline static const Kernel<T>& potential();
+  inline static const Kernel<T>& gradient();
 };
 
-//#ifdef PVFMM_QUAD_T
-//template<> const Kernel<QuadReal_t>& LaplaceKernel<QuadReal_t>::potn_ker(){ return laplace_potn_q; };
-//template<> const Kernel<QuadReal_t>& LaplaceKernel<QuadReal_t>::grad_ker(){ return laplace_grad_q; };
-//#endif
-
-template<> const Kernel<double>& LaplaceKernel<double>::potn_ker(){ return laplace_potn_d; };
-template<> const Kernel<double>& LaplaceKernel<double>::grad_ker(){ return laplace_grad_d; };
-
-template<> const Kernel<float>& LaplaceKernel<float>::potn_ker(){ return laplace_potn_f; };
-template<> const Kernel<float>& LaplaceKernel<float>::grad_ker(){ return laplace_grad_f; };
-
-////////////////////////////////////////////////////////////////////////////////
-////////                   STOKES KERNEL                             ////////
-////////////////////////////////////////////////////////////////////////////////
-
-/**
- * \brief Green's function for the Stokes's equation. Kernel tensor
- * dimension = 3x3.
- */
-template <class T>
-void stokes_vel(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-template <class T>
-void stokes_sym_dip(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-template <class T>
-void stokes_press(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-template <class T>
-void stokes_stress(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-template <class T>
-void stokes_grad(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-
-#ifndef __MIC__
-#ifdef USE_SSE
-template <>
-void stokes_vel<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr);
-
-template <>
-void stokes_press<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr);
-
-template <>
-void stokes_stress<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr);
-
-template <>
-void stokes_grad<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr);
-#endif
-#endif
-
-const Kernel<double> ker_stokes_vel   =BuildKernel<double, stokes_vel, stokes_sym_dip>("stokes_vel"   , 3, std::pair<int,int>(3,3));
-
-const Kernel<double> ker_stokes_press =BuildKernel<double, stokes_press              >("stokes_press" , 3, std::pair<int,int>(3,1));
-
-const Kernel<double> ker_stokes_stress=BuildKernel<double, stokes_stress             >("stokes_stress", 3, std::pair<int,int>(3,9));
-
-const Kernel<double> ker_stokes_grad  =BuildKernel<double, stokes_grad               >("stokes_grad"  , 3, std::pair<int,int>(3,9));
-
-////////////////////////////////////////////////////////////////////////////////
-////////                  BIOT-SAVART KERNEL                            ////////
-////////////////////////////////////////////////////////////////////////////////
-
-template <class T>
-void biot_savart(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-const Kernel<double> ker_biot_savart=BuildKernel<double, biot_savart>("biot_savart", 3, std::pair<int,int>(3,3));
-
-////////////////////////////////////////////////////////////////////////////////
-////////                   HELMHOLTZ KERNEL                             ////////
-////////////////////////////////////////////////////////////////////////////////
-
-/**
- * \brief Green's function for the Helmholtz's equation. Kernel tensor
- * dimension = 2x2.
- */
-template <class T>
-void helmholtz_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-template <class T>
-void helmholtz_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr);
-
-
+template<class T>
+struct StokesKernel{
+  inline static const Kernel<T>& velocity();
+  inline static const Kernel<T>& pressure();
+  inline static const Kernel<T>& stress  ();
+  inline static const Kernel<T>& vel_grad();
+};
 
-const Kernel<double> ker_helmholtz     =BuildKernel<double, helmholtz_poten>("helmholtz"     , 3, std::pair<int,int>(2,2));
+template<class T>
+struct BiotSavartKernel{
+  inline static const Kernel<T>& potential();
+};
 
-const Kernel<double> ker_helmholtz_grad=BuildKernel<double, helmholtz_grad >("helmholtz_grad", 3, std::pair<int,int>(2,6));
+template<class T>
+struct HelmholtzKernel{
+  inline static const Kernel<T>& potential();
+};
 
 }//end namespace
 #ifdef __INTEL_OFFLOAD

include/kernel.txx

@@ -866,7 +866,7 @@ void laplace_poten_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value,
   #undef SRC_BLK
 }
 
-template <class T, int newton_iter>
+template <class T, int newton_iter=0>
 void laplace_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
   #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
         generic_kernel<Real_t, 1, 1, laplace_poten_uKernel<Real_t,Vec_t, rinv_intrin##nwtn<Vec_t,Real_t> > > \
@@ -977,7 +977,7 @@ void laplace_dbl_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, M
   #undef SRC_BLK
 }
 
-template <class T, int newton_iter>
+template <class T, int newton_iter=0>
 void laplace_dbl_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
   #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
         generic_kernel<Real_t, 4, 1, laplace_dbl_uKernel<Real_t,Vec_t, rinv_intrin##nwtn<Vec_t,Real_t> > > \
@@ -1089,7 +1089,7 @@ void laplace_grad_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value,
   #undef SRC_BLK
 }
 
-template <class T, int newton_iter>
+template <class T, int newton_iter=0>
 void laplace_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
   #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
         generic_kernel<Real_t, 1, 3, laplace_grad_uKernel<Real_t,Vec_t, rinv_intrin##nwtn<Vec_t,Real_t> > > \
@@ -1134,6 +1134,31 @@ void laplace_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cn
 }
 
 
+template<class T> const Kernel<T>& LaplaceKernel<T>::potential(){
+  static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,1>, laplace_dbl_poten<T,1> >("laplace"     , 3, std::pair<int,int>(1,1));
+  return potn_ker;
+}
+template<class T> const Kernel<T>& LaplaceKernel<T>::gradient(){
+  static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,1>, laplace_dbl_poten<T,1> >("laplace"     , 3, std::pair<int,int>(1,1));
+  static Kernel<T> grad_ker=BuildKernel<T, laplace_grad <T,1>                         >("laplace_grad", 3, std::pair<int,int>(1,3),
+      &potn_ker, &potn_ker, NULL, &potn_ker, &potn_ker, NULL, &potn_ker, NULL);
+  return grad_ker;
+}
+
+template<> const Kernel<double>& LaplaceKernel<double>::potential(){
+  typedef double T;
+  static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,2>, laplace_dbl_poten<T,2> >("laplace"     , 3, std::pair<int,int>(1,1));
+  return potn_ker;
+}
+template<> const Kernel<double>& LaplaceKernel<double>::gradient(){
+  typedef double T;
+  static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,2>, laplace_dbl_poten<T,2> >("laplace"     , 3, std::pair<int,int>(1,1));
+  static Kernel<T> grad_ker=BuildKernel<T, laplace_grad <T,2>                         >("laplace_grad", 3, std::pair<int,int>(1,3),
+      &potn_ker, &potn_ker, NULL, &potn_ker, &potn_ker, NULL, &potn_ker, NULL);
+  return grad_ker;
+}
+
+
 ////////////////////////////////////////////////////////////////////////////////
 ////////                     STOKES KERNEL                              ////////
 ////////////////////////////////////////////////////////////////////////////////
@@ -2176,6 +2201,23 @@ void stokes_grad<double>(double* r_src, int src_cnt, double* v_src_, int dof, do
 #endif
 #endif
 
+template<class T> const Kernel<T>& StokesKernel<T>::velocity(){
+  static Kernel<T> ker=BuildKernel<T, stokes_vel, stokes_sym_dip>("stokes_vel"   , 3, std::pair<int,int>(3,3));
+  return ker;
+}
+template<class T> const Kernel<T>& StokesKernel<T>::pressure(){
+  static Kernel<T> ker=BuildKernel<T, stokes_press              >("stokes_press" , 3, std::pair<int,int>(3,1));
+  return ker;
+}
+template<class T> const Kernel<T>& StokesKernel<T>::stress(){
+  static Kernel<T> ker=BuildKernel<T, stokes_stress             >("stokes_stress", 3, std::pair<int,int>(3,9));
+  return ker;
+}
+template<class T> const Kernel<T>& StokesKernel<T>::vel_grad(){
+  static Kernel<T> ker=BuildKernel<T, stokes_grad               >("stokes_grad"  , 3, std::pair<int,int>(3,9));
+  return ker;
+}
+
 
 ////////////////////////////////////////////////////////////////////////////////
 ////////                  BIOT-SAVART KERNEL                            ////////
@@ -2217,6 +2259,11 @@ void biot_savart(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cn
   }
 }
 
+template<class T> const Kernel<T>& BiotSavartKernel<T>::potential(){
+  static Kernel<T> ker=BuildKernel<T, biot_savart>("biot_savart", 3, std::pair<int,int>(3,3));
+  return ker;
+}
+
 
 ////////////////////////////////////////////////////////////////////////////////
 ////////                   HELMHOLTZ KERNEL                             ////////
@@ -2244,7 +2291,7 @@ void helmholtz_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg
         if (R!=0){
           R = sqrt(R);
           T invR=1.0/R;
-          T G[2]={cos(mu*R)*invR, sin(mu*R)*invR};
+          T G[2]={(T)cos(mu*R)*invR, (T)sin(mu*R)*invR};
           p[0] += v_src[(s*dof+i)*2+0]*G[0] - v_src[(s*dof+i)*2+1]*G[1];
           p[1] += v_src[(s*dof+i)*2+0]*G[1] + v_src[(s*dof+i)*2+1]*G[0];
         }
@@ -2255,9 +2302,9 @@ void helmholtz_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg
   }
 }
 
-template <class T>
-void helmholtz_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr){
-  //TODO Implement this.
+template<class T> const Kernel<T>& HelmholtzKernel<T>::potential(){
+  static Kernel<T> ker=BuildKernel<T, helmholtz_poten>("helmholtz"     , 3, std::pair<int,int>(2,2));
+  return ker;
 }
 
 }//end namespace