diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index adb858a119921812cf0c9b7084a2cedfb0397856..50fb8e5887266235ac7864fee91986217e06be96 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -441,16 +441,16 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
     # {{{ fmm-based execution
 
     @memoize_method
-    def expansion_wrangler_code_container(self, base_kernel, out_kernels):
+    def expansion_wrangler_code_container(self, fmm_kernel, out_kernels):
         mpole_expn_class = \
-                self.expansion_factory.get_multipole_expansion_class(base_kernel)
+                self.expansion_factory.get_multipole_expansion_class(fmm_kernel)
         local_expn_class = \
-                self.expansion_factory.get_local_expansion_class(base_kernel)
+                self.expansion_factory.get_local_expansion_class(fmm_kernel)
 
         from functools import partial
-        fmm_mpole_factory = partial(mpole_expn_class, base_kernel)
-        fmm_local_factory = partial(local_expn_class, base_kernel)
-        qbx_local_factory = partial(local_expn_class, base_kernel)
+        fmm_mpole_factory = partial(mpole_expn_class, fmm_kernel)
+        fmm_local_factory = partial(local_expn_class, fmm_kernel)
+        qbx_local_factory = partial(local_expn_class, fmm_kernel)
 
         if self.fmm_backend == "sumpy":
             from pytential.qbx.fmm import \
@@ -514,67 +514,23 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
         strengths = (evaluate(insn.density).with_queue(queue)
                 * self.weights_and_area_elements())
 
-        # {{{ get expansion wrangler
-
-        base_kernel = None
-        out_kernels = []
-
-        from sumpy.kernel import AxisTargetDerivativeRemover
-        for knl in insn.kernels:
-            candidate_base_kernel = AxisTargetDerivativeRemover()(knl)
-
-            if base_kernel is None:
-                base_kernel = candidate_base_kernel
-            else:
-                assert base_kernel == candidate_base_kernel
-
         out_kernels = tuple(knl for knl in insn.kernels)
-
-        if base_kernel.is_complex_valued or strengths.dtype.kind == "c":
-            value_dtype = self.complex_dtype
-        else:
-            value_dtype = self.real_dtype
-
-        # {{{ build extra_kwargs dictionaries
-
-        # This contains things like the Helmholtz parameter k or
-        # the normal directions for double layers.
-
-        def reorder_sources(source_array):
-            if isinstance(source_array, cl.array.Array):
-                return (source_array
-                        .with_queue(queue)
-                        [geo_data.tree().user_source_ids]
-                        .with_queue(None))
-            else:
-                return source_array
-
-        kernel_extra_kwargs = {}
-        source_extra_kwargs = {}
-
-        from sumpy.tools import gather_arguments, gather_source_arguments
-        from pytools.obj_array import with_object_array_or_scalar
-        for func, var_dict in [
-                (gather_arguments, kernel_extra_kwargs),
-                (gather_source_arguments, source_extra_kwargs),
-                ]:
-            for arg in func(out_kernels):
-                var_dict[arg.name] = with_object_array_or_scalar(
-                        reorder_sources,
-                        evaluate(insn.kernel_arguments[arg.name]))
-
-        # }}}
+        fmm_kernel = self.get_fmm_kernel(out_kernels)
+        output_and_expansion_dtype = (
+                self.get_fmm_output_and_expansion_dtype(fmm_kernel, strengths))
+        kernel_extra_kwargs, source_extra_kwargs = (
+                self.get_fmm_expansion_wrangler_extra_kwargs(
+                    queue, out_kernels, geo_data.tree().user_source_ids,
+                    insn.kernel_arguments, evaluate))
 
         wrangler = self.expansion_wrangler_code_container(
-                base_kernel, out_kernels).get_wrangler(
-                        queue, geo_data, value_dtype,
+                fmm_kernel, out_kernels).get_wrangler(
+                        queue, geo_data, output_and_expansion_dtype,
                         self.qbx_order,
                         self.fmm_level_to_order,
                         source_extra_kwargs=source_extra_kwargs,
                         kernel_extra_kwargs=kernel_extra_kwargs)
 
