diff --git a/doc/discretization.rst b/doc/discretization.rst
index 3229e267a0e94ac4542c2cdc5402cf4a1c1f3d75..ad2eb738a04b14f72044a72f07120787a007f535 100644
--- a/doc/discretization.rst
+++ b/doc/discretization.rst
@@ -15,6 +15,11 @@ and you can start computing.
 
 .. automodule:: pytential.qbx
 
+Unregularized discretization
+-------
+
+.. automodule:: pytential.unregularized
+
 Sources
 -------
 
diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index 2212b0057477f84c9cac14838d0f003e18d9eaaa..0c4bfdf366d9adb96eea1cdfe298e9345afac07e 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -589,22 +589,6 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
                     for knl in kernels],
                 value_dtypes=value_dtype)
 
-    @memoize_method
-    def get_p2p(self, kernels):
-        # needs to be separate method for caching
-
-        from pytools import any
-        if any(knl.is_complex_valued for knl in kernels):
-            value_dtype = self.density_discr.complex_dtype
-        else:
-            value_dtype = self.density_discr.real_dtype
-
-        from sumpy.p2p import P2P
-        p2p = P2P(self.cl_context,
-                    kernels, exclude_self=False, value_dtypes=value_dtype)
-
-        return p2p
-
     @memoize_method
     def get_qbx_target_numberer(self, dtype):
         assert dtype == np.int32
diff --git a/pytential/source.py b/pytential/source.py
index f4dccda518392c1468b181792b2b1df093377c1a..6d78da5723b10dc30ae39569aba9628c422c4b7a 100644
--- a/pytential/source.py
+++ b/pytential/source.py
@@ -190,6 +190,22 @@ class LayerPotentialSourceBase(PotentialSource):
     def complex_dtype(self):
         return self.density_discr.complex_dtype
 
+    @memoize_method
+    def get_p2p(self, kernels):
+        # needs to be separate method for caching
+
+        from pytools import any
+        if any(knl.is_complex_valued for knl in kernels):
+            value_dtype = self.density_discr.complex_dtype
+        else:
+            value_dtype = self.density_discr.real_dtype
+
+        from sumpy.p2p import P2P
+        p2p = P2P(self.cl_context,
+                  kernels, exclude_self=False, value_dtypes=value_dtype)
+
+        return p2p
+
     # {{{ weights and area elements
 
     @memoize_method
diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py
index 3eef6c4eda28200e8528250dbf6882af9bf433ee..47319841fa52fc7a6085438fb87384c406f6ede4 100644
--- a/pytential/symbolic/mappers.py
+++ b/pytential/symbolic/mappers.py
@@ -340,6 +340,33 @@ class DerivativeBinder(DerivativeBinderBase, IdentityMapper):
 # }}}
 
 
+# {{{ Unregularized preprocessor
+
+class UnregularizedPreprocessor(IdentityMapper):
+
+    def __init__(self, source_name, places):
+        self.source_name = source_name
+        self.places = places
+
+    def map_int_g(self, expr):
+        if expr.qbx_forced_limit in (-1, 1):
+            raise ValueError(
+                    "Unregularized evaluation does not support one-sided limits")
+
+        expr = expr.copy(
+                qbx_forced_limit=None,
+                kernel=expr.kernel,
+                density=self.rec(expr.density),
+                kernel_arguments=dict(
+                    (name, self.rec(arg_expr))
+                    for name, arg_expr in expr.kernel_arguments.items()
+                    ))
+
+        return expr
+
+# }}}
+
+
 # {{{ QBX preprocessor
 
 class QBXPreprocessor(IdentityMapper):
