diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index df9c0468..e528e63e 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -35,9 +35,10 @@
     cast,
     overload,
 )
+from warnings import warn
 
 import numpy as np
-from typing_extensions import Self, override
+from typing_extensions import override
 
 import loopy as lp
 import pymbolic.primitives as prim
@@ -107,6 +108,18 @@
     :show-inheritance:
     :members: mapper_method
 
+.. [Pozrikidis1992] C. Pozrikidis,
+    *Boundary Integral and Singularity Methods for Linearized Viscous Flow*,
+    Cambridge University Press, 1992.
+
+.. [Hsiao2008] G. C. Hsiao, W. L. Wendland,
+    *Boundary Integral Equations*,
+    Springer, 2008.
+
+.. [Kress2013] R. Kress,
+    *Linear Integral Equations*,
+    Springer Science & Business Media, 2013.
+
 Derivatives
 -----------
 
@@ -483,6 +496,13 @@ def is_complex_valued(self) -> bool:
 # {{{ PDE kernels
 
 class LaplaceKernel(ExpressionKernel):
+    r"""A kernel for the Laplace equation (see e.g. Theorem 6.2 from [Kress2013]_).
+
+    .. math::
+
+        \Delta K(\mathbf{x}, \mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}).
+    """
+
     mapper_method: ClassVar[str] = "map_laplace_kernel"
 
     def __init__(self, dim: int) -> None:
@@ -531,6 +551,13 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 
 
 class BiharmonicKernel(ExpressionKernel):
+    r"""A kernel for the biharmonic equation.
+
+    .. math::
+
+        \Delta^2 K(\mathbf{x}, \mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}).
+    """
+
     mapper_method: ClassVar[str] = "map_biharmonic_kernel"
 
     def __init__(self, dim: int) -> None:
@@ -585,7 +612,13 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 
 @dataclass(frozen=True, repr=False)
 class HelmholtzKernel(ExpressionKernel):
-    """
+    r"""A kernel for the Helmholtz equation (see e.g. Example 12.14 in [Kress2013]_).
+
+    .. math::
+
+        \Delta K(\mathbf{x}, \mathbf{y}) + k^2 K(\mathbf{x}, \mathbf{y})
+            = \delta(\mathbf{x} - \mathbf{y}).
+
     .. autoattribute:: helmholtz_k_name
     .. autoattribute:: allow_evanescent
     """
@@ -673,7 +706,13 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 
 @dataclass(frozen=True, repr=False)
 class YukawaKernel(ExpressionKernel):
-    """
+    r"""A kernel for the Yukawa equation.
+
+    .. math::
+
+        \Delta K(\mathbf{x}, \mathbf{y}) - \lambda^2 K(\mathbf{x}, \mathbf{y})
+            = \delta(\mathbf{x} - \mathbf{y}).
+
     .. autoattribute:: yukawa_lambda_name
     """
 
@@ -764,7 +803,15 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 
 @dataclass(frozen=True, repr=False)
 class ElasticityKernel(ExpressionKernel):
-    """
+    r"""A kernel for the linear elasticity (Navier-Cauchy) equations
+    (see e.g. Section 2.2 in [Hsiao2008]_).
+
+    .. math::
+
+        \mu \Delta K_{ij}(\mathbf{x}, \mathbf{y})
+        + \frac{\mu}{1 - 2 \nu} \nabla (\nabla \cdot K_{ij}(\mathbf{x}, \mathbf{y}))
+        = \delta_{ij} \delta(\mathbf{x} - \mathbf{y}).
+
     .. autoattribute:: icomp
     .. autoattribute:: jcomp
     .. autoattribute:: viscosity_mu
@@ -779,59 +826,54 @@ class ElasticityKernel(ExpressionKernel):
 
     viscosity_mu: float | SpatialConstant
     r"""The argument name to use for the dynamic viscosity :math:`\mu` when
-    generating functions to evaluate this kernel. Can also be a numeric value.
+    generating functions to evaluate this kernel.
     """
     poisson_ratio: float | SpatialConstant
     r"""The argument name to use for Poisson's ratio :math:`\nu` when generating
-    functions to evaluate this kernel. Can also be a numeric value.
+    functions to evaluate this kernel.
     """
 