-        # }}}
-
         from pytential.qbx.geometry import target_state
         if (geo_data.user_target_to_center().with_queue(queue)
                 == target_state.FAILED).any().get():
diff --git a/pytential/source.py b/pytential/source.py
index 6d78da5723b10dc30ae39569aba9628c422c4b7a..75d0db70ee6babe762f9589f67fc8da877fd2379 100644
--- a/pytential/source.py
+++ b/pytential/source.py
@@ -206,6 +206,60 @@ class LayerPotentialSourceBase(PotentialSource):
 
         return p2p
 
+    # {{{ fmm setup helpers
+
+    def get_fmm_kernel(self, kernels):
+        fmm_kernel = None
+
+        from sumpy.kernel import AxisTargetDerivativeRemover
+        for knl in kernels:
+            candidate_fmm_kernel = AxisTargetDerivativeRemover()(knl)
+
+            if fmm_kernel is None:
+                fmm_kernel = candidate_fmm_kernel
+            else:
+                assert fmm_kernel == candidate_fmm_kernel
+
+        return fmm_kernel
+
+    def get_fmm_output_and_expansion_dtype(self, base_kernel, strengths):
+        if base_kernel.is_complex_valued or strengths.dtype.kind == "c":
+            return self.complex_dtype
+        else:
+            return self.real_dtype
+
+    def get_fmm_expansion_wrangler_extra_kwargs(
+            self, queue, out_kernels, tree_user_source_ids, arguments, evaluator):
+        # This contains things like the Helmholtz parameter k or
+        # the normal directions for double layers.
+
+        def reorder_sources(source_array):
+            if isinstance(source_array, cl.array.Array):
+                return (source_array
+                        .with_queue(queue)
+                        [tree_user_source_ids]
+                        .with_queue(None))
+            else:
+                return source_array
+
+        kernel_extra_kwargs = {}
+        source_extra_kwargs = {}
+
+        from sumpy.tools import gather_arguments, gather_source_arguments
+        from pytools.obj_array import with_object_array_or_scalar
+        for func, var_dict in [
+                (gather_arguments, kernel_extra_kwargs),
+                (gather_source_arguments, source_extra_kwargs),
+                ]:
+            for arg in func(out_kernels):
+                var_dict[arg.name] = with_object_array_or_scalar(
+                        reorder_sources,
+                        evaluator(arguments[arg.name]))
+
+        return kernel_extra_kwargs, source_extra_kwargs
+
+    # }}}
+
     # {{{ weights and area elements
 
     @memoize_method
diff --git a/pytential/unregularized.py b/pytential/unregularized.py
index 4ea94bf1094fd3c820b786e3f814d43e67a7e347..fcd5298e798c41bc9853c2d64491d0e84fce2c3c 100644
--- a/pytential/unregularized.py
+++ b/pytential/unregularized.py
@@ -27,7 +27,15 @@ THE SOFTWARE.
 
 import six
 
+import numpy as np
+import loopy as lp
+
+from boxtree.tools import DeviceDataRecord
 from pytential.source import LayerPotentialSourceBase
+from pytools import memoize_method
+
+import pyopencl as cl
+import pyopencl.array  # noqa
 
 import logging
 logger = logging.getLogger(__name__)
@@ -43,17 +51,43 @@ __doc__ = """
 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.
+
+    .. attribute:: fmm_level_to_order
     """
 
     def __init__(self, density_discr,
+            fmm_order=False,
+            fmm_level_to_order=None,
+            expansion_factory=None,
             # begin undocumented arguments
             # FIXME default debug=False once everything works
             debug=True):
         """
+        :arg fmm_order: `False` for direct calculation.
         """
         self.density_discr = density_discr
         self.debug = debug
 