diff --git a/pytential/unregularized.py b/pytential/unregularized.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ea94bf1094fd3c820b786e3f814d43e67a7e347
--- /dev/null
+++ b/pytential/unregularized.py
@@ -0,0 +1,136 @@
+# -*- coding: utf-8 -*-
+from __future__ import division, absolute_import
+
+__copyright__ = """
+Copyright (C) 2017 Matt Wala
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import six
+
+from pytential.source import LayerPotentialSourceBase
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+__doc__ = """
+.. autoclass:: UnregularizedLayerPotentialSource
+"""
+
+
+# {{{ (panel-based) unregularized layer potential source
+
+class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
+    """A source discretization for a layer potential discretized with a Nyström
+    method that uses panel-based quadrature and does not modify the kernel.
+    """
+
+    def __init__(self, density_discr,
+            # begin undocumented arguments
+            # FIXME default debug=False once everything works
+            debug=True):
+        """
+        """
+        self.density_discr = density_discr
+        self.debug = debug
+
+    @property
+    def fine_density_discr(self):
+        return self.density_discr
+
+    def resampler(self, queue, f):
+        return f
+
+    def with_refinement(self):
+        raise NotImplementedError
+
+    def copy(
+            self,
+            density_discr=None,
+            debug=None
+            ):
+        return type(self)(
+                density_discr=density_discr or self.density_discr,
+                debug=debug if debug is not None else self.debug)
+
+    def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate):
+        from pytools.obj_array import with_object_array_or_scalar
+
+        def evaluate_wrapper(expr):
+            value = evaluate(expr)
+            return with_object_array_or_scalar(lambda x: x, value)
+
+        func = self.exec_compute_potential_insn_direct
+        return func(queue, insn, bound_expr, evaluate_wrapper)
+
+    def op_group_features(self, expr):
+        from sumpy.kernel import AxisTargetDerivativeRemover
+        result = (
+                expr.source, expr.density,
+                AxisTargetDerivativeRemover()(expr.kernel),
+                )
+
+        return result
+
+    def preprocess_optemplate(self, name, discretizations, expr):
+        """
+        :arg name: The symbolic name for *self*, which the preprocessor
+            should use to find which expressions it is allowed to modify.
+        """
+        from pytential.symbolic.mappers import UnregularizedPreprocessor
+        return UnregularizedPreprocessor(name, discretizations)(expr)
+
+    def exec_compute_potential_insn_direct(self, queue, insn, bound_expr, evaluate):
+        kernel_args = {}
+
+        for arg_name, arg_expr in six.iteritems(insn.kernel_arguments):
+            kernel_args[arg_name] = evaluate(arg_expr)
+
+        strengths = (evaluate(insn.density).with_queue(queue)
+                * self.weights_and_area_elements())
+
+        result = []
+        p2p = None
+
+        for o in insn.outputs:
+            target_discr = bound_expr.get_discretization(o.target_name)
+
+            if p2p is None:
+                p2p = self.get_p2p(insn.kernels)
+
+            evt, output_for_each_kernel = p2p(queue,
+                    target_discr.nodes(), self.fine_density_discr.nodes(),
+                    [strengths], **kernel_args)
+
+            result.append((o.name, output_for_each_kernel[o.kernel_index]))
+
+        return result, []
+
+# }}}
+
+
+__all__ = (
+        UnregularizedLayerPotentialSource,
+        )
+
+# vim: fdm=marker
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index 0f3d35bb9394359433c376394303121884456a63..994c96735834eb6c76f36956731a1e84644c1f12 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -1273,6 +1273,52 @@ def test_off_surface_eval_vs_direct(ctx_getter,  do_plot=False):
 # }}}
 
 
+# {{{ unregularized tests
+
+
+def test_unregularized_with_ones_kernel(ctx_getter):
+    cl_ctx = ctx_getter()
+    queue = cl.CommandQueue(cl_ctx)
+
+    nelements = 10
+    order = 8
+
+    mesh = make_curve_mesh(partial(ellipse, 1),
+            np.linspace(0, 1, nelements+1),
+            order)
+
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
+
+    discr = Discretization(cl_ctx, mesh,
+            InterpolatoryQuadratureSimplexGroupFactory(order))
+
+    from pytential.unregularized import UnregularizedLayerPotentialSource
+    lpot_src = UnregularizedLayerPotentialSource(discr)
+
+    from sumpy.kernel import one_kernel_2d
+
+    expr = sym.IntG(one_kernel_2d, sym.var("sigma"), qbx_forced_limit=None)
+
+    from pytential.target import PointsTarget
+    op_self = bind(lpot_src, expr)
+    op_nonself = bind((lpot_src, PointsTarget(np.zeros((2, 1), dtype=float))), expr)
+
+    with cl.CommandQueue(cl_ctx) as queue:
+        sigma = cl.array.zeros(queue, discr.nnodes, dtype=float)
+        sigma.fill(1)
+        sigma.finish()
+
+        result_self = op_self(queue, sigma=sigma)
+        result_nonself = op_nonself(queue, sigma=sigma)
+
+    assert np.allclose(result_self.get(), 2 * np.pi)
+    assert np.allclose(result_nonself.get(), 2 * np.pi)
+
+# }}}
+
+
 # You can test individual routines by typing
 # $ python test_layer_pot.py 'test_routine()'