diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 3f745b7086dec81bd1bdd0e8d39e5368e8e43781..a06e97023329d468306f911b7d4dd39ec550a049 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -222,7 +222,7 @@ class ExpressionKernel(Kernel):
             raise ValueError("dist_vec length does not match expected dimension")
 
         expr = expr.subs([
-            ("d_%d" % i, dist_vec_i)
+            ("d%d" % i, dist_vec_i)
             for i, dist_vec_i in enumerate(dist_vec)
             ])
 
@@ -427,6 +427,84 @@ class StokesletKernel(ExpressionKernel):
 
     mapper_method = "map_stokeslet_kernel"
 
+
+class StressletKernel(ExpressionKernel):
+    init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name",
+        "stresslet_vector_name")
+
+    def __init__(self, dim=None, icomp=None, jcomp=None, viscosity_mu_name="mu",
+                        stresslet_vector_name="stresslet_vec"):
+        """
+        :arg viscosity_mu_name: The argument name to use for
+                dynamic viscosity :math:`\mu` the then generating functions to
+                evaluate this kernel.
+        """
+        mu = var(viscosity_mu_name)
+         #Remove mu?  (Not needed inside kernel?)
+
+        if dim == 2:
+            d = make_sym_vector("d", dim)
+            n = make_sym_vector(stresslet_vector_name, dim)
+            r = pymbolic_real_norm_2(d)
+            expr = (
+                sum(n[axis]*d[axis] for axis in range(dim))
+                * 
+                d[icomp]*d[jcomp]/r**4
+                )
+            scaling = 1/(var("pi"))
+
+        elif dim is None:
+            expr = None
+            scaling = None
+        else:
+            raise RuntimeError("unsupported dimensionality")
+
+        self.viscosity_mu_name = viscosity_mu_name
+        self.stresslet_vector_name = stresslet_vector_name
+        self.icomp = icomp
+        self.jcomp = jcomp
+
+        ExpressionKernel.__init__(
+                self,
+                dim,
+                expression=expr,
+                scaling=scaling,
+                is_complex_valued=False)
+
+    def __getinitargs__(self):
+        return (self._dim, self.icomp, self.jcomp, self.viscosity_mu_name,
+                      self.stresslet_vector_name)
+
+    def update_persistent_hash(self, key_hash, key_builder):
+        key_hash.update(type(self).__name__.encode())
+        key_builder.rec(key_hash, (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name, self.stresslet_vector_name))
+
+    def __repr__(self):
+        return "StressletKnl%dD_%d%d[%s]" % (self.dim, self.icomp, self.jcomp,
+                self.stresslet_vector_name)
+
+    def get_args(self):
+        return [
+                KernelArgument(
+                    loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64),
+                    ),
+                KernelArgument(
+                        lp.GlobalArg(self.stresslet_vector_name,
+                        None,
+                        shape=(self.dim, "nsources"),
+                        dim_tags="sep,C"),
+                    )
+                    ]
+
+    def transform_to_code(self, expr):
+        from sumpy.codegen import VectorComponentRewriter
+        vcr = VectorComponentRewriter([self.stresslet_vector_name])
+        from pymbolic.primitives import Variable
+        return _VectorIndexAdder(self.stresslet_vector_name, (Variable("isrc"),))(
+                vcr(expr))
+
+    mapper_method = "map_stresslet_kernel"
+
 # }}}
 
 
@@ -638,6 +716,8 @@ class KernelIdentityMapper(KernelMapper):
 
     map_laplace_kernel = map_expression_kernel
     map_helmholtz_kernel = map_expression_kernel
+    map_stokeslet_kernel = map_expression_kernel
+    map_stresslet_kernel = map_expression_kernel
 
     def map_axis_target_derivative(self, kernel):
         return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel))
@@ -669,6 +749,8 @@ class DerivativeCounter(KernelCombineMapper):
 
     map_laplace_kernel = map_expression_kernel
     map_helmholtz_kernel = map_expression_kernel
+    map_stokeslet_kernel = map_expression_kernel
+    map_stresslet_kernel = map_expression_kernel
 
     def map_axis_target_derivative(self, kernel):
         return 1 + self.rec(kernel.inner_kernel)
@@ -707,6 +789,29 @@ class KernelDimensionSetter(KernelIdentityMapper):
                 helmholtz_k_name=kernel.helmholtz_k_name,
                 allow_evanescent=kernel.allow_evanescent)
 
+    def map_stokeslet_kernel(self, kernel):
+        if kernel._dim is not None and self.dim != kernel.dim:
+            raise RuntimeError("cannot set kernel dimension to new value (%d)"
+                    "different from existing one (%d)"
+                    % (self.dim, kernel.dim))
+
+        return StressletKernel(self.dim,
+                kernel.icomp,
+                kernel.jcomp,
+                viscosity_mu_name=kernel.viscosity_mu_name)
+
+    def map_stresslet_kernel(self, kernel):
+        if kernel._dim is not None and self.dim != kernel.dim:
+            raise RuntimeError("cannot set kernel dimension to new value (%d)"
+                    "different from existing one (%d)"
+                    % (self.dim, kernel.dim))
+
+        return StressletKernel(self.dim,
+                kernel.icomp,
+                kernel.jcomp,
+                viscosity_mu_name=kernel.viscosity_mu_name,
+                stresslet_vector_name=kernel.stresslet_vector_name) 
+
 # }}}
 
 
diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py
index b2ab466dc477b4503f4b8f428ac82952b54041be..338ced1491702b11d9411a01ad673f27db4b5ddb 100644
--- a/sumpy/symbolic.py
+++ b/sumpy/symbolic.py
@@ -31,6 +31,7 @@ import sympy as sp
 import numpy as np
 from pymbolic.mapper import IdentityMapper as IdentityMapperBase
 from pymbolic.sympy_interface import PymbolicToSympyMapper
+import pymbolic.primitives as prim
 
 import logging
 logger = logging.getLogger(__name__)
@@ -210,4 +211,11 @@ class PymbolicToSympyMapperWithSymbols(PymbolicToSympyMapper):
         else:
             return PymbolicToSympyMapper.map_variable(self, expr)
 
+    def map_subscript(self, expr):
+        if isinstance(expr.aggregate, prim.Variable) and isinstance(expr.index, int):
+            return sp.Symbol("%s%d" % (expr.aggregate.name, expr.index))
+        else:
+            raise RuntimeError("do not know how to translate '%s' to sympy"
+                    % expr)
+
 # vim: fdm=marker