+        if fmm_order is not False and fmm_level_to_order is not None:
+            raise TypeError("may not specify both fmm_order and fmm_level_to_order")
+
+        if fmm_level_to_order is None:
+            if fmm_order is not False:
+                def fmm_level_to_order(level):
+                    return fmm_order
+            else:
+                fmm_level_to_order = False
+
+        self.density_discr = density_discr
+        self.fmm_level_to_order = fmm_level_to_order
+
+        if expansion_factory is None:
+            from sumpy.expansion import DefaultExpansionFactory
+            expansion_factory = DefaultExpansionFactory()
+        self.expansion_factory = expansion_factory
+
+        self.debug = debug
+
     @property
     def fine_density_discr(self):
         return self.density_discr
@@ -67,9 +101,12 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
     def copy(
             self,
             density_discr=None,
-            debug=None
+            fmm_level_to_order=None,
+            debug=None,
             ):
         return type(self)(
+                fmm_level_to_order=(
+                    fmm_level_to_order or self.fmm_level_to_order),
                 density_discr=density_discr or self.density_discr,
                 debug=debug if debug is not None else self.debug)
 
@@ -80,7 +117,11 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
             value = evaluate(expr)
             return with_object_array_or_scalar(lambda x: x, value)
 
-        func = self.exec_compute_potential_insn_direct
+        if self.fmm_level_to_order is False:
+            func = self.exec_compute_potential_insn_direct
+        else:
+            func = self.exec_compute_potential_insn_fmm
+
         return func(queue, insn, bound_expr, evaluate_wrapper)
 
     def op_group_features(self, expr):
@@ -126,6 +167,259 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
 
         return result, []
 
+    # {{{ fmm-based execution
+
+    @memoize_method
+    def expansion_wrangler_code_container(self, fmm_kernel, out_kernels):
+        mpole_expn_class = \
+                self.expansion_factory.get_multipole_expansion_class(fmm_kernel)
+        local_expn_class = \
+                self.expansion_factory.get_local_expansion_class(fmm_kernel)
+
+        from functools import partial
+        fmm_mpole_factory = partial(mpole_expn_class, fmm_kernel)
+        fmm_local_factory = partial(local_expn_class, fmm_kernel)
+
+        from sumpy.fmm import SumpyExpansionWranglerCodeContainer
+        return SumpyExpansionWranglerCodeContainer(
+                self.cl_context,
+                fmm_mpole_factory,
+                fmm_local_factory,
+                out_kernels)
+
+    @property
+    @memoize_method
+    def fmm_geometry_code_container(self):
+        return _FMMGeometryCodeContainer(
+                self.cl_context, self.ambient_dim, self.debug)
+
+    def fmm_geometry_data(self, targets):
+        return _FMMGeometryData(
+                self,
+                self.fmm_geometry_code_container,
+                targets,
+                self.debug)
+
+    def exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate):
+
+        # {{{ gather unique target discretizations used
+
+        target_name_to_index = {}
+        targets = []
+
+        for o in insn.outputs:
+            assert o.qbx_forced_limit not in (-1, 1)
+
+            if o.target_name in target_name_to_index:
+                continue
+
+            target_name_to_index[o.target_name] = len(targets)
+            targets.append(bound_expr.places[o.target_name])
+
+        targets = tuple(targets)
+
+        # }}}
+
+        # {{{ get wrangler
+
+        geo_data = self.fmm_geometry_data(targets)
+
+        strengths = (evaluate(insn.density).with_queue(queue)
+                * self.weights_and_area_elements())
+
+        out_kernels = tuple(knl for knl in insn.kernels)
+        fmm_kernel = self.get_fmm_kernel(out_kernels)
+        output_and_expansion_dtype = (
+                self.get_fmm_output_and_expansion_dtype(fmm_kernel, strengths))
+        kernel_extra_kwargs, source_extra_kwargs = (
+                self.get_fmm_expansion_wrangler_extra_kwargs(
+                    queue, out_kernels, geo_data.tree().user_source_ids,
+                    insn.kernel_arguments, evaluate))
+
+        wrangler = self.expansion_wrangler_code_container(
+                fmm_kernel, out_kernels).get_wrangler(
+                    queue,
+                    geo_data.tree(),
+                    output_and_expansion_dtype,
+                    self.fmm_level_to_order,
+                    source_extra_kwargs=source_extra_kwargs,
+                    kernel_extra_kwargs=kernel_extra_kwargs)
+
+        # }}}
+
+        from boxtree.fmm import drive_fmm
+        all_potentials_on_every_tgt = drive_fmm(
+                geo_data.traversal(), wrangler, strengths)
+
+        # {{{ postprocess fmm
+
+        result = []
+
+        for o in insn.outputs:
+            target_index = target_name_to_index[o.target_name]
+            target_slice = slice(*geo_data.target_info().target_discr_starts[
+                    target_index:target_index+2])
+
+            result.append(
+                    (o.name,
+                        all_potentials_on_every_tgt[o.kernel_index][target_slice]))
+
+        # }}}
+
+        return result, []
+
+    # }}}
+
+# }}}
+
+
+# {{{ fmm tools
+
+class _FMMGeometryCodeContainer(object):
+
+    def __init__(self, cl_context, ambient_dim, debug):
+        self.cl_context = cl_context
+        self.ambient_dim = ambient_dim
+        self.debug = debug
+
+    @memoize_method
+    def copy_targets_kernel(self):
+        knl = lp.make_kernel(
+            """{[dim,i]:
+                0<=dim<ndims and
+                0<=i<npoints}""",
+            """
+                targets[dim, i] = points[dim, i]
+                """,
+            default_offset=lp.auto, name="copy_targets")
+
+        knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
+
+        knl = lp.split_iname(knl, "i", 128, inner_tag="l.0", outer_tag="g.0")
+        knl = lp.tag_array_axes(knl, "points", "sep, C")
+
+        knl = lp.tag_array_axes(knl, "targets", "stride:auto, stride:1")
+        return lp.tag_inames(knl, dict(dim="ilp"))
+
+    @property
+    @memoize_method
+    def build_tree(self):
+        from boxtree import TreeBuilder
+        return TreeBuilder(self.cl_context)
+
+    @property
+    @memoize_method
+    def build_traversal(self):
+        from boxtree.traversal import FMMTraversalBuilder
+        return FMMTraversalBuilder(self.cl_context)
+
+
+class _TargetInfo(DeviceDataRecord):
+    """
+    .. attribute:: targets
+
+        Shape: ``[dim,ntargets]``
+
+    .. attribute:: target_discr_starts
+
+        Shape: ``[ndiscrs+1]``
+
+    .. attribute:: ntargets
+    """
+
+
+class _FMMGeometryData(object):
+
+    def __init__(self, lpot_source, code_getter, target_discrs, debug=True):
+        self.lpot_source = lpot_source
+        self.code_getter = code_getter
+        self.target_discrs = target_discrs
+        self.debug = debug
+
+    @property
+    def cl_context(self):
+        return self.lpot_source.cl_context
+
+    @property
+    def coord_dtype(self):
+        return self.lpot_source.fine_density_discr.nodes().dtype
+
+    @property
+    def ambient_dim(self):
+        return self.lpot_source.fine_density_discr.ambient_dim
+
+    @memoize_method
+    def traversal(self):
+        with cl.CommandQueue(self.cl_context) as queue:
+            trav, _ = self.code_getter.build_traversal(queue, self.tree(),
+                    debug=self.debug)
+
+            return trav
+
+    @memoize_method
+    def tree(self):
+        """Build and return a :class:`boxtree.tree.Tree`
+        for this source with these targets.
+
+        |cached|
+        """
+
+        code_getter = self.code_getter
+        lpot_src = self.lpot_source
+        target_info = self.target_info()
+
+        with cl.CommandQueue(self.cl_context) as queue:
+            nsources = lpot_src.fine_density_discr.nnodes
+            nparticles = nsources + target_info.ntargets
+
+            refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32)
+            refine_weights[:nsources] = 1
+            refine_weights.finish()
+
+            MAX_LEAF_REFINE_WEIGHT = 32  # noqa
+
+            tree, _ = code_getter.build_tree(queue,
+                    particles=lpot_src.fine_density_discr.nodes(),
+                    targets=target_info.targets,
+                    max_leaf_refine_weight=MAX_LEAF_REFINE_WEIGHT,
+                    refine_weights=refine_weights,
+                    debug=self.debug,
+                    kind="adaptive")
+
+            return tree
+
+    @memoize_method
+    def target_info(self):
+        code_getter = self.code_getter
+        lpot_src = self.lpot_source
+        target_discrs = self.target_discrs
+
+        with cl.CommandQueue(self.cl_context) as queue:
+            ntargets = 0
+            target_discr_starts = []
+
+            for target_discr in target_discrs:
+                target_discr_starts.append(ntargets)
+                ntargets += target_discr.nnodes
+
+            target_discr_starts.append(ntargets)
+
+            targets = cl.array.empty(
+                    self.cl_context,
+                    (lpot_src.ambient_dim, ntargets),
+                    self.coord_dtype)
+
+            for start, target_discr in zip(target_discr_starts, target_discrs):
+                code_getter.copy_targets_kernel()(
+                        queue,
+                        targets=targets[:, start:start+target_discr.nnodes],
+                        points=target_discr.nodes())
+
+            return _TargetInfo(
+                    targets=targets,
+                    target_discr_starts=target_discr_starts,
+                    ntargets=ntargets).with_queue(None)
+
 # }}}
 
 
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index 280a918ca3399a3caa8427e1334d6d69babd0699..31e9cc5607ddd1e379ed4e43851b38dec013ae79 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -470,6 +470,50 @@ def test_unregularized_with_ones_kernel(ctx_getter):
     assert np.allclose(result_self.get(), 2 * np.pi)
     assert np.allclose(result_nonself.get(), 2 * np.pi)
 
+
+def test_unregularized_off_surface_fmm_vs_direct(ctx_getter):
+    cl_ctx = ctx_getter()
+    queue = cl.CommandQueue(cl_ctx)
+
+    nelements = 300
+    target_order = 8
+    fmm_order = 4
+
+    mesh = make_curve_mesh(WobblyCircle.random(8, seed=30),
+                np.linspace(0, 1, nelements+1),
+                target_order)
+
+    from pytential.unregularized import UnregularizedLayerPotentialSource
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
+
+    density_discr = Discretization(
+            cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
+    direct = UnregularizedLayerPotentialSource(
+            density_discr,
+            fmm_order=False,
+            )
+    fmm = direct.copy(fmm_level_to_order=lambda _: fmm_order)
+
+    sigma = density_discr.zeros(queue) + 1
+
+    fplot = FieldPlotter(np.zeros(2), extent=5, npoints=100)
+    from pytential.target import PointsTarget
+    ptarget = PointsTarget(fplot.points)
+    from sumpy.kernel import LaplaceKernel
+
+    op = sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=None)
+
+    direct_fld_in_vol = bind((direct, ptarget), op)(queue, sigma=sigma)
+    fmm_fld_in_vol = bind((fmm, ptarget), op)(queue, sigma=sigma)
+
+    err = cl.clmath.fabs(fmm_fld_in_vol - direct_fld_in_vol)
+
+    linf_err = cl.array.max(err).get()
+    print("l_inf error:", linf_err)
+    assert linf_err < 5e-3
+
 # }}}