diff --git a/source/pysph/base/kernels.h b/source/pysph/base/kernels.h
index 41b0b4b..ef5f97d 100644
--- a/source/pysph/base/kernels.h
+++ b/source/pysph/base/kernels.h
@@ -1,7 +1,23 @@
 // C declarations for the SPH smoothing kernels
+#ifndef CL_KERNEL_H
+#define CL_KERNEL_H
 
 #define PI 2.0f*acos(0.0f)
 
+enum KernelType {
+    CUBIC_SPLINE = 1,
+    GAUSSIAN = 2,
+    QUINTIC_SPLINE = 3,
+    WENDLAND_QUINTIC_SPLINE = 4,
+    HARMONIC = 5,
+    M6_SPLINE = 6,
+    W8 = 7,
+    W10 = 8,
+    REPULSIVE = 9,
+    POLY6 = 10
+};
+
+
 REAL cubic_spline_fac(unsigned int dim, REAL h)
 {
   REAL fac = 0.0;
@@ -72,3 +88,37 @@ void cubic_spline_gradient(REAL4 pa, REAL4 pb, REAL4* grad,
   grad[0].s2 = val * rab.s2;
 
 }
+
+REAL kernel_function(REAL4 pa, REAL4 pb, unsigned int dim, int kernel_type)
+{
+    REAL w;
+
+    switch (kernel_type) {
+        case CUBIC_SPLINE:
+            w = cubic_spline_function(pa, pb, dim);
+            break;
+        // FIXME: implement the rest!
+        default:
+            w = cubic_spline_function(pa, pb, dim);
+            break;
+    } // switch (kernel_type).
+
+    return w;
+}
+
+void kernel_gradient(REAL4 pa, REAL4 pb, REAL4* grad, 
+                     unsigned int dim, int kernel_type)
+{
+    switch (kernel_type) {
+        case CUBIC_SPLINE:
+            cubic_spline_gradient(pa, pb, grad, dim);
+            break;
+        // FIXME: implement the rest!
+        default:
+            cubic_spline_gradient(pa, pb, grad, dim);
+            break;
+    } // switch (kernel_type).
+}
+
+#endif // CL_KERNEL_H
+
diff --git a/source/pysph/base/kernels.pyx b/source/pysph/base/kernels.pyx
index debb0f7..f6290b3 100755
--- a/source/pysph/base/kernels.pyx
+++ b/source/pysph/base/kernels.pyx
@@ -25,6 +25,25 @@ cdef:
     double SQRT_1_PI = 1.0/sqrt(PI)
     double infty = numpy.inf
 
+
+##############################################################################
+# These are simple ints to define the different types of kernels.  This
+# is what is returned by the `KernelBase.get_type` method. The only
+# reason we do this is to unify the numbers used between Cython and
+# OpenCL.
+cdef:
+    int CUBIC_SPLINE = 1
+    int GAUSSIAN = 2
+    int QUINTIC_SPLINE = 3
+    int WENDLAND_QUINTIC_SPLINE = 4
+    int HARMONIC = 5
+    int M6_SPLINE = 6
+    int W8 = 7
+    int W10 = 8
+    int REPULSIVE = 9
+    int POLY6 = 10
+
+
 cdef inline double h_dim(double h, int dim):
     if dim == 1:
         return 1/h
@@ -300,6 +319,9 @@ cdef class Poly6Kernel(KernelBase):
     cpdef double radius(self):
         return 1.0
 
+    def get_type(self):
+        return POLY6 
+
 ##############################################################################
 # DummyKernel` class.
 ##############################################################################
@@ -424,7 +446,8 @@ cdef class CubicSplineKernel(KernelBase):
         return 2
 
     def get_type(self):
-        return 1
+        return CUBIC_SPLINE
+
 ##############################################################################
 
 
@@ -526,6 +549,8 @@ cdef class QuinticSplineKernel(KernelBase):
     cpdef double radius(self):
         return 3
 
+    def get_type(self):
+        return QUINTIC_SPLINE 
 
 ##############################################################################
 #`WendlandQuinticSplineKernel`
@@ -605,6 +630,9 @@ cdef class WendlandQuinticSplineKernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return WENDLAND_QUINTIC_SPLINE 
+
 ##############################################################################
 #`HarmonicKernel`
 ##############################################################################
@@ -697,6 +725,10 @@ cdef class HarmonicKernel(KernelBase):
 
     cpdef double radius(self):
         return 2
+
+    def get_type(self):
+        return HARMONIC
+
 #############################################################################
 
 ##############################################################################
@@ -769,6 +801,10 @@ cdef class M6SplineKernel(KernelBase):
 
     cpdef double radius(self):
         return 2
+
+    def get_type(self):
+        return M6_SPLINE
+
 #############################################################################
 
 ##############################################################################
@@ -826,6 +862,9 @@ cdef class GaussianKernel(KernelBase):
 
     cpdef double radius(self):
         return 3.0
+
+    def get_type(self):
+        return GAUSSIAN 
 ##############################################################################
 
 
@@ -903,6 +942,9 @@ cdef class W8Kernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return W8 
+
 ##############################################################################
 
 
@@ -978,6 +1020,9 @@ cdef class W10Kernel(KernelBase):
     cpdef double radius(self):
         return 2
 
+    def get_type(self):
+        return W10 
+
 
 ################################################################################
 # `RepulsiveBoundaryKernel` class.
@@ -1056,5 +1101,7 @@ cdef class RepulsiveBoundaryKernel(KernelBase):
             raise ValueError
         else: return (1.0/h)
 
+    def get_type(self):
+        return REPULSIVE 
 
 
diff --git a/source/pysph/solver/cl_common.h b/source/pysph/solver/cl_common.h
index a106eb5..2823604 100644
--- a/source/pysph/solver/cl_common.h
+++ b/source/pysph/solver/cl_common.h
@@ -58,4 +58,4 @@ unsigned int get_gid(unsigned int work_dim)
 
 }
 
-#endif CL_COMMON
+#endif // CL_COMMON
diff --git a/source/pysph/sph/funcs/density_funcs.cl b/source/pysph/sph/funcs/density_funcs.cl
index aa51bcb..95d705a 100644
--- a/source/pysph/sph/funcs/density_funcs.cl
+++ b/source/pysph/sph/funcs/density_funcs.cl
@@ -24,19 +24,7 @@ __kernel void SPHRho(int const kernel_type, int const dim, int const nbrs,
     {
       REAL4 pb = (REAL4)( s_x[i], s_y[i], s_z[i], s_h[i] );
       REAL mb = s_m[i];
-      
-      switch (kernel_type)
-	{
-	case 1:
-	  w = cubic_spline_function(pa, pb, dim);
-	  break;
-	  
-	default:
-	  w = cubic_spline_function(pa, pb, dim);
-	  break;
-	      
-	} // switch
-	  
+      w = kernel_function(pa, pb, dim, kernel_type); 
       wmb += w*mb;
   	  	  
     } // for i