-    def __new__(cls,
-                dim: int, icomp: int, jcomp: int,
-                viscosity_mu: float | str | SpatialConstant = "mu",
-                poisson_ratio: float | str | SpatialConstant = "nu",
-        ) -> Self:
-        if poisson_ratio == 0.5:  # noqa: RUF069
-            return super().__new__(StokesletKernel)
-        else:
-            return super().__new__(cls)
-
     def __init__(self,
-                 dim: int, icomp: int, jcomp: int,
+                 dim: int,
+                 icomp: int,
+                 jcomp: int,
                  viscosity_mu: float | str | SpatialConstant = "mu",
                  poisson_ratio: float | str | SpatialConstant = "nu") -> None:
         if isinstance(viscosity_mu, str):
             mu = SpatialConstant(viscosity_mu)
         else:
+            warn("Passing a non-str 'viscosity_mu' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
             mu = viscosity_mu
 
         if isinstance(poisson_ratio, str):
             nu = SpatialConstant(poisson_ratio)
         else:
+            warn("Passing a non-str 'poisson_ratio' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
+            if poisson_ratio == 0.5:  # noqa: RUF069
+                raise ValueError(
+                    f"{type(self)} cannot accept 'poisson_ratio' of 0.5 "
+                    "(use StokesletKernel instead)")
+
             nu = poisson_ratio
 
+        d = make_sym_vector("d", dim)
+        r = pymbolic_real_norm_2(d)
+        delta_ij = 1 if icomp == jcomp else 0
+
         if dim == 2:
-            d = make_sym_vector("d", dim)
-            r = pymbolic_real_norm_2(d)
             # See (Berger and Karageorghis 2001)
-            expr = (
-                -var("log")(r)*((3 - 4 * nu) if icomp == jcomp else 0)
-                +
-                d[icomp]*d[jcomp]/r**2
-                )
+            expr = -var("log")(r)*(3 - 4 * nu)*delta_ij + d[icomp]*d[jcomp]/r**2
             scaling = -1/(8*var("pi")*(1 - nu)*mu)
-
         elif dim == 3:
-            d = make_sym_vector("d", dim)
-            r = pymbolic_real_norm_2(d)
             # Kelvin solution
-            expr = (
-                (1/r)*((3 - 4*nu) if icomp == jcomp else 0)
-                +
-                d[icomp]*d[jcomp]/r**3
-                )
+            expr = (1/r)*(3 - 4*nu)*delta_ij + d[icomp]*d[jcomp]/r**3
             scaling = -1/(16*var("pi")*(1 - nu)*mu)
-
         else:
             raise NotImplementedError(f"unsupported dimension: '{dim}'")
 
@@ -842,23 +884,33 @@ def __init__(self,
         object.__setattr__(self, "viscosity_mu", mu)
         object.__setattr__(self, "poisson_ratio", nu)
 
+    @override
+    def __str__(self) -> str:
+        return (
+            f"ElasticityKnl{self.dim}D_{self.icomp}{self.jcomp}"
+            f"({self.viscosity_mu}, {self.poisson_ratio})")
+
     @override
     def __reduce__(self):
+        # TODO: remove this once we stop allowing non-str values for viscosity_mu
+        viscosity_mu = (
+            self.viscosity_mu.name
+            if isinstance(self.viscosity_mu, SpatialConstant)
+            else self.viscosity_mu)
+        poisson_ratio = (
+            self.poisson_ratio.name
+            if isinstance(self.poisson_ratio, SpatialConstant)
+            else self.poisson_ratio)
+
         return (
             type(self),
-            (self.dim, self.icomp, self.jcomp, self.viscosity_mu, self.poisson_ratio))
+            (self.dim, self.icomp, self.jcomp, viscosity_mu, poisson_ratio))
 
     @property
     @override
     def is_complex_valued(self) -> bool:
         return False
 
-    @override
-    def __str__(self) -> str:
-        return (
-            f"ElasticityKnl{self.dim}D_{self.icomp}{self.jcomp}"
-            f"({self.viscosity_mu}, {self.poisson_ratio})")
-
     @memoize_method
     @override
     def get_args(self) -> Sequence[KernelArgument]:
@@ -892,48 +944,137 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 
 
 @dataclass(frozen=True, repr=False)
-class StokesletKernel(ElasticityKernel):
-    """
+class StokesletKernel(ExpressionKernel):
+    r"""A kernel for the Stokes equations (see e.g. Chapter 2 in [Pozrikidis1992]_).
+
+    .. math::
+
+        \begin{cases}
+        -\mu \Delta K_{ij}(\mathbf{x}, \mathbf{y})
+            + \nabla_i P_j(\mathbf{x}, \mathbf{y})
+            = \delta_{ij} \delta(\mathbf{x} - \mathbf{y}), \\
+        \nabla_i K_{ij}(\mathbf{x}, \mathbf{y}) = 0, \\
+        \end{cases}
+
+    where pressure kernel :math:`P_j = \partial_j K` is the derivative of the
+    Laplace kernel. This kernel is often called the Stokeslet or the Oseen-Burgers
+    tensor and it represents the velocity field.
+
     .. autoattribute:: icomp
     .. autoattribute:: jcomp
     .. autoattribute:: viscosity_mu
     """
 
-    def __new__(cls,
-                dim: int,
-                icomp: int,
-                jcomp: int,
-                viscosity_mu: float | str | SpatialConstant = "mu",
-                poisson_ratio: float | str | SpatialConstant | None = None,
-            ) -> Self:
-        return object.__new__(cls)
+    mapper_method: ClassVar[str] = "map_stokeslet_kernel"
+
+    icomp: int
+    """Component index for the kernel."""
+    jcomp: int
+    """Component index for the kernel."""
+
+    viscosity_mu: float | SpatialConstant
+    r"""The argument name to use for the dynamic viscosity :math:`\mu` when
+    generating functions to evaluate this kernel.
+    """
 
     def __init__(self,
-                dim: int,
-                icomp: int,
-                jcomp: int,
-                viscosity_mu: float | str | SpatialConstant = "mu",
-                poisson_ratio: float | str | SpatialConstant | None = None) -> None:
-        if poisson_ratio is None:
-            poisson_ratio = 0.5
-
-        if poisson_ratio != 0.5:  # noqa: RUF069
-            raise ValueError(
-                "'StokesletKernel' must have a Poisson ratio of 0.5: "
-                f"got '{poisson_ratio}'")
+                 dim: int,
+                 icomp: int,
+                 jcomp: int,
+                 viscosity_mu: float | str | SpatialConstant = "mu",
+                 poisson_ratio: float | str | SpatialConstant | None = None) -> None:
+        if poisson_ratio is not None:
+            warn(f"Passing 'poisson_ratio' to {type(self)} is deprecated. This "
+                 "argument will be removed in Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
+            if poisson_ratio != 0.5:  # noqa: RUF069
+                raise ValueError("'poisson_ratio' must be 0.5 (if provided)")
+
+        if isinstance(viscosity_mu, str):
+            mu = SpatialConstant(viscosity_mu)
+        else:
+            warn("Passing a non-str 'viscosity_mu' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
+            mu = viscosity_mu
+
+        d = make_sym_vector("d", dim)
+        r = pymbolic_real_norm_2(d)
+        delta_ij = 1 if icomp == jcomp else 0
 
-        super().__init__(dim, icomp, jcomp, viscosity_mu, poisson_ratio)
+        if dim == 2:
+            expr = -var("log")(r)*delta_ij + d[icomp]*d[jcomp]/r**2
+            scaling = -1/(4*var("pi")*mu)
+        elif dim == 3:
+            expr = (1/r)*delta_ij + d[icomp]*d[jcomp]/r**3
+            scaling = -1/(8*var("pi")*mu)
+        else:
+            raise NotImplementedError(f"unsupported dimension: '{dim}'")
+
+        super().__init__(dim, expression=expr, global_scaling_const=scaling)
+        object.__setattr__(self, "icomp", icomp)
+        object.__setattr__(self, "jcomp", jcomp)
+        object.__setattr__(self, "viscosity_mu", mu)
 
     @override
     def __str__(self) -> str:
         return (
             f"StokesletKnl{self.dim}D_{self.icomp}{self.jcomp}"
-            f"({self.viscosity_mu}, {self.poisson_ratio})")
+            f"({self.viscosity_mu})")
+
+    @override
+    def __reduce__(self):
+        # TODO: remove this once we stop allowing non-str values for viscosity_mu
+        viscosity_mu = (
+            self.viscosity_mu.name
+            if isinstance(self.viscosity_mu, SpatialConstant)
+            else self.viscosity_mu)
+
+        return (
+            type(self),
+            (self.dim, self.icomp, self.jcomp, viscosity_mu))
+
+    @property
+    @override
+    def is_complex_valued(self) -> bool:
+        return False
+
+    @memoize_method
+    @override
+    def get_args(self) -> Sequence[KernelArgument]:
+        # TODO: remove this once we stop allowing non-str values for viscosity_mu
+        if isinstance(self.viscosity_mu, SpatialConstant):
+            return [
+                KernelArgument(
+                    loopy_arg=lp.ValueArg(self.viscosity_mu.name, np.float64),
+                )
+            ]
+        else:
+            return []
+
+    @override
+    def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
+        from sumpy.expansion.diff_op import laplacian, make_identity_diff_op
+
+        w = make_identity_diff_op(self.dim)
+        return laplacian(laplacian(w))
 
 
 @dataclass(frozen=True, repr=False)
 class StressletKernel(ExpressionKernel):
-    """
+    r"""A kernel for the Stokes equations (see e.g. Chapter 2 in [Pozrikidis1992]_).
+
+    .. math::
+
+        K_{ijk}(\mathbf{x}, \mathbf{y}) =
+            -P_j \delta_{ik}
+            + \mu (\partial_k K_{ij} + \partial_i K_{kj})
+
+    where the two-index :math:`K_{ij}` is the :class:`~sumpy.kernel.StokesletKernel`.
+    This kernel is often called the Stresslet and it represents the stress tensor.
+
     .. autoattribute:: icomp
     .. autoattribute:: jcomp
     .. autoattribute:: kcomp
@@ -947,36 +1088,37 @@ class StressletKernel(ExpressionKernel):
     """Component index for the kernel."""
     kcomp: int
     """Component index for the kernel."""
+
     viscosity_mu: float | SpatialConstant
     r"""The argument name to use for the dynamic viscosity :math:`\mu` when
-    generating functions to evaluate this kernel. Can also be a numeric value.
+    generating functions to evaluate this kernel.
     """
 
     def __init__(self,
-                 dim: int, icomp: int, jcomp: int, kcomp: int,
+                 dim: int,
+                 icomp: int,
+                 jcomp: int,
+                 kcomp: int,
                  viscosity_mu: float | str | SpatialConstant = "mu") -> None:
         # mu is unused but kept for consistency with the Stokeslet.
         if isinstance(viscosity_mu, str):
             mu = SpatialConstant(viscosity_mu)
         else:
+            warn("Passing a non-str 'viscosity_mu' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
             mu = viscosity_mu
 
+        d = make_sym_vector("d", dim)
+        r = pymbolic_real_norm_2(d)
+
         if dim == 2:
-            d = make_sym_vector("d", dim)
-            r = pymbolic_real_norm_2(d)
-            expr = (
-                d[icomp]*d[jcomp]*d[kcomp]/r**4
-                )
+            expr = d[icomp]*d[jcomp]*d[kcomp]/r**4
             scaling = 1/(var("pi"))
-
         elif dim == 3:
-            d = make_sym_vector("d", dim)
-            r = pymbolic_real_norm_2(d)
-            expr = (
-                d[icomp]*d[jcomp]*d[kcomp]/r**5
-                )
+            expr = d[icomp]*d[jcomp]*d[kcomp]/r**5
             scaling = 3/(4*var("pi"))
-
         else:
             raise NotImplementedError(f"unsupported dimension: '{dim}'")
 
@@ -989,10 +1131,15 @@ def __init__(self,
 
     @override
     def __reduce__(self) -> tuple[object, ...]:
+        # TODO: remove this once we stop allowing non-str values for viscosity_mu
+        viscosity_mu = (
+            self.viscosity_mu.name
+            if isinstance(self.viscosity_mu, SpatialConstant)
+            else self.viscosity_mu)
+
         return (
             self.__class__,
-            (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu),
-        )
+            (self.dim, self.icomp, self.jcomp, self.kcomp, viscosity_mu))
 
     @property
     @override
@@ -1017,6 +1164,7 @@ def get_args(self) -> Sequence[KernelArgument]:
     @override
     def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
         from sumpy.expansion.diff_op import laplacian, make_identity_diff_op
+
         w = make_identity_diff_op(self.dim)
         return laplacian(laplacian(w))
 
@@ -1024,9 +1172,9 @@ def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
 @dataclass(frozen=True, repr=False)
 class LineOfCompressionKernel(ExpressionKernel):
     """A kernel for the line of compression or dilatation of constant strength
-    along the axis "axis" from zero to negative infinity.
+    along an axis from zero to negative infinity.
 
-    This is used for the explicit solution to half-space Elasticity problem.
+    This is used for the explicit solution to half-space linear elasticity problem.
     See [Mindlin1936]_ for details.
 
     .. [Mindlin1936] R. D. Mindlin (1936).
@@ -1043,13 +1191,14 @@ class LineOfCompressionKernel(ExpressionKernel):
 
     axis: int
     """Axis number (defaulting to 2 for the z axis)."""
+
     viscosity_mu: float | SpatialConstant
     r"""The argument name to use for the dynamic viscosity :math:`\mu` when
-    generating functions to evaluate this kernel. Can also be a numeric value.
+    generating functions to evaluate this kernel.
     """
     poisson_ratio: float | SpatialConstant
     r"""The argument name to use for Poisson's ratio :math:`\nu` when
-    generating functions to evaluate this kernel. Can also be a numeric value.
+    generating functions to evaluate this kernel.
     """
 
     def __init__(self,
@@ -1061,15 +1210,28 @@ def __init__(self,
         if isinstance(viscosity_mu, str):
             mu = SpatialConstant(viscosity_mu)
         else:
+            warn("Passing a non-str 'viscosity_mu' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
             mu = viscosity_mu
+
         if isinstance(poisson_ratio, str):
             nu = SpatialConstant(poisson_ratio)
         else:
+            warn("Passing a non-str 'poisson_ratio' is deprecated. This will "
+                 "raise a ValueError starting with Q1 2027.",
+                 DeprecationWarning, stacklevel=2)
+
+            if poisson_ratio == 0.5:  # noqa: RUF069
+                raise ValueError(f"{type(self)} cannot accept 'poisson_ratio' of 0.5")
+
             nu = poisson_ratio
 
         if dim == 3:
             d = make_sym_vector("d", dim)
             r = pymbolic_real_norm_2(d)
+
             # Kelvin solution
             expr = d[axis] * var("log")(r + d[axis]) - r
             scaling = (1 - 2*nu)/(4*var("pi")*mu)
@@ -1082,11 +1244,27 @@ def __init__(self,
         object.__setattr__(self, "viscosity_mu", mu)
         object.__setattr__(self, "poisson_ratio", nu)
 
+    @override
+    def __str__(self) -> str:
+        return (
+            f"LineOfCompressionKnl{self.dim}D_{self.axis}"
+            f"({self.viscosity_mu}, {self.poisson_ratio})")
+
     @override
     def __reduce__(self) -> tuple[object, ...]:
+        # TODO: remove this once we stop allowing non-str values for viscosity_mu
+        viscosity_mu = (
+            self.viscosity_mu.name
+            if isinstance(self.viscosity_mu, SpatialConstant)
+            else self.viscosity_mu)
+        poisson_ratio = (
+            self.poisson_ratio.name
+            if isinstance(self.poisson_ratio, SpatialConstant)
+            else self.poisson_ratio)
+
         return (
             self.__class__,
-            (self.dim, self.axis, self.viscosity_mu, self.poisson_ratio),
+            (self.dim, self.axis, viscosity_mu, poisson_ratio),
         )
 
     @property
@@ -1094,10 +1272,6 @@ def __reduce__(self) -> tuple[object, ...]:
     def is_complex_valued(self) -> bool:
         return False
 
-    @override
-    def __str__(self) -> str:
-        return f"LineOfCompressionKnl{self.dim}D_{self.axis}"
-
     @memoize_method
     @override
     def get_args(self) -> Sequence[KernelArgument]:
@@ -1115,6 +1289,7 @@ def get_args(self) -> Sequence[KernelArgument]:
     @override
     def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
         from sumpy.expansion.diff_op import laplacian, make_identity_diff_op
+
         w = make_identity_diff_op(self.dim)
         return laplacian(w)
 
@@ -1589,6 +1764,7 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> Kernel:
     map_yukawa_kernel = map_expression_kernel
     map_elasticity_kernel = map_expression_kernel
     map_line_of_compression_kernel = map_expression_kernel
+    map_stokeslet_kernel = map_expression_kernel
     map_stresslet_kernel = map_expression_kernel
 
     def map_axis_target_derivative(self, kernel: AxisTargetDerivative) -> Kernel:
@@ -1660,7 +1836,9 @@ def map_expression_kernel(self, kernel: ExpressionKernel) -> int:
     map_biharmonic_kernel = map_expression_kernel
     map_helmholtz_kernel = map_expression_kernel
     map_yukawa_kernel = map_expression_kernel
+    map_elasticity_kernel = map_expression_kernel
     map_line_of_compression_kernel = map_expression_kernel
+    map_stokeslet_kernel = map_expression_kernel
     map_stresslet_kernel = map_expression_kernel
 
     @override
diff --git a/sumpy/test/test_misc.py b/sumpy/test/test_misc.py
index ce85f30b..9be2a919 100644
--- a/sumpy/test/test_misc.py
+++ b/sumpy/test/test_misc.py
@@ -497,22 +497,20 @@ def test_as_scalar_pde_elasticity():
 # }}}
 
 
-# {{{ test_elasticity_new
+# {{{ test_elasticity_pickle
 
-def test_elasticity_new():
+def test_elasticity_pickle():
     from pickle import dumps, loads
-    stokes_knl = StokesletKernel(3, 0, 1, "mu1", 0.5)
-    stokes_knl2 = ElasticityKernel(3, 0, 1, "mu1", 0.5)
-    elasticity_knl = ElasticityKernel(3, 0, 1, "mu1", "nu")
-    elasticity_helper_knl = LineOfCompressionKernel(3, 0, "mu1", "nu")
+    stokes_knl = StokesletKernel(
+        3, 0, 1, viscosity_mu="mu1")
+    elasticity_knl = ElasticityKernel(
+        3, 0, 1, viscosity_mu="mu1", poisson_ratio="nu1")
+    elasticity_helper_knl = LineOfCompressionKernel(
+        3, 0, viscosity_mu="mu1", poisson_ratio="nu1")
 
-    assert isinstance(stokes_knl2, StokesletKernel)
-    assert stokes_knl == stokes_knl2
     assert loads(dumps(stokes_knl)) == stokes_knl
-
-    for knl in [elasticity_knl, elasticity_helper_knl]:
-        assert not isinstance(knl, StokesletKernel)
-        assert loads(dumps(knl)) == knl
+    assert loads(dumps(elasticity_knl)) == elasticity_knl
+    assert loads(dumps(elasticity_helper_knl)) == elasticity_helper_knl
 
 # }